% Fig21.m
%
% This is a script devised to compare the gain matrices obtained various 
% sampling matrix.  The comparison is done by comparing the distribution of
% the elements of log10(|G|).
%
% It then plots the G matrix for 2 towers. In both cases, we use a full
% set of Haar wavelets on a 64 x 64 mesh.
%
% Jaideep Ray
%
% 2/20/2013
% -------------------------------------------------------------------------

close all ; clear ;
% -------------------------------------------------------------------------

%% Set paths

SPARSELABLIBPATH1 = '/home/jaray/Projects/ascr-graph/lib/SparseLab2.1-Core/';
Wavelab850Path    = '../../../../lib/Wavelab850/Orthogonal/' ;

%% do not change below.
addpath('../../sparse-inversion/') ;
addpath('../../msrf-selection/') ; % for wavelet transforms
addpath('../../incoherence/')    ; % for incoherent samplers
addpath(Wavelab850Path); % for MakeONFilter and wavelet utils
addpath(genpath(SPARSELABLIBPATH1));

%% 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 parameters

obsDir  = '../../data/tower.obs/' ;  % To get tower names
hDir    = '../../H-matrices/' ; % to get H matrices
nlightsFile  = '../../data/nlights-0.1x0.1/nightlights.dat' ;
nlightsLats  = '../../data/nlights-0.1x0.1/nightlights_lat.dat' ;
nlightsLongs = '../../data/nlights-0.1x0.1/nightlights_long.dat' ;
USBndryFile  = '../../data/US_coords.dat' ; % To chop-off canada in US nightlights
edgarDir = '../../data/edgar2005.1x1/' ;

% ---- Basic info about problem
info.nUSCells = 816 ;
info.nObsPerPeriod = 8*8 ;
info.nCompSamples = 10 ;

% ---- Time period over which info is required
buffer  = 2 ; % buffer periods before and after the duration of interest
periods = [21 22] ; % Want to tackle Days 241:272 i.e. 32 days
% periods = [3 43] ; % Want to tackle Days 241:248 i.e. 8 days

% ---- Sampling matrices
samplers = [ {'Gaussian'}, {'Circulant'}, {'Hadamard'}, {'sbHadamard'}, {'Noiselet'} ] ;

% ---- Plotting related
lc      = [{'-k'}, {'-b'}, {'-r'}, {'-g'}, {'-m'}];
bc      = [{'k'}, {'b'}, {'r'}, {'g'}, {'m'}];
sc1      = [{'ok'}, {'or'}, {'og'}, {'om'}, {'ob'},] ;
sc2      = [{'+k'}, {'+r'}, {'+g'}, {'+m'}, {'+b'},] ;
marker  = 'dw' ;
ms      = 8 ; % marker size for tower locations
mfc     = 'w'; % marger fill color
lw      = 2  ; % linewidth
fs      = 16 ; % fontsize
lfs     = 14 ; % legend fontsize

logtransform = true ;
visualize   = false ;
levelsToZero = 0 ;
% ---- End, user-settable parameters --------------------------------------

%% ---- Compute the Gain matrices using H and full wavelet bases
useNlights = true ;
[Gtowers, GCS, Phi64] = computeGainUsingH_v2(obsDir, hDir, info, periods, ...
                           buffer, useNlights, nlightsFile,nlightsLats, ...
                           nlightsLongs, USBndryFile, edgarDir, ...
                           levelsToZero, visualize) ;
% ---- Dump out incoherence of towers
signalSize = size(Phi64, 1) ;
for itower = 1 : length(Gtowers)
    % Find the max absolute element of G, for each tower
    tmp = Gtowers{itower} ;
    maxGscaled = sqrt(signalSize) * max( abs( nonzeros(tmp) ) ) ;
    fprintf(' Incoherence of tower %d = %f; max = %d\n', itower, maxGscaled, ...
             sqrt(signalSize)) ;
end
fprintf(' Incoherence of R PhiClower48 = %f\n', sqrt(signalSize) * max( abs( nonzeros(GCS) ) )) ;

                                      
%% ---- Compute Gain matrices using conventional samplers
ncols = size(Phi64, 1) ;
% Just generate a gain matrix that representative of 1 tower
nrows = info.nObsPerPeriod * (periods(2) - periods(1) + 1 + 2*buffer) ;
Gsampler = cell(0) ;
for count = 1 : length(samplers)
    
    % ---- Compute gain for this sampler
    mtype = samplers{count} ;
    PsiM  = constructSamplingMatrix(nrows, ncols, mtype) ;
    Gsampler{count} = PsiM * Phi64 ;

    % ---- Print info
    fprintf(' %s - Incoherence of PsiM * Phi64 = %f\n', mtype, ...
        sqrt(signalSize) * max( abs( nonzeros(Gsampler{count}) ) )) ;
end
  
%% ---- Compare the PDF of G_{ij}, for G_sampler
figNo = 0 ;
for count = 1 : length(samplers)
    
    % ---- Plot the PDF from the sampler. Logtransform, if necessary
    gs = abs( nonzeros(Gsampler{count}) ) ;
    if (logtransform == true)
        gs = log10( gs ) ;
    end
    if (count < 3 )
        [f, x] = ksdensity(gs) ;
        plot(x, f, lc{count}, 'LineWidth', lw) ;  hold on ;    
    else
        [f, x] = hist(gs, 30) ;
        f      = f / sum(f) ; % Probability mass function
        bar(x, f, 0.05, bc{count}) ; hold on;
    end
end

%% ---- Plot the G_{ij} for G_towers. But only 2 of them
towerNumbers = [2, 21] ;
for itower = 1 : length(towerNumbers)
    
    % Take out the G matrix for this tower
    tmp = Gtowers{ towerNumbers(itower)} ;
    
    % Take out the non-zeros. Log-transform, if necessary
    g = abs( nonzeros(tmp) ) ;
    if (logtransform == true)
        g = log10(g) ;
    end
    
    % Make the kernel density estimate and plot
    [f, x] = ksdensity(g) ;
    plot(x, f, sc1{itower}, 'MarkerSize', ms) ;  hold on ;
end

%% ---- Decorate plot
title(' Distribution of A_{i,j} from A_{\Psi} and A_H', 'FontSize', fs) ;
xlim([-7, 0]) ;
ylabel('Density', 'FontSize', fs) ;
if (logtransform == true)
    xlabel('log_{10}(|A_{i,j}|)', 'FontSize', fs) ;
else
    xlabel('(|A_{i,j}|)', 'FontSize', fs) ;
end
% plots = [samplers, {'Tower 2; full basis set'}, {'Tower 21; full basis set'}, ....
%                    {'Tower 2; reduced basis set'}, {'Tower 21; reduced basis set'},];
plots = [samplers, {'Tower 2'}, {'Tower 21'}] ;
h = legend(plots, 'Location', 'Best');
set(h, 'FontSize', lfs); legend boxoff ; legend boxoff ;
set(gca, 'FontSize', fs);

% Save to disk
print('-dpdf', '-painters', 'Fig21.pdf') ;

