% Fig24.m
%
% The aim of this script is to compute the total emissions error and
% correlation between true and sparsified emissions at different threshold
% levels. We model emissions with Haars and sparsify.
%
% This is a script that reads in 8-day-averaged, 1-deg mesh Vulcan emissions in
% ../../data/vulcan.8day.1x1/EightDayAvSources*dat. It then does
% wavelet transform on it,  so that we can see the wavelet weights. We then
% loop through a number of thresholds for wavelet weights (below which we
% set them to zero) and see how the errors rise.
%
%
% Jaideep Ray, 07/01/2013
% -------------------------------------------------------------------------

close all; clear ;

%% ---- User-settable paths to toolboxes

Wavelab850Path    = '../../../../lib/Wavelab850/Orthogonal/' ;

%% ---- Do not change path below
addpath('../../sparse-inversion/') ; % For getEmissionsFromFile
addpath('../../msrf-selection/')   ; % For createEmptyWCS(), computeWaveletTransform() 
addpath(Wavelab850Path); % for MakeONFilter and wavelet utils

%% Legal splash screen
fprintf('Copyright 2014 Sandia Corporation. Under terms of Contract DE-AC04-94AL85000\n') ;
fprintf('there is a non-exclusive license for the use of this work by or on behalf of the U.S.\n');
fprintf('Government.\n');


%% ==== User-settable variables

vDir      = '../../data/vulcan.8day.1x1/' ;% Directory with 1-deg, 8-day vulcan
periods   = [1 45]                               ; % Which periods are we intested in
dyadicSq  = [64, 64]  ; % The dyadic square grid to do wavelet business on
waveletName  = 'Haar' ; % [help MakeONFilter for possible names]
fs = 16 ; % Font-size
lw = 2  ; % Line-width
ofile = 'Fig24.pdf' ;
% ==== End, user-settable variables 

%% Read in all EightDayAverage Vulcan sources and get an yearly average on
%  a power-of-2 square grid
[emissions, lats2D, longs2D] = getEmissionsFromFile(vDir, periods) ;
annualAv   = mean(emissions, 2) ; % Annuanl average.
annualAv   = reshape(annualAv, size(lats2D, 1), size(lats2D,2));
nr         = size(annualAv, 1) ; 
nc         = size(annualAv, 2) ;
dyadicGrid = zeros(dyadicSq(1), dyadicSq(2)) ;
dyadicGrid(1:nr, 1:nc) = annualAv ; 

%% Data structs to save data in 
nlevels       = log2(size(dyadicGrid, 1)) ;
nwavelets     = length(waveletName) ;

%% Do a wavelet transform of the average source, and zero-out coeffs below
%  a certain thresold. Reconstruct and check accuracy of reconstruction

% Threshold levels of interest
logFrac = -6 : 0.5 : -0.0 ;
data = zeros( length(logFrac), 4) ;

% Loop over threshold rations
for ithresh = 1 : length(logFrac)
    
    % Print what we are doing
    frac = 10^logFrac(ithresh) ;
    fprintf(' Doing log10(threshold) = %f\n', frac) ;
        
    % Do a wavelet transform of the average source
    wcs = createEmptyWCS(dyadicGrid, waveletName, 2) ;
    wcs = computeWaveletTransform(dyadicGrid, wcs) ;
    
    % Loop over all levels. Extract wavelets, remove very small ones and
    % stick back in
    nz = 0 ; 
    for ilevel = 1 : nlevels        
        wcv = extractWaveletCoeffs(wcs, ilevel) ;
        awcv = abs(wcv) ; eps = frac * max(awcv) ; inds = awcv <= eps ;
        wcv(inds) = 0.0 ;
        wcs = insertWaveletCoeffs(wcs, wcv, ilevel)  ;
        nz  = nnz(wcv) + nz ; % Running sum of non-zero wavelets
    end
    dyadicGridRecons = inverseWaveletTransform( wcs, nlevels ) ;
    annualAvRecons   = dyadicGridRecons(1:nr, 1:nc) ;
    
    % Metrics of reconstruction accuracy
    sparsity         = 1.0 - double(nz) / 4.0^nlevels ;
    % emRatio          = sum( annualAvRecons(:) ) / sum( annualAv(:) ) ;
    reconsErr          = norm(annualAvRecons(:) - annualAv(:)) / norm( annualAv(:) ) ;
    rho              = corr(annualAvRecons(:), annualAv(:)) ;
    data(ithresh, :) = [frac, sparsity, reconsErr, rho] ;
    fprintf('Threshold = %e, sparsity = %f, reconsErr = %f, rho = %f\n', ...
            frac, sparsity, reconsErr, rho) ;
end

%% Plot and save to file
[ax, h1, h2] = plotyy(data(:, 1), data(:, 2), data(:, 1), data(:, 4)) ; hold on ;
set(h1, 'Linewidth', lw, 'Linestyle', '--', 'Marker', 'x') ;
set(h2, 'Linewidth', lw, 'Linestyle', '-', 'Marker', 'o') ;
set(ax, 'Xscale', 'log') ;
set(ax(1), 'ylim', [0, 1.1], 'ytick', [0 0.2, 0.4, 0.6, 0.8, 1.0]) ;
set(ax(2), 'ylim', [0, 1.1], 'ytick', [0 0.2, 0.4, 0.6, 0.8, 1.0]) ;
xlabel('Threshold fraction \alpha', 'Fontsize', fs) ;
ylabel(ax(1), 'Sparsity & reconstruction error', 'Fontsize', fs) ;
ylabel(ax(2), 'Pearson correlation', 'Fontsize', fs) ;
title('Reconstruction fidelity versus sparsity', 'Fontsize', fs) ;
plot(ax(1), data(:,1), data(:, 3), '-bo', 'Linewidth', lw) ;

% Legend and save
h = legend('Sparsity','Reconstruction error', 'Correlation', 'Location', ...
           'West') ; legend boxoff ;
set(h, 'Fontsize', fs) ;       
print('-dpdf', '-painters', ofile) ;
