% created by Nils Hase 2016 (nilshase@math.uni-bremen.de)

% Script to run the Tikhonov based inversion methods for the Barnett
% scenario

%% Which method to use
% These can be calculated
run_tp_l2 = 1;          % L2
run_tp_l1 = 1;          % L1
run_tp_l2_pos = 1;      % L2 POS
run_tp_l1_pos = 1;      % L1 POS
run_tp_l1_dic = 1;      % L1 DIC
run_tp_l1_dic_pos = 1;  % L1 DIC POS

% Calculate or just load results from 'inverse_methods_solutions.mat'
load_saved = 0; 
filename = 'barnett_scenario_solutions.mat';

% For these methods from Miller et al. (2014) only the result is available
% (loaded from 'inverse_methods_solutions.mat')
show_st_inv = 1;        % standard inversion
show_lm = 1;            % Lagrange multiplier method
show_tr_inv = 1;        % transform inversion
show_gibbs = 1;         % Gibbs sampler


%% Load footprints, measurements and setup data
load('synthetic_methane_data.mat')

% Which measurements to use?
stat_select = [sel1;sel2;sel3;sel4;sel5;sel6;sel7;sel8;sel9;sel10];

A = A(stat_select,:);
R = R(stat_select,stat_select);
y = y(stat_select);

% Plot EDGAR fluxes
map_result(f_EDGAR,'EDGAR fluxes');

% Add emissions from the Barnett shale formation that was active during
% 2007/2008 
% Leakage of p% -> p*438.66 mol/s = f(Barnett).*area(Barnett)
% Leakage of 1.5% -> 658 mol/s 
ff = zeros(size(f_EDGAR));
total_barnett = 1.5*438.66;
b_index = get_index(-99,-97,32,33);
ff(b_index) = [2,2.5,6,6,2,1.5];
ff(b_index) = ff(b_index)* total_barnett/sum(ff(b_index).*area(b_index));
f_barnett = f_EDGAR + ff;

% Recalculate noisy synthetic with added Barnett emissions
y = y + A*ff;

% Plot EDGAR + Barnett emissions
map_result(f_barnett,'EDGAR + Barnett');

%% Setup initial values and parameters
delta = 1;                                              % discrepancy is normalized
L = sparse(1:size(R,1),1:size(R,1),1./sqrt(diag(R)));   % weighting matrix

alpha = 10.^(5:-0.1:-3);                                % regularization parameter
x_a = zeros(size(f_EDGAR));                             % apriori estimate of surface fluxes (will not work if different from 0)

show_barnett('Barnett',f_barnett,[],b_index);


%% Tikhonov Regularization with l2-penalty
if run_tp_l2 == 1
    methodname = 'TP L2';
    if load_saved
        load(filename,'f_l2','alpha_l2'); f_opt = f_l2; alpha_opt = alpha_l2; runtime = 0;
    else
        tic; [f_opt,alpha_opt] = TP_CGLS_fast(L*A,L*y,alpha,delta,0); runtime = toc;
    end
    f_l2 = f_opt; alpha_l2 = alpha_opt;
    show_barnett(methodname,f_opt,alpha_opt,b_index);
end

%% Tikhonov Regularization with l1-penalty
if run_tp_l1 == 1
    methodname = 'TP L1';
    if load_saved
        load(filename,'f_l1','alpha_l1'); f_opt = f_l1; alpha_opt = alpha_l1; runtime = 0;
    else
        tic; [f_opt,alpha_opt] = FISTA_L1(L*A,L*y,x_a,alpha,0); runtime = toc;
    end
    f_l1 = f_opt; alpha_l1 = alpha_opt;
    show_barnett(methodname,f_opt,alpha_opt,b_index);
end

%% Tikhonov Regularization with l2-penalty & positivity constraint
if run_tp_l2_pos == 1
    methodname = 'TP L2 POS';
    if load_saved
        load(filename,'f_l2_pos','alpha_l2_pos'); f_opt = f_l2_pos; alpha_opt = alpha_l2_pos; runtime = 0;
    else
        tic; [f_opt,alpha_opt] = FISTA_L2_pos(L*A,L*y,x_a,alpha,0); runtime = toc;
    end
    f_l2_pos = f_opt; alpha_l2_pos = alpha_opt;
    show_barnett(methodname,f_opt,alpha_opt,b_index);
end
%% Tikhonov Regularization with l1-penalty & positivity constraint
if run_tp_l1_pos == 1
    methodname = 'TP L1 POS';
    if load_saved
        load(filename,'f_l1_pos','alpha_l1_pos'); f_opt = f_l1_pos; alpha_opt = alpha_l1_pos; runtime = 0;
    else
        tic; [f_opt,alpha_opt] = FISTA_L1_pos(L*A,L*y,x_a,alpha,0); runtime = toc;
    end
    f_l1_pos = f_opt; alpha_l1_pos = alpha_opt;
    show_barnett(methodname,f_opt,alpha_opt,b_index);
end

%% Tikhonov Regularization with l1-penalry on a stemp-mean-dictionary
if run_tp_l1_dic == 1
    methodname = 'TP L1 DIC';
    if load_saved
        load(filename,'f_l1_dic','alpha_l1_dic'); f_opt = f_l1_dic; alpha_opt = alpha_l1_dic; runtime = 0;
    else
        tic; [f_opt,alpha_opt] = FISTA_DIC(L*A,L*y,alpha,0); runtime = toc;
    end
    f_l1_dic = f_opt; alpha_l1_dic = alpha_opt;
    show_barnett(methodname,f_opt,alpha_opt,b_index);
end

%% Tikhonov Regularization with l1-penalty on a stemp-mean-dictionary with positivity constraint
if run_tp_l1_dic_pos == 1
    methodname = 'TP L1 DIC POS';
    if load_saved
        load(filename,'f_l1_dic_pos','alpha_l1_dic_pos'); f_opt = f_l1_dic_pos; alpha_opt = alpha_l1_dic_pos; runtime = 0;
    else
        tic; [f_opt,alpha_opt] = FISTA_DIC_POS(L*A,L*y,alpha,0); runtime = toc;
    end
    f_l1_dic_pos = f_opt; alpha_l1_dic_pos = alpha_opt;
    show_barnett(methodname,f_opt,alpha_opt,b_index);
end

%% Standard inversion
if show_st_inv == 1
    methodname = 'standard inversion';
    load(filename,'f_st_inv'); f_opt = f_st_inv;
    show_barnett(methodname,f_opt,[],b_index);
end

%% Lagrange multiplier
if show_st_inv == 1
    methodname = 'Lagrange multiplier';
    load(filename,'f_lm'); f_opt = f_lm;
    show_barnett(methodname,f_opt,[],b_index);
end

%% Transform inversion
if show_st_inv == 1
    methodname = 'transform inversion';
    load(filename,'f_tr_inv'); f_opt = f_tr_inv;
    show_barnett(methodname,f_opt,[],b_index);
end

%% Gibbs sampler
if show_st_inv == 1
    methodname = 'Gibbs sampler';
    load(filename,'f_gibbs'); f_opt = f_gibbs;
    show_barnett(methodname,f_opt,[],b_index);
end