% SubdomainLcurveFigure_v1.m

% Lea-Sophie Hoeyns, October 2023 

% This script calculates for different subdomains in '.shp' file format
% the corresponding subdomain L-curve by calculating the cost function at 
% this positions. The output subdomain L-curve figures are created afterwards. 

%clearvars
close all
tic 

%% DATA PATHS 
issmscriptpath='/albedo/work/projects/landice/ISSM/ISSM_installations/ISSM_v4p22/trunk/';
functionpath='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Functions';
subdomainpath='~/1_Models/WAIS_FrictionOptimization_v2/Data/subdomains/WAIS/';
addpath(genpath(issmscriptpath),functionpath,subdomainpath)
% File names:
lcurvefiles={'/albedo/scratch/user/lhoeyns/1_ISSM_execution/Zwischenspeicher/Convergence_test/Run_6/WAIS.dragInv_HO_m1_Ncuas_1DTemp_Lcurve_05-Jan-2024.mat';...
        '/albedo/scratch/user/lhoeyns/1_ISSM_execution/Zwischenspeicher/OldLcurveRange/Run_1/WAIS.dragInv_HO_m3_Ncuas_1DTemp_Lcurve_17-Jan-2024.mat'};
% Load subdomains from exp or shp file 
files=dir(sprintf('~/1_Models/WAIS_FrictionOptimization_v2/Data/subdomains/WAIS/*.shp'));
subdomains={files(:).name};

%% PARAMETERS
% Model descriptions:
descriptext={'m=1, Budd, N_{CUAS}','m=3, Budd, N_{CUAS}'};  % cell array of strings

% Ticks and limits:
jreglims=[1e-9,1e-2];                 % [1x2] unitless
jobslims=[1e-6,1e-0];                 % [1x2] unitless

% Aesthetic settings
% Colors:
newblue='#5EA3D8';                    % valid color string
orange='#ED553B';                     % valid color string
lcurvecolors={newblue,orange};        % {1xN} valid color strings
totcurvecolor=gray;                   % valid color string
datacurvecolor=newblue;               % valid color string
regcurvecolor=orange;                 % valid color string
backplanecolor=[.75,.75,.75];         % valid color string or RGB

% Fontsizes:    
axesfontsize=25;                     % points
titlefontsize=25;                    % points
legendfontsize=42;                   % points
Fontsize=51;                         % points
LabelFontsize=19;                    % points
modelmarkersize=15;                  % points
curvelinesize=7;                     % points
cornerlinesize=21;                   % points
bestmarkersize=32;                   % points
% Linewidths:
LabelLinewidth=7;                    % points
Linewidth=5;                         % points
axeslinewidth=5;                     % points
% Styles:
modelmarker='o';                     % valid marker string
curvelinestyle='-';                  % valid linestyle string
cornerlinestyle='-';                 % valid linestyle string
bestmarker='d';                      % valid marker string

% Other parameter: 
numlambda=25;                        % integer
outlierthreshold=4;%2;               % integer

% Output resolution:
resolution=300;                      % dpi


%% PREPARATION 
numlcurves=length(lcurvefiles);

% Load basic model:
load(lcurvefiles{1},'BasicModel','MyInversionParameters','*best*','Lambda')

% Rename model structure:
md=BasicModel;

% Unpack inversion parameters:
unpack(MyInversionParameters)

% Compute mesh properties:
% Extract 2D mesh:
try
    x2d=md.mesh.x2d;
    y2d=md.mesh.y2d;
    elements2d=md.mesh.elements2d;
    oldnumlayers=md.mesh.numberoflayers;
catch
    x2d=md.mesh.x;
    y2d=md.mesh.y;
    elements2d=md.mesh.elements;
    oldnumlayers=1;
end
numelements2d=size(elements2d,1);
numvertices2d=length(x2d);
% Compute edge lengths:
a=sqrt((x2d(elements2d(:,1))-x2d(elements2d(:,2))).^2+(y2d(elements2d(:,1))-y2d(elements2d(:,2))).^2);
b=sqrt((x2d(elements2d(:,2))-x2d(elements2d(:,3))).^2+(y2d(elements2d(:,2))-y2d(elements2d(:,3))).^2);
c=sqrt((x2d(elements2d(:,3))-x2d(elements2d(:,1))).^2+(y2d(elements2d(:,3))-y2d(elements2d(:,1))).^2);
% Compute area:
s=(a+b+c)/2;
Area=sqrt(s.*(s-a).*(s-b).*(s-c));

% Compute gradient of linear basis functions:
BasisGradX_elements=[(y2d(elements2d(:,2))-y2d(elements2d(:,3)))./(x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))...
    +x2d(elements2d(:,3)).*y2d(elements2d(:,1))+x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))),...
    (y2d(elements2d(:,3))-y2d(elements2d(:,1)))./(x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))...
    +x2d(elements2d(:,1)).*y2d(elements2d(:,2))+x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))),...
    (y2d(elements2d(:,1))-y2d(elements2d(:,2)))./(x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))...
    +x2d(elements2d(:,2)).*y2d(elements2d(:,3))+x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1)))];
BasisGradY_elements=[-(x2d(elements2d(:,2))-x2d(elements2d(:,3)))./(x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))...
    +x2d(elements2d(:,3)).*y2d(elements2d(:,1))+x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))),...
    -(x2d(elements2d(:,3))-x2d(elements2d(:,1)))./(x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))...
    +x2d(elements2d(:,1)).*y2d(elements2d(:,2))+x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))),...
    -(x2d(elements2d(:,1))-x2d(elements2d(:,2)))./(x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))...
    +x2d(elements2d(:,2)).*y2d(elements2d(:,3))+x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1)))];

%% COLLECT DATA 
InversionParameters=cell(numlcurves,1);
FittingParameters=cell(numlcurves,1);
LCurveData=struct('LogLambda_points',[],'LogRegCost_points',[],'LogDataCost_points',[],'IsOutlier',[],'LogRegCost_curve',[],'LogDataCost_curve',[],...
    'LogDataCurvature_curve',[],'LogRegCurvature_curve',[],'TotalCurvature_curve',[],...
    'bestlambdaind_maxallowable',[],'bestlambdaind_minallowable',[],'bestlambdaind',[]);

%% CALCULATE L-curve 
for i=1:length(subdomains)
    subdomains(i)
    savename=['Lcurve_',subdomains{i},'.png'];
    % Determine elements or nodes inside the subdomains 
    in_node = ContourToMesh(elements2d,x2d,y2d,subdomains{i},'node',1);

    % Set Cosfunction (data and regu costs depending on in_nod vector)    
    % Pre-allocate cost values:
    DataCost=zeros(numlambda,1);
    RegularizationCost=zeros(numlambda,1);

    % Loop through models
    for thismodel=1:numlcurves
        % Load the basic file:
        load(lcurvefiles{thismodel},'Log*Cost_bestcurve','*Curvature_bestcurve','Lambda','DataCost',...
            'RegularizationCost','IsOutlier','bestlambda*','bestmodelind','MyInversionParameters','LCurveFittingParameters');%'numlambdasubsamples','lambdainterpmethod','maxlambdawavelength','numlambdawavelengths')
       
        % Unpack parameters:
        unpack(MyInversionParameters)
        unpack(LCurveFittingParameters)

        % Record parameters:
        InversionParameters{thismodel}=MyInversionParameters;
        FittingParameters{thismodel}=LCurveFittingParameters;

        % Loop through lambda values:
        for jj=1:numlambda
            % Communicate:
            disp('...')
            disp(['...lambda=',num2str(Lambda(jj))])
        
            % Construct this file name:
            slashnum=strfind(lcurvefiles{thismodel},'/');
            thisbase=10^floor(log10(Lambda(jj)));
            thisname=[lcurvefiles{thismodel}(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(jj)/thisbase)),'.mat'];
        
            % Load this model:
            load([lcurvefiles{thismodel}(1:slashnum(end)),thisname],'md')

            % Get model velocities:
            ModelVx=md.results.StressbalanceSolution.Vx(end-numvertices2d+1:end);
            ModelVy=md.results.StressbalanceSolution.Vy(end-numvertices2d+1:end);
            ModelVel=md.results.StressbalanceSolution.Vel(end-numvertices2d+1:end);
        
            % Get observed velocities:
            ObsVx=md.inversion.vx_obs(end-numvertices2d+1:end);
            ObsVy=md.inversion.vy_obs(end-numvertices2d+1:end);
            ObsVel=md.inversion.vel_obs(end-numvertices2d+1:end);
        
            % Get drag coefficient:
            ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
        
            % Get cost function coefficients:
            for ii=1:length(costfunctions)
                if costfunctions(ii)==101
                    ObsCostCoefficient=md.inversion.cost_functions_coefficients(end-numvertices2d+1:end,ii);
                elseif costfunctions(ii)==501
                    RegCostCoefficient=md.inversion.cost_functions_coefficients(end-numvertices2d+1:end,ii);
                else
                    error('Sorry, I have only implemented cost functions 101 and 501')
                end
            end
        
            % Interpolate reg coefficient to elements:
            RegCostCoefficient_elements=mean(RegCostCoefficient(elements2d),2);
            
            % Compute obs misfit:
            ObsMisfit=.5*((ModelVx-ObsVx).^2+(ModelVy-ObsVy).^2)/secondsperyear^2;
        
            % Compute reg misfit:
            % Compute gradients:
            DragCoeffGradX_elements=sum(BasisGradX_elements.*ModelDragCoeff(elements2d),2);
            DragCoeffGradY_elements=sum(BasisGradY_elements.*ModelDragCoeff(elements2d),2);
        
            % Compute misfit from squared magnitude of gradient:
            RegMisfit_elements=.5*(DragCoeffGradX_elements.^2+DragCoeffGradY_elements.^2);
        
            % Determine which elements are in our region:
            in_element=mean(in_node(elements2d),2);  % can have fractional values, that's okay
        
            % Integrate over area:
            J0=sum(mean(ObsMisfit(elements2d).*ObsCostCoefficient(elements2d),2).*Area.*in_element);   % interpolate product of obs misfit and coefficient to the elements, rather than interpolate individually and then multiply
            R=sum(RegMisfit_elements.*Area.*RegCostCoefficient_elements.*in_element);
        
            % Record and display values:
            DataCost(jj)=J0;
            RegularizationCost(jj)=R/Lambda(jj);   
        end 

        %% PREPARE L-curve:
        % Identify outlier models:
        % preparation:
        % if thismodel==2
        %     pickOutlier_m3_Nocp; 
        %     pickOutlier_m3_Ncuas;
        % end
        % IsOutlier=isnan(DataCost);
        done=0;
        numiterations_outlierid=1;
        % loop:
        while done==0
            % record last outlier list:
            IsOutlier_last=IsOutlier;
            % isolateremaining good points:
            TheseDataCost=DataCost(~IsOutlier_last);
            TheseRegCost=RegularizationCost(~IsOutlier_last);
            % check remaining good points for new outliers:
            IsOutlier(~IsOutlier_last)=[TheseDataCost(1)>outlierthreshold*TheseDataCost(2);TheseDataCost(2:end-1)>outlierthreshold*sqrt(TheseDataCost(1:end-2).*TheseDataCost(3:end));...
                TheseDataCost(end)>outlierthreshold*TheseDataCost(end-1)]|...
                [TheseRegCost(1)>outlierthreshold*TheseRegCost(2);TheseRegCost(2:end-1)>outlierthreshold*sqrt(TheseRegCost(1:end-2).*TheseRegCost(3:end));...
                TheseRegCost(end)>outlierthreshold*TheseRegCost(end-1)];
            % check for change in outlier list:
            if isequal(IsOutlier,IsOutlier_last)
                % break from loop:
                done=1;
            elseif numiterations_outlierid>=numlambda
                error('Unable to identify outlier models')
            else
                % advance counter:
                numiterations_outlierid=numiterations_outlierid+1;
            end
        end

        % Convert points to log scale:
        LogDataCost_points=log(DataCost);
        LogRegCost_points=log(RegularizationCost);
        LogLambda_points=log(Lambda);
        
        % Compute lambda subsampling:
        LogLambda_path=linspace(LogLambda_points(find(~IsOutlier,1,'first')),LogLambda_points(find(~IsOutlier,1,'last')),numlambdasubsamples)';
        
        % Interpolate log(variables) onto the new lambdasubsampling:
        LogDataCost_path=interp1(LogLambda_points(~IsOutlier),LogDataCost_points(~IsOutlier),LogLambda_path,lambdainterpmethod);
        LogRegCost_path=interp1(LogLambda_points(~IsOutlier),LogRegCost_points(~IsOutlier),LogLambda_path,lambdainterpmethod);
        
        % Compute dlambda:
        dloglambda_path=LogLambda_path(2)-LogLambda_path(1);
        
        % Compute minimum lambda wavelength:
        minlambdawavelength=.5*(LogLambda_points(2)-LogLambda_points(1));
        
        % Generate wavelength sampling:
        LambdaWavelengths=linspace(minlambdawavelength,maxlambdawavelength,numlambdawavelengths);
        
        % Compute characteristic scales:
        if strcmp(curvenormstyle,'quad')
            % Estimate scales from quadratic best-fits:
            A=[ones(sum(~IsOutlier),1),LogLambda_points(~IsOutlier),LogLambda_points(~IsOutlier).^2];
            equation_data=(A'*A)\A'*LogDataCost_points(~IsOutlier);
            equation_reg=(A'*A)\A'*LogRegCost_points(~IsOutlier);
            Fit_data=A*equation_data;
            Fit_reg=A*equation_reg;
            Misfit_data=LogDataCost_points(~IsOutlier)-Fit_data;
            Misfit_reg=LogRegCost_points(~IsOutlier)-Fit_reg;
            curvescale=sqrt(equation_data(3)^2+equation_reg(3)^2);
            scatterscale=sqrt(mean(Misfit_data.^2+Misfit_reg.^2));
        elseif strcmp(curvenormstyle,'max')
            curvescale=1;
            scatterscale=1;
        else
            error('Parameter "curvenormstyle" must be a string equal to "quad" or "max".')
        end
        
        % Pre-allocate:
        LogDataCost_curves=zeros(numlambdasubsamples,numlambdawavelengths);
        LogRegCost_curves=zeros(numlambdasubsamples,numlambdawavelengths);
        CurvatureVariance_curves=zeros(1,numlambdawavelengths);
        ScatterVariance_curves=zeros(1,numlambdawavelengths);
            
        % Loop through potential smoothings:
        for ii=1:numlambdawavelengths
            % Check wavelength:
            if LambdaWavelengths(ii)>0
                % Smooth the variables:
                LogDataCost_curves(:,ii)=intuitive_lowpass(LogDataCost_path,LambdaWavelengths(ii)/dloglambda_path,'quasilocal');
                LogRegCost_curves(:,ii)=intuitive_lowpass(LogRegCost_path,LambdaWavelengths(ii)/dloglambda_path,'quasilocal');
            else
                % No smoothing required:
                LogDataCost_curves(:,ii)=LogDataCost_path;
                LogRegCost_curves(:,ii)=LogRegCost_path;
            end
            % Compute curvature of this curve:
            ThisLogDataCurvature=(LogDataCost_curves(1:end-2,ii)+LogDataCost_curves(3:end,ii)-2*LogDataCost_curves(2:end-1,ii))/(dloglambda_path^2);
            ThisLogRegCurvature=(LogRegCost_curves(1:end-2,ii)+LogRegCost_curves(3:end,ii)-2*LogRegCost_curves(2:end-1,ii))/(dloglambda_path^2);
            % Compute curvature variance:
            if strcmp(curvevarstyle,'var')
                CurvatureVariance_curves(ii)=mean(ThisLogDataCurvature.^2+ThisLogRegCurvature.^2)/(curvescale^2);
            elseif strcmp(curvevarstyle,'rms')
                CurvatureVariance_curves(ii)=sqrt(mean(ThisLogDataCurvature.^2+ThisLogRegCurvature.^2))/curvescale;
            else
                error('Parameter "curvevarstyle" must be a string equal to "var" or "rms".')
            end
            % Compute scatter about this curve:
            ThisDataScatter=LogDataCost_points(~IsOutlier)-interp1(LogLambda_path,LogDataCost_curves(:,ii),LogLambda_points(~IsOutlier),lambdainterpmethod);
            ThisRegScatter=LogRegCost_points(~IsOutlier)-interp1(LogLambda_path,LogRegCost_curves(:,ii),LogLambda_points(~IsOutlier),lambdainterpmethod);
            % Compute scatter variance:
            if strcmp(curvevarstyle,'var')
                ScatterVariance_curves(ii)=mean(ThisDataScatter.^2+ThisRegScatter.^2)/(scatterscale^2);
            else % rms
                ScatterVariance_curves(ii)=sqrt(mean(ThisDataScatter.^2+ThisRegScatter.^2))/scatterscale;
            end
        end
        
        % Finish normalizing curvature variances:
        if strcmp(curvenormstyle,'max')
            CurvatureVariance_curves=CurvatureVariance_curves/max(CurvatureVariance_curves);
            ScatterVariance_curves=ScatterVariance_curves/max(ScatterVariance_curves);
        end
        
        % Determine best wavelength:
        TotalVariance_curves=ScatterVariance_curves+CurvatureVariance_curves;
        bestwavelengthind=find(TotalVariance_curves==min(TotalVariance_curves),1,'first');
        bestwavelength=LambdaWavelengths(bestwavelengthind);
        
        % Isolate best curve:
        LogDataCost_bestcurve=LogDataCost_curves(:,bestwavelengthind);
        LogRegCost_bestcurve=LogRegCost_curves(:,bestwavelengthind);
        
        % Compute gradient and curvature of best curve:
        LogDataGradient_bestcurve=(LogDataCost_bestcurve(2:end)-LogDataCost_bestcurve(1:end-1))/dloglambda_path;
        LogRegGradient_bestcurve=(LogRegCost_bestcurve(2:end)-LogRegCost_bestcurve(1:end-1))/dloglambda_path;
        LogDataCurvature_bestcurve=(LogDataCost_bestcurve(1:end-2)+LogDataCost_bestcurve(3:end)-2*LogDataCost_bestcurve(2:end-1))/(dloglambda_path^2);
        LogRegCurvature_bestcurve=(LogRegCost_bestcurve(1:end-2)+LogRegCost_bestcurve(3:end)-2*LogRegCost_bestcurve(2:end-1))/(dloglambda_path^2);
        TotalCurvature_bestcurve=LogDataCurvature_bestcurve+LogRegCurvature_bestcurve;
        
        % Indexing note:  curvature has size numlambdasubsamples-2.
        % ie, it is missing the first and last elements.  The "find"
        % commands below have been modified so that the resulting index
        % points to the correct position in the original list.
        
        % Identify best lambda:
        maxcurve=max(TotalCurvature_bestcurve(ceil(.5*bestwavelength/dloglambda_path):floor(numlambdasubsamples-.5*bestwavelength/dloglambda_path)));
        bestlambdaind=1+find(TotalCurvature_bestcurve==maxcurve,1,'first');
        bestlambda=exp(LogLambda_path(bestlambdaind));
        bestmodelind=find(abs(LogLambda_points-log(bestlambda))==min(abs(LogLambda_points-log(bestlambda))),1,'first');
        
        % Compute max/min allowable lambdas:
        bestlambdaind_minallowable=find(TotalCurvature_bestcurve(1:bestlambdaind)<.5*maxcurve,1,'last');
        if isempty(bestlambdaind_minallowable)
            bestlambdaind_minallowable=1;
        end
        bestlambdaind_maxallowable=bestlambdaind-1+find(TotalCurvature_bestcurve(bestlambdaind-1:end)<.5*maxcurve,1,'first');
        if isempty(bestlambdaind_maxallowable)
            bestlambdaind_maxallowable=numlambdasubsamples;
        end
        bestlambda_minallowable=exp(LogLambda_path(bestlambdaind_minallowable));
        bestlambda_maxallowable=exp(LogLambda_path(bestlambdaind_maxallowable));
        
        % Compute uncertainty:
        bestlambda_loguncertainty=.5*bestwavelength;
        
        % Display results:
        disp(['Min acceptable lambda=',num2str(bestlambda_minallowable)])
        disp(['Best lambda=',num2str(bestlambda)])
        disp(['Max acceptable lambda=',num2str(bestlambda_maxallowable)])
        disp(['Lambda uncertainty ratio=',num2str(exp(bestlambda_loguncertainty))])
        disp(['Drag coeff scale=',num2str(dragcoeffstdev)])
        disp(['Velocity scale=',num2str(velrms)])
        disp(['Best model Jobs=',num2str(DataCost(bestmodelind))])
        disp(['Best model Jreg=',num2str(RegularizationCost(bestmodelind))])
        
        %% MAKE L-curve FIGURE
        if thismodel==1
            FigH = figure('Position', get(0, 'Screensize'));
        end
        % Plot cost vs cost, log scale:
        % Plot model costs:
        hgoodpoints=plot(RegularizationCost(~IsOutlier),DataCost(~IsOutlier),'Marker',modelmarker,'MarkerFaceColor',lcurvecolors{thismodel},'Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',modelmarkersize);
        hold on
        if sum(IsOutlier)>0
            houtliers=plot(RegularizationCost(IsOutlier),DataCost(IsOutlier),'Marker',modelmarker,'MarkerFaceColor','none','Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',modelmarkersize,'LineWidth',Linewidth);
        end
        % Plot best curve:
        hcurve=plot(exp(LogRegCost_bestcurve),exp(LogDataCost_bestcurve),'Color',lcurvecolors{thismodel},'LineStyle',curvelinestyle,'LineWidth',curvelinesize);
        % Plot corner region:
        hcorner(thismodel)=plot(exp(LogRegCost_bestcurve(bestlambdaind_minallowable:bestlambdaind_maxallowable)),...
            exp(LogDataCost_bestcurve(bestlambdaind_minallowable:bestlambdaind_maxallowable)),'Marker','none','Color',lcurvecolors{thismodel},'LineStyle',curvelinestyle,'LineWidth',cornerlinesize); 

        % Plot and label best value:
        angle1=atan2(LogDataCost_bestcurve(bestlambdaind+1)-LogDataCost_bestcurve(bestlambdaind),...
            LogRegCost_bestcurve(bestlambdaind+1)-LogRegCost_bestcurve(bestlambdaind));
        plot(exp(LogRegCost_bestcurve(bestlambdaind)),exp(LogDataCost_bestcurve(bestlambdaind)),'Color',lcurvecolors{thismodel},'MarkerFaceColor','w','Marker',bestmarker,'MarkerSize',bestmarkersize,'LineWidth',LabelLinewidth)
        lastpow10=10^floor(log10(bestlambda));
%         plot(exp(LogRegCost_bestcurve(bestlambdaind)-[0,.2*cos(angle1+pi/2)]),exp(LogDataCost_bestcurve(bestlambdaind)-[0,.2*sin(angle1+pi/2)]),'Color',lcurvecolors{thismodel},'LineWidth',Linewidth)
%         text(exp(LogRegCost_bestcurve(bestlambdaind)+.2),exp(LogDataCost_bestcurve(bestlambdaind)+.2),{['\lambda_{best}=',num2str(.1*lastpow10*round(10*bestlambda/lastpow10))]},...
%             'HorizontalAlignment','left','VerticalAlignment','bottom','Color',lcurvecolors{thismodel},'FontSize',Fontsize) %'Best value:';

        % Set up axes:
        xlim(jreglims.*[.999,1.001])
        ylim(jobslims.*[.999,1.001])
        set(gca,'TickDir','in')%'out')
        set(gca,'XScale','log')
        set(gca,'YScale','log')
        set(gca,'PlotBoxAspectRatio',[(log(jreglims(2))-log(jreglims(1)))/(log(jobslims(2))-log(jobslims(1))),1,1])
 %       xlabel('Regularization cost J_{reg} (unitless)','FontSize',titlefontsize)
 %       ylabel('Data cost J_{obs} (unitless)','FontSize',titlefontsize)
        set(gca,'FontSize',axesfontsize,'LineWidth',axeslinewidth)
        text(2e-9,5e-1,['(',char(96+1+i),')'],'color','k','FontSize',Fontsize,'HorizontalAlignment','left') 
        hold on 
    end 
    % Make legend:
    legend(hcorner,descriptext,'location','NorthEast','FontSize',legendfontsize)

    % Save figure:
    saveas(FigH,savename,'png');
    exportgraphics(FigH,savename,'resolution',resolution)

end

disp('Done!')
toc