% LcurveComparisonFigure_N_m_v1

% Lea-Sophie Höyns, 01/16/2024

% This script makes a figure comparing the results of different effective 
% pressure fields.  Makes a table with stats about the comparison, and makes 
% a figure showing comparative L-curves. 

% Major part of this script is taken from the script InversionComparisonFigure_v1.m 
% by M. Wolovick, 12/9/2022 (https://doi.org/10.5281/zenodo.7798650).

% a. Different L-curves (m,Weertman; m,Nop; m,Ncuas); b. Cost curvature for 
% m, Weertman; c. Cost curvature for m, Budd using Nop; d. Cost curvature 
% for m, Budd using Ncuas

clearvars
tic
close all

%% PATHS
issmscriptpath='/albedo/work/projects/landice/ISSM/ISSM_installations/ISSM_v4p22/trunk/';
functionpath='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Functions';
addpath(genpath(issmscriptpath),functionpath)
% File names:
lcurvefiles={'/albedo/work/projects/p_modsem23/Lea/WAIS/OldLcurveRange/Run_3/WAIS.dragInv_HO_m3_noN_1DTemp_Lcurve_20-Jan-2024.mat';... % cell array of strings, include trailing ".mat", first one is reference
        '/albedo/work/projects/p_modsem23/Lea/WAIS/OldLcurveRange/Run_2/WAIS.dragInv_HO_m3_Nocp_1DTemp_Lcurve_19-Jan-2024.mat';...
        '/albedo/work/projects/p_modsem23/Lea/WAIS/OldLcurveRange/Run_1/WAIS.dragInv_HO_m3_Ncuas_1DTemp_Lcurve_17-Jan-2024.mat'};
%lcurvefiles={'/albedo/work/projects/p_modsem23/Lea/WAIS/Convergence_test/Run_9/WAIS.dragInv_HO_m1_noN_1DTemp_Lcurve_08-Jan-2024.mat';... 
%        '/albedo/work/projects/p_modsem23/Lea/WAIS/Convergence_test/Run_7/WAIS.dragInv_HO_m1_Nocp_1DTemp_Lcurve_06-Jan-2024.mat';...
%        '/albedo/work/projects/p_modsem23/Lea/WAIS/Convergence_test/Run_6/WAIS.dragInv_HO_m1_Ncuas_1DTemp_Lcurve_05-Jan-2024.mat'};
figname_lcurves='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/LcurveComparisonFigure_N_m3_v1_final.png'; 

%% PARAMETERS 
% N descriptions: %'WeertmanReference',
ndescriptext={'Weertman','Budd, N_{op}','Budd, N_{CUAS}'};  % cell array of strings
descriptext={'m=3, Weertman','m=3, Budd, N_{op}','m=3, Budd, N_{CUAS}'};  % cell array of strings
%descriptext={'m=1, Weertman','m=1, Budd, N_{op}','m=1, Budd, N_{CUAS}'};  % cell array of strings

% Ticks and limits:
jreglims=[1e-4,1e-1];                 % [1x2] unitless
jobslims=[1e-3,1];                    % [1x2] unitless
%jreglims=[1e-5,0.02]                   % [1x2] unitless
%jobslims=[1e-3,1];                     % [1x2] unitless

% Aesthetic settings
% Colors:
orange='#ED553B';                             % valid color string
yellow='#F2B134';                             % valid color string
turquoise='#4FB99F';                          % valid color string
blue='#068587';                               % valid color string
darkblue='#112F41';                           % valid color string
newblue='#5EA3D8';                            % valid color string
gray='#8C8B8C';                               % valid color string
lcurvecolors={yellow,darkblue,turquoise};     % {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=22;                              % points
titlefontsize=22;                             % points
legendfontsize=22;                            % points
axesfontsize2=18;                             % points
titlefontsize2=18;                            % points
legendfontsize2=10;                           % points
modelmarkersize=9;                            % points
curvelinesize=3;                              % points
cornerlinesize=8;                             % points
bestmarkersize=18;                            % points
legendfontsize3=17;
% Linewidths: 
axeslinewidth=2;                              % points
% Styles:
modelmarker='o';                              % valid marker string
curvelinestyle='-';                           % valid linestyle string
cornerlinestyle='-';                          % valid linestyle string
bestmarker='d';                               % valid marker string

% Threshold for analyzing fast-flow areas:
velanalysisthrehsold=15;               % m/yr (smoothing weighting btwn 1/2 and 2x this value)
minvelstabilizer=1e-3;                 % m/yr
mintaustabilizer=1;                    % Pa

% Output size and resolution:
TileSpacing='compact';%'tight'
Padding='compact';
resolution=300;                        % dpi

%% PREPARATION

% Check number of N:
numlcurves=length(lcurvefiles);
if length(ndescriptext)~=numlcurves
    error('Length of "ndescriptext" did not match length of "lcurvefiles".')
end

% 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 hydraulic diameter:
HydraulicDiameter_elements=2*Area./s;
% Compute element centers:
CentX=mean(x2d(elements2d),2);
CentY=mean(y2d(elements2d),2);
% Find mesh limits:
xlims=[min(x2d),max(x2d)];
ylims=[min(y2d),max(y2d)];
% Extract masks on 2D mesh and truncate:
OceanLevelSet_vertices=max(-1,min(1,md.mask.ocean_levelset(1:numvertices2d)));
IceLevelSet_vertices=max(-1,min(1,md.mask.ice_levelset(1:numvertices2d)));
% Compute ocean and ice area fractions:
OceanAreaFraction_elements=(1-mean(OceanLevelSet_vertices(elements2d),2))/2;
IceAreaFraction_elements=(1-mean(IceLevelSet_vertices(elements2d),2))/2;
% Compute area totals:
totalarea=sum(Area);
totalicearea=sum(Area.*IceAreaFraction_elements);
totalgroundedarea=sum(Area.*IceAreaFraction_elements.*(1-OceanAreaFraction_elements));
% 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)))];

% Generate the active mask:
ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0&md.mesh.vertexonboundary(1:numvertices2d)==0;
ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);

% Unpack the good data mask:
GoodData_obs_vertices=md.miscellaneous.dummy(1:numvertices2d);
GoodDataWeight_elements=mean(GoodData_obs_vertices(elements2d).^2,2);

% Extract observed velocity:
U_obs_vertices=md.inversion.vx_obs(end-numvertices2d+1:end);
V_obs_vertices=md.inversion.vy_obs(end-numvertices2d+1:end);
VelMag_obs_vertices=md.inversion.vel_obs(end-numvertices2d+1:end);

% Compute weighting function for fast-flowing areas:
FastFlowWeight_vertices=ones(numelements2d,1);
FastFlowWeight_vertices(VelMag_obs_vertices<.5*velanalysisthrehsold)=0;
FastFlowWeight_vertices(VelMag_obs_vertices>=.5*velanalysisthrehsold&VelMag_obs_vertices<2*velanalysisthrehsold)=...
    .5*(1+sin(.5*pi*log2(VelMag_obs_vertices(VelMag_obs_vertices>=.5*velanalysisthrehsold&VelMag_obs_vertices<2*velanalysisthrehsold)/velanalysisthrehsold)));
FastFlowWeight_elements=mean(FastFlowWeight_vertices(elements2d),2);

%% COLLECT DATA 

% Pre-allocate:
EffectivePressure_vertices=ones(numvertices2d,numlcurves);
EffectivePressure_elements=ones(numelements2d,numlcurves);
DragCoeff_vertices=zeros(numvertices2d,numlcurves);
DragCoeff_elements=zeros(numelements2d,numlcurves);
Drag_vertices=zeros(numvertices2d,numlcurves);
Drag_elements=zeros(numelements2d,numlcurves);
VelMag_vertices=zeros(numvertices2d,numlcurves);   
VelMag_elements=zeros(numelements2d,numlcurves);  
VelMagMisfit_vertices=zeros(numvertices2d,numlcurves);   
VelMagMisfit_elements=zeros(numelements2d,numlcurves);  
ObsMisfit_vertices=zeros(numvertices2d,numlcurves);   
ObsMisfit_elements=zeros(numelements2d,numlcurves);   
RegMisfit_elements=zeros(numelements2d,numlcurves);  % no interpolation from vertices to elements required for reg term
BestLambda=zeros(numlcurves,1);
MaxLambda=zeros(numlcurves,1);
MinLambda=zeros(numlcurves,1);
DeltaLambdaRatio=zeros(numlcurves,1);
DragCoeffScale=zeros(numlcurves,1);
VelScale=zeros(numlcurves,1);
InversionParameters=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',[]);

% Loop through models:
for thismodel=1:numlcurves

    % Check if we're loading remote files:
    if exist('serverloc','var')
        % Parse filename:
        slashnum=strfind(lcurvefiles{thismodel},'/');
        remotepath=lcurvefiles{thismodel}(1:slashnum(end));
        lcurvename=lcurvefiles{thismodel}(slashnum(end)+1:end);
        % Fetch the file:
        if ~exist([localpath,lcurvename],'file')
            ssh2_struct=scp_simple_get(serverloc,serverusername,serverpassword,lcurvename,localpath,remotepath);
        end
        % Load the basic file:
        load([localpath,lcurvename],'Log*Cost_bestcurve','*Curvature_bestcurve','Lambda','DataCost','RegularizationCost','IsOutlier','bestlambda*','bestmodelind','MyInversionParameters')
    else
        % Load the basic file:
        load(lcurvefiles{thismodel},'Log*Cost_bestcurve','*Curvature_bestcurve','Lambda','DataCost','RegularizationCost','IsOutlier','bestlambda*','bestmodelind','MyInversionParameters')
    end
    
    % Unpack parameters:
    unpack(MyInversionParameters)
    
    % Record parameters:
    InversionParameters{thismodel}=MyInversionParameters;

    % Construct this file name:
    slashnum=strfind(lcurvefiles{thismodel},'/');
    thisbase=10^floor(log10(Lambda(bestmodelind)));
    thisname=[lcurvefiles{thismodel}(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(bestmodelind)/thisbase)),'.mat'];

    % Check if we're loading remote files:
    if exist('serverloc','var')
        % Fetch the file:
        if ~exist([localpath,thisname],'file')
            ssh2_struct=scp_simple_get(serverloc,serverusername,serverpassword,thisname,localpath,remotepath);
        end
        % Load the file:
        load([localpath,thisname],'md')
    else
        % Load this model:
        load([lcurvefiles{thismodel}(1:slashnum(end)),thisname],'md')
    end

    % Assign simple things:
    BestLambda(thismodel)=bestlambda;
    MaxLambda(thismodel)=bestlambda_maxallowable;
    MinLambda(thismodel)=bestlambda_minallowable;
    DeltaLambdaRatio(thismodel)=exp(bestlambda_loguncertainty);
    DragCoeffScale(thismodel)=dragcoeffstdev;
    
    % Assign velocity scale:
     if strcmp(velabsscale,'mean')
         VelScale(thismodel)=velmean;
    elseif strcmp(velabsscale,'median')
        VelScale(thismodel)=velmedian;
    elseif strcmp(velabsscale,'stdev')
        VelScale(thismodel)=velstdev;
    elseif strcmp(velabsscale,'rms')
        VelScale(thismodel)=velrms;
     else
        error('I do not recognize the value of "velabsscale".')
    end
    
    % Assign lcurve data:
    LCurveData(thismodel).LogLambda_points=log(Lambda);
    LCurveData(thismodel).LogRegCost_points=log(RegularizationCost);
    LCurveData(thismodel).LogDataCost_points=log(DataCost);
    LCurveData(thismodel).IsOutlier=IsOutlier;
    LCurveData(thismodel).LogRegCost_curve=LogRegCost_bestcurve;
    LCurveData(thismodel).LogDataCost_curve=LogDataCost_bestcurve;
    LCurveData(thismodel).LogDataCurvature_curve=LogDataCurvature_bestcurve;
    LCurveData(thismodel).LogRegCurvature_curve=LogRegCurvature_bestcurve;
    LCurveData(thismodel).TotalCurvature_curve=TotalCurvature_bestcurve;
    LCurveData(thismodel).bestlambdaind=bestlambdaind;
    LCurveData(thismodel).bestlambdaind_maxallowable=bestlambdaind_maxallowable;
    LCurveData(thismodel).bestlambdaind_minallowable=bestlambdaind_minallowable;

    % Check mesh size:
    if numlayers>1 && ~strcmp(flowequation,'SSA')
        if md.mesh.numberofvertices2d~=numvertices2d || md.mesh.numberofelements2d~=numelements2d
            error(['Mesh number ',num2str(thismodel),' does not match reference mesh.'])
        end
    else
        if md.mesh.numberofvertices~=numvertices2d || md.mesh.numberofelements~=numelements2d
            error(['Mesh number ',num2str(thismodel),' does not match reference mesh.'])
        end
    end

    % Extract variables we need:
    if numlayers>1 && ~strcmp(flowequation,'SSA')
        ModelSurfU=md.results.StressbalanceSolution.Vx(md.mesh.vertexonsurface==1);
        ModelSurfV=md.results.StressbalanceSolution.Vy(md.mesh.vertexonsurface==1);
        ModelSurfVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonsurface==1);
        ModelBasalVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonbase==1);
        DragCoeff_vertices(:,thismodel)=md.results.StressbalanceSolution.FrictionCoefficient(md.mesh.vertexonbase==1);
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices(:,thismodel)=md.friction.effective_pressure(md.mesh.vertexonbase==1);
        else
            EffectivePressure_vertices(:,thismodel)=ones(numvertices2d,1);
        end
    else
        ModelBasalVel=md.results.StressbalanceSolution.Vel;
        ModelSurfVel=ModelBasalVel;
        ModelSurfU=md.results.StressbalanceSolution.Vx;
        ModelSurfV=md.results.StressbalanceSolution.Vy;
        DragCoeff_vertices(:,thismodel)=md.results.StressbalanceSolution.FrictionCoefficient;
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices(:,thismodel)=md.friction.effective_pressure;
        else
            EffectivePressure_vertices(:,thismodel)=ones(numvertices2d,1);
        end
    end
    
    % Interpolate cost function coefficients to elements:
    AreReg=costfunctions>500;
    CostFunctionCoeffs_elements=zeros(numelements2d,length(costfunctions));
    for jj=1:length(costfunctions)
        CostFunctionCoeffs_elements(:,jj)=mean(reshape(md.inversion.cost_functions_coefficients(elements2d,jj),[numelements2d,3]),2);
        if AreReg(jj)
            CostFunctionCoeffs_elements(:,jj)=CostFunctionCoeffs_elements(:,jj)/Lambda(bestmodelind);
        end
    end
    
    % Compute regularization term:
    % Compute gradients:
    DragCoeffGradX_elements=sum(BasisGradX_elements.*reshape(DragCoeff_vertices(elements2d,thismodel),[numelements2d,3]),2);
    DragCoeffGradY_elements=sum(BasisGradY_elements.*reshape(DragCoeff_vertices(elements2d,thismodel),[numelements2d,3]),2);
    % Compute (squared) magnitude and apply coefficients:
    RegMisfit_elements(:,thismodel)=.5*(DragCoeffGradX_elements.^2+DragCoeffGradY_elements.^2).*CostFunctionCoeffs_elements(:,AreReg);

    % Replace zeros:
    DragCoeff_vertices(:,thismodel)=max(mindragcoeff,DragCoeff_vertices(:,thismodel));

    % Compute drag:
    Drag_vertices(:,thismodel)=(DragCoeff_vertices(:,thismodel).^2).*((ModelBasalVel/secondsperyear).^(1/m));
    if ~strcmp(nsource,'none')
        Drag_vertices(:,thismodel)=Drag_vertices(:,thismodel).*(EffectivePressure_vertices(:,thismodel).^r);
    end
    % Set drag to zero under floating shelf and in ice-free areas:
    Drag_vertices(OceanLevelSet_vertices<0|IceLevelSet_vertices>0,thismodel)=0;

    % Compute velocity magnitude misfit:
    VelMagMisfit_vertices(:,thismodel)=ModelSurfVel-VelMag_obs_vertices;
    
    % Compute observational misfit:
    ObsMisfit_vertices(:,thismodel)=.5*((ModelSurfU-U_obs_vertices).^2+(ModelSurfV-V_obs_vertices).^2)/(secondsperyear^2);
    
    % Assign (sliding) velocity:
    VelMag_vertices(:,thismodel)=ModelBasalVel;

    % Interpolate to elements:
    EffectivePressure_elements(:,thismodel)=mean(reshape(EffectivePressure_vertices(elements2d,thismodel),[numelements2d,3]),2);
    DragCoeff_elements(:,thismodel)=mean(reshape(DragCoeff_vertices(elements2d,thismodel),[numelements2d,3]),2);
    Drag_elements(:,thismodel)=mean(reshape(Drag_vertices(elements2d,thismodel),[numelements2d,3]),2);
    VelMag_elements(:,thismodel)=mean(reshape(VelMag_vertices(elements2d,thismodel),[numelements2d,3]),2);
    VelMagMisfit_elements(:,thismodel)=mean(reshape(VelMagMisfit_vertices(elements2d,thismodel),[numelements2d,3]),2);
    ObsMisfit_elements(:,thismodel)=mean(reshape(ObsMisfit_vertices(elements2d,thismodel),[numelements2d,3]),2).*CostFunctionCoeffs_elements(:,~AreReg);  % apply coefficients as well

end

% Convert to ln(C) and ln(N):
LnC_elements=2*log(DragCoeff_elements);
% Lnk_elements=log(DragCoeff_elements);
% LnC_elements=Lnk_elements;
LnN_elements=log(EffectivePressure_elements);
LnU_elements=log(max(VelMag_elements,minvelstabilizer));
LnTau_elements=log(max(Drag_elements,mintaustabilizer));

%% COMPUTE STATS

% Compute total area for stats computation:
totalstatsarea=[sum(Area.*ActiveMask_elements.*GoodDataWeight_elements),sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements)];

% Pre-allocate:
meanc=zeros(numlcurves,2);
rmsc=zeros(numlcurves,2);
meanlnc=zeros(numlcurves,2);
rmslnc=zeros(numlcurves,2);
meann=zeros(numlcurves,2);
rmsn=zeros(numlcurves,2);
meanlnn=zeros(numlcurves,2);
rmslnn=zeros(numlcurves,2);
meantau=zeros(numlcurves,2);
rmstau=zeros(numlcurves,2);
meanlntau=zeros(numlcurves,2);
rmslntau=zeros(numlcurves,2);
meanlnu=zeros(numlcurves,2);
rmslnu=zeros(numlcurves,2);
rlnclnc=zeros(numlcurves,2);
rnc=zeros(numlcurves,2);
rlnnlnc=zeros(numlcurves,2);
rtautau=zeros(numlcurves,2);
rntau=zeros(numlcurves,2);
rmslnclnc=zeros(numlcurves,2);
rmstautau=zeros(numlcurves,2);
rmsvelmisfit=zeros(numlcurves,2);
jobs=zeros(numlcurves,2);
jreg=zeros(numlcurves,2);
jobsvelequiv=zeros(numlcurves,2);

% for all stats variables, first column is whole domain and second column
% is fast-flow domain.

% Loop through models:
for ii=1:numlcurves
    
    % Compute basic means and variances:
    % linear C:  (corrected from k to C)
    meanc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(DragCoeff_elements(:,ii).^2))/totalstatsarea(1);
    rmsc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(((DragCoeff_elements(:,ii).^2)-meanc(ii,1)).^2))/totalstatsarea(1));
    meanc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(DragCoeff_elements(:,ii).^2))/totalstatsarea(2);
    rmsc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(((DragCoeff_elements(:,ii).^2)-meanc(ii,2)).^2))/totalstatsarea(2));
    % lnC:
    meanlnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnC_elements(:,ii))/totalstatsarea(1);
    rmslnc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,1)).^2))/totalstatsarea(1));
    meanlnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnC_elements(:,ii))/totalstatsarea(2);
    rmslnc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,2)).^2))/totalstatsarea(2));
    % linear N:
    meann(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*EffectivePressure_elements(:,ii))/totalstatsarea(1);
    rmsn(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((EffectivePressure_elements(:,ii)-meann(ii,1)).^2))/totalstatsarea(1));
    meann(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*EffectivePressure_elements(:,ii))/totalstatsarea(2);
    rmsn(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((EffectivePressure_elements(:,ii)-meann(ii,2)).^2))/totalstatsarea(2));
    % lnN:
    meanlnn(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnN_elements(:,ii))/totalstatsarea(1);
    rmslnn(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnN_elements(:,ii)-meanlnn(ii,1)).^2))/totalstatsarea(1));
    meanlnn(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnN_elements(:,ii))/totalstatsarea(2);
    rmslnn(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnN_elements(:,ii)-meanlnn(ii,2)).^2))/totalstatsarea(2));
    % linear tau:
    meantau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*Drag_elements(:,ii))/totalstatsarea(1);
    rmstau(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-meantau(ii,1)).^2))/totalstatsarea(1));
    meantau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*Drag_elements(:,ii))/totalstatsarea(2);
    rmstau(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-meantau(ii,2)).^2))/totalstatsarea(2));
    % ln tau:
    meanlntau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnTau_elements(:,ii))/totalstatsarea(1);
    rmslntau(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnTau_elements(:,ii)-meanlntau(ii,1)).^2))/totalstatsarea(1));
    meanlntau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnTau_elements(:,ii))/totalstatsarea(2);
    rmslntau(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnTau_elements(:,ii)-meanlntau(ii,2)).^2))/totalstatsarea(2));
    % ln U:
    meanlnu(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnU_elements(:,ii))/totalstatsarea(1);
    rmslnu(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnU_elements(:,ii)-meanlnu(ii,1)).^2))/totalstatsarea(1));
    meanlnu(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnU_elements(:,ii))/totalstatsarea(2);
    rmslnu(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnU_elements(:,ii)-meanlnu(ii,2)).^2))/totalstatsarea(2));
    
    % Compute correlations:
    % R(N,Cw): (corrected from k to C)
    rnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,1)).*((DragCoeff_elements(:,1).^2)-meanc(1,1)))/(totalstatsarea(1)*rmsn(ii,1)*rmsc(1,1));
    rnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,2)).*((DragCoeff_elements(:,1).^2)-meanc(1,2)))/(totalstatsarea(2)*rmsn(ii,2)*rmsc(1,2));
    % R(ln(N),ln(Cw)):
    rlnnlnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(LnN_elements(:,ii)-meanlnn(ii,1)).*(LnC_elements(:,1)-meanlnc(1,1)))/(totalstatsarea(1)*rmslnn(ii,1)*rmslnc(1,1));
    rlnnlnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(LnN_elements(:,ii)-meanlnn(ii,2)).*(LnC_elements(:,1)-meanlnc(1,2)))/(totalstatsarea(2)*rmslnn(ii,2)*rmslnc(1,2));
    % R(ln(C),ln(Cw)):
    rlnclnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(LnC_elements(:,ii)-meanlnc(ii,1)).*(LnC_elements(:,1)-meanlnc(1,1)))/(totalstatsarea(1)*rmslnc(ii,1)*rmslnc(1,1));
    rlnclnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(LnC_elements(:,ii)-meanlnc(ii,2)).*(LnC_elements(:,1)-meanlnc(1,2)))/(totalstatsarea(2)*rmslnc(ii,2)*rmslnc(1,2));
    % R(tau,tau_w):
    rtautau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(Drag_elements(:,ii)-meantau(ii,1)).*(Drag_elements(:,1)-meantau(1,1)))/(totalstatsarea(1)*rmstau(ii,1)*rmstau(1,1));
    rtautau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(Drag_elements(:,ii)-meantau(ii,2)).*(Drag_elements(:,1)-meantau(1,2)))/(totalstatsarea(2)*rmstau(ii,2)*rmstau(1,2));
    % R(N,tau_w):
    rntau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,1)).*(Drag_elements(:,1)-meantau(1,1)))/(totalstatsarea(1)*rmsn(ii,1)*rmstau(1,1));
    rntau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,2)).*(Drag_elements(:,1)-meantau(1,2)))/(totalstatsarea(2)*rmsn(ii,2)*rmstau(1,2));
    
    % Compute RMS differences:
    rmslnclnc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,1)-(LnC_elements(:,1)-meanlnc(1,1))).^2))/totalstatsarea(1));
    rmslnclnc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,2)-(LnC_elements(:,1)-meanlnc(1,2))).^2))/totalstatsarea(2));
    rmstautau(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-Drag_elements(:,1)).^2))/totalstatsarea(1));
    rmstautau(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-Drag_elements(:,1)).^2))/totalstatsarea(2));
    
    % Compute RMS vel misfit:
    rmsvelmisfit(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(VelMagMisfit_elements(:,ii).^2))/totalstatsarea(1));
    rmsvelmisfit(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(VelMagMisfit_elements(:,ii).^2))/totalstatsarea(2));
    
    % Compute cost terms:
    jobs(ii,1)=sum(Area.*ObsMisfit_elements(:,ii));
    jobs(ii,2)=sum(Area.*FastFlowWeight_elements.*ObsMisfit_elements(:,ii))*(totalstatsarea(1)/totalstatsarea(2));
    jreg(ii,1)=sum(Area.*RegMisfit_elements(:,ii));
    jreg(ii,2)=sum(Area.*FastFlowWeight_elements.*RegMisfit_elements(:,ii))*(totalstatsarea(1)/totalstatsarea(2));
    
    % Compute velocity equivalent of obs misfit:
    jobsvelequiv(ii,:)=sqrt(jobs(ii,:))*VelScale(ii);
    
end

% Compute squared quantities:
varlnc=rmslnc.^2;
varlnn=rmslnn.^2;
varlntau=rmslntau.^2;
varlnu=rmslnu.^2;
r2nc=rnc.^2;
r2lnnlnc=rlnnlnc.^2;
r2lnclnc=rlnclnc.^2;
r2tautau=rtautau.^2;
r2ntau=rntau.^2;

% compute variance ratios:
varratiolnclnc=varlnc./repmat(varlnc(1,:),[numlcurves,1]);
varratiolnclntau=varlnc./varlntau;
varratiolnclnu=varlnc./varlnu;
varratiojobsjobs=jobs./repmat(jobs(1,:),[numlcurves,1]);
varratiototal=varratiolnclnc.*varratiojobsjobs;

%% COMMUNICATE STATS

% Display stats, whole domain:
string1='RMS(ln(C))=';
string2='var(ln(C))=';
string3='Variance ratio (ln(C)/ln(C1))=';
string4='R(N,C1)=';
string5='R2(N,C1)=';
string6='R(ln(N),ln(C1))=';
string7='R2(ln(N),ln(C1))=';
string8='R(ln(C),ln(C1))=';
string9='R2(ln(C),ln(C1))=';
string10='R(tau,tau1)=';
string11='R2(tau,tau1)=';
string12='R(N,tau1)=';
string13='R2(N,tau1)=';
string14='RMS(ln(C)-ln(C1))=';
string15='RMS(tau-tau1)=';
string16='RMS(u_model-u_obs)=';
string17='Jobs=';
string18='Jreg=';
string19='RMS(ln(N))=';
string20='var(ln(N))=';
string21='Jobs, equiv du=';
string22='RMS(ln(tau))=';
string23='var(ln(tau))=';
string24='RMS(ln(Ub))=';
string25='var(ln(Ub))=';
string26='Variance ratio (ln(C)/ln(tau))=';
string27='Variance ratio (ln(C)/ln(u))=';
string28='Variance ratio (Jobs/Jobs(1))=';
string29='Variance ratio total (Jobsratio*lnCratio)=';
for ii=1:numlcurves
    string1=[string1,num2str(.01*round(100*rmslnc(ii,1))),', '];
    string2=[string2,num2str(.01*round(100*varlnc(ii,1))),', '];
    string3=[string3,num2str(.01*round(100*varratiolnclnc(ii,1))),', '];
    string4=[string4,num2str(.01*round(100*rnc(ii,1))),', '];
    string5=[string5,num2str(.01*round(100*r2nc(ii,1))),', '];
    string6=[string6,num2str(.01*round(100*rlnnlnc(ii,1))),', '];
    string7=[string7,num2str(.01*round(100*r2lnnlnc(ii,1))),', '];
    string8=[string8,num2str(.01*round(100*rlnclnc(ii,1))),', '];
    string9=[string9,num2str(.01*round(100*r2lnclnc(ii,1))),', '];
    string10=[string10,num2str(.01*round(100*rtautau(ii,1))),', '];
    string11=[string11,num2str(.01*round(100*r2tautau(ii,1))),', '];
    string12=[string12,num2str(.01*round(100*rntau(ii,1))),', '];
    string13=[string13,num2str(.01*round(100*r2ntau(ii,1))),', '];
    string14=[string14,num2str(.01*round(100*rmslnclnc(ii,1))),', '];
    string15=[string15,num2str(.01*round(100*rmstautau(ii,1)/1000)),', '];
    string16=[string16,num2str(.01*round(100*rmsvelmisfit(ii,1))),', '];
    string17=[string17,num2str(.01*round(100*jobs(ii,1)*10000)),', '];
    string18=[string18,num2str(.01*round(100*jreg(ii,1)*10000)),', '];
    string19=[string19,num2str(.01*round(100*rmslnn(ii,1))),', '];
    string20=[string20,num2str(.01*round(100*varlnn(ii,1))),', '];
    string21=[string21,num2str(.01*round(100*jobsvelequiv(ii,1))),', '];
    string22=[string22,num2str(.01*round(100*rmslntau(ii,1))),', '];
    string23=[string23,num2str(.01*round(100*varlntau(ii,1))),', '];
    string24=[string24,num2str(.01*round(100*rmslnu(ii,1))),', '];
    string25=[string25,num2str(.01*round(100*varlnu(ii,1))),', '];
    string26=[string26,num2str(.01*round(100*varratiolnclntau(ii,1))),', '];
    string27=[string27,num2str(.01*round(100*varratiolnclnu(ii,1))),', '];
    string28=[string28,num2str(.01*round(100*varratiojobsjobs(ii,1))),', '];
    string29=[string29,num2str(.01*round(100*varratiototal(ii,1))),', '];
end
disp('...')
disp('Whole Domain Stats:')
disp(string1(1:end-1))
disp(string2(1:end-1))
disp(string3(1:end-1))
disp(string4(1:end-1))
disp(string5(1:end-1))
disp(string6(1:end-1))
disp(string7(1:end-1))
disp(string8(1:end-1))
disp(string9(1:end-1))
disp(string10(1:end-1))
disp(string11(1:end-1))
disp(string12(1:end-1))
disp(string13(1:end-1))
disp(string14(1:end-1))
disp([string15(1:end-1),' kPa'])
disp([string16(1:end-1),' m/yr'])
disp([string17(1:end-1),' x10^-4'])
disp([string18(1:end-1),' x10^-4'])
disp(string19(1:end-1))
disp(string20(1:end-1))
disp([string21(1:end-1),' m/yr'])
disp(string22(1:end-1))
disp(string23(1:end-1))
disp(string24(1:end-1))
disp(string25(1:end-1))
disp(string26(1:end-1))
disp(string27(1:end-1))
disp(string28(1:end-1))
disp(string29(1:end-1))

% Display stats, fast flow:
string1='RMS(ln(C))=';
string2='var(ln(C))=';
string3='Variance ratio (ln(C)/ln(C1))=';
string4='R(N,C1)=';
string5='R2(N,C1)=';
string6='R(ln(N),ln(C1))=';
string7='R2(ln(N),ln(C1))=';
string8='R(ln(C),ln(C1))=';
string9='R2(ln(C),ln(C1))=';
string10='R(tau,tau1)=';
string11='R2(tau,tau1)=';
string12='R(N,tau1)=';
string13='R2(N,tau1)=';
string14='RMS(ln(C)-ln(C1))=';
string15='RMS(tau-tau1)=';
string16='RMS(u_model-u_obs)=';
string17='Jobs=';
string18='Jreg=';
string19='RMS(ln(N))=';
string20='var(ln(N))=';
string21='Jobs, equiv du=';
string22='RMS(ln(tau))=';
string23='var(ln(tau))=';
string24='RMS(ln(Ub))=';
string25='var(ln(Ub))=';
string26='Variance ratio (ln(C)/ln(tau))=';
string27='Variance ratio (ln(C)/ln(u))=';
string28='Variance ratio (Jobs/Jobs(1))=';
string29='Variance ratio total (Jobsratio*lnCratio)=';
for ii=1:numlcurves
    string1=[string1,num2str(.01*round(100*rmslnc(ii,2))),', '];
    string2=[string2,num2str(.01*round(100*varlnc(ii,2))),', '];
    string3=[string3,num2str(.01*round(100*varratiolnclnc(ii,2))),', '];
    string4=[string4,num2str(.01*round(100*rnc(ii,2))),', '];
    string5=[string5,num2str(.01*round(100*r2nc(ii,2))),', '];
    string6=[string6,num2str(.01*round(100*rlnnlnc(ii,2))),', '];
    string7=[string7,num2str(.01*round(100*r2lnnlnc(ii,2))),', '];
    string8=[string8,num2str(.01*round(100*rlnclnc(ii,2))),', '];
    string9=[string9,num2str(.01*round(100*r2lnclnc(ii,2))),', '];
    string10=[string10,num2str(.01*round(100*rtautau(ii,2))),', '];
    string11=[string11,num2str(.01*round(100*r2tautau(ii,2))),', '];
    string12=[string12,num2str(.01*round(100*rntau(ii,2))),', '];
    string13=[string13,num2str(.01*round(100*r2ntau(ii,2))),', '];
    string14=[string14,num2str(.01*round(100*rmslnclnc(ii,2))),', '];
    string15=[string15,num2str(.01*round(100*rmstautau(ii,2)/1000)),', '];
    string16=[string16,num2str(.01*round(100*rmsvelmisfit(ii,2))),', '];
    string17=[string17,num2str(.01*round(100*jobs(ii,2)*10000)),', '];
    string18=[string18,num2str(.01*round(100*jreg(ii,2)*10000)),', '];
    string19=[string19,num2str(.01*round(100*rmslnn(ii,2))),', '];
    string20=[string20,num2str(.01*round(100*varlnn(ii,2))),', '];
    string21=[string21,num2str(.01*round(100*jobsvelequiv(ii,2))),', '];
    string22=[string22,num2str(.01*round(100*rmslntau(ii,2))),', '];
    string23=[string23,num2str(.01*round(100*varlntau(ii,2))),', '];
    string24=[string24,num2str(.01*round(100*rmslnu(ii,2))),', '];
    string25=[string25,num2str(.01*round(100*varlnu(ii,2))),', '];
    string26=[string26,num2str(.01*round(100*varratiolnclntau(ii,2))),', '];
    string27=[string27,num2str(.01*round(100*varratiolnclnu(ii,2))),', '];
    string28=[string28,num2str(.01*round(100*varratiojobsjobs(ii,2))),', '];
    string29=[string29,num2str(.01*round(100*varratiototal(ii,2))),', '];
end
disp('...')
disp('Fast-Flow Stats:')
disp(string1(1:end-1))
disp(string2(1:end-1))
disp(string3(1:end-1))
disp(string4(1:end-1))
disp(string5(1:end-1))
disp(string6(1:end-1))
disp(string7(1:end-1))
disp(string8(1:end-1))
disp(string9(1:end-1))
disp(string10(1:end-1))
disp(string11(1:end-1))
disp(string12(1:end-1))
disp(string13(1:end-1))
disp(string14(1:end-1))
disp([string15(1:end-1),' kPa'])
disp([string16(1:end-1),' m/yr'])
disp([string17(1:end-1),' x10^-4'])
disp([string18(1:end-1),' x10^-4'])
disp(string19(1:end-1))
disp(string20(1:end-1))
disp([string21(1:end-1),' m/yr'])
disp(string22(1:end-1))
disp(string23(1:end-1))
disp(string24(1:end-1))
disp(string25(1:end-1))
disp(string26(1:end-1))
disp(string27(1:end-1))
disp(string28(1:end-1))
disp(string29(1:end-1))
disp('...')

%% MAKE L-Curves FIGURE 
FigH = figure('Position', get(0, 'Screensize'));
tiledlayout(3,3,'TileSpacing',TileSpacing,'Padding',Padding) 
% Plot cost vs cost, log scale:
% Make subplot:
nexttile([3,2])
% pre-allocate handles:
thesehandles=zeros(numlcurves,1);
% loop through L-curves:
for thismodel=1:numlcurves
    % Plot model costs:
    plot(exp(LCurveData(thismodel).LogRegCost_points(~LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogDataCost_points(~LCurveData(thismodel).IsOutlier)),...
        'Marker',modelmarker,'MarkerFaceColor',lcurvecolors{thismodel},'Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',modelmarkersize)
    hold on
    plot(exp(LCurveData(thismodel).LogRegCost_points(LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)),...
        'Marker',modelmarker,'MarkerFaceColor','none','Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',modelmarkersize)
    % Plot smooth costs:
    plot(exp(LCurveData(thismodel).LogRegCost_curve),exp(LCurveData(thismodel).LogDataCost_curve),...
        'Marker','none','Color',lcurvecolors{thismodel},'LineStyle',curvelinestyle,'LineWidth',curvelinesize);
    % Plot corner region:
    thesehandles(thismodel)=plot(exp(LCurveData(thismodel).LogRegCost_curve(LCurveData(thismodel).bestlambdaind_minallowable:LCurveData(thismodel).bestlambdaind_maxallowable)),...
        exp(LCurveData(thismodel).LogDataCost_curve(LCurveData(thismodel).bestlambdaind_minallowable:LCurveData(thismodel).bestlambdaind_maxallowable)),...
        'Marker','none','Color',lcurvecolors{thismodel},'LineStyle',cornerlinestyle,'LineWidth',cornerlinesize);
    % Plot best lambda value:
    plot(exp(LCurveData(thismodel).LogRegCost_curve(LCurveData(thismodel).bestlambdaind)),exp(LCurveData(thismodel).LogDataCost_curve(LCurveData(thismodel).bestlambdaind)),...
        'Marker',bestmarker,'MarkerFaceColor','w','Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',bestmarkersize)
end
% Note outliers:
if m==3
for thismodel=1:numlcurves
    if thismodel==1
        if sum(LCurveData(thismodel).IsOutlier)>0
            numoutliers=sum(LCurveData(thismodel).IsOutlier);
            minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
            minoutlierpow10=floor(log10(minoutlierjobs));
            t=text(jreglims(2)*.9,jobslims(2),[num2str(sum(LCurveData(thismodel).IsOutlier)),' outlier for ',ndescriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'}'],...
		  'HorizontalAlignment','right','VerticalAlignment','top','Color',lcurvecolors{thismodel},'FontSize',legendfontsize3);
            get(t,'Position')
        end
    elseif thismodel==2
        if sum(LCurveData(thismodel).IsOutlier)>0
            numoutliers=sum(LCurveData(thismodel).IsOutlier);
            minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
            minoutlierpow10=floor(log10(minoutlierjobs));
            t=text(jreglims(2)*.9,jobslims(2),[num2str(sum(LCurveData(thismodel).IsOutlier)),' outlier for ',ndescriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'}'],...
                'Position',[0.0013 0.6 0],'Color',lcurvecolors{thismodel},'FontSize',legendfontsize3);
            get(t,'Position')
        end
    elseif thismodel==3
        if sum(LCurveData(thismodel).IsOutlier)>0
            numoutliers=sum(LCurveData(thismodel).IsOutlier);
            minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
            minoutlierpow10=floor(log10(minoutlierjobs));
            t=text(jreglims(2)*.9,jobslims(2),[num2str(sum(LCurveData(thismodel).IsOutlier)),' outliers for ',ndescriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'}'],...
                'Position',[0.09 -10 0],'Color',lcurvecolors{thismodel},'FontSize',legendfontsize3);
        end
    end 
end
else
% Note outliers:
    for thismodel=1:numlcurves
        if thismodel==1
            if sum(LCurveData(thismodel).IsOutlier)>0
                numoutliers=sum(LCurveData(thismodel).IsOutlier);
                minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
                minoutlierpow10=floor(log10(minoutlierjobs));
                t=text(jreglims(2)*.9,jobslims(2),[num2str(sum(LCurveData(thismodel).IsOutlier)),' outlier for ',ndescriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'}'],...
		'HorizontalAlignment','right','VerticalAlignment','top','Color',lcurvecolors{thismodel},'FontSize',legendfontsize3);
                get(t,'Position')
            end
        elseif thismodel==2
            if sum(LCurveData(thismodel).IsOutlier)>0
                numoutliers=sum(LCurveData(thismodel).IsOutlier);
                minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
                minoutlierpow10=floor(log10(minoutlierjobs));
                t=text(jreglims(2)*.9,jobslims(2),[num2str(sum(LCurveData(thismodel).IsOutlier)),' outlier for ',ndescriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'}'],...
                    'Position',[0.00017 0.55 0],'Color',lcurvecolors{thismodel},'FontSize',legendfontsize3);
                get(t,'Position') %[0.06 0.7 0] 0.0045
            end
        elseif thismodel==3
            if sum(LCurveData(thismodel).IsOutlier)>0
                numoutliers=sum(LCurveData(thismodel).IsOutlier);
                minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
                minoutlierpow10=floor(log10(minoutlierjobs));
                t=text(jreglims(2)*.9,jobslims(2),[num2str(sum(LCurveData(thismodel).IsOutlier)),' outliers for ',ndescriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'}'],...
                    'Position',[0.00012 0.35 0],'Color',lcurvecolors{thismodel},'FontSize',legendfontsize3); %[0.09 -10 0]
            end %0.0038 0.00015
        end 
    end
end
% Set up axes etc:
xlim(jreglims.*[.999,1.001])
ylim(jobslims.*[.999,1.001])
set(gca,'XScale','log')
set(gca,'YScale','log')
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize,'LineWidth',axeslinewidth)
set(gca,'PlotBoxAspectRatio',[(log(jreglims(2))-log(jreglims(1)))/(log(jobslims(2))-log(jobslims(1))),1,1])
xlabel('J_{reg} (unitless)')
ylabel('J_{obs} (unitless)')
legend(thesehandles,descriptext,'location','southwest','FontSize',legendfontsize)
title('\rm (a) L-Curve Plots')

% Plot cost curvatures:
for thismodel=1:numlcurves
    % Make subplot:
    nexttile
    % Compute lambda subsampling:
    LogLambda_path=linspace(LCurveData(thismodel).LogLambda_points(find(~LCurveData(thismodel).IsOutlier,1,'first')),...
        LCurveData(thismodel).LogLambda_points(find(~LCurveData(thismodel).IsOutlier,1,'last')),length(LCurveData(thismodel).LogDataCost_curve))';
    % Identify max curvature:
    maxcurve=LCurveData(thismodel).TotalCurvature_curve(LCurveData(thismodel).bestlambdaind-1);
    % Plot best curve curvatures:
    plot(exp(LogLambda_path(2:end-1)),LCurveData(thismodel).LogDataCurvature_curve,'Color',datacurvecolor,'LineWidth',3)
    hold on
    plot(exp(LogLambda_path(2:end-1)),LCurveData(thismodel).LogRegCurvature_curve,'Color',regcurvecolor,'LineWidth',3)
    plot(exp(LogLambda_path(2:end-1)),LCurveData(thismodel).TotalCurvature_curve,'Color',totcurvecolor,'LineWidth',3)
    % Show best value and uncertainty:
    plot(BestLambda(thismodel),maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
    text(BestLambda(thismodel),maxcurve*1.025,{'Best';'value'},'HorizontalAlignment','center','VerticalAlignment','bottom')
    plot(BestLambda(thismodel)*[1/DeltaLambdaRatio(thismodel),DeltaLambdaRatio(thismodel)],maxcurve*[1,1],'k')
    plot((BestLambda(thismodel)/DeltaLambdaRatio(thismodel))*[1,1],maxcurve*[.975,1.025],'k')
    plot(BestLambda(thismodel)*DeltaLambdaRatio(thismodel)*[1,1],maxcurve*[.975,1.025],'k')
    % Show max/min lambda and uncertainties:
    %plot(BestLambda(thismodel)*[1,1],maxcurve*[.5,1],'--k')
    %plot([LCurveData(thismodel).bestlambda_minallowable,LCurveData(thismodel).bestlambda_maxallowable],.5*maxcurve*[1,1],'--k')
    plot(MinLambda(thismodel),.5*maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
    plot(MinLambda(thismodel)*[1/DeltaLambdaRatio(thismodel),DeltaLambdaRatio(thismodel)],.5*maxcurve*[1,1],'k')
    plot((MinLambda(thismodel)/DeltaLambdaRatio(thismodel))*[1,1],.5*maxcurve*[.975,1.025],'k')
    plot(MinLambda(thismodel)*DeltaLambdaRatio(thismodel)*[1,1],.5*maxcurve*[.975,1.025],'k')
    plot(MinLambda(thismodel)*[1,.5],.5*maxcurve+[0,.01],'k')
    text(MinLambda(thismodel)*.5,.5*maxcurve+.01,{'Min';'value'},'HorizontalAlignment','right','VerticalAlignment','bottom','Color','k')
    plot(MaxLambda(thismodel),.5*maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
    plot(MaxLambda(thismodel)*[1/DeltaLambdaRatio(thismodel),DeltaLambdaRatio(thismodel)],.5*maxcurve*[1,1],'k')
    plot((MaxLambda(thismodel)/DeltaLambdaRatio(thismodel))*[1,1],.5*maxcurve*[.975,1.025],'k')
    plot(MaxLambda(thismodel)*DeltaLambdaRatio(thismodel)*[1,1],.5*maxcurve*[.975,1.025],'k')
    plot(MaxLambda(thismodel)*[1,2],.5*maxcurve+[0,.01],'k')
    text(MaxLambda(thismodel)*2,.5*maxcurve+.01,{'Max';'value'},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
    % Determine axes limits:
    thesexlims=10.^[floor(log10(exp(LCurveData(thismodel).LogLambda_points(1)))),ceil(log10(exp(LCurveData(thismodel).LogLambda_points(end))))];
    theseylims=maxcurve*[-.25,1.25];
    % Plot zero line:
    plot(thesexlims,[0,0],'--k')
    % Set up axes etc:
    xlim(thesexlims)
    ylim(theseylims)
    set(gca,'XScale','log')
    set(gca,'FontSize',axesfontsize2,'LineWidth',axeslinewidth)
    set(gca,'TickDir','out')
    if thismodel==3
        set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
        xlabel('\lambda value (unitless)')
    else 
        set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
        set(gca,'XTickLabel',[])
    end 
    %ylabel('$\frac{d^{2}(ln(J))}{d(ln(\lambda))^{2}}$','interpreter','latex')
    ylabel('Cost curvature','FontSize',axesfontsize2)
    title(['\rm (',char(96+1+thismodel),') ',ndescriptext{thismodel}])
    if log(BestLambda(thismodel))>mean([LCurveData(thismodel).LogLambda_points(1),LCurveData(thismodel).LogLambda_points(end)])
        locationstring='NorthWest';
    else
        locationstring='NorthEast';
    end
    legend('J_{obs}','J_{reg}','Total','location',locationstring,'FontSize',legendfontsize2)
end


%% SAVE FIGURE 
saveas(FigH,figname_lcurves,'png');
exportgraphics(FigH,figname_lcurves,'resolution',resolution)

% Final display:
disp('Done!')
toc
