%% Prepare inputs for Optimization
% Use this script to:
%   1. generate/prepare data for optimization.
%   2. estimate min-max ranges (constraints) for optimization.
%
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));
clear
close all
clc

%% Settings: Output
% Set output file name and range of the output time series. The gravity
% time vector will be cutted to this range and used as reference for
% output.
output_file = fullfile('..','..','Data','Optimization','Optimization_Inputs.tsf');
time_start = datenum(2015,06,19);
time_end = datenum(2016,6,30);

% Set names of output channels (keep in accordance with *.out_column
% parameter)
channels = {'Gravity','Atmo','GHE','Discharge','Rain','ET0'};
units =    {'nm/s^2', 'hPa', 'nm/s^2','mm','mm','mm'};

% Initial comment (additianal rows will be added in 'Load data' section)
c = 1;comment(c) = {'All input time series linearly interpolated to Gravity time vector'};

%% Setting: Inputs
input_file.grav = fullfile('..','..','Data','Grav','iGrav006_SG030_residuals.tsf');
input_column.grav = 1;
output_column.grav = 1;
% Barometric pressure series
input_file.atmo = fullfile('..','..','Data','Grav','iGrav006_SG030_residuals.tsf');
input_column.atmo = 3;
output_column.atmo = 2;
% Global hydrological effect. Set all input time series/models that will be
% loaded and than used to compute average effect (+/- amplitude factor). It
% is not necessary to set the data column as fixed mGlobe output is
% expected
input_file.ghe = fullfile('..','..','Data','Grav','GHE');
input_column.ghe = {'CLM','MOS','NOAH025','VIC'};
output_column.ghe = 3; 
% Rain
input_file.rain = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
input_column.rain = 1;
output_column.rain = 5;
% Input evapotranspiration time series. Use hourly ET0. Cumulative
% sum will be computed here
input_file.et0 = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
input_column.et0 = 4;
output_column.et0 = 6;

% Input Streamflow time series. Use discharge not baseflow (will be
% determined within optimization). 
input_file.disch = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
input_column.disch = 5; 
output_column.disch = 4;

%% Load data
in = {'grav','atmo','ghe','rain','et0','disch'};
for i = 1:length(in)
    if i == 1
        [temp_time,temp_data] = loadtsf(input_file.(in{i}));
        time = temp_time;
        data = temp_data(:,input_column.(in{i}));
        data(time<time_start,:) = [];
        time(time<time_start,:) = [];
        data(time>time_end,:) = [];
        time(time>time_end,:) = [];
    elseif strcmp(in{i},'ghe')
        for j = 1:length(input_column.ghe)
            [temp_time,temp_data] = loadtsf(fullfile(input_file.ghe,['iGrav_',char(input_column.ghe(j)),'_24h_gre1_ant1_dem1_thr010.tsf']));
            % Interpolate and subtract mean value so all time series vary around 
            % zero and can be used for regression analysis
            ghe.data(:,j) = demean(interp1(temp_time,temp_data(:,1),time)); % mGlobe output allways in first column
        end
        data(:,output_column.(in{i})) = mean(ghe.data,2);
        % Min-max range
        for j = 1:length(input_column.ghe)
            ghe.reg_par(j,:) = regress(ghe.data(:,j),data(:,output_column.(in{i})));
            ghe.cor(j,1) = mmcorr(ghe.data(:,j),data(:,output_column.(in{i})));
            comment(c) = {sprintf('GHE input models used for mean series: %s. Corr. coefficient = %.3f',input_column.ghe{j},ghe.cor(j,1))};c = c + 1;
        end
        temp = sprintf('GHE: min reg. parameter = %.3f',ghe.reg_par(ghe.reg_par==min(ghe.reg_par)));
        disp(temp);comment(c) = {temp};c = c + 1;
        temp = sprintf('GHE: max reg. parameter = %.3f',ghe.reg_par(ghe.reg_par==max(ghe.reg_par)));
        disp(temp);comment(c) = {temp};c = c + 1;
    elseif strcmp(in{i},'disch') || strcmp(in{i},'et0') || strcmp(in{i},'rain')
        [temp_time,temp_data] = loadtsf(input_file.(in{i}));
        temp_time = temp_time - 0.5/24;
        % Cut the interval first, then compute cumulative sum (so it starts at
        % zero) and after these steps interpolate to gravity/reference time vector
        temp_data(temp_time<time_start,:) = [];
        temp_time(temp_time<time_start,:) = [];
        temp_data = cumsum(temp_data(:,input_column.(in{i})));
        data(:,output_column.(in{i})) = interp1(temp_time,temp_data,time);
    else
        [temp_time,temp_data] = loadtsf(input_file.(in{i}));
        data(:,output_column.(in{i})) = temp_data(time >= time_start & time <= time_end,input_column.(in{i}));
    end
    comment(c) = {sprintf('%s: %s, column %d',in{i},input_file.(in{i}),output_column.(in{i}))};c = c + 1;
end

%% Subtract first value and Show outputs
figure('Units','normalized','Position',[0.1 0.1 0.6 0.8],'PaperPositionMode','auto','Name','Output time series');
for i = 1:length(channels)
    subplot(5,2,i)
    plot(time,data(:,i));
    mm_setaxes(gca,'title',channels{i},'ylabel',units{i},'xlim',[time(1),time(end)],...
        'xtick',sort([datenum(2015,1:3:12,1),datenum(2016,1:3:12,1)]),'dateformat','dd.mm.yyyy');
end

%% Write output
% Prepare header: channels and units
for i = 1:length(channels)
    header(i,1:4) = {'Wettzell','seeComment',channels{i},units{i}};
end
writetsf([datevec(time),data],header,output_file,3,comment);

rmpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
rmpath(fullfile('..','MatlabLibrary','aux_fce'));