% Fig02.m
%
% The aim of this script is to generate figures that show Haars give the
% sparsest representation of Vulcan.
%
% 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 sparsity of wavelet coefficients on a level-by-level basis
%   (fraction non-zero on each level)
%   - The average magnitude of wavelet coeffs on each level
%
% This is done for Haars, Daubechies4, Daubechies8, Symmlet4, Symmlet6 in 2
% plots with multiple lines.
%
% Jaideep Ray, 12/17/2012
% -------------------------------------------------------------------------

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
frac      = 1.0e-3    ; % small wavelet coeffs are frac*max(abs(coeff)) ;
waveletName  = {'Haar', 'Daubechies', 'Daubechies', 'Symmlet', 'Symmlet'} ; % [help MakeONFilter for possible names]
waveletOrder = [2, 4, 8, 4, 6] ; % [help MakeONFilter for possible wavelet orders]
ls = [{'-k'}, {'-r'}, {'-b'}, {'-m'}, {'-g'}] ; % Line styles to use when plotting
ms = [{'--k'}, {'--r'}, {'--b'}, {'--m'}, {'--g'}] ; % Line styles to use when plotting
fs = 16 ; % Font-size
lw = 2  ; % Line-width

% ==== 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) ;
nonZeroCoeffs = zeros(nlevels, nwavelets) ;
meanCoeff     = nonZeroCoeffs ;
sdCoeff       = nonZeroCoeffs ;

%% Do a wavelet transform of the average source

for iwavelet = 1 : nwavelets
    
    % Print what we are doing
    fprintf(' Doing %s of order %d\n', waveletName{iwavelet}, ...
            waveletOrder(iwavelet)) ;
        
    %% Do a wavelet transform of the average source
    wcs = createEmptyWCS(dyadicGrid, waveletName(iwavelet), waveletOrder(iwavelet)) ;
    wcs = computeWaveletTransform(dyadicGrid, wcs) ;
    
    %% Loop over all levels. Extract wavelets, remove very small ones, and
    %  count the number of non-zeros, mean and average of the non-zero wavelet
    %  coefficients
    for ilevel = 1 : nlevels
        % Extract wavelet coeffs
        wcv = extractWaveletCoeffs(wcs, ilevel) ;
        
        % Find the small ones and zero them out. Small is defined as a fraction
        % of the max value
        awcv      = abs(wcv) ;
        eps       = frac * max(awcv) ;
        inds      = awcv < eps ;
        wcv(inds) = 0.0 ;
        fprintf(' Doing level %d, max = %f, min = %f, nonzeros = %d\n', ...
            ilevel, max(wcv), min(wcv), sum(inds));
        
        % Collect stats on them
        tmp = nonzeros(wcv) ; % Extract all non-zeros from the vector
        nonZeroCoeffs(ilevel, iwavelet) = length(tmp) / (4^ilevel - 4^(ilevel-1));
        meanCoeff(ilevel, iwavelet)     = mean(tmp)   ;
        sdCoeff(ilevel, iwavelet)       = std(tmp)    ;
    end
    fprintf('\n') ;
end

%% Plot the results. Figure 1 is sparsity versus level for all wavelets.
%  Figure 2 has 2 y-axes. On the left, mean non-zero wavelet coeff versus
%  levels; on the right y-axes, sd of non-zero wavelet coeffs.
f = figure(1) ;
for iwavelet = 1 : nwavelets
    plot(1:nlevels, nonZeroCoeffs(:,iwavelet), ls{iwavelet}, 'LineWidth', lw)  ;
    hold on ;
end
hold off;
c = get(f, 'children') ;  set(c, 'FontSize', fs) ; 
axis([1, 6, 0.2, 1.2]) ;
set(gca, 'XTick', 1:nlevels) ;
set(gca, 'XTickLabel', 1:nlevels) ;
xlabel('Scale', 'FontSize', fs) ; 
ylabel('Fraction of non-zeros', 'FontSize', fs) ;
title('Sparsity of various wavelet choices', 'FontSize', fs) ;
legend('Haar', 'Daubechies 4', 'Daubechies 8', 'Symlet 4', 'Symlet 6', ...
       'Location', 'NorthEast') ; legend boxoff ;
print('-dpdf', '-painters', 'Fig02a.pdf') ;

f = figure(2) ;
for iwavelet = 1 : nwavelets
    plot(1:nlevels, meanCoeff(:,iwavelet), ls{iwavelet}, 'LineWidth', lw) ;
    hold on ;  
end
c = get(f, 'children') ;  set(c, 'FontSize', fs) ; 
axis([1, 6, -1, 2]) ;
legend('Haar', 'Daubechies 4', 'Daubechies 8', 'Symlet 4', 'Symlet 6', ...
       'Location', 'NorthEast') ; legend boxoff ;
   
for iwavelet = 1 : nwavelets
    plot(1:nlevels, sdCoeff(:,iwavelet), ms{iwavelet}, 'LineWidth', lw) ;
    hold on ;  
end  
hold off;
set(gca, 'XTick', 1:nlevels) ;
set(gca, 'XTickLabel', 1:nlevels) ;
xlabel('Scale', 'FontSize', fs) ; 
ylabel('Mean and standard deviation', 'FontSize', fs) ;
title('Statistics of non-zero wavelet coefficients', 'FontSize', fs) ;
legend('Haar', 'Daubechies 4', 'Daubechies 8', 'Symlet 4', 'Symlet 6', ...
       'Location', 'NorthEast') ; legend boxoff ;
print('-dpdf', '-painters', 'Fig02b.pdf') ;
