% Fig23.m
%
% The aim of this script is to create plots of sparsity as a function of
% order for Haars, daubechies, symmlet and coiflet, at level 4 and 6.
%
% 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. This will 
% motivate sparse modeling.
%
% 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
frac      = 1.0e-3    ; % Threshold for determining smallness of wavelet coeffs
waveletFamily   = {'Haar','Daubechies','Symmlet','Coiflet'} ; % [help MakeONFilter for possible names]
waveletOrder{1} = 1 ; % [help MakeONFilter for possible wavelet orders]
waveletOrder{2} = 4:2:20 ; 
waveletOrder{3} = 4:2:10 ;
waveletOrder{4} = 1:1:5 ;
levels = [4, 6] ;
ls = [{'-ok'}, {'-or'}, {'-ob'}, {'-om'}, {'-og'}] ; % Line styles to use when plotting
fs = 16 ; % Font-size
lw = 2  ; % Line-width
ofile = 'Fig23.pdf';
% ofile = 'Fig23.jpg';

% ==== 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(waveletFamily) ;

%% Do a wavelet transform of the average source
figure() ;

% Loop over families of wavelets
for iwavelet = 1 : nwavelets
    
    % Find wavelet family and orders
    waveletName = waveletFamily{iwavelet};
    myOrders = waveletOrder{iwavelet} ;
    sparsity = zeros(length(myOrders),2) ;
    
    % Loop over orders
    for iorder = 1 : length(myOrders)
        order = myOrders(iorder) ;
        fprintf(' Doing %s of order %d\n', waveletName, order) ;
        
        % Do a wavelet transform of the average source
        wcs = createEmptyWCS(dyadicGrid, waveletName, order) ;
        wcs = computeWaveletTransform(dyadicGrid, wcs) ;
        
        for l = 1 : length(levels)
            % Extract wavelets from level level, remove very small ones, and
            % count the number of non-zeros, mean and average of the non-zero wavelet
            % coefficients
            wcv = extractWaveletCoeffs(wcs, levels(l)) ;
            
            % Absolute value for thresholding. Set all wavelet coeffs below
            % threshold to zero
            awcv = abs(wcv) ;  eps  = frac * max(awcv) ; inds = awcv < eps ;
            wcv(inds) = 0.0 ;
            
            tmp = nonzeros(wcv) ; % Extract all non-zeros from the vector
            sparsity(iorder, l) = 1 - length(tmp) / (4^levels(l) - 4^(levels(l)-1));
            fprintf(' Doing level = %d, order %d, sparsity = %f\n', levels(l), ...
                order, sparsity(iorder, l)) ;
        end
    end
    ax1 = subplot(1,2,1) ;
    plot(ax1, myOrders, sparsity(:, 1), ls{iwavelet}, 'LineWidth', lw) ; hold on ;
    ax2 = subplot(1, 2, 2) ;
    plot(ax2, myOrders, sparsity(:, 2), ls{iwavelet}, 'LineWidth', lw) ; hold on ;
    fprintf('\n') ;
end

% Put in the title for each plot
title(ax1, 's = 4', 'FontSize', fs-2) ;
title(ax2, 's = 6', 'FontSize', fs-2) ;

% Ditto x, y labels
xlabel(ax1, 'Order', 'FontSize', fs-2) ;
xlabel(ax2, 'Order', 'FontSize', fs-2) ;
ylabel(ax1, 'Sparsity', 'FontSize', fs-2) ;
ylabel(ax2, 'Sparsity', 'FontSize', fs-2) ;

% Legend
legend(ax1, 'Haar','Daubechies','Symlet','Coiflet', 'Location', 'Best') ; legend(ax1,'boxoff') ;
legend(ax2, 'Haar','Daubechies','Symlet','Coiflet', 'Location', 'Best') ; legend(ax2, 'boxoff');

%% Save to file
print('-dpdf', '-painters', ofile) ;
% print('-djpeg', '-painters', ofile) ;


