clear;
close all;
clc;

catchments_list_local =   ["leogra"; "astico"; "seugne"; "umpqua"];

tableErr1 = nan(numel(catchments_list_local),3);
tableErr2 = nan(numel(catchments_list_local),3);

for hh_local=1:numel(catchments_list_local)
    name_catchment = catchments_list_local(hh_local);
    disp("Processing catchment: "+upper(name_catchment));
    
    %Import and process DTM data
    load("results_"+name_catchment+".mat","network_Lengths","areaKm2","hillslope_Means");
    emp = hillslope_Means;
    init_cond = emp(1);

    % Computing analitical approximations of the mean
    % Linearized with empirical init. cond.
    konst = init_cond^2;
    handLinMean = - network_Lengths + sqrt((network_Lengths.^2)+konst*ones(size(network_Lengths)));
    % Linearized with approx. (sqrt area) init. cond.
    konst = areaKm2;
    handLinMeanAppr = - network_Lengths + sqrt((network_Lengths.^2)+konst*ones(size(network_Lengths)));
    % Horton
    horton = 0.5*areaKm2./network_Lengths;

    % Specify the interval to be considered (max val for x=L/sqrt(A))
    maxval = 3.1;
    % Computing errors
    errHorton = computeErr(emp,horton,network_Lengths/sqrt(areaKm2),maxval);
    errModel = computeErr(emp,handLinMean,network_Lengths/sqrt(areaKm2),maxval);
    errModelAppr = computeErr(emp,handLinMeanAppr,network_Lengths/sqrt(areaKm2),maxval);
    tableErr2(hh_local,1) = 100*errHorton;
    tableErr2(hh_local,2) = 100*errModel;
    tableErr2(hh_local,3) = 100*errModelAppr;

end

disp("ERRORI SU L/sqrt{A} IN [0, min(3,maxObserved)]");
disp(round(tableErr1,0));
disp("ERRORI SU L/sqrt{A} IN [0, 3]");
disp(round(tableErr2,0));

catchAveraged1 = sum(tableErr1,1)/4;
catchAveraged2 = sum(tableErr2,1)/4;

disp("MEDIE INTERCATCHMENTS SU L/sqrt{A} IN [0, min(3,maxObserved)]");
disp(catchAveraged1);
disp("MEDIE INTERCATCHMENTS SU L/sqrt{A} IN [0, 3]");
disp(catchAveraged2);

disp("AutoComputeMeansPreliminary COMPLETED his TASK.");






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

function err = computeErrAss(emp,model,xvals,maxval,dx)

    err = 0;
    cou = 0;
    for jj=1:numel(emp)
        if (xvals(jj)<maxval)
            cou = cou + 1;
            err = err + (emp(jj)-model(jj))^2;
        else
            break;
        end
    end
    err = sqrt(err/cou)*dx;

end

function err = computeErr(emp,model,xvals,maxval)

    err = 0;
    cou = 0;
    for jj=1:numel(emp)
        if (xvals(jj)<maxval)
            cou = cou + 1;
            err = err + abs((emp(jj)-model(jj))/emp(jj));
        else
            break;
        end
    end
    err = err/cou;

end

% function err = computeErr(emp,model,xvals,maxval)
% 
%     err = 0;
%     cou = 0;
%     for jj=1:numel(emp)
%         if (xvals(jj)<maxval)
%             cou = cou + 1;
%             err = err + ((emp(jj)-model(jj))/emp(jj))^2;
%         else
%             break;
%         end
%     end
%     err = sqrt(err/cou);
% 
% end


% function mean = computeMeanOfPdfOld(Sol,dir,vals)
% 
%     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)
%                 mean = vals*Sol;
%             else
%                 disp("   ERROR: incompatible inputs:");
%             end
%         else
%             if (size(vals,1)==size(Sol,2))&&(size(vals,2)==1)
%                 mean = Sol*vals;
%             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

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

function [mean, var] = computeVarianceOfPdf(Sol,dir,vals)

    if ((isnumeric(dir)) && ((dir==1) || (dir==2)))
        mean = computeMeanOfPdf(Sol,dir,vals,vals);
        meanSq = computeMeanOfPdf(Sol,dir,vals,vals.^2);
        var = meanSq - (mean.^2);
        % funV = (vals-mean).^2;
        % var = computeMeanOfPdf(Sol,dir,vals,funV);
    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.");
        var = nan(2,1);
    else
        disp("ERROR at computeMeanOfPdf: invalid input dir.");
    end

end

function [mean, var, tmp, skew] = computeSkewnessOfPdf(Sol,dir,vals)

    if ((isnumeric(dir)) && ((dir==1) || (dir==2)))
        mean = computeMeanOfPdf(Sol,dir,vals,vals);
        meanSq = computeMeanOfPdf(Sol,dir,vals,vals.^2);
        meanTh = computeMeanOfPdf(Sol,dir,vals,vals.^3);
        var = meanSq - (mean.^2);
        tmp = meanTh + 2*(mean.^3) - 3*mean.*meanSq;
        skew = tmp ./ (var.^(3/2));
    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.");
        var = nan(2,1);
    else
        disp("ERROR at computeMeanOfPdf: invalid input dir.");
    end

end

function w = getIdealAlphaAreaSol(T,a,tau_bc,w0_bc)

    %Want to return the solution at T onto tau values tau_bc
    %Initialization
    tau = tau_bc;
    step = tau_bc(1);
    maxind = numel(tau);
    % Tvect = T*ones(size(tau));
    % %Actual computations (through definition)
    % linComb = tau + a*Tvect;
    % coeff1 = linComb./tau;
    % linComb = linComb + Tvect;
    % frac = tau./linComb;
    % frac = frac.^(1/(1+a));
    % frac1 = frac.^a;
    % scale = tau.*(frac.^(-1));
    % coeff = coeff1.*frac1;
    % w(1:maxind) = coeff .* bcFunction(scale, w0_bc);
    w = nan(size(tau));
    for kk=1:maxind
        l = tau_bc(kk);
        indIC = l * (1 + (1+a)*T/l)^(1/(1+a));
        indIC = floor(indIC/step);
        if ((indIC>0)&&(indIC<=maxind))
            w(kk) = w0_bc(indIC) * (a*T+l)/l*(l/(l+(1+a)*T))^(a/(1+a));
        else
            w(kk) = 0;
        end
    end
    % hold on;
    % plot(tau_bc,w0_bc);
    % plot(tau_bc,w);
    % hold off;

end