%% Assess the performance of optimization
% Use lysimeter time series as "should" values and try to estimate actual
% values based on reference/penman-monteith proxy.
% Results will be printed as figures
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));
close all
clear
clc

%% Main settings
% Set input files for Evapotranspiration and Lysimeter. These time series
% will be used for optimization. Set the 'time_shift' parameter to adjust
% heterogeneous time stamp convention (some time series point ot previous
% 24 hours, other to current hour/day...)
% Reference/potential evapotranspiration (grass) 
evap.in_file = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf'); 
evap.in_column = 4;
% Lysimeter ET/Rain = as iGrav signal to assess the method accucary
lysi.in_file = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf'); 
lysi.in_column = 3;

% Set time period analysed (all points outside this interval will not be
% used)
start = [2015 07 24];
stop =  [2016 06 28];
% Set gravimeter precission (will be used to account for offset due to
% noise = the first point of the gravity time series will thus not
% determine the level of the time series)
response = 0.478;
grav_prec = 3; % nm/s^2

% Optimization parameters: [evaporation, offset]
param_spec.ini_val = [-1,  1];
param_spec.max_val = [ 0,  grav_prec/response];
param_spec.min_val = [-1.0, -grav_prec/response];

% Set length of the moving window in days
wind_length = [9,11,13]; % scalar or vector

% Plot settings
font_size = 11;

%% Load data
% Evapotranspiration (proxy): hourly
[evap.time,temp] = loadtsf(evap.in_file);
temp(evap.time<datenum(start),:) = [];
evap.time(evap.time<datenum(start),:) = [];
evap.data = cumsum(temp(:,evap.in_column));
% Evapotranspiration: day
[evap.time_day,evap.data_day] = data2daily(evap.time,temp(:,evap.in_column),2,0);
evap.time_day = datenum(evap.time_day) + 12/24;
% Cut data outside required time interval
evap.data_day = evap.data_day(evap.time_day>datenum(start) & evap.time_day<datenum(stop));
evap.time_day = evap.time_day(evap.time_day>datenum(start) & evap.time_day<datenum(stop));

% Lysimeter (~gravimeter signal)
[lysi.time,temp] = loadtsf(lysi.in_file);
temp(lysi.time<datenum(start),:) = [];
lysi.time(lysi.time<datenum(start),:) = [];
lysi.data = temp(:,lysi.in_column);
lysi.data = lysi.data(lysi.time>datenum(start) & lysi.time<datenum(stop));
lysi.time = lysi.time(lysi.time>datenum(start) & lysi.time<datenum(stop));
[lysi.time_day,lysi.data_day] = data2daily(lysi.time,lysi.data,2,0);
lysi.time_day = datenum(lysi.time_day) + 12/24; % add 12 hours to point to the noon. Independent of 'time_shift' parameter
    
for s = 1:length(wind_length)
    %% Declare output variables
    % Optimized parameters = evaporation multiplier and offset
    param = zeros(length(lysi.time),length(param_spec.ini_val))*NaN;
    % Min slope = only gravity time series (window) with slope smaller then
    % this value will be used
    pmin = 0;

    %% Main computation
    % pmin = polyfit([0,wind_length*2],[0,-grav_prec/response],1);
    % Run loop movinw the window along the whole time series (except edges)
    for i = wind_length(s):length(lysi.time)-wind_length(s)
        % Show current status of computation (percentage)
        disp(i/length(lysi.time)*100);
        % Use lysimeter data as input instead of gravity
        in_grav = -cumsum(lysi.data(lysi.time>lysi.time(i)-wind_length(s)/2 & lysi.time<=lysi.time(i)+wind_length(s)/2));
        in_grav = in_grav - in_grav(1);
        in_time = lysi.time(lysi.time>lysi.time(i)-wind_length(s)/2 & lysi.time<=lysi.time(i)+wind_length(s)/2);

        % Get P, EV and runoff time series for identical time period. Always
        % subtract first value so all series start at the same point (the
        % estimated offset is limited to a certain range)
        in_evap = evap.data(evap.time>lysi.time(i)-wind_length(s)/2 & evap.time<=lysi.time(i)+wind_length(s)/2);
        in_evap = in_evap - in_evap(1);

        % Call optimization function only if input gravity time series does not
        % contain NaNs and all input time series have same length
        try 
            % Declare bias/offset variable
            in_bias = ones(length(in_evap),1);
            % Try to compute. Will result in an error if dimensions are not
            % correct = missing data is present
            temp = in_evap + in_grav;
            % Optimization will be carried out only if error_flag is equal 0
            error_flag = 0;
        catch
            error_flag = 1;
        end
        % Check for errors and NaNs
        if error_flag == 0 && (sum(~isnan(in_grav)) == length(in_grav)) 
            % Optimization. Use function only if gravity decreased during the
            % selected interval = only negative (below threshold) slope is
            % accepted 
            p = polyfit(in_time - in_time(1),in_grav,1);
            if p(1) <= pmin
                % Use -in_disc as simple adding will be performed (unlike in
                % case of in_evap where optimized parameter can have only
                % negative values)
                param(i,:) = Optimization_iGrav_TWS_daily_function(in_grav,[in_evap.*0,in_evap.*0],[in_evap,in_bias],param_spec);
            end
        end
        clc
    end

    %% Aggregate to daily values
    param_hour(:,s) = param(:,1);
    [time_out,param_out(:,s)] = data2daily(lysi.time,param_hour(:,s),1,0);
    time_out = datenum(time_out)+0.5;
    et_mean(:,s) = -param_out(:,s).*evap.data_day;

    %% Show the results
    figure('Position',[150,100,1700,900],'paperpositionmode','auto')
    aL1 = axes('units','normalized','Position',[0.1,0.50,0.8,0.45],'Tag','axesL1');	
    aL2 = axes('units','normalized','Position',[0.1,0.08,0.8,0.33],'Tag','axesL2','YDir','reverse');	

    % Set xtick label
    xt = sort([datenum(2015,1:1:12,1),datenum(2016,1:1:12,1)]);
    % Evapotranspiration comparison
    bar(aL1,evap.time_day,evap.data_day,1,'k','EdgeColor','none');hold(aL1,'on');
    bar(aL1,lysi.time_day,lysi.data_day,1,'b','edgecolor','none')
    bar(aL1,time_out,et_mean(:,s),1,'facecolor','none','EdgeColor',[1,0,0]);
    mm_setaxes(aL1,'xlim',[datenum(start),datenum(stop)],'xticklabel',[],'xtick',xt,...
        'ylabel','mm/day','legend',{'reference','lysi','estim'},'fontsize',font_size,...
        'title',sprintf('Reference vs Lysimeter vs Estimation (window = %2d)',wind_length(s)),...
        'ylim',[0 6]);

    % Errors
    er(:,s) = lysi.data_day-et_mean(:,s);
    [er_std(:,s),~,er_min(:,s),er_max(:,s),~] = mm_statnan(er(:,s));
    er_ref(:,s) = evap.data_day-et_mean(:,s);
    [er_ref_std(:,s),~,er_ref_min(:,s),er_ref_max(:,s),~] = mm_statnan(er_ref(:,s));
    bar(aL2,lysi.time_day,lysi.data_day,1,'k','EdgeColor','k');hold(aL2,'on');
    bar(aL2,lysi.time_day,abs(er(:,s)),1,'r','EdgeColor','r');
    mm_setaxes(aL2,'xlim',[datenum(start),datenum(stop)],'xtick',xt,'ylim',[0 6],...
        'ylabel','mm/day','legend',{'lysi','error'},'fontsize',font_size,'dateformat','dd.mm.yyyy','title',...
        sprintf('Error SD = %.2f, min = %.2f, max = %.2f mm',er_std(:,s),er_min(:,s),er_max(:,s)));
    % Add quantiles (use try just to prevent error related to missing
    % floating (GFZ) licence)
    try
        qi = [0.68 0.95 0.99];
        qo = quantile(abs(er(:,s)),qi); % the median of x
        col = {'r','g','b'};
        for i = 1:length(qi)
            plot([datenum(start),datenum(stop)]',[qo(i);qo(i)],col{i});
            text(datenum(start)+wind_length(s)+10,5.5-(i*2)/4,sprintf('Quantile of abs(error): %02d%%=%.2f mm',round(qi(i)*100),qo(i)),'color',col{i});
        end
    end
    
    %% Remove parameters needed for next loop run
    clear er qi qo col in_grav in_time param in_bias error_flag p pmin
end
rmpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
rmpath(fullfile('..','MatlabLibrary','aux_fce'));