clear;close all;clc;

addpath('/Users/xudo627/donghui/CODE/Setup-E3SM-Mac/matlab-scripts-for-mosart/');
addpath('/Users/xudo627/donghui/mylib/m/');
datadir = '/Volumes/LaCie/PROJECTS/Pantanal/';

xv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xv');
yv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yv');
xc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xc');
yc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yc');

forcing = 'CRUNCEPv7';%'GSWP3v1';%'CRUNCEPv7';
%exs     = {'default','geo','man2','sub','sur','hydraulic','run3','hydraulic2'};
%labels  = {'Observation','Default','Inundation','Geometry','Manning','Subsurface','Surface','Hydraulic','Hydrology','Overall'};
exs     = {'default','hydraulic','run5','hydraulic2'};
labels  = {'Observation','Default','Hydraulic experiment','Hydrology experiment','Two-steps calibration'};
labels2 = {'Default','Hydraulic','Hydrology','Calibration'};
[lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);

cols = getPyPlot_cMap('Set1',9);
yr1 = 1979;
yr2 = 2009;
[tmon] = get_time_number(yr1,yr2);
max_mon = NaN(length(labels),yr2-yr1+1);

figure; set(gcf,'Position',[10 10 1400 800]);

ax(1) = axes('Position',[0.050 0.525 0.425 0.40]);
gsim_yr = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
[~,max_mon(1,:)] = max(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]));
plot(gsim_yr,'k--','LineWidth',3); hold on; grid on;
set(gca,'xtick',1:12, ...
        'xticklabel',{'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'}, ...
        'xlim',[1 12],'FontSize',13);
ylabel('Discharge [m^{3}/s]','FontSize',16,'FontWeight','bold');
k = 1;

metrics = NaN(4,length(exs));
for i = 1 : length(exs)
    tag = ['outputs_' exs{i} '_CLM' forcing];
    if strcmp(exs{i},'default')
        ntot = 10;
    elseif strcmp(exs{i},'geo') || strcmp(exs{i},'sub')
        ntot = 100;
    end
    load(['data/' tag '.mat'],'q','R2','KGE','NSE','ind2','qdraiall','qrunoffall');
    
    if strcmp(exs{i},'hydraulic2')
        qcal = q;
        ical = ind2;
    end
%     if strcmp(exs{i},'geo')
%         rng(1);
%         rwid = 2.6 + (20.2-2.6)*rand(ntot,1);   % rwid [2.6 - 20.2];
%         rdep = 0.12 + (0.63-0.12)*rand(ntot,1); % rdep [0.12 0.63];
%         dist = sqrt((rwid - 7.2).^2 + (rdep - 0.26).^2);
%         ind3 = find(dist == min(dist));
%         q_yr0= nanmean(reshape(q(ind3,1:end-12)',[12,yr2-yr1+1]),2);
%         [~,max_mon(k+1,:)] = max(reshape(q(ind3,1:end-12)',[12,yr2-yr1+1]));
%         plot(q_yr0,'-','Color',cols(k,:),'LineWidth',3); 
%         [~,~,metrics(1,k+1)] = LSE(q(ind3,1:end-12)',mu_gsim(yr >= yr1 & yr <= yr2));
%         metrics(1,k+1) = sqrt(metrics(1,k+1));
%         metrics(2,k+1) = nanmean(q(ind3,1:end-12))/nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
%         metrics(3,k+1) = nanstd(q(ind3,1:end-12))/nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
%         [~,~,metrics(4,k+1)] =  estimate_evaluation_metric(mu_gsim(yr >= yr1 & yr <= yr2),q(ind3,1:end-12)');
%         k = k + 1;
%     end
    if strcmp(exs{i},'default')
        q_yr = nanmean(reshape(q(1,1:end-12)',[12,yr2-yr1+1]),2);
        [~,max_mon(k+1,:)] = max(reshape(q(1,1:end-12)',[12,yr2-yr1+1]));
        [~,~,metrics(1,k+1)] = LSE(q(1,1:end-12)',mu_gsim(yr >= yr1 & yr <= yr2));
        metrics(1,k+1) = sqrt(metrics(1,k+1));
        metrics(2,k+1) = nanmean(q(1,1:end-12))/nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
        metrics(3,k+1) = nanstd(q(1,1:end-12))/nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
        metrics(4,k+1) =  nanmean(nanmean(qdraiall((1-1)*191+1:1*191,1:end-12),2)./nanmean(qrunoffall((1-1)*191+1:1*191,1:end-12),2));
        %estimate_evaluation_metric(mu_gsim(yr >= yr1 & yr <= yr2),q(1,1:end-12)');
    else
        q_yr = nanmean(reshape(q(ind2,1:end-12)',[12,yr2-yr1+1]),2);
        [~,max_mon(k+1,:)] = max(reshape(q(ind2,1:end-12)',[12,yr2-yr1+1]));
        [~,~,metrics(1,k+1)] = LSE(q(ind2,1:end-12)',mu_gsim(yr >= yr1 & yr <= yr2));
        metrics(1,k+1) = sqrt(metrics(1,k+1));
        metrics(2,k+1) = nanmean(q(ind2,1:end-12))/nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
        metrics(3,k+1) = nanstd(q(ind2,1:end-12))/nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
        metrics(4,k+1) =  nanmean(nanmean(qdraiall((ind2-1)*191+1:ind2*191,1:end-12),2)./nanmean(qrunoffall((ind2-1)*191+1:ind2*191,1:end-12),2));
        %estimate_evaluation_metric(mu_gsim(yr >= yr1 & yr <= yr2),q(ind2,1:end-12)');
    end
    
    plot(q_yr,'-','Color',cols(k,:),'LineWidth',3); 
    k = k + 1;
end
leg1 = legend(labels,'FontSize',15,'FontWeight','bold');

%metrics = metrics(1:end-1,:);

ax(2) = axes('Position',[0.55 0.525 0.425 0.40]);
X = categorical({'\rho','AR','rSD','BFI'});
X = reordercats(X,{'\rho','AR','rSD','BFI'});
b = bar(X,metrics(:,2:end),0.9); grid on;
set(gca,'FontSize',13);
ylabel('Metrics [-]','FontSize',16,'FontWeight','bold');
ax(2).XAxis.FontSize = 15;ax(2).XAxis.FontWeight = 'bold';
for i = 1 : length(labels2)
    b(i).FaceColor = cols(i,:);
end
% Add values on the bar
dataSample = metrics(:,2:end);
dataSample = round(dataSample,2);
nModel = size(metrics(:,2:end),2);
nCat = size(metrics(:,2:end),1);
xPosAmpl = 0.28;%0.3682626-0.3298725*exp(-0.407004*(nModel-1)); % position amplitude
xPosInc = 2*xPosAmpl/(nModel-1);
modelNames = [];
for idxModel=1:nModel
    bar_xPos = 1:nCat;
    if nModel~=1
        bar_xPos = bar_xPos-xPosAmpl+(idxModel-1)*xPosInc;
    end
    y = dataSample(:,idxModel);
    y(y < 0) = 0;
    text(bar_xPos,y,num2str(dataSample(:,idxModel)), ...
         'vert','bottom','horiz','center','FontSize',13,'fontname','times'); 
end
%legend(modelNames)

ax(3) = axes('Position',[0.050 0.075 0.925 0.375]);
p(1) = plot(tmon,mu_gsim(yr >= yr1 & yr <= yr2),'k--','LineWidth',3); hold on; grid on;
p(2) = plot(tmon,qcal(ical,1:end-12),'-','Color',cols(4,:),'LineWidth',3);
xlim([tmon(1) tmon(end)]);
datetick('x','mmmyy','keeplimits');
set(gca,'FontSize',13);
ylabel('Discharge [m^{3}/s]','FontSize',16,'FontWeight','bold');
leg = legend(p,{'Observation','Two-steps calibration'});
leg.FontSize = 18; leg.FontWeight = 'bold';
ax(3).XAxis.TickValues = tmon(1:12:end);
ax(3).XAxis.TickLabelRotation = 45;

% grid on;
% set(gca,'xtick',1:12,'xticklabel',{'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
% set(gca,'FontSize',13); xlim([0.5 12.5]);
% ylabel('Counts','FontSize',15,'FontWeight','bold');

add_title(ax(1),'(a). Streamflow seasonality',18);
add_title(ax(2),'(b). Evaluation metrics',18);
add_title(ax(3),'(c). Streamflow monthly time series',18);

figure;
plot(nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1])),'k--','LineWidth',3); hold on;
plot(nanmean(reshape(qcal(ical,1:end-12)',[12,yr2-yr1+1])),'-','Color',cols(4,:),'LineWidth',3); hold on;