clear; close all; clc; % ------------------------------------------------------------------------- % CODE DESCRIPTION: % This program automatically perform the calibration of the hyperbolic % parameter "a" for all of the catchments. % Please, make sure that the program "p1_empiricalHillslopeLengthDistributions.m" % has been runned for all of the catchments before starting the calibration. % The user is not requested to specify/set the value of any % variable of this code before running it. % Results of the calibration will be stored just in one location for all the catchment grouped together, % in a folder named "CalibrationOfHyperbolicParam". % Results will be stored in a file named "aOptimNoDiff.mat"; in particular, % the last column of the variable "a_optimals" contained in that file will % corresponds to the calibrated values of "a" (one row correspond to a % certain catchment) over the range of the channel length L observed in the field (or an approximation of it, where % there is no good information about that range). % ------------------------------------------------------------------------- catchments_list = ["astico"; "borsoia"; "coweeta12"; "coweeta33"; "coweeta40"; "fernow14"; "fernow16"; "fernow37"; "focobon"; "leogra"; "montecalvello"; "poverty25"; "poverty35"; "turboloEst"; "turboloOvest"; "southForkPotts"; "rietholzback"; "fuciade"; "seugne"; "hubbardBrook"; "valfredda"; "cordevole"]; thresh_lower_lims = [200; 1000; 1000; 1000; 1000; 50; 50; 200; 250; 1000; 1000; 5000; 4000; 300; 300; 500; 10000; 21500000; 3463; 31000000; 500000; 1000000]; thresh_upper_lims = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 22000000; 0; 31341400; 501400; 2000000]; %DTM resolutions in m dxDTMs = [50; 5; 1; 1; 1; 3; 3; 3; 20; 25; 5; 1; 0.9; 5; 5; 3; 2; 0.5; 26.87; 1; 1; 10]; %DTM area in m^2 AreasDTMs = 1000*[136470; 13790; 120; 310; 410; 140; 160; 370; 8420; 67410; 3750; 250; 340; 630; 460; 730; 3120; 5600; 911100; 31340; 5270; 832290]; %DTM area in pixels AreasDTMsPx = [54587; 551461; 124989; 306933; 407453; 15300; 17612; 41109; 21043; 107850; 149888; 249506; 426424; 25343; 18367; 81102; 780311; 22412400; 1262280; 31341400; 5268450; 8322850]; %DTM area in m^2 (accurate) AreasDTMsM2Accurate = AreasDTMsPx.*dxDTMs.*dxDTMs; aInterval = [0.4; 1.8]; num_sub_intervals = 3; num_drainDensity_vals = 3; limitDrainVal_local = nan(size(catchments_list)); limitDrainVal_local(6) = 3; limitDrainVal_local(7) = 3; limitDrainVal_local(8) = 3; limitDrainVal_local(16) = 2; limitDrainVal_local(17) = 4; global anotherLimToDd; saveplots = true; appendToUniqueFile = true; uniqueFileName = "lookingForOptimalHyperbolicsParam.pdf"; global savefolder; savefolder = "CalibrationOfHyperbolicParam"; global readMatFiles; %Set this to true if you have already run this code once and want to skip %the computation of the preliminary variables readMatFiles = false; global areaPx; global dx; a_optimals = nan(size(dxDTMs,1),num_drainDensity_vals); D_densities = nan(size(dxDTMs,1),num_drainDensity_vals); maxhill = nan(size(dxDTMs)); Lmax = nan(size(dxDTMs)); for hh=1:numel(catchments_list)-5 name_catchment = catchments_list(hh); areaPx = AreasDTMsPx(hh); dx = dxDTMs(hh); anotherLimToDd = limitDrainVal_local(hh); disp("Operating on catchment: "+upper(name_catchment)); disp("Looking for the optimal value of the hyperbolic parameter"); if (hh==1) [a_optimals(hh,:), D_densities(hh,:), maxhill(hh), Lmax(hh)] = cmpForSingleBasin(name_catchment,false,uniqueFileName,aInterval,num_sub_intervals,saveplots,num_drainDensity_vals); else [a_optimals(hh,:), D_densities(hh,:), maxhill(hh), Lmax(hh)] = cmpForSingleBasin(name_catchment,appendToUniqueFile,uniqueFileName,aInterval,num_sub_intervals,saveplots,num_drainDensity_vals); end close all; clc; end disp("Saving the computed values..."); save("aOptimNoDiff.mat","catchments_list","AreasDTMs","AreasDTMsPx","AreasDTMsM2Accurate","dxDTMs","Lmax","maxhill","a_optimals","D_densities"); disp(" "); disp("AutoFindOptimalHyperbolics COMPLETED his TASK."); function [a_opt, ddensi, maxtau, Tmax] = cmpForSingleBasin(name_catchment,appendToUniqueFile,uniqueFileName,aBounds,num_sub_int,saveplots,num_drainDensity_vals) %-------------------------------------------------------------------------- %INPUT SPECIFICATIONS (TO BE SET BY USER) saveformat = "png"; approx_color = [[0 0.4470 0.7410]; [0.8500 0.3250 0.0980]; [0.9290 0.6940 0.1250]; [0.4940 0.1840 0.5560]; [0.4940 0.1840 0.5560]; [0.4660 0.6740 0.1880]; [0.4660 0.6740 0.1880]; [0.3010 0.7450 0.9330]; [0.3010 0.7450 0.9330]; [0.4940 0.1840 0.5560]; [0.4660 0.6740 0.1880]; [0.3010 0.7450 0.9330]; [0.6350 0.0780 0.1840]; [1 0 0]; [0 0 0]]; %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- %PRELIMINARY AND IMPORT OPERATIONS: %SETTING WORKING DIRECTORY TO READ OUTPUT FILES restore_path = cd(); cd("Results_"+name_catchment); cd("Output_files"); %INITIALIZATION OF VARIABLES load("Workspace_Final.mat","network_Lengths","distributions_support","distributions_w0"); tau = distributions_support(:,1); maxtau = max(tau); T = network_Lengths; % support = distributions_support; % alpha = distributions_alpha; % w0 = 2*distributions_w0; w0 = distributions_w0; % w = distributions_w; load("totpixel.mat"); %RESTORING INITIAL WORKING DIRECTORY cd(restore_path); %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- %RETRIEVING INFORMATION ABOUT BOUNDARY CONDITION (AT T=0) %Readed minimal network's data is T=1 w0_bc_loc = w0(:,1); global hillStep; hillStep = tau(1); %Rearranging the boundaries to hold at T=0 (instead of T=1) Tmax = max(T); paddInit = (hillStep:hillStep:1)'; numPadd = numel(paddInit); tau_bc_loc = [paddInit; tau+ones(size(tau))]; w0_bc_loc = [zeros(numPadd-1); 1/totpixel; w0_bc_loc*(totpixel-1)/totpixel]; %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- %CORE OPERATIONS: %COMPUTING SOLUTIONS FOR MEAN ALPHA COEFFICIENTS maskExp = isnan(w0); w0(maskExp) = 0; global w0_glob; w0_glob = w0; global leg; leg = []; global tau_bc; tau_bc = tau_bc_loc; global w0_bc; w0_bc = w0_bc_loc; global Tglob; Tglob = T; global stepsForErrors; stepsForErrors = nan(size(T)); stepsForErrors(1) = 1; stepsForErrors(2:numel(T)-1) = abs(T(3:numel(T))-T(1:numel(T)-2))./2; global numValidSteps; numValidSteps = numel(T) - sum(isnan(stepsForErrors),'all'); global areaPx; global dx; global anotherLimToDd; maxDd = min(anotherLimToDd,5); L_corresp_to_maxDd = maxDd*areaPx*dx/1000; for bl=1:numel(T) if (T(bl)>L_corresp_to_maxDd) break; end end bl=bl-1; limit = min(bl,numValidSteps); global local_lim; a_opt = zeros(1,num_drainDensity_vals); ddensi = zeros(1,num_drainDensity_vals); for jhj=1:num_drainDensity_vals disp("ROUND "+jhj); disp(""); local_lim = floor(jhj*limit/(num_drainDensity_vals)); % global sumValidSteps; % sumValidSteps = sum(stepsForErrors(1:numValidSteps),'all'); % global reducedStepsForErrors; % reducedStepsForErrors = stepsForErrors(1:numValidSteps); global sumValidSteps; sumValidSteps = sum(stepsForErrors(1:local_lim),'all'); global reducedStepsForErrors; reducedStepsForErrors = stepsForErrors(1:local_lim); figPareto=figure; hold on; %Determining optimal a (within the given interval) tau_prov = zeros(num_sub_int,1); fit_prov = zeros(num_sub_int,1); step = (aBounds(2)-aBounds(1))/num_sub_int; lbstart = aBounds(1); for ii=1:num_sub_int lb = lbstart+(ii-1)*step; ub = lb +step; disp(" Attempting optimization on subinterval: "+ii+" out of "+num_sub_int); [tau_prov(ii), fit_prov(ii)] = fminbnd(@(a) evalAChoice(a),lb,ub); end fit_min = min(fit_prov); mask = fit_min*ones(size(fit_prov)) == fit_prov; a_opt_tmp = tau_prov(mask); a_opt(jhj) = a_opt_tmp(1); scatter(a_opt(jhj), fit_min, 36, 'r', "x","LineWidth",2); legend("a_{opt}="+a_opt(jhj),'Location','best'); %legend(leg,'Location','bestoutside'); %xlabel("Overall error w.r.t empirical [-]"); %ylabel("Preservation of unitary integral [-]"); % xlabel("\epsilon_1 = 100 (2T_{max})^{ -1} \Sigma_{T,\tau} |\omega_0- \omega_{0,dtm}|"); % ylabel("\epsilon_2 = T_{max}^{ -1} \Sigma_{T} |1- \Sigma_{\tau} \omega_{0} |"); xlabel("a"); ylabel("mean error"); title(" "+name_catchment+": mono-objective optimization"); consMaxDd = T(local_lim)*1000/(areaPx*dx); ddensi(jhj) = consMaxDd; subtitle("Considered max drain density: "+consMaxDd+" [km^{-1}]"); if (saveplots) global savefolder; mkdir(savefolder); cd(savefolder); % mkdir("images"); % cd("images"); % saveas(figPareto,name_catchment+"_objectives"+jhj+"."+saveformat,saveformat); % cd ..; if ((appendToUniqueFile)||(jhj>1)) exportgraphics(figPareto,uniqueFileName,"ContentType","vector","Append",true); else exportgraphics(figPareto,uniqueFileName,"ContentType","vector","Append",false); end cd ..; end end %-------------------------------------------------------------------------- end %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- %AUXILIARY FUNCTIONS function fit = evalAChoice(a) % global leg global w0_glob global tau_bc global w0_bc global Tglob global numValidSteps w0_newExp1 = nan(size(w0_glob)); nrows=size(w0_newExp1,1); fine = min(numValidSteps,size(w0_newExp1,2)); for kk=1:fine tmp = getIdealAlphaAreaSol(Tglob(kk),a,tau_bc,w0_bc); w0_newExp1(:,kk) = tmp(1:nrows); end diff = abs(w0_newExp1 - w0_glob); [err_newExp1, intTo1_newExp1, errTot_newExp1] = generalizedComputeErrorParameters(diff,w0_newExp1); % leg = [leg; "a="+a]; % scatter(errTot_newExp1, intTo1_newExp1, 36, [0 0.4470 0.7410], "o","LineWidth",1); scatter(a, errTot_newExp1, 36, [0 0.4470 0.7410], "o","LineWidth",1); fit = errTot_newExp1; end %Functions to compute quality parameters of analitical solutions function [err, intTo1, errTot] = generalizedComputeErrorParameters(diff,w0) global reducedStepsForErrors; global local_lim; global sumValidSteps; global hillStep; err = nan(size(diff,2),1); intTo1_tmp = nan(size(diff,2),1); for kk=1:local_lim mask=~isnan(diff(:,kk)); err(kk)=hillStep*sum(diff(mask,kk),"all")*100/2; %mask=~isnan(w0(:,kk)); intTo1_tmp(kk)=abs(1-hillStep*sum(w0(mask,kk),"all")); end % mask2 = ~isnan(stepsForErrors); % validSteps = stepsForErrors(mask2); % intTo1 = (validSteps*intTo1_tmp(mask2))/size(diff,2); % errTot = (validSteps*err(mask2))/size(diff,2); intTo1 = (reducedStepsForErrors*intTo1_tmp(1:local_lim))/sumValidSteps; errTot = (reducedStepsForErrors*err(1:local_lim))/sumValidSteps; 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