% InversionResMisfitFigure_v1

% Mike Wolovick, 11/29/2022

% This script makes a figure showing how drag coefficient differences
% change as a function of resolution and regularization.  All runs are
% compared with the highest-resolution mesh at the same lambda value.  

% Comparison metrics:  RMS(log(C)), R2(log(C)), RMS(drag), R2(drag)


clearvars
tic
close all

%% Parameters:

% File names and paths:
BasicModelFiles={'/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';};
figname='/home/ollie/mwolovic/InversionResMisfitFigure_v5.png';

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

% Font sizes:
axesfontsize=13;                       % points
titlefontsize=14;                          % points
colorbarfontsize=13;                   % points

% Subplot positioning controls:
vertaxesbuffer=.06;            % unitless
titlebuffer=.05;               % unitless
horzaxesbuffer=.05;            % unitless
horzcolorbarbuffer=.05;       % unitless
horznotextbuffer=.00;          % unitless

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


%% Preparation:

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

% Pre-allocate:
BestLambda=zeros(1,numlcurves);
MaxLambda=zeros(1,numlcurves);
MinLambda=zeros(1,numlcurves);
DeltaLambdaRatio=zeros(1,numlcurves);
DragCoeffScale=zeros(1,numlcurves);
VelScale=zeros(1,numlcurves);
MeshSize_fastflow=zeros(1,numlcurves);
MeshSize_grounded=zeros(1,numlcurves);

% Communicate:
disp('Looping through L-curves')

% Loop through L-Curves:
for ii=1:numlcurves
    
    % Communicate:
    disp('...')
    disp(['... L-curve number ',num2str(ii),'/',num2str(numlcurves)])
    
    %% Load file:
    load(BasicModelFiles{ii},'BasicModel','MyInversionParameters','*best*','Lambda','hassolved')
    
    % Check hassolved:
    if hassolved~=1
        error('L-curve was not solved')
    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;
    
    % 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 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);
    % Scale weighting function for grounded and ice fractions:
    FastFlowWeight_elements=FastFlowWeight_elements.*IceAreaFraction_elements.*(1-OceanAreaFraction_elements);
    
    %% Compute median fast-flow grid size:
    % sort by element size:
    [SortedSize,SortedInd]=sort(HydraulicDiameter_elements(FastFlowWeight_elements>1e-3&IceAreaFraction_elements>1e-3&OceanAreaFraction_elements<1-1e-3));
    % generate cumulative sorted area:
    ThisArea=Area(FastFlowWeight_elements>1e-3&IceAreaFraction_elements>1e-3&OceanAreaFraction_elements<1-1e-3).*...
        FastFlowWeight_elements(FastFlowWeight_elements>1e-3&IceAreaFraction_elements>1e-3&OceanAreaFraction_elements<1-1e-3);
    CumulativeSortedArea=cumsum(ThisArea(SortedInd));
    % Interpolate 50th percentile:
    MeshSize_fastflow(ii)=interp1(CumulativeSortedArea/CumulativeSortedArea(end),SortedSize,.5); 
    
    % Compute median grounded grid size:
    % sort by element size:
    [SortedSize,SortedInd]=sort(HydraulicDiameter_elements(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5));
    % generate cumulative sorted area:
    ThisArea=Area(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5);
    CumulativeSortedArea=cumsum(ThisArea(SortedInd));
    % Interpolate 50th percentile:
    MeshSize_grounded(ii)=interp1(CumulativeSortedArea/CumulativeSortedArea(end),SortedSize,.5); 
    
    % Check if this is the first mesh:
    if ii==1
        % Record lambda sampling:
        Lambda_all=Lambda;
        numlambda=length(Lambda_all);
        % Pre-allocate things:
        RMS_C_fastflow=zeros(numlambda,numlcurves);
        RMS_C_grounded=zeros(numlambda,numlcurves);
        R_C_fastflow=ones(numlambda,numlcurves);
        R_C_grounded=ones(numlambda,numlcurves);
        RMS_Drag_fastflow=zeros(numlambda,numlcurves);
        RMS_Drag_grounded=zeros(numlambda,numlcurves);
        R_Drag_fastflow=ones(numlambda,numlcurves);
        R_Drag_grounded=ones(numlambda,numlcurves);
        LogC_ref=zeros(numelements2d,numlambda);
        Drag_ref=zeros(numelements2d,numlambda);
        % Record first mesh values:
        x2d_ref=x2d;
        y2d_ref=y2d;
        elements2d_ref=elements2d;
        numvertices2d_ref=numvertices2d;
        numelements2d_ref=numelements2d;
        FastFlowWeight_ref_elements=FastFlowWeight_elements;
        IceAreaFraction_ref_elements=IceAreaFraction_elements;
        OceanAreaFraction_ref_elements=OceanAreaFraction_elements;
        Area_ref=Area;
    else
        % Check lambda sampling:
        if ~isequal(Lambda,Lambda_all)
            error('Lambda sampling for this mesh did not match reference lambda sampling')
        end
    end
    
    % Compute area totals:
    if ii==1
        fastflowtotarea=sum(Area_ref.*FastFlowWeight_ref_elements);
        groundedtotarea=sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements);
    end
    
    % Communicate:
    disp('...looping through lambda samples')
    
    % Loop through lambda values:
    for jj=1:numlambda
        
        % Communicate:
        disp(['... ...lambda sample number ',num2str(jj),'/',num2str(numlambda)])
        
        % Construct this file name:
        slashnum=strfind(BasicModelFiles{ii},'/');
        thisbase=10^floor(log10(Lambda(jj)));
        thisname=[BasicModelFiles{ii}(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(jj)/thisbase))];
        
        % Load this model:
        load([BasicModelFiles{ii}(1:slashnum(end)),thisname,'.mat'],'md','hassolved')
        
        % Check hassolved:
        if hassolved~=1
            error('This model was not solved')
        end
        
        % Get this drag coefficient:
        ThisLogC_vertices=log(max(MyInversionParameters.mindragcoeff,md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d)).^2);
        
        % Compute this basal drag:
        ThisDrag_vertices=(md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d).^2).*((md.results.StressbalanceSolution.Vel(1:numvertices2d)/MyInversionParameters.secondsperyear).^(1/MyInversionParameters.m));
        if ~strcmp(MyInversionParameters.nsource,'none')
            ThisDrag_vertices=ThisDrag_vertices.*(md.friction.effective_pressure.^MyInversionParameters.r);
        end
        
        % Check if this is the first mesh:
        if ii==1
            % Average to elements for reference drag values:
            LogC_ref(:,jj)=mean(ThisLogC_vertices(elements2d),2);
            Drag_ref(:,jj)=mean(ThisDrag_vertices(elements2d),2);
        else
            % Interpolate onto original mesh:
            ThisLogC_origvertices=InterpFromMeshToMesh2d(elements2d,x2d,y2d,ThisLogC_vertices,x2d_ref,y2d_ref);
            ThisDrag_origvertices=InterpFromMeshToMesh2d(elements2d,x2d,y2d,ThisDrag_vertices,x2d_ref,y2d_ref);
            % Average onto elements:
            ThisLogC_origelements=mean(ThisLogC_origvertices(elements2d_ref),2);
            ThisDrag_origelements=mean(ThisDrag_origvertices(elements2d_ref),2);
            % Compute RMS misfits:
            RMS_C_fastflow(jj,ii)=sqrt(sum(Area_ref.*FastFlowWeight_ref_elements.*((ThisLogC_origelements-LogC_ref(:,jj)).^2))/fastflowtotarea);
            RMS_C_grounded(jj,ii)=sqrt(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*((ThisLogC_origelements-LogC_ref(:,jj)).^2))/groundedtotarea);
            RMS_Drag_fastflow(jj,ii)=sqrt(sum(Area_ref.*FastFlowWeight_ref_elements.*((ThisDrag_origelements-Drag_ref(:,jj)).^2))/fastflowtotarea);
            RMS_Drag_grounded(jj,ii)=sqrt(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*((ThisDrag_origelements-Drag_ref(:,jj)).^2))/groundedtotarea);
            % Compute correlation coefficients:
            % compute individual means:
            thisfastflowmeanc=sum(Area_ref.*FastFlowWeight_ref_elements.*ThisLogC_origelements)/fastflowtotarea;
            thisgroundedmeanc=sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*ThisLogC_origelements)/groundedtotarea;
            thisfastflowrefmeanc=sum(Area_ref.*FastFlowWeight_ref_elements.*LogC_ref(:,jj))/fastflowtotarea;
            thisgroundedrefmeanc=sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*LogC_ref(:,jj))/groundedtotarea;
            thisfastflowmeandrag=sum(Area_ref.*FastFlowWeight_ref_elements.*ThisDrag_origelements)/fastflowtotarea;
            thisgroundedmeandrag=sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*ThisDrag_origelements)/groundedtotarea;
            thisfastflowrefmeandrag=sum(Area_ref.*FastFlowWeight_ref_elements.*Drag_ref(:,jj))/fastflowtotarea;
            thisgroundedrefmeandrag=sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*Drag_ref(:,jj))/groundedtotarea;
            % compute R values:
            R_C_fastflow(jj,ii)=(sum(Area_ref.*FastFlowWeight_ref_elements.*(ThisLogC_origelements-thisfastflowmeanc).*(LogC_ref(:,jj)-thisfastflowrefmeanc))/fastflowtotarea)/...
                (sqrt(sum(Area_ref.*FastFlowWeight_ref_elements.*((ThisLogC_origelements-thisfastflowmeanc).^2))/fastflowtotarea)*sqrt(sum(Area_ref.*FastFlowWeight_ref_elements.*((LogC_ref(:,jj)-thisfastflowrefmeanc).^2))/fastflowtotarea));
            R_C_grounded(jj,ii)=(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*(ThisLogC_origelements-thisgroundedmeanc).*(LogC_ref(:,jj)-thisgroundedrefmeanc))/groundedtotarea)/...
                (sqrt(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*((ThisLogC_origelements-thisgroundedmeanc).^2))/groundedtotarea)*sqrt(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*((LogC_ref(:,jj)-thisgroundedrefmeanc).^2))/groundedtotarea));
            R_Drag_fastflow(jj,ii)=(sum(Area_ref.*FastFlowWeight_ref_elements.*(ThisDrag_origelements-thisfastflowmeandrag).*(Drag_ref(:,jj)-thisfastflowrefmeandrag))/fastflowtotarea)/...
                (sqrt(sum(Area_ref.*FastFlowWeight_ref_elements.*((ThisDrag_origelements-thisfastflowmeandrag).^2))/fastflowtotarea)*sqrt(sum(Area_ref.*FastFlowWeight_ref_elements.*((Drag_ref(:,jj)-thisfastflowrefmeandrag).^2))/fastflowtotarea));
            R_Drag_grounded(jj,ii)=(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*(ThisDrag_origelements-thisgroundedmeandrag).*(Drag_ref(:,jj)-thisgroundedrefmeandrag))/groundedtotarea)/...
                (sqrt(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*((ThisDrag_origelements-thisgroundedmeandrag).^2))/groundedtotarea)*sqrt(sum(Area_ref.*IceAreaFraction_ref_elements.*OceanAreaFraction_ref_elements.*((Drag_ref(:,jj)-thisgroundedrefmeandrag).^2))/groundedtotarea));
        end
        
    end
    
end

% Communicate:
disp('...')
disp('Done with looping through models!')

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

% Compute best-fits:
A=[ones(numlcurves,1),MeshSize_fastflow'/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'));

% Compute squared correlation coefficients:
R2_C_fastflow=R_C_fastflow.^2;
R2_C_grounded=R_C_grounded.^2;
R2_Drag_fastflow=R_Drag_fastflow.^2;
R2_Drag_grounded=R_Drag_grounded.^2;

% Make edges of lambda grid:
dloglambda=log(Lambda(2)/Lambda(1));
Lambda_edges=repmat([exp(log(Lambda(1))-.5*dloglambda);sqrt(Lambda(1:end-1).*Lambda(2:end));exp(log(Lambda(end))+.5*dloglambda)],[1,numlcurves]);

% Adjust lambda grid for changing scales:
LambdaRescaleRatio_edges=[.5*(LambdaRescaleRatio(1:end-1)+LambdaRescaleRatio(2:end)),LambdaRescaleRatio(end)];
Lambda_edges=Lambda_edges.*repmat(LambdaRescaleRatio_edges,[numlambda+1,1]);

% Make edges of mesh size grid:
MeshSize_fastflow_edges=repmat([.5*(MeshSize_fastflow(1:end-1)+MeshSize_fastflow(2:end)),MeshSize_fastflow(end)+.5*(MeshSize_fastflow(end)-MeshSize_fastflow(end-1))],[numlambda+1,1]);

% % Create patch face indexing:
Faces=zeros((numlcurves-1)*numlambda,4);
Faces(1:numlambda,:)=[linspace(1,numlambda,numlambda)',linspace(2,numlambda+1,numlambda)',linspace(numlambda+3,2*numlambda+2,numlambda)',linspace(numlambda+2,2*numlambda+1,numlambda)'];
for ii=2:numlcurves-1
    Faces((ii-1)*numlambda+1:ii*numlambda,:)=repmat(Faces((ii-1)*numlambda,:)+1,[numlambda,1])+repmat(linspace(1,numlambda,numlambda)',[1,4]);
end

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

% Make figure:
figure(1)

% Fast-flow RMS(C):
% make subplot:
subplot('Position',Boxes{1,1})
% Make a patch:
patch('Vertices',[MeshSize_fastflow_edges(:)/1e3,Lambda_edges(:)],'Faces',Faces,'EdgeColor','none','FaceColor','flat','FaceVertexCData',reshape(RMS_C_fastflow(:,2:end),[],1))
% Plot best-fits:
hold on
DummyX=linspace(0,max(MeshSize_fastflow_edges(:)/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',1,'Color','k')
plot(DummyX,exp(equation_max(1)+equation_max(2)*DummyX),'LineWidth',1,'Color','k')
% Set up the rest:
set(gca,'YScale','log')
set(gca,'TickDir','out')
ylim([min(Lambda_edges(:)),max(Lambda_edges(:))])
xlim([MeshSize_fastflow_edges(1),MeshSize_fastflow_edges(end)]/1e3)
%xlabel('Median mesh size (km)','FontSize',axesfontsize)
set(gca,'XTickLabel',[])
ylabel('Regularization \lambda (unitless)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
caxis([0,max(RMS_C_fastflow(:))])
hc=colorbar;
set(get(hc,'YLabel'),'String','ln(C) units','FontSize',colorbarfontsize)
set(hc,'FontSize',colorbarfontsize)
title({'a) RMS misfit of ln(C)';'misfit computed w.r.t. Mesh #1 at that \lambda value'},'FontSize',titlefontsize)

% Fast-flow R2(C):
% make subplot:
subplot('Position',Boxes{2,1})
% Make a patch:
patch('Vertices',[MeshSize_fastflow_edges(:)/1e3,Lambda_edges(:)],'Faces',Faces,'EdgeColor','none','FaceColor','flat','FaceVertexCData',reshape(R2_C_fastflow(:,2:end),[],1))
% Plot best-fits:
hold on
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',1,'Color','k')
plot(DummyX,exp(equation_max(1)+equation_max(2)*DummyX),'LineWidth',1,'Color','k')
% Set up the rest:
set(gca,'YScale','log')
set(gca,'TickDir','out')
ylim([min(Lambda_edges(:)),max(Lambda_edges(:))])
xlim([MeshSize_fastflow_edges(1),MeshSize_fastflow_edges(end)]/1e3)
xlabel('Median mesh size (km)','FontSize',axesfontsize)
ylabel('Regularization \lambda (unitless)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
caxis([0,1])
hc=colorbar;
set(get(hc,'YLabel'),'String','R^2 (unitless)','FontSize',colorbarfontsize)
set(hc,'FontSize',colorbarfontsize)
title({'b) Correlation of ln(C),';'correlation computed w.r.t. Mesh #1 at that \lambda value'},'FontSize',titlefontsize)

% Fast-flow RMS(drag):
% make subplot:
subplot('Position',Boxes{1,2})
% Make a patch:
patch('Vertices',[MeshSize_fastflow_edges(:)/1e3,Lambda_edges(:)],'Faces',Faces,'EdgeColor','none','FaceColor','flat','FaceVertexCData',reshape(RMS_Drag_fastflow(:,2:end)/1e3,[],1))
% Plot best-fits:
hold on
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',1,'Color','k')
plot(DummyX,exp(equation_max(1)+equation_max(2)*DummyX),'LineWidth',1,'Color','k')
% Set up the rest:
set(gca,'YScale','log')
set(gca,'TickDir','out')
ylim([min(Lambda_edges(:)),max(Lambda_edges(:))])
xlim([MeshSize_fastflow_edges(1),MeshSize_fastflow_edges(end)]/1e3)
%xlabel('Median mesh size (km)','FontSize',axesfontsize)
%ylabel('Regularization \lambda (unitless)','FontSize',axesfontsize)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'FontSize',axesfontsize)
caxis([0,max(RMS_Drag_fastflow(:))/1e3])
hc=colorbar;
set(get(hc,'YLabel'),'String','kPa','FontSize',colorbarfontsize)
set(hc,'FontSize',colorbarfontsize)
title({'c) RMS misfit of basal drag';'misfit computed w.r.t. Mesh #1 at that \lambda value'},'FontSize',titlefontsize)

% Fast-flow R2(drag):
% make subplot:
subplot('Position',Boxes{2,2})
% Make a patch:
patch('Vertices',[MeshSize_fastflow_edges(:)/1e3,Lambda_edges(:)],'Faces',Faces,'EdgeColor','none','FaceColor','flat','FaceVertexCData',reshape(R2_Drag_fastflow(:,2:end),[],1))
% Plot best-fits:
hold on
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',1,'Color','k')
plot(DummyX,exp(equation_max(1)+equation_max(2)*DummyX),'LineWidth',1,'Color','k')
% Set up the rest:
set(gca,'YScale','log')
set(gca,'TickDir','out')
ylim([min(Lambda_edges(:)),max(Lambda_edges(:))])
xlim([MeshSize_fastflow_edges(1),MeshSize_fastflow_edges(end)]/1e3)
xlabel('Median mesh size (km)','FontSize',axesfontsize)
%ylabel('Regularization \lambda (unitless)','FontSize',axesfontsize)
set(gca,'YTickLabel',[])
set(gca,'FontSize',axesfontsize)
caxis([0,1])
hc=colorbar;
set(get(hc,'YLabel'),'String','R^2 (unitless)','FontSize',colorbarfontsize)
set(hc,'FontSize',colorbarfontsize)
title({'d) Correlation of basal drag,';'correlation computed w.r.t. Mesh #1 at that \lambda value'},'FontSize',titlefontsize)

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

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