clear;
close all;
clc;

%--------------------------------------------------------------------------
% This code computes the empirical and modeled mean hillslope lengths and
% plot them.
% Make sure to have run the code "p1_empiricalHillslopeLengthDistributions.m"
% (available at the same deposit of this code) for every catchment before
% running this code.
% Empirical mean will be stored in a file "smoothEmpMean.mat", contained in
% the same folder of the distributions that has been computed by the program
% "p1_empiricalHillslopeLengthDistributions.m".
%--------------------------------------------------------------------------


saveplots = false;

catchments_list_local =   ["borsoia"; "coweeta12"; "coweeta33"; "coweeta40"; "fernow14"; "fernow16"; "fernow37"; "focobon"; "montecalvello"; "poverty25"; "poverty35"; "turboloEst"; "turboloOvest"; "southForkPotts"; "rietholzback"];
catchments_maxDd_obse = [0.56; 6.07; 3.40; 5.10; 2.96; 2.42; 2.32; 2.24; 3.69; 5.11; 5.32; 7.78; 3.50; 2.00; 3.51];

maxnumToPlot = numel(catchments_list_local);

global smoothSignals;
global keepExtremes;
global mergeBins;
smoothSignals = true;
keepExtremes = true;
mergeBins = true;

labSz = 14;

fig1 = figure();
t=tiledlayout(3,5,'TileSpacing','tight');
hold on;

fig2 = figure();
t2=tiledlayout(3,5,'TileSpacing','tight');
hold on;

colors = [0 0.4470 0.7410; 0.8500 0.3250 0.0980; 0.9290 0.6940 0.1250; 0.4940 0.1840 0.5560; 0.4660 0.6740 0.1880];

leg = ["Empirical $\qquad$"; "$\langle\ell\rangle=\sqrt{L^2+\langle\ell_{_0}\rangle^2}-L$ $\qquad$"; "$\langle\ell\rangle=\sqrt{L^2+A}-L$ $\qquad$"; "$\langle\ell\rangle=1/(2D_d)$ $\qquad$"];
legflagdone = false;

letters = ["a"; "b"; "c"; "d"; "e"; "f"; "g"; "h"; "i"; "l"; "m"; "n"; "o"; "p"; "q"; "r"; "s"];

xboundtot = 3;

for hh_local=1:maxnumToPlot
    name_catchment = catchments_list_local(hh_local);
    disp("Processing catchment: "+upper(name_catchment));
    maxDd = catchments_maxDd_obse(hh_local);
    
    %Import and process DTM data
    orig_path = cd("Results_"+name_catchment);
    cd("Output_files");
    load("FinalResults.mat","network_Lengths","dx","totpixel","distributions_support","distributions_w0");
    xbound = maxDd*sqrt(totpixel)*dx/1000;
    smoothEmpMean = computeMeanOfPdf(distributions_w0,1,distributions_support(:,1)',distributions_support(:,1)');
    smoothEmpMeanKm = smoothEmpMean*dx/1000;
    init_cond = smoothEmpMean(1) + 1;
    init_cond_Km = init_cond*dx/1000;
    save("smoothEmpMean.mat","smoothEmpMean","smoothEmpMeanKm","init_cond","init_cond_Km");
    cd(orig_path);

    figure(fig1);
    hold on;
    nexttile;
    hold on;
    % Find out the proper indices at which to plot the markers
    inds = [1];
    numMarks = 15;
    for klkh=1:numel(network_Lengths)
        for mm1=1:numMarks
            if ((network_Lengths(klkh)>=3/numMarks*mm1*sqrt(totpixel)-sqrt(2)) && (network_Lengths(klkh)<3/numMarks*mm1*sqrt(totpixel)))
                inds = [inds; klkh];
            end
        end
    end
    % Linearized mean with empirical init. cond.
    konst = init_cond^2;
    handLinMean = - network_Lengths + sqrt((network_Lengths.^2)+konst*ones(size(network_Lengths)));
    handLinMean = handLinMean*dx/1000;
    p2 = plot(network_Lengths/sqrt(totpixel),handLinMean/(init_cond*dx/1000),'Color',colors(2,:),'LineWidth',2.2,'LineStyle',':');
    % Linearized mean with approx. (sqrt area) init. cond.
    konst = totpixel;
    handLinMeanAppr = - network_Lengths + sqrt((network_Lengths.^2)+konst*ones(size(network_Lengths)));
    handLinMeanAppr = handLinMeanAppr*dx/1000;
    p3 = plot(network_Lengths/sqrt(totpixel),handLinMeanAppr/(init_cond*dx/1000),'Color',colors(3,:),'LineWidth',2.2);
    % Horton
    p4 = plot(network_Lengths/sqrt(totpixel),0.5*totpixel*dx/1000./network_Lengths/(init_cond*dx/1000),'Color',colors(4,:),'LineStyle','--');
    % Empirical mean
    p1 = plot(network_Lengths/sqrt(totpixel),smoothEmpMeanKm/(init_cond*dx/1000),'Color',colors(1,:),'LineStyle',"none",'Marker','o','MarkerSize',4,'MarkerIndices',inds,'MarkerFaceColor',[1 1 1]);
    % Feasible region (observed range)
    ybound = 1.5;
    if (xbound<xboundtot)
        X = [xbound, xboundtot, xboundtot, xbound]';
        Y = [0, 0, ybound, ybound]';
        f1 = fill(X,Y,'k','FaceAlpha',0.05,'EdgeAlpha',0.05);
        if (~legflagdone)
            % Legend
            figure(fig1);
            hold on;
            bblg = legend([p1,p2,p3,p4],leg,'Orientation','horizontal','Interpreter','latex','Box','off','FontSize',14);
            bblg.Layout.Tile = 'north';
        end
    end
    % Finishing
    xlim([0, xboundtot]);
    ylim([0, ybound]);
    text(0.9,0.9,letters(hh_local),'HorizontalAlignment','center','Units','normalized','FontSize',13);
    if (hh_local==6)
        ylabel("$\langle\ell\rangle /\langle\ell_{_0}\rangle\quad [-]$","Interpreter","latex",'FontSize',labSz);
    elseif (hh_local==13)
        xlabel("$L/\sqrt{A}\quad [-]$","Interpreter","latex",'FontSize',labSz);
    end

    figure(fig2);
    hold on;
    nexttile;
    hold on;
    % COMPUTING THE FUNCTIONS OF PERCENTAGE DERIVATIVES:
    % Empirical percentage derivative
    derivEmp = computeDerivative(network_Lengths*dx/1000,smoothEmpMeanKm);
    percEmp = -derivEmp.*network_Lengths*dx/1000./smoothEmpMeanKm;
    [modifiedX, percEmp] = movingAverage(percEmp,network_Lengths*dx/1000);
    % Linearized percentage derivative with empirical init. cond.
    derivHandLin = - handLinMean ./ (network_Lengths*dx/1000 + handLinMean);
    percHandLin = -derivHandLin.*network_Lengths*dx/1000./handLinMean;
    % Linearized percentage derivative with approx. init. cond.
    derivHandLinAppr = - handLinMeanAppr ./ (network_Lengths*dx/1000 + handLinMeanAppr);
    percHandLinAppr = -derivHandLinAppr.*network_Lengths*dx/1000./handLinMeanAppr;
    % PLOTTING:
    p6 = plot(network_Lengths/sqrt(totpixel),percHandLin,'Color',colors(2,:),'LineWidth',2.2,'LineStyle',':');
    p7 = plot(network_Lengths/sqrt(totpixel),percHandLinAppr,'Color',colors(3,:),'LineWidth',2.2);
    p8 = plot(network_Lengths/sqrt(totpixel),ones(size(network_Lengths)),'Color',colors(4,:),'LineStyle','--');
    % Empirical percentage derivative
    p5 = plot(network_Lengths/sqrt(totpixel),percEmp,'Color',colors(1,:),'LineStyle',"none",'Marker','o','MarkerSize',4,'MarkerIndices',inds,'MarkerFaceColor',[1 1 1]);
    % Feasible region (observed range)
    ybound = 1.8;
    if (xbound<xboundtot)
        X = [xbound, xboundtot, xboundtot, xbound]';
        Y = [0, 0, ybound, ybound]';
        f2 = fill(X,Y,'k','FaceAlpha',0.05,'EdgeAlpha',0.05);
        if (~legflagdone)
            % Legend
            figure(fig2);
            hold on;
            bblg = legend([p5,p6,p7,p8],leg,'Orientation','horizontal','Interpreter','latex','Box','off','FontSize',14);
            bblg.Layout.Tile = 'north';
            legflagdone = true;
        end
    end
    % Finishing
    xlim([0, xboundtot]);
    ylim([0, ybound]);
    text(0.9,0.9,letters(hh_local),'HorizontalAlignment','center','Units','normalized','FontSize',13);
    if (hh_local==6)
        ylabel("$-\frac{L}{\langle\ell\rangle}\frac{d\langle\ell\rangle}{dL}\quad [-]$","Interpreter","latex",'FontSize',labSz+3);
    elseif (hh_local==13)
        xlabel("$L/\sqrt{A}\quad [-]$","Interpreter","latex",'FontSize',labSz);
    end

end


if (saveplots)
    exportgraphics(fig1,"Means.pdf","ContentType","vector","Append",false);
    exportgraphics(fig2,"PercentDeriv.pdf","ContentType","vector","Append",false);
end


disp("p3_computeMeansAndPlot.m COMPLETED his TASK.");






%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% AUXILIARY FUNCTIONS:


function [modifiedX, smoothed] = movingAverage(vect,vectx)

    global smoothSignals;
    global keepExtremes;
    global mergeBins;


    modifiedX = vectx;
    if (smoothSignals)

        if (mergeBins)
            smoothed = nan(size(vect));
            n = numel(vect);
            width = floor(n/40);
            for jj=1:n
                local_width = min([jj-1,width,n-jj]);
                extract = vect(jj-local_width:jj+local_width);
                mask = isnan(extract);
                extract(mask) = 0;
                tmp = sum(extract,"all");
                smoothed(jj) = tmp/(numel(extract)-sum(mask,"all"));
            end
            if (~keepExtremes)
                smoothed(1:max(1,floor(width/5))) = nan;
            end
        else
            n = numel(vect);
            numBins = 16;
            width = n/numBins;
            modifiedX = zeros(numBins,1);
            smoothed = zeros(numBins,1);
            for jj=1:numBins
                local_pos = floor((jj-1)*width) + 1 : floor(jj*width);
                extract = vect(local_pos);
                mask = isnan(extract);
                extract(mask) = 0;
                tmp = sum(extract,"all");
                smoothed(jj) = sum(tmp,"all")/numel(local_pos);
                modifiedX(jj) = vectx(floor((jj-0.5)*width));
            end
        end

    else
        smoothed = vect;
    end

end

function deriv = computeDerivative(xvals,vect)

    num = size(vect);
    deriv = nan(num);
    num = max(num,[],"all");
    deriv(1:num-1) = (vect(2:num)-vect(1:num-1))./(xvals(2:num)-xvals(1:num-1));

end

function mean = computeMeanOfPdf(Sol,dir,vals,funVals)

    if ((isnumeric(dir)) && ((dir==1) || (dir==2)))
        mask = isnan(Sol);
        Sol(mask) = 0;
        if (dir==1)
            if (size(vals,2)==size(Sol,1))&&(size(vals,1)==1)
                steps = vals(2:numel(vals)) - vals(1:numel(vals)-1);
                steps = [steps, steps(numel(steps))];
                tmp = funVals.*steps;
                mean = tmp*Sol;
            else
                disp("   ERROR: incompatible inputs:");
            end
        else
            if (size(vals,1)==size(Sol,2))&&(size(vals,2)==1)
                steps = vals(2:numel(vals)) - vals(1:numel(vals)-1);
                steps = [steps; steps(numel(steps))];
                tmp = funVals.*steps;
                mean = Sol*tmp;
            else
                disp("   ERROR: incompatible inputs:");
            end
        end
    elseif ((isstring(dir) && strcmpi(dir,"all")) || (ischar(dir) && sum(dir=='all',"all")==3))
        disp("ERROR at computeMeanOfPdf: still don't know how to deal with multivariate pdfs.");
        mean = nan(2,1);
    else
        disp("ERROR at computeMeanOfPdf: invalid input dir.");
    end

end