% make_L1SPIRiT_plots.m -- make plots in Figs. 1, 3, 6, 7, and 8 in MRM
% paper and Figs. S2, S4 in the supplement 

% Copyright 2014-07-22 Daniel Weller, University of Michigan. All rights
% reserved.

% update path to image files
imfiles = { ...
'../data/real1.mat', ...
'../data/real2.mat', ...
% '/y/wellerds/MRI_data/MCSURE_v1.0/real1.mat', ...
% '/y/wellerds/MRI_data/MCSURE_v1.0/real2.mat' ...
};
imfile = 1; % for real image #1
imfile = 2; % for real image #2

% update paths to SPIRiT stuff (version 0.1 or 0.3)
if ~exist('SPIRiT', 'class')
%	addpath('../SPIRiT_v0.1/');
%	addpath('../ESPIRiT/');
%	addpath('../ESPIRiT/SPIRiT_code/');
	addpath('~/l/dat/code/matlab/ESPIRiT/');
	addpath('~/l/dat/code/matlab/ESPIRiT/SPIRiT_code/');
end
if ~exist('zpad.m', 'file')
%	addpath('../SPIRiT_v0.1/utils/');
%	addpath('../ESPIRiT/utils/');
	addpath('~/l/dat/code/matlab/ESPIRiT/utils/');
end
if ~exist('ticker.m', 'file')
	error('code requires IRT to be in path');
end


if ~exist('figures', 'var')
%     figures = 1; % just make Fig. 1
%     figures = 3; % just make Fig. 3
%     figures = 6; % just make Fig. 6
%     figures = 7; % just make Fig. 7
%     figures = 8; % just make Fig. 8
%     figures = 10; % just make Fig. S2
%     figures = 12; % just make Fig. S4
    
    figures = 1:12; % make all figures
end

switch imfile
    case 1
        amp_factor_R_figs310 = [6,5]; 
        amp_factor_R_figs612 = [48,2;48,3;48,4;48,5;48,6];
        amp_factor_R_fig7 = [48,4];
        amp_factor_R_fig8 = [48,4];
        
        window_image = [0,1];
        window_diff = [0,0.2];
        irect = [69,97,72,64];
        
        lambda = 1e4;
    case 2
        amp_factor_R_figs310 = [4,5]; 
        amp_factor_R_figs612 = [8,2;8,3;8,4;8,5;8,6];
        amp_factor_R_fig7 = [8,4];
        amp_factor_R_fig8 = [8,4];

        window_image = [0,1];
        window_diff = [0,0.2];
        irect = [129,29,84,63];
        
        lambda = 1e2;
    otherwise
        error('Bad image #');
end

rng('default');

%% start pool
if any(figures == 3) || any(figures == 10) || any(figures == 8)
    matlabpool('open','local');
end

%% load 8-channel coil sensitivities and brain slice from file 
load(imfiles{imfile},'ims','noisecov','mask*');

% re-cast as doubles (seems to be useful for proper BayesDESIGN/CG behavior
ims = double(ims);
noisecov = double(noisecov);

noisecov = noisecov./numel(mask); % use unitary Fourier transform in L1 SPIRiT code

% enable/disable sensitivity normalization
smap_norm = 1;

% get reference image
[Nx,Ny,NCha] = size(ims);
k_ref_init = fft2c(ims);
im_ref_init = sos(ims)./smap_norm;

mag = 100*max(size(im_ref_init)./irect(3:4));

noisecov_init = noisecov;
clear ims;

%% get amplitude factors
amp_factors = zeros(0,2);
if any(figures == 3) || any(figures == 10)
    amp_factors = union(amp_factors,amp_factor_R_figs310,'rows');
end
if any(figures == 6) || any(figures == 12)
    amp_factors = union(amp_factors,amp_factor_R_figs612,'rows');
end
if any(figures == 7)
    amp_factors = union(amp_factors,amp_factor_R_fig7,'rows');
end
if any(figures == 8)
    amp_factors = union(amp_factors,amp_factor_R_fig8,'rows');
end

%% start simulating 
WSUREs = cell([size(amp_factors,1),1]);
Y_SPIRiT_ref_normsqds = cell([size(amp_factors,1),1]);
wavWeights = cell([size(amp_factors,1),1]);
OracleMSEs = cell([size(amp_factors,1),1]);
WMSEs = cell([size(amp_factors,1),1]);
k_refs = cell([size(amp_factors,1),1]);
Ks = cell([size(amp_factors,1),1]);
noisecovs = cell([size(amp_factors,1),1]);
wavWeight_OMSEs = cell([size(amp_factors,1),1]);
wavWeight_WSUREs = cell([size(amp_factors,1),1]);
wavWeight_lcurves = cell([size(amp_factors,1),1]);
Y_SPIRiTs = cell([size(amp_factors,1),1]);
im_SPIRiTs = cell([size(amp_factors,1),1]);
im_SPIRiT_diffs = cell([size(amp_factors,1),1]);
Y_SPIRiT_OMSEs = cell([size(amp_factors,1),1]);
im_SPIRiT_OMSEs = cell([size(amp_factors,1),1]);
im_SPIRiT_OMSE_diffs = cell([size(amp_factors,1),1]);
Y_SPIRiT_WSUREs = cell([size(amp_factors,1),1]);
im_SPIRiT_WSUREs = cell([size(amp_factors,1),1]);
im_SPIRiT_WSURE_diffs = cell([size(amp_factors,1),1]);
Y_SPIRiT_lcurves = cell([size(amp_factors,1),1]);
im_SPIRiT_lcurves = cell([size(amp_factors,1),1]);
im_SPIRiT_lcurve_diffs = cell([size(amp_factors,1),1]);

rseed = rng();

% repeat over range of amplification factors
for iamp = 1:size(amp_factors,1)
    rng(rseed);
    amp_factor = amp_factors(iamp,1);
    R = amp_factors(iamp,2);
    
    %% amplify noise (optional)
    if amp_factor > 1
        noise = sqrtm((amp_factor-1).*noisecov_init/2)*complex(randn([NCha,Nx*Ny]),randn([NCha,Nx*Ny]));
        k_full = k_ref_init + reshape(noise.',[Nx,Ny,NCha]);
        noisecov = amp_factor.*noisecov_init;
    else
        noisecov = noisecov_init;
        k_full = k_ref_init;
    end
    im_full = sos(ifft2c(k_full))./smap_norm;
    noisecovs{iamp} = noisecov;

    %% undersample k-space
    mask_name = sprintf('mask_pdrand_crop_%dx',R);
    if exist(mask_name,'var')
        eval(['K = ' mask_name ';']);
    else
        error('No suitable mask available.');
    end
    fprintf(1,'Poisson-disc mask: reduction factor = %g.\n',numel(K)/sum(K(:))); drawnow;

    % ACS lines for SPIRiT calibration
    NACSx = 24; NACSy = 24;
    ACSFirstx = floor(Nx/2) + 1 - floor(NACSx/2);
    ACSFirsty = floor(Ny/2) + 1 - floor(NACSy/2);
    K(ACSFirstx:ACSFirstx+NACSx-1,ACSFirsty:ACSFirsty+NACSy-1) = true;
    k_ACS = k_full(ACSFirstx:ACSFirstx+NACSx-1,ACSFirsty:ACSFirsty+NACSy-1,:);
    Ks{iamp} = K;

    k_acq = k_full;
    k_acq(~K(:,:,ones(1,NCha))) = 0;

    %% normalize the data
    dcomp = zeros(size(K),class(im_ref_init)) + (sum(K(:))-NACSx*NACSy)/(numel(K)-NACSx*NACSy);
    dcomp(ACSFirstx:ACSFirstx+NACSx-1,ACSFirsty:ACSFirsty+NACSy-1) = 1;
    k_comp = k_acq.*dcomp(:,:,ones(1,NCha));
    scale_fctr = norm(k_comp(:))/sqrt(NCha)/20; % from demo_l1_spirit_pocs.m
    k_ref = k_ref_init./scale_fctr;
    im_ref = im_ref_init./scale_fctr;
    k_full = k_full./scale_fctr;
    im_full = im_full./scale_fctr;
    k_acq = k_acq./scale_fctr;
    k_ACS = k_ACS./scale_fctr;
    noisecov = noisecov./scale_fctr^2;
    k_refs{iamp} = k_ref;

    %% perform SPIRiT calibration
    kSize = [5,5];
    CalibTyk = 0.01;
    Kernel = zeros([kSize,NCha,NCha]);

    [AtA,~] = corrMatrix(k_ACS,kSize);
    for n=1:NCha
        Kernel(:,:,:,n) = calibrate(AtA,kSize,NCha,n,CalibTyk);
    end
    SPIRiTop = SPIRiT(Kernel, 'fft',[Nx,Ny]);

    %% (L1)SPIRiT options
    nSPIRiTiters = 25;

    %% do SPIRiT recon
    epsilon_scale_factor = norm(col(k_acq))/sqrt(numel(k_acq));
    Y_SPIRiT = cgSPIRiT(k_acq,SPIRiTop,nSPIRiTiters,0,k_acq);
    im_SPIRiT = sos(ifft2c(Y_SPIRiT))./smap_norm;
    im_SPIRiT_diff = im_SPIRiT-im_ref;
    Y_SPIRiTs{iamp} = Y_SPIRiT;
    im_SPIRiTs{iamp} = im_SPIRiT;
    im_SPIRiT_diffs{iamp} = im_SPIRiT_diff;

    %% do SPIRiT recon of reference k-space
    k_acq_ref = k_ref;
    k_acq_ref(~K(:,:,ones(1,NCha))) = 0;
    Y_SPIRiT_ref = cgSPIRiT(k_acq_ref,SPIRiTop,nSPIRiTiters,0,k_acq_ref);
    Y_SPIRiT_ref_normsqd = Y_SPIRiT_ref;
    Y_SPIRiT_ref_normsqd(K(:,:,ones(1,NCha))) = 0;
    Y_SPIRiT_ref_normsqd = abs(Y_SPIRiT_ref_normsqd(:)'*Y_SPIRiT_ref_normsqd(:));
    Y_SPIRiT_ref_normsqds{iamp} = Y_SPIRiT_ref_normsqd;
    im_SPIRiT_ref = sos(ifft2c(Y_SPIRiT_ref))./smap_norm;
    im_SPIRiT_ref_diff = im_SPIRiT_ref-im_ref;

    %% generate noise for WSURE estimation
    if (any(figures == 3) || any(figures == 10) || any(figures == 8)) && (any(all(bsxfun(@eq,amp_factor_R_figs310,amp_factors(iamp,:)),2)) || any(all(bsxfun(@eq,amp_factor_R_fig8,amp_factors(iamp,:)),2)))
        Nbs = 20;
        bs = complex(zeros([Nx,Ny,NCha,Nbs]));
        Y_SPIRiT_bs = complex(zeros([Nx,Ny,NCha,Nbs]));
        Kcha = K(:,:,ones(1,NCha));
        parfor ib = 1:Nbs
            breal = rand([Nx,Ny,NCha]);
            breal(breal >= 0.5) = 1;
            breal(breal < 0.5) = -1;
            bimag = rand([Nx,Ny,NCha]);
            bimag(bimag >= 0.5) = 1;
            bimag(bimag < 0.5) = -1;
            b = complex(breal,bimag)./sqrt(2);
            b(~Kcha) = 0;
            bs(:,:,:,ib) = b;

            % compute SPIRiT(Lambda*b)
            SPIRiTop_b = SPIRiT(Kernel, 'fft',[Nx,Ny]);
            Y_SPIRiT_b = reshape((noisecov*reshape(b,[Nx*Ny,NCha]).').',[Nx,Ny,NCha]);
            Y_SPIRiT_b = cgSPIRiT(Y_SPIRiT_b,SPIRiTop_b,nSPIRiTiters,0,Y_SPIRiT_b);
            Y_SPIRiT_bs(:,:,:,ib) = Y_SPIRiT_b;
        end
        b = bs(:,:,:,1);
        Y_SPIRiT_b = Y_SPIRiT_bs(:,:,:,1);
    else
        Kcha = K(:,:,ones(1,NCha));
        breal = rand([Nx,Ny,NCha]);
        breal(breal >= 0.5) = 1;
        breal(breal < 0.5) = -1;
        bimag = rand([Nx,Ny,NCha]);
        bimag(bimag >= 0.5) = 1;
        bimag(bimag < 0.5) = -1;
        b = complex(breal,bimag)./sqrt(2);
        b(~Kcha) = 0;

        % compute SPIRiT(Lambda*b)
        Y_SPIRiT_b = reshape((noisecov*reshape(b,[Nx*Ny,NCha]).').',[Nx,Ny,NCha]);
        Y_SPIRiT_b = cgSPIRiT(Y_SPIRiT_b,SPIRiTop,nSPIRiTiters,0,Y_SPIRiT_b);
    end

    %% compare WSURE, WMSE, OMSE values for Fig. 3
    if (any(figures == 3) || any(figures == 10) || any(figures == 8)) && (any(all(bsxfun(@eq,amp_factor_R_figs310,amp_factors(iamp,:)),2)) || any(all(bsxfun(@eq,amp_factor_R_fig8,amp_factors(iamp,:)),2)))
        if (any(figures == 3) || any(figures == 10)) && any(all(bsxfun(@eq,amp_factor_R_figs310,amp_factors(iamp,:)),2))
            epsilons = logspace(-5,-2,4); % try a few to get an idea of what is small enough
        else
            epsilons = 1e-4; % just do one
        end
        wavWeights1 = union(logspace(-4,-1,13),logspace(-3,-2,9),'sorted');
        WSUREs1 = zeros(length(wavWeights1),length(epsilons),Nbs);
        WMSEs1 = zeros(length(wavWeights1),1);
        OracleMSEs1 = zeros(length(wavWeights1),1);
        fprintf(1,'Comparing WSURE, WMSE, OMSE values for various wavWeights, epsilons...\n'); drawnow;
        for iw = 1:length(wavWeights1)
            wavWeight = wavWeights1(iw);
            for ie = 1:length(epsilons)
                epsilon = epsilons(ie);
                if ie == 1
                    [WSUREs1(iw,ie,1),Y_L1SPIRiT] = L1SPIRiT_WSURE(k_acq,Y_SPIRiT,bs(:,:,:,1),Y_SPIRiT_bs(:,:,:,1),SPIRiTop,nSPIRiTiters,wavWeight,epsilon*epsilon_scale_factor);
                    parfor ib = 2:Nbs
                        SPIRiTop_b = SPIRiT(Kernel, 'fft',[Nx,Ny]);
                        WSUREs1(iw,ie,ib) = L1SPIRiT_WSURE(k_acq,Y_SPIRiT,bs(:,:,:,ib),Y_SPIRiT_bs(:,:,:,ib),SPIRiTop_b,nSPIRiTiters,wavWeight,epsilon*epsilon_scale_factor,Y_L1SPIRiT);
                    end
                else
                    parfor ib = 1:Nbs
                        SPIRiTop_b = SPIRiT(Kernel, 'fft',[Nx,Ny]);
                        WSUREs1(iw,ie,ib) = L1SPIRiT_WSURE(k_acq,Y_SPIRiT,bs(:,:,:,ib),Y_SPIRiT_bs(:,:,:,ib),SPIRiTop_b,nSPIRiTiters,wavWeight,epsilon*epsilon_scale_factor,Y_L1SPIRiT);
                    end
                end
                fprintf(1,' %9g',mean(WSUREs1(iw,ie,:))+Y_SPIRiT_ref_normsqd); drawnow;
            end
            OMSE = Y_L1SPIRiT - k_ref;
            OMSE(K(:,:,ones(1,NCha))) = 0;
            OMSE = abs(OMSE(:)'*OMSE(:));
            OracleMSEs1(iw) = OMSE;
            WMSE = Y_L1SPIRiT - Y_SPIRiT_ref;
            WMSE(K(:,:,ones(1,NCha))) = 0;
            WMSE = abs(WMSE(:)'*WMSE(:));
            WMSEs1(iw) = WMSE;
            fprintf(1,' %9g %9g\n',WMSEs1(iw),OracleMSEs1(iw)); drawnow;
        end
        WSUREs{iamp} = WSUREs1;
        wavWeights{iamp} = wavWeights1;
        OracleMSEs{iamp} = OracleMSEs1;
        WMSEs{iamp} = WMSEs1;
    end

    %% done with multiple-noise part
    if (any(figures == 3) || any(figures == 10) || any(figures == 8)) && (any(all(bsxfun(@eq,amp_factor_R_figs310,amp_factors(iamp,:)),2)) || any(all(bsxfun(@eq,amp_factor_R_fig8,amp_factors(iamp,:)),2)))
        clear bs Y_SPIRiT_bs;
    end

    %% do brute-force search for optimal wavWeight that minimizes true MSE
    if (any(figures == 6) || any(figures == 7) || any(figures == 12)) && (any(all(bsxfun(@eq,amp_factor_R_figs612,amp_factors(iamp,:)),2)) || any(all(bsxfun(@eq,amp_factor_R_fig7,amp_factors(iamp,:)),2)))
        tSPIRiT_OMSE = tic();
        fprintf(1,'Using fminsearch() to find Oracle WMSE-optimal wavWeight...\n'); drawnow;
        [logwavWeight_opt,OMSE_opt,~,fminsearch_output] = fminsearch(@(logwavWeight) sum(sum(abs(reshape(~K(:,:,ones(1,NCha)).*(pocsSPIRiT(k_acq,SPIRiTop,nSPIRiTiters,k_acq,10.^logwavWeight,0) - k_ref),[Nx*Ny,NCha])).^2)),log10(1e-2),optimset('Display','iter','FunValCheck','on'));
        wavWeight_OMSE = 10.^logwavWeight_opt;
        tSPIRiT_OMSE = toc(tSPIRiT_OMSE);
        SPIRiTevals_OMSE = fminsearch_output.funcCount;
        fprintf(1,'Oracle WMSE search: min(OMSE) = %g at gamma = %g.\n',OMSE_opt,wavWeight_OMSE); drawnow;
        Y_SPIRiT_OMSE = pocsSPIRiT(k_acq,SPIRiTop,nSPIRiTiters,k_acq,wavWeight_OMSE,0);
        im_SPIRiT_OMSE = sos(ifft2c(Y_SPIRiT_OMSE))./smap_norm;
        im_SPIRiT_OMSE_diff = im_SPIRiT_OMSE-im_ref;
        wavWeight_OMSEs{iamp} = wavWeight_OMSE;
        Y_SPIRiT_OMSEs{iamp} = Y_SPIRiT_OMSE;
        im_SPIRiT_OMSEs{iamp} = im_SPIRiT_OMSE;
        im_SPIRiT_OMSE_diffs{iamp} = im_SPIRiT_OMSE_diff;
    end

    %% do brute-force search for optimal wavWeight that minimizes WSURE
    if (any(figures == 6) || any(figures == 7) || any(figures == 12)) && (any(all(bsxfun(@eq,amp_factor_R_figs612,amp_factors(iamp,:)),2)) || any(all(bsxfun(@eq,amp_factor_R_fig7,amp_factors(iamp,:)),2)))
        tSPIRiT_WSURE = tic();
        epsilon = 1e-4; % base on earlier evaluations for what works best
        NwavWeights = 9;
        logwavWeight_range = [-8,2];
        Nlevels = 2;
        WSUREs_search = zeros(Nlevels,NwavWeights);
        logwavWeights = zeros(Nlevels,NwavWeights+2);
        fprintf(1,'Using coarse-to-fine parameter sweep to find WSURE-optimal wavWeight...\n'); drawnow;
        for ilvl = 1:Nlevels
            logwavWeights(ilvl,:) = linspace(logwavWeight_range(1),logwavWeight_range(2),NwavWeights+2);
            for iw = 1:NwavWeights
                [WSUREs_search(ilvl,iw),Y_L1SPIRiT] = L1SPIRiT_WSURE(k_acq,Y_SPIRiT,b,Y_SPIRiT_b,SPIRiTop,nSPIRiTiters,10.^logwavWeights(ilvl,iw+1),epsilon*epsilon_scale_factor);
                fprintf(1,'WSURE search level %d: WSURE(%g) = %g + ||GMx||^2.\n',ilvl,10.^logwavWeights(ilvl,iw+1),WSUREs_search(ilvl,iw)); drawnow;
            end
            [WSURE_opt,iw] = min(WSUREs_search(ilvl,:));
            wavWeight_WSURE = 10.^logwavWeights(ilvl,iw+1);
            fprintf(1,'WSURE search level %d: min(WSURE) = %g + ||GMx||^2 at wavWeight = %g.\n',ilvl,WSURE_opt,wavWeight_WSURE); drawnow;
            logwavWeight_range = logwavWeights(ilvl,[iw,iw+2]);
        end
        tSPIRiT_WSURE = toc(tSPIRiT_WSURE);
        SPIRiTevals_WSURE = 2*NwavWeights*Nlevels;
        Y_SPIRiT_WSURE = pocsSPIRiT(k_acq,SPIRiTop,nSPIRiTiters,k_acq,wavWeight_WSURE,0);
        im_SPIRiT_WSURE = sos(ifft2c(Y_SPIRiT_WSURE))./smap_norm;
        im_SPIRiT_WSURE_diff = im_SPIRiT_WSURE-im_ref;
        wavWeight_WSUREs{iamp} = wavWeight_WSURE;
        Y_SPIRiT_WSUREs{iamp} = Y_SPIRiT_WSURE;
        im_SPIRiT_WSUREs{iamp} = im_SPIRiT_WSURE;
        im_SPIRiT_WSURE_diffs{iamp} = im_SPIRiT_WSURE_diff;
    end

    %% do "L-curve" type of analysis
    if (any(figures == 6) || any(figures == 7) || any(figures == 12)) && (any(all(bsxfun(@eq,amp_factor_R_figs612,amp_factors(iamp,:)),2)) || any(all(bsxfun(@eq,amp_factor_R_fig7,amp_factors(iamp,:)),2)))
        tSPIRiT_Lcurve = tic();
        NwavWeights = 9;
        logwavWeight_range = [-8,2];
        Nlevels = 2;
        L1objs = zeros(Nlevels,NwavWeights);
        SPIRiTobjs = zeros(Nlevels,NwavWeights);
        logwavWeights = zeros(Nlevels,NwavWeights+2);
        fprintf(1,'Using coarse-to-fine parameter sweep to do l-curve analysis...\n'); drawnow;
        [~,SPIRiTobj_min,L1obj_min] = pocsSPIRiT_wrapper(k_acq,SPIRiTop,nSPIRiTiters,k_acq,10.^logwavWeight_range(1),0);
        [~,SPIRiTobj_max,L1obj_max] = pocsSPIRiT_wrapper(k_acq,SPIRiTop,nSPIRiTiters,k_acq,10.^logwavWeight_range(2),0);
        for ilvl = 1:Nlevels
            logwavWeights(ilvl,:) = linspace(logwavWeight_range(1),logwavWeight_range(2),NwavWeights+2);
            for iw = 1:NwavWeights
                [~,SPIRiTobjs(ilvl,iw),L1objs(ilvl,iw)] = pocsSPIRiT_wrapper(k_acq,SPIRiTop,nSPIRiTiters,k_acq,10.^logwavWeights(ilvl,iw+1),0);
                fprintf(1,'L-curve search level %d: SPIRiTobj(%g) = %g, L1obj(%g) = %g.\n',ilvl,10.^logwavWeights(ilvl,iw+1),SPIRiTobjs(ilvl,iw),10.^logwavWeights(ilvl,iw+1),L1objs(ilvl,iw)); drawnow;
            end
            L1deriv_post = (log10([L1objs(ilvl,2:end),L1obj_max(ilvl)]) - log10(L1objs(ilvl,:)));
            L1deriv_pre = (log10(L1objs(ilvl,:)) - log10([L1obj_min(ilvl),L1objs(ilvl,1:end-1)]));
            SPIRiTderiv_post = log10([SPIRiTobjs(ilvl,2:end),SPIRiTobj_max(ilvl)]) - log10(SPIRiTobjs(ilvl,:));
            SPIRiTderiv_pre = log10(SPIRiTobjs(ilvl,:)) - log10([SPIRiTobj_min(ilvl),SPIRiTobjs(ilvl,1:end-1)]);

            deriv = (L1deriv_post + L1deriv_pre)./(SPIRiTderiv_post + SPIRiTderiv_pre);
            deriv2 = 4.*(L1deriv_post-L1deriv_pre)./(SPIRiTderiv_post + SPIRiTderiv_pre).^2;

            firstvalidderiv = find(deriv <= 0 & ~isnan(deriv),1,'first');
            if isempty(firstvalidderiv), firstvalidderiv = length(deriv)+1; end
            deriv2 = deriv2(firstvalidderiv:end);
            firstinvalidderiv = find(deriv(firstvalidderiv:end) > 0 | isnan(deriv(firstvalidderiv:end)),1,'first');
            if ~isempty(firstinvalidderiv)
                deriv2 = deriv2(1:firstinvalidderiv-1);
                if firstinvalidderiv > 1
    %                 deriv2(firstinvalidderiv-1) = -L1deriv_pre(firstinvalidderiv-1)./SPIRiTderiv_pre(firstinvalidderiv-1).^2;
                    deriv2 = deriv2(1:firstinvalidderiv-2);
                end
            end
            if isempty(deriv2)
                fprintf(2,'L-curve search level %d: no valid part of pareto curve found.\n',ilvl); drawnow;
                deriv2_max = NaN; iw = floor(NwavWeights/2)+1;
%                 break;
            else
                [deriv2_max,iw] = max(deriv2);
                iw = iw + (firstvalidderiv-1);
                if deriv2_max <= 0
                    fprintf(1,'L-curve search level %d: pareto curve seems to be concave. Continuing...\n',ilvl); drawnow;
                end
            end

            wavWeight_lcurve = 10.^logwavWeights(ilvl,iw+1);
            fprintf(1,'L-curve search level %d: max(deriv2) = %g at wavWeight = %g.\n',ilvl,deriv2_max,wavWeight_lcurve); drawnow;
            if iw == 1, L1obj_min(ilvl+1) = L1obj_min(ilvl); else L1obj_min(ilvl+1) = L1objs(ilvl,iw-1); end
            if iw == NwavWeights, L1obj_max(ilvl+1) = L1obj_max(ilvl); else L1obj_max(ilvl+1) = L1objs(ilvl,iw+1); end
            if iw == 1, SPIRiTobj_min(ilvl+1) = SPIRiTobj_min(ilvl); else SPIRiTobj_min(ilvl+1) = SPIRiTobjs(ilvl,iw-1); end
            if iw == NwavWeights, SPIRiTobj_max(ilvl+1) = SPIRiTobj_max(ilvl); else SPIRiTobj_max(ilvl+1) = SPIRiTobjs(ilvl,iw+1); end
            logwavWeight_range = logwavWeights(ilvl,[iw,iw+2]);
        end
        tSPIRiT_Lcurve = toc(tSPIRiT_Lcurve);
        SPIRiTevals_Lcurve = NwavWeights*Nlevels+2;
        Y_SPIRiT_lcurve = pocsSPIRiT(k_acq,SPIRiTop,nSPIRiTiters,k_acq,wavWeight_lcurve,0);
        im_SPIRiT_lcurve = sos(ifft2c(Y_SPIRiT_lcurve))./smap_norm;
        im_SPIRiT_lcurve_diff = im_SPIRiT_lcurve-im_ref;
        wavWeight_lcurves{iamp} = wavWeight_lcurve;
        Y_SPIRiT_lcurves{iamp} = Y_SPIRiT_lcurve;
        im_SPIRiT_lcurves{iamp} = im_SPIRiT_lcurve;
        im_SPIRiT_lcurve_diffs{iamp} = im_SPIRiT_lcurve_diff;
    end

%% end simulations
end

%% cleanup

if any(figures == 3) || any(figures == 10) || any(figures == 8)
    matlabpool('close');
end

%% do plots

%% generate plot of image, inset -- Fig. 1
if any(figures == 1)
    if ~exist('im_ref','var'), error('Need to run any other experiment to get normalized reference image'); end
    figure; imshow(squeeze(im_ref),window_image,'Border','tight'); hold on;
    hrect = rectangle('Position',irect([2,1,4,3]),'LineWidth',2,'EdgeColor','w');

    figure; imshow(squeeze(im_ref(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_image,'Border','tight','InitialMagnification',mag);
end

%% WSURE vs. true and approx. WMSE plot -- Fig. 3

if any(figures == 3)
    iamp = find(all(bsxfun(@eq,amp_factors,amp_factor_R_figs310(1,:)),2),1);
    
    eps_index = 2; % epsilon = 1e-4
    nstds = 5;
    WSURE_val = WSUREs{iamp}(:,eps_index,1) + Y_SPIRiT_ref_normsqds{iamp};
    WSURE_mean = mean(WSUREs{iamp}(:,eps_index,:),3) + Y_SPIRiT_ref_normsqds{iamp};
    WSURE_std = std(WSUREs{iamp}(:,eps_index,:),0,3);
    WSURE_lower = WSURE_mean - nstds.*WSURE_std;
    WSURE_upper = WSURE_mean + nstds.*WSURE_std;

    % make plot
    figure; 
%     subplot(2,1,1);
    plot(log10(wavWeights{iamp}),10*log10(OracleMSEs{iamp}),'-o',log10(wavWeights{iamp}),10*log10(WMSEs{iamp}),'-s'); hold all;
    errorbar(log10(wavWeights{iamp}),10*log10(WSURE_val),10*log10(WSURE_val./WSURE_lower),10*log10(WSURE_upper./WSURE_val));
    hcs = get(gca,'Children');
    legend(hcs([3,1,2]),'True WMSE (Eq. [3])','WSURE (\epsilon = 10^{-4})','Approx. WMSE (Eq. [4])','Location','NorthWest');
    xlim([floor(log10(min(wavWeights{iamp}))),ceil(log10(max(wavWeights{iamp})))]);
    xticks = floor(log10(min(wavWeights{iamp}))):ceil(log10(max(wavWeights{iamp})));
    xtickstrs = arrayfun(@(x) sprintf('%g',10.^x),xticks,'UniformOutput',false);
    set(gca,'Box','off','XTick',xticks,'XTickLabel',xtickstrs);
    xlabel('Regularization parameter (\gamma)');
    ylabel('Squared error (dB)');
end

%% WSURE vs. true and approx. WMSE plot for violation at R = 6 -- Fig. 8

if any(figures == 8)
    iamp = find(all(bsxfun(@eq,amp_factors,amp_factor_R_fig8(1,:)),2),1);
    
    eps_index = 1; % for violation curve only
    nstds = 5;
    WSURE_val = WSUREs{iamp}(:,eps_index,1) + Y_SPIRiT_ref_normsqds{iamp};
    WSURE_mean = mean(WSUREs{iamp}(:,eps_index,:),3) + Y_SPIRiT_ref_normsqds{iamp};
    WSURE_std = std(WSUREs{iamp}(:,eps_index,:),0,3);
    WSURE_lower = WSURE_mean - nstds.*WSURE_std;
    WSURE_upper = WSURE_mean + nstds.*WSURE_std;

    % make plot
    figure; 
%     subplot(2,1,1);
    plot(log10(wavWeights{iamp}),10*log10(OracleMSEs{iamp}),'-o',log10(wavWeights{iamp}),10*log10(WMSEs{iamp}),'-s'); hold all;
    errorbar(log10(wavWeights{iamp}),10*log10(WSURE_val),10*log10(WSURE_val./WSURE_lower),10*log10(WSURE_upper./WSURE_val));
    hcs = get(gca,'Children');
    legend(hcs([3,1,2]),'True WMSE (Eq. [3])','WSURE (\epsilon = 10^{-4})','Approx. WMSE (Eq. [4])','Location','NorthWest');
    xlim([floor(log10(min(wavWeights{iamp}))),ceil(log10(max(wavWeights{iamp})))]);
    xticks = floor(log10(min(wavWeights{iamp}))):ceil(log10(max(wavWeights{iamp})));
    xtickstrs = arrayfun(@(x) sprintf('%g',10.^x),xticks,'UniformOutput',false);
    set(gca,'Box','off','XTick',xticks,'XTickLabel',xtickstrs);
    xlabel('Regularization parameter (\gamma)');
    ylabel('Squared error (dB)');
end

%% WSURE vs. true and approx. WMSE plot for all epsilons -- Fig. S2

if any(figures == 10)
    iamp = find(all(bsxfun(@eq,amp_factors,amp_factor_R_figs310(1,:)),2),1);

    % compute WSURE stuff
    WSURE_vals = WSUREs{iamp}(:,:,1) + Y_SPIRiT_ref_normsqds{iamp};
    WSURE_means = mean(WSUREs{iamp},3) + Y_SPIRiT_ref_normsqds{iamp};
    WSURE_stds = std(WSUREs{iamp},0,3);
    WSURE_lowers = WSURE_means - WSURE_stds;
    WSURE_uppers = WSURE_means + WSURE_stds;

    % make plot
    figure; 
%     subplot(2,1,1); 
    plot(log10(wavWeights{iamp}),10*log10(OracleMSEs{iamp}),'-o',log10(wavWeights{iamp}),10*log10(WMSEs{iamp}),'-s',log10(wavWeights{iamp}),10*log10(WSURE_vals),'--');
    hcs = get(gca,'Children');
    legendstrs = {'True WMSE (Eq. [3])','WSURE (\epsilon = 10^{-5})','WSURE (\epsilon = 10^{-4})','WSURE (\epsilon = 0.001)','WSURE (\epsilon = 0.01)','Approx. WMSE (Eq. [4])'};
    legend(hcs([6,4:-1:1,5]),legendstrs,'Location','NorthWest');
    xlim([floor(log10(min(wavWeights{iamp}))),ceil(log10(max(wavWeights{iamp})))]);
    xticks = floor(log10(min(wavWeights{iamp}))):ceil(log10(max(wavWeights{iamp})));
    xtickstrs = arrayfun(@(x) sprintf('%g',10.^x),xticks,'UniformOutput',false);
    set(gca,'Box','off','XTick',xticks,'XTickLabel',xtickstrs);
    xlabel('Regularization parameter (\gamma)');
    ylabel('Squared error (dB)');
end

%% compute trend plot for WSURE-optimal and WMSE-optimal DESIGN recons vs. noise level -- Fig. 6 or Fig. S4

if any(figures == 6) || any(figures == 12)
    [~,iamps] = intersect(amp_factors,amp_factor_R_figs612,'rows');
    
    SNRs = zeros(size(iamps));
    accels = zeros(size(iamps));
    gammas_SPIRiT_OMSE = zeros(size(iamps));
    gammas_SPIRiT_Lcurve = zeros(size(iamps));
    gammas_SPIRiT_WSURE = zeros(size(iamps));
    WMSEs_SPIRiT = zeros(size(iamps));
    WMSEs_SPIRiT_OMSE = zeros(size(iamps));
    WMSEs_SPIRiT_Lcurve = zeros(size(iamps));
    WMSEs_SPIRiT_WSURE = zeros(size(iamps));

    for ii = 1:length(iamps)
        iamp = iamps(ii);
        accels(ii) = numel(Ks{iamp})/sum(Ks{iamp}(:));
        SNRs(ii) = norm(k_refs{iamp}(:))^2/(numel(Ks{iamp})*sum(abs(svd(noisecovs{iamp}))));
        gammas_SPIRiT_OMSE(ii) = wavWeight_OMSEs{iamp};
        gammas_SPIRiT_Lcurve(ii) = wavWeight_lcurves{iamp};
        gammas_SPIRiT_WSURE(ii) = wavWeight_WSUREs{iamp};
        KCha = Ks{iamp}(:,:,:,ones(1,size(noisecovs{iamp},1)));
        WMSEs_SPIRiT(ii) = norm((~KCha(:)).*(Y_SPIRiTs{iamp}(:)-k_refs{iamp}(:)))^2;
        WMSEs_SPIRiT_OMSE(ii) = norm((~KCha(:)).*(Y_SPIRiT_OMSEs{iamp}(:)-k_refs{iamp}(:)))^2;
        WMSEs_SPIRiT_Lcurve(ii) = norm((~KCha(:)).*(Y_SPIRiT_lcurves{iamp}(:)-k_refs{iamp}(:)))^2;
        WMSEs_SPIRiT_WSURE(ii) = norm((~KCha(:)).*(Y_SPIRiT_WSUREs{iamp}(:)-k_refs{iamp}(:)))^2;
    end

    sortedvecs = sortrows([accels(:),gammas_SPIRiT_OMSE(:),gammas_SPIRiT_Lcurve(:),gammas_SPIRiT_WSURE(:),WMSEs_SPIRiT(:),WMSEs_SPIRiT_OMSE(:),WMSEs_SPIRiT_Lcurve(:),WMSEs_SPIRiT_WSURE(:)],1);

    figure; subplot(2,1,2); plot(sortedvecs(:,1),sortedvecs(:,2),'-x',sortedvecs(:,1),sortedvecs(:,4),'-o',sortedvecs(:,1),sortedvecs(:,3),'-*');
    xlabel('Undersampling factor (R)');
    ylabel('Regularization parameter (\gamma)');
    hcs = get(gca,'Children');
    set(hcs(3),'MarkerSize',14);
    xlim(sortedvecs([1,end],1));
    set(gca,'Box','off');

    subplot(2,1,1); plot(sortedvecs(:,1),10*log10(sortedvecs(:,6)),'-x',sortedvecs(:,1),10*log10(sortedvecs(:,8)),'-o',sortedvecs(:,1),10*log10(sortedvecs(:,7)),'-*',sortedvecs(:,1),10*log10(sortedvecs(:,5)),'-s');
    xlabel('Undersampling factor (R)');
    ylabel('WMSE (dB)');
    hcs = get(gca,'Children');
    set(hcs(4),'MarkerSize',14);
    legend(hcs([1,2,3,4]),'SPIRiT','Lcurve-opt. L_1-SPIRiT','WSURE-opt. L_1-SPIRiT','WMSE-opt. L_1-SPIRiT','Location','NorthWest');
    xlim(sortedvecs([1,end],1));
    set(gca,'Box','off');
end

%% show reconstructed images -- Fig. 7

if any(figures == 7)
    iamp = find(all(bsxfun(@eq,amp_factors,amp_factor_R_fig7(1,:)),2),1);
    
    % SPIRiT alone
    figure; imshow(squeeze(im_SPIRiTs{iamp}),window_image,'Border','tight'); hold on;
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiTs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_image); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);
    figure; imshow(squeeze(im_SPIRiT_diffs{iamp}),window_diff,'Border','tight'); hold on;
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_diffs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_diff); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);

    % oracle
    figure; imshow(squeeze(im_SPIRiT_OMSEs{iamp}),window_image,'Border','tight'); hold on;
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_OMSEs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_image); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);
    figure; imshow(squeeze(im_SPIRiT_OMSE_diffs{iamp}),window_diff,'Border','tight'); hold on;
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_OMSE_diffs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_diff); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);

    % lcurve
    figure; imshow(squeeze(im_SPIRiT_lcurves{iamp}),window_image,'Border','tight'); hold on;
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_lcurves{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_image); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);
    figure; imshow(squeeze(im_SPIRiT_lcurve_diffs{iamp}),window_diff,'Border','tight'); hold on;
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_lcurve_diffs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_diff); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);

    % WSURE
    figure; imshow(squeeze(im_SPIRiT_WSUREs{iamp}),window_image,'Border','tight'); hold on;
    hcb = colorbar('Location','South');
    set(hcb,'XAxisLocation','bottom','XTick',window_image);
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_WSUREs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_image); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);
    figure; imshow(squeeze(im_SPIRiT_WSURE_diffs{iamp}),window_diff,'Border','tight'); hold on;
    hcb = colorbar('Location','South');
    set(hcb,'XAxisLocation','bottom','XTick',window_diff);
    hin = axes('Position',[0.5,0,0.5,0.5]);
    imagesc(squeeze(im_SPIRiT_WSURE_diffs{iamp}(irect(1) + (0:irect(3)-1),irect(2) + (0:irect(4)-1))),window_diff); colormap gray; axis image;
    set(hin,'XTick',[],'YTick',[],'Box','on','LineWidth',2);
end
