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

% script to compare the estimates in local, regional and total measure

%% Load reconstructions

load('inverse_methods_solutions.mat')
load('synthetic_methane_data.mat','area','land_grids')

methods = {'st inv','tr inv','lm','gibbs','l2','l2 pos','l1','l1 pos','l1 dic','l1 dic pos'};
f = [f_st_inv,f_tr_inv,f_lm,f_gibbs,f_l2,f_l2_pos,f_l1,f_l1_pos,f_l1_dic,f_l1_dic_pos];

select = [1,2,3,4,5,6,7,8,9,10];
methods = methods(select);
f = f(:,select);

f_EDGAR = f_EDGAR.*area;
for k = 1:size(f,2)
    f(:,k) = f(:,k).*area;
end

N = length(select);

%% Calculate different measures

% 1.) Pixel by pixel comparison scaled by norm of solution
error_1 = zeros(size(select));
error_2 = zeros(size(select));
rel_error_1 = zeros(size(select));
rel_error_2 = zeros(size(select));

for k = 1:N
    error_1(k) = sum(abs(f_EDGAR - f(:,k)));
    error_2(k) = norm(f_EDGAR - f(:,k),2);
    rel_error_1(k) = error_1(k)/sum(abs(f_EDGAR));
    rel_error_2(k) = error_2(k)/norm(f_EDGAR,2);    
end
local_error = error_1;
rel_local_error = rel_error_1;

% 2) Integrated region pixel by pixel
reg_error_1 = zeros(1,N);
rel_reg_error_1 = zeros(1,N);

for k = 1:N
    reg_error_1(k) = sum(abs(local_int_f(f_EDGAR - f(:,k),fspecial('disk',2),0)));
    reg_error_2(k) = norm(local_int_f(f_EDGAR - f(:,k),fspecial('disk',2),0));
    rel_reg_error_1(k) = reg_error_1(k)/sum(local_int_f(f_EDGAR,fspecial('disk',2),0));
    rel_reg_error_2(k) = reg_error_2(k)/norm(local_int_f(f_EDGAR,fspecial('disk',2),0));
end


% 3) Total budget
budget_error = zeros(size(select));
rel_budget_error = zeros(size(select));

for k = 1:N
    budget_error(k) = sum(f_EDGAR) - sum(f(:,k));
    rel_budget_error(k) = budget_error(k)/sum(f_EDGAR);
end

%% Summarize measures

table(methods',budget_error',round(rel_budget_error',3),reg_error_1',round(rel_reg_error_1',3),local_error',round(rel_local_error',3), ...
      'VariableName',{'methods','total','rel_total','regional','rel_regional','local','rel_local'});
  
  
%% Compare results in each region

[regions,region_name] = get_regions(1);

regional_budget_opt = zeros(length(regions),N);
regional_budget_sol = zeros(length(regions),1);
for k = 1:length(regions)
    regional_budget_sol(k) = sum(f_EDGAR(regions{k}));
    for kk = 1:N
        regional_budget_opt(k,kk) = sum(f(regions{k},kk));
    end
end


% This one is breakyaxis (and some tricks for correct display)
figure('position',[500,480,625,280]);
hp = bar([regional_budget_sol,regional_budget_opt]);
set(gca,'XTick',1:length(regions),'XTickLabel',region_name);
set(gca,'YTick',[5000,10000,35000,40000,45000]);
xlim([0.5,length(regions)+0.5])
ylim([0,4.3e4])
h = breakyaxis([1.3,3.2]*10^4);
set(h.lowAxes,'ytick',[5000,10000],'yticklabel',[5,10]);
set(h.highAxes,'ytick',[35000,40000],'yticklabel',[35,40]);
hh = text(-0.05,1.03,'$k\ mol\ s^{-1}$','units','normalized');
set(hh,'interpreter','latex','fontsize',12);

hh = gcf;
uistack(hh.Children(3),'top');
axes(hh.Children(1));
hh = legend({'EDGAR',methods{1:end}});
set(hh,'color','w','position',[0.7020,0.4368,0.1967,0.4728]);
