% InversionComparison_v1

% Mike Wolovick, 11/28/2022

% This script makes a figure showing how the best-value regularization
% lambda value converges as a function of mesh resolution

% v2:  Extra subplots, now including:  
% a) all L-curves
% b) lambda values as a function of dx (the old plot)
% c) misfit with best dx at best lambda
% d) correlation with best dx at best lambda

clearvars
tic
close all

%% Parameters:

% File names and paths:
% lcurvefiles={'/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh2.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh3.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh4.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh5.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh6.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh7.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh8.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh9.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh10.mat';};
% localpath='/home/mvolovic/Dropbox/FRISModeling/LCurves/';
% issmscriptpath='~/Dropbox/ISSM/';
% myscriptpath='~/Dropbox/MATLAB/';
% figname='/home/mvolovic/Dropbox/FRISModeling/InversionResConvergenceFigure_v3.png';
lcurvefiles={'/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh2.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh3.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh4.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh5.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh6.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh7.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh8.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh9.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh10.mat';};
localpath='/home/mvolovic/Dropbox/FRISModeling/LCurves/';
issmscriptpath='/home/ollie/mwolovic/ISSM_installations/trunk-jpl/';
myscriptpath='/home/ollie/mwolovic/MATLAB/';
figname='/home/ollie/mwolovic/InversionResConvergenceFigure_v5.png';

% Server control:
%serverloc='albedo0.dmawi.de';  % can be commented
serverusername='mwolovic';
serverpassword='';  % sorry for saving this in plain text.  
% Maybe if your graphical connection wasn't so damn slow I wouldn't need to run matlab from my home PC, AWI.

% Threshold for analyzing fast-flow areas:
velanalysisthreshold=15;               % m/yr (smoothing weighting btwn 1/2 and 2x this value)

% Axes limits and ticks:
jreglims=[1e-5,1e-1];                  % [1x2] unitless
jobslims=[1e-4,1e-1];                  % [1x2] unitless
cmisfitlims=[0,1.75];                   % [1x2] unitless
taumisfitlims=[0,40];                  % [1x2] kPa
r2lims=[.4,1];                         % [1x2] unitless
dxlims=[0,10];                         % [1x2] km

% Other aesthetic settings:
dxcmap='parula';                       % valid colormap string
flipdxcmap=1;                          % logical
modelmarker='o';                       % valid marker string
cornerlinestyle='-';                   % valid linestyle string
wholecurvelinestyle='-';               % valid linestyle string
bestlambdamarker='d';                  % valid marker string
maxlambdamarker='s';                   % valid marker string
minlambdamarker='^';                   % valid marker string
bestlambdacolor='k';                   % valid color string or RGB
minlambdacolor='r';                    % valid color string or RGB
maxlambdacolor='b';                    % valid color string or RGB
lnccolor='b';                          % valid color string or RGB
taubcolor='r';                         % valid color string or RGB
groundedlinestyle='--';                % valid linestyle string
fastlinestyle='-';                     % valid linestyle string
axesfontsize=13;                       % points
titlefontsize=14;                          % points
colorbarfontsize=13;                   % points

% Subplot controls:
vertaxesbuffer=.05;                    % unitless
titlebuffer=.04;                       % unitless
horzaxesbuffer=.08;                    % unitless
horznotextbuffer=.01;                  % unitless

% Output size and resolution:
pagesize=[18,12];                      % [1x2] inches 
resolution=300;                        % dpi


%% Preparation:

% Compute subplot dimensions:
% Hardcoded subplot arrangement:
numvertplots=2;
numhorzplots=2;
% Compute available space:
availablehorzspace=1-(numhorzplots+1)*(horzaxesbuffer);
availablevertspace=1-numvertplots*(titlebuffer+vertaxesbuffer);
% Compute boxes:
Boxes=cell(numvertplots,numhorzplots);
for ii=1:numvertplots
    for jj=1:numhorzplots
        Boxes{numvertplots-ii+1,jj}=[horzaxesbuffer+(jj-1)*(horzaxesbuffer+availablehorzspace/numhorzplots),...
            vertaxesbuffer+(ii-1)*(vertaxesbuffer+titlebuffer+availablevertspace/numvertplots),...
            availablehorzspace/numhorzplots,availablevertspace/numvertplots];
    end
end

% Add the ISSM script path:
addpath(genpath(issmscriptpath))

% Compute number of L-curves:
numlcurves=length(lcurvefiles);

% Pre-allocate:
BestLambda=zeros(numlcurves,1);
MaxLambda=zeros(numlcurves,1);
MinLambda=zeros(numlcurves,1);
DeltaLambdaRatio=zeros(numlcurves,1);
DragCoeffScale=zeros(numlcurves,1);
VelScale=zeros(numlcurves,1);
MeshSize=zeros(numlcurves,1);
GroundedRMSlnc=zeros(numlcurves,1);
FastRMSlnc=zeros(numlcurves,1);
GroundedRMStaub=zeros(numlcurves,1);
FastRMStaub=zeros(numlcurves,1);
GroundedR2lnc=ones(numlcurves,1);
FastR2lnc=ones(numlcurves,1);
GroundedR2taub=ones(numlcurves,1);
FastR2taub=ones(numlcurves,1);
LCurveData=struct('LogRegCost_points',[],'LogDataCost_points',[],'LogRegCost_curve',[],'LogDataCost_curve',[],...
    'bestlambdaind_maxallowable',[],'bestlambdaind_minallowable',[],'bestlambdaind',[]);

% Loop through model files:
for ii=1:numlcurves

    % Check if we're loading remote files:
    if exist('serverloc','var')
        % Parse filename:
        slashnum=strfind(lcurvefiles{ii},'/');
        remotepath=lcurvefiles{ii}(1:slashnum(end));
        lcurvename=lcurvefiles{ii}(slashnum(end)+1:end);
        % Fetch the basic 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],'BasicModel','Log*Cost_bestcurve','Lambda','DataCost','RegularizationCost','bestlambda*','bestmodelind','MyInversionParameters')
    else
        % Load the file:
        load(lcurvefiles{ii},'BasicModel','Log*Cost_bestcurve','Lambda','DataCost','RegularizationCost','bestlambda*','bestmodelind','MyInversionParameters')
    end

    % Construct the file name for the best model:
    slashnum=strfind(lcurvefiles{ii},'/');
    thisbase=10^floor(log10(Lambda(bestmodelind)));
    thisname=[lcurvefiles{ii}(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 best model file:
        if ~exist([localpath,thisname],'file')
            ssh2_struct=scp_simple_get(serverloc,serverusername,serverpassword,thisname,localpath,remotepath);
        end
        % Load the best model:
        load([localpath,thisname],'md')
    else
        % Load the best model:
        load([lcurvefiles{ii}(1:slashnum(end)),thisname],'md')
    end

    
    % Assign simple things:
    BestLambda(ii)=bestlambda;
    MaxLambda(ii)=bestlambda_maxallowable;
    MinLambda(ii)=bestlambda_minallowable;
    DeltaLambdaRatio(ii)=exp(bestlambda_loguncertainty);
    DragCoeffScale(ii)=MyInversionParameters.dragcoeffstdev;
    VelScale(ii)=MyInversionParameters.velrms;

    % Assign lcurve data:
    LCurveData(ii).LogRegCost_points=log(RegularizationCost);
    LCurveData(ii).LogDataCost_points=log(DataCost);
    LCurveData(ii).LogRegCost_curve=LogRegCost_bestcurve;
    LCurveData(ii).LogDataCost_curve=LogDataCost_bestcurve;
    LCurveData(ii).bestlambdaind=bestlambdaind;
    LCurveData(ii).bestlambdaind_maxallowable=bestlambdaind_maxallowable;
    LCurveData(ii).bestlambdaind_minallowable=bestlambdaind_minallowable;
    
    % Get 2D mesh:
    try
        x2d=BasicModel.mesh.x2d;
        y2d=BasicModel.mesh.y2d;
        elements2d=BasicModel.mesh.elements2d;
    catch
        x2d=BasicModel.mesh.x;
        y2d=BasicModel.mesh.y;
        elements2d=BasicModel.mesh.elements;
    end
    numelements2d=size(elements2d,1);
    numvertices2d=length(x2d);
    
    % Compute mesh size metrics:
    % 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;
    
    % Extract masks on 2D mesh and truncate:
    OceanLevelSet_vertices=max(-1,min(1,BasicModel.mask.ocean_levelset(1:numvertices2d)));
    IceLevelSet_vertices=max(-1,min(1,BasicModel.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 grounded ice fraction:
    GroundedIceAreaFraction_elements=IceAreaFraction_elements.*(1-OceanAreaFraction_elements);
    
    % Compute fast-flow weighting function:
    % Extract velocity observations:
    VelMag_vertices=BasicModel.inversion.vel_obs(end-numvertices2d+1:end);
    % Compute weighting function for fast-flowing areas:
    FastFlowWeight_vertices=ones(numvertices2d,1);
    FastFlowWeight_vertices(VelMag_vertices<.5*velanalysisthreshold)=0;
    FastFlowWeight_vertices(VelMag_vertices>=.5*velanalysisthreshold&VelMag_vertices<2*velanalysisthreshold)=...
        .5*(1+sin(.5*pi*log2(VelMag_vertices(VelMag_vertices>=.5*velanalysisthreshold&VelMag_vertices<2*velanalysisthreshold)/velanalysisthreshold)));
    % Average weighting function to the elements:
    FastFlowWeight_elements=mean(FastFlowWeight_vertices(elements2d),2);
    
    % Compute median fast-flow grid size:
    % sort by element size:
    [SortedSize,SortedInd]=sort(HydraulicDiameter_elements(FastFlowWeight_elements>1e-3&GroundedIceAreaFraction_elements>1e-3));
    % generate cumulative sorted area:
    ThisArea=Area(FastFlowWeight_elements>1e-3&GroundedIceAreaFraction_elements>1e-3).*...
        FastFlowWeight_elements(FastFlowWeight_elements>1e-3&GroundedIceAreaFraction_elements>1e-3).*GroundedIceAreaFraction_elements(FastFlowWeight_elements>1e-3&GroundedIceAreaFraction_elements>1e-3);
    CumulativeSortedArea=cumsum(ThisArea(SortedInd));
    % Interpolate 50th percentile:
    MeshSize(ii)=interp1(CumulativeSortedArea/CumulativeSortedArea(end),SortedSize,.5); 


    % Unpack parameters:
    unpack(MyInversionParameters)

    % Extract variables we need:
    if numlayers>1 && ~strcmp(flowequation,'SSA')
        ModelSurfVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonsurface==1);
        ModelBasalVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonbase==1);
        DragCoeff_vertices=md.results.StressbalanceSolution.FrictionCoefficient(md.mesh.vertexonbase==1);
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices=md.friction.effective_pressure(md.mesh.vertexonbase==1);
        end
    else
        ModelBasalVel=md.results.StressbalanceSolution.Vel;
        ModelSurfVel=ModelBasalVel;
        DragCoeff_vertices=md.results.StressbalanceSolution.FrictionCoefficient;
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices=md.friction.effective_pressure;
        end
    end

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

    % Check if this is the first model:
    if ii==1
        % Record reference model:
        x2d_ref=x2d;
        y2d_ref=y2d;
        elements2d_ref=elements2d;
        DragCoeff_ref=mean(2*log(max(mindragcoeff,DragCoeff_vertices(elements2d))),2);  % this variable represents ln(C) on elements
        Drag_ref=mean(Drag_vertices(elements2d),2); % interpolate to elements
        FastFlowWeight_ref=FastFlowWeight_elements;
        GroundedIceAreaFraction_ref=GroundedIceAreaFraction_elements;
        Area_ref=Area;
        % Compute area totals:
        groundedicearea=sum(Area.*GroundedIceAreaFraction_elements);
        fasticearea=sum(Area.*GroundedIceAreaFraction_elements.*FastFlowWeight_elements);
        % Compute means:
        refmeangroundedc=sum(Area.*GroundedIceAreaFraction_elements.*DragCoeff_ref)/groundedicearea;
        refmeanfastc=sum(Area.*GroundedIceAreaFraction_elements.*FastFlowWeight_elements.*DragCoeff_ref)/fasticearea;
        refmeangroundedtau=sum(Area.*GroundedIceAreaFraction_elements.*Drag_ref)/groundedicearea;
        refmeanfasttau=sum(Area.*GroundedIceAreaFraction_elements.*FastFlowWeight_elements.*Drag_ref)/fasticearea;
        % Compute standard deviations:
        refstdgroundedc=sqrt(sum(Area.*GroundedIceAreaFraction_elements.*((DragCoeff_ref-refmeangroundedc).^2))/groundedicearea);
        refstdfastc=sqrt(sum(Area.*GroundedIceAreaFraction_elements.*FastFlowWeight_elements.*((DragCoeff_ref-refmeanfastc).^2))/fasticearea);
        refstdgroundedtau=sqrt(sum(Area.*GroundedIceAreaFraction_elements.*((Drag_ref-refmeangroundedtau).^2))/groundedicearea);
        refstdfasttau=sqrt(sum(Area.*GroundedIceAreaFraction_elements.*FastFlowWeight_elements.*((Drag_ref-refmeanfasttau).^2))/fasticearea);
    else
        % Interpolate to reference mesh vertices:
        ThisDragCoeff=InterpFromMeshToMesh2d(elements2d,x2d,y2d,2*log(max(mindragcoeff,DragCoeff_vertices)),x2d_ref,y2d_ref);
        ThisDrag=InterpFromMeshToMesh2d(elements2d,x2d,y2d,Drag_vertices,x2d_ref,y2d_ref);
        % Interpolate to elements:
        ThisDragCoeff_elements=mean(ThisDragCoeff(elements2d_ref),2);
        ThisDrag_elements=mean(ThisDrag(elements2d_ref),2);
        % Compute means:
        thismeangroundedc=sum(Area_ref.*GroundedIceAreaFraction_ref.*ThisDragCoeff_elements)/groundedicearea;
        thismeanfastc=sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*ThisDragCoeff_elements)/fasticearea;
        thismeangroundedtau=sum(Area_ref.*GroundedIceAreaFraction_ref.*ThisDrag_elements)/groundedicearea;
        thismeanfasttau=sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*ThisDrag_elements)/fasticearea;
        % Compute standard deviations:
        thisstdgroundedc=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*((ThisDragCoeff_elements-thismeangroundedc).^2))/groundedicearea);
        thisstdfastc=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*((ThisDragCoeff_elements-thismeanfastc).^2))/fasticearea);
        thisstdgroundedtau=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*((ThisDrag_elements-thismeangroundedtau).^2))/groundedicearea);
        thisstdfasttau=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*((ThisDrag_elements-thismeanfasttau).^2))/fasticearea);
        % Compute RMS misfits:
        GroundedRMSlnc(ii)=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*((ThisDragCoeff_elements-DragCoeff_ref).^2))/groundedicearea);
        FastRMSlnc(ii)=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*((ThisDragCoeff_elements-DragCoeff_ref).^2))/fasticearea);
        GroundedRMStaub(ii)=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*((ThisDrag_elements-Drag_ref).^2))/groundedicearea);
        FastRMStaub(ii)=sqrt(sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*((ThisDrag_elements-Drag_ref).^2))/fasticearea);
        % Compute R2 correlations:
        GroundedR2lnc(ii)=(sum(Area_ref.*GroundedIceAreaFraction_ref.*(ThisDragCoeff_elements-thismeangroundedc).*(DragCoeff_ref-refmeangroundedc))/(groundedicearea*thisstdgroundedc*refstdgroundedc))^2;
        FastR2lnc(ii)=(sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*(ThisDragCoeff_elements-thismeanfastc).*(DragCoeff_ref-refmeanfastc))/(fasticearea*thisstdfastc*refstdfastc))^2;
        GroundedR2taub(ii)=(sum(Area_ref.*GroundedIceAreaFraction_ref.*(ThisDrag_elements-thismeangroundedtau).*(Drag_ref-refmeangroundedtau))/(groundedicearea*thisstdgroundedtau*refstdgroundedtau))^2;
        FastR2taub(ii)=(sum(Area_ref.*GroundedIceAreaFraction_ref.*FastFlowWeight_ref.*(ThisDrag_elements-thismeanfasttau).*(Drag_ref-refmeanfasttau))/(fasticearea*thisstdfasttau*refstdfasttau))^2;
    end
    
end

% Correct the best lambda values for changing scales:
BestLambda=BestLambda.*((DragCoeffScale(1)./DragCoeffScale).^2).*((VelScale./VelScale(1)).^2);
MaxLambda=MaxLambda.*((DragCoeffScale(1)./DragCoeffScale).^2).*((VelScale./VelScale(1)).^2);
MinLambda=MinLambda.*((DragCoeffScale(1)./DragCoeffScale).^2).*((VelScale./VelScale(1)).^2);

% Compute best-fits:
A=[ones(numlcurves,1),MeshSize/1e3];
W=eye(numlcurves).*repmat(1./(log(DeltaLambdaRatio).^2),[1,numlcurves]);
equation_best=(A'*W*A)\(A'*W*log(BestLambda));
equation_max=(A'*W*A)\(A'*W*log(MaxLambda));
equation_min=(A'*W*A)\(A'*W*log(MinLambda));


%% Figure:
close
% Make figure:
figure(1)

% Make a cmap:
thiscmap=colormap(dxcmap);
if flipdxcmap
    thiscmap=flipud(thiscmap);
end

% Make all L-curves:
% make subplot:
subplot('Position',Boxes{1,1})
% loop through models:
for ii=numlcurves:-1:1
    % get this color:
    thiscolor=interp1(linspace(0,1,size(thiscmap,1))',thiscmap,max(0,min(1,(MeshSize(ii)/1000-dxlims(1))/(dxlims(2)-dxlims(1)))));
    % plot smooth curve (whole thing):
    plot((DragCoeffScale(ii)/DragCoeffScale(1))*exp(LCurveData(ii).LogRegCost_curve),(VelScale(ii)/VelScale(1))*exp(LCurveData(ii).LogDataCost_curve),...
        'LineStyle',wholecurvelinestyle,'Color',thiscolor,'LineWidth',2)
    hold on
    % plot smooth curve (corner region):
    plot((DragCoeffScale(ii)/DragCoeffScale(1))*exp(LCurveData(ii).LogRegCost_curve(LCurveData(ii).bestlambdaind_minallowable:LCurveData(ii).bestlambdaind_maxallowable)),...
        (VelScale(ii)/VelScale(1))*exp(LCurveData(ii).LogDataCost_curve(LCurveData(ii).bestlambdaind_minallowable:LCurveData(ii).bestlambdaind_maxallowable)),...
        'LineStyle',cornerlinestyle,'Color',thiscolor,'LineWidth',2)
    % plot model points:
%     plot((DragCoeffScale(1)/DragCoeffScale(ii))*exp(LCurveData(ii).LogRegCost_points),(VelScale(1)/VelScale(ii))*exp(LCurveData(ii).LogDataCost_points),...
%         'Marker',modelmarker,'Color',thiscolor,'MarkerFaceColor',thiscolor)
    % Plot best/min/max lambda:
%     plot((DragCoeffScale(1)/DragCoeffScale(ii))*exp(LCurveData(ii).LogRegCost_curve(LCurveData(ii).bestlambdaind_minallowable)),...
%         (VelScale(1)/VelScale(ii))*exp(LCurveData(ii).LogDataCost_curve(LCurveData(ii).bestlambdaind_minallowable)),...
%         'LineStyle','none','Color',thiscolor,'Marker',minlambdamarker,'MarkerFaceColor',thiscolor)
    plot((DragCoeffScale(ii)/DragCoeffScale(1))*exp(LCurveData(ii).LogRegCost_curve(LCurveData(ii).bestlambdaind)),...
        (VelScale(ii)/VelScale(1))*exp(LCurveData(ii).LogDataCost_curve(LCurveData(ii).bestlambdaind)),...
        'LineStyle','none','Color',thiscolor,'Marker',bestlambdamarker,'MarkerFaceColor',thiscolor)
%     plot((DragCoeffScale(1)/DragCoeffScale(ii))*exp(LCurveData(ii).LogRegCost_curve(LCurveData(ii).bestlambdaind_maxallowable)),...
%         (VelScale(1)/VelScale(ii))*exp(LCurveData(ii).LogDataCost_curve(LCurveData(ii).bestlambdaind_maxallowable)),...
%         'LineStyle','none','Color',thiscolor,'Marker',maxlambdamarker,'MarkerFaceColor',thiscolor)
end
% Set up axes etc:
caxis(dxlims)
colormap(thiscmap)
hc=colorbar;
set(get(hc,'Ylabel'),'String','Median mesh size (km)','FontSize',colorbarfontsize)
set(hc,'FontSize',colorbarfontsize)
set(gca,'XScale','log')
set(gca,'YScale','log')
xlim(jreglims)
ylim(jobslims)
set(gca,'PlotBoxAspectRatio',[(log(jreglims(2))-log(jreglims(1)))/(log(jobslims(2))-log(jobslims(1))),1,1])
xlabel('Regularization cost J_{reg} (unitless)','FontSize',axesfontsize)
ylabel('Data cost J_{obs} (unitless)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
title('a) All L-curves','FontSize',titlefontsize)


% Plot best-fits:
% make subplot:
subplot('Position',Boxes{1,2})
% Plot data points with error bars:
errorbar(MeshSize/1e3,BestLambda,BestLambda.*(1-1./DeltaLambdaRatio),BestLambda.*(DeltaLambdaRatio-1),'LineStyle','none','Marker','s','Color','k','MarkerFaceColor','k')
hold on
errorbar(MeshSize/1e3,MinLambda,MinLambda.*(1-1./DeltaLambdaRatio),MinLambda.*(DeltaLambdaRatio-1),'LineStyle','none','Marker','o','Color','r','MarkerFaceColor','r')
errorbar(MeshSize/1e3,MaxLambda,MaxLambda.*(1-1./DeltaLambdaRatio),MaxLambda.*(DeltaLambdaRatio-1),'LineStyle','none','Marker','d','Color','b','MarkerFaceColor','b')
% Plot best-fits:
DummyX=linspace(0,max(MeshSize/1e3),10);
plot(DummyX,exp(equation_best(1)+equation_best(2)*DummyX),'LineWidth',2,'Color','k')
plot(DummyX,exp(equation_min(1)+equation_min(2)*DummyX),'LineWidth',2,'Color','r')
plot(DummyX,exp(equation_max(1)+equation_max(2)*DummyX),'LineWidth',2,'Color','b')
% Label best-fits:
text(MeshSize(9)/1e3,.8*MaxLambda(9)/DeltaLambdaRatio(9),{'Best-fit to best values:';['\lambda=',num2str(.01*round(100*exp(equation_best(1)))),'*exp({\Delta}x/',num2str(.01*round(100*1/equation_best(2))),' km)']},...
    'HorizontalAlignment','center','VerticalAlignment','top','FontSize',axesfontsize)
text(MeshSize(9)/1e3,.5*MinLambda(9)/DeltaLambdaRatio(9),{'Best-fit to min values:';['\lambda=',num2str(.01*round(100*exp(equation_min(1)))),'*exp({\Delta}x/',num2str(.01*round(100*1/equation_min(2))),' km)']},...
    'HorizontalAlignment','center','VerticalAlignment','top','FontSize',axesfontsize)
text(MeshSize(3)/1e3+.25,1.5*MaxLambda(3)*DeltaLambdaRatio(3),{'Best-fit to max values:';['\lambda=',num2str(.01*round(100*exp(equation_max(1)))),'*exp({\Delta}x/',num2str(.01*round(100*1/equation_max(2))),' km)']},...
    'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',axesfontsize)
% Set up everything else:
xlim(dxlims)
set(gca,'YScale','log')
ylim([min(lambdamin,.9*min(MinLambda./DeltaLambdaRatio)),max(lambdamax,1.1*max(MaxLambda.*DeltaLambdaRatio))])
xlabel('Median mesh size (km)','FontSize',axesfontsize)
ylabel('Regularization \lambda (unitless)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
title('b) Mesh size dependence of regularization','FontSize',titlefontsize)
legend('Best \lambda values','Min \lambda values','Max \lambda values','location','SouthEast','FontSize',axesfontsize)


% Plot RMS misfits:
% make subplot:
subplot('Position',Boxes{2,1})
% make first y-axis:
yyaxis left
% plot RMS(ln(c))
plot(MeshSize(2:end)/1e3,GroundedRMSlnc(2:end),'Color',lnccolor,'LineStyle',groundedlinestyle)
hold on
plot(MeshSize(2:end)/1e3,FastRMSlnc(2:end),'Color',lnccolor,'LineStyle',fastlinestyle)
% set up axes etc:
xlim(dxlims)
ylim(cmisfitlims)
set(gca,'YColor',lnccolor)
set(gca,'XColor','k')
xlabel('Median mesh size (km)','FontSize',axesfontsize)
ylabel('RMS(ln(C)) (unitless)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
title('c) RMS misfit w.r.t. Mesh #1 at \lambda_{best}','FontSize',titlefontsize)
% make second y-axis:
yyaxis right
% plot RMS(taub)
plot(MeshSize(2:end)/1e3,GroundedRMStaub(2:end)/1e3,'Color',taubcolor,'LineStyle',groundedlinestyle)
hold on
plot(MeshSize(2:end)/1e3,FastRMStaub(2:end)/1e3,'Color',taubcolor,'LineStyle',fastlinestyle)
% set up axes etc:
xlim(dxlims)
ylim(taumisfitlims)
set(gca,'YColor',taubcolor)
set(gca,'XColor','k')
xlabel('Median mesh size (km)','FontSize',axesfontsize)
ylabel('RMS(\tau_b) (kPa)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
% Make legend:
h1=plot(-9999*[1,1],-9999*[1,1],'Color','k','Marker','none','LineStyle',groundedlinestyle);
h2=plot(-9999*[1,1],-9999*[1,1],'Color','k','Marker','none','LineStyle',fastlinestyle);
legend([h1;h2],'Whole grounded domain','Grounded fast-flow','location','SouthEast','FontSize',axesfontsize)


% Plot correlations:
% make subplot:
subplot('Position',Boxes{2,2})
% plot R2(ln(c))
plot(MeshSize(2:end)/1e3,GroundedR2lnc(2:end),'Color',lnccolor,'LineStyle',groundedlinestyle)
hold on
plot(MeshSize(2:end)/1e3,FastR2lnc(2:end),'Color',lnccolor,'LineStyle',fastlinestyle)
% plot R2(taub)
plot(MeshSize(2:end)/1e3,GroundedR2taub(2:end),'Color',taubcolor,'LineStyle',groundedlinestyle)
plot(MeshSize(2:end)/1e3,FastR2taub(2:end),'Color',taubcolor,'LineStyle',fastlinestyle)
% set up axes etc:
xlim(dxlims)
ylim(r2lims)
xlabel('Median mesh size (km)','FontSize',axesfontsize)
ylabel('R^2 (unitless)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
title('d) Correlation coefficient w.r.t. Mesh #1 at \lambda_{best}','FontSize',titlefontsize)
% Make legend:
legend('Grounded domain, ln(C)','Fast-flow, ln(C)','Grounded domain, \tau_b','Fast-flow, \tau_b','location','SouthWest','FontSize',axesfontsize)


% Export figure:
set(gcf,'PaperSize',pagesize)
set(gcf,'PaperPosition',[0,0,pagesize])
print('-dpng',figname,['-r',num2str(resolution)])

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