function [ bases, latlon, weights ] = ...
    constructBasesFromNlightEmissionModel(nlightsFile, nlightsLatsFile, ...
                    nlightsLongsFile, USBndryFile, edgarDir, levelsToZero, visualize)
% constructBasesFromNlightEmissionsModel Construct MsRF, bases on Haars, fron
% nightlights, at 1 deg X 1 deg resolution
%
%   constructBasesFromNlightEmissionModel(nlightsFile, nlightsLatsFile, ...
%         nlightsLongsFile, USBndryFile, edgarDir, levelsToZero, visualize)
%
%   contructs a MsRF, based on Haar wavelets, for ffCO2 emissions at 1 deg
%   X 1 deg resolution. It does so using nightlights. We create an ffCO2
%   emissions field by Em = c * Nlights, where c is calculated by matching
%   the US-integrated Nlights with US-integrated EDGAR 2005. The wavlets
%   retained in the MsRF can model ffCO2 emissions over US - and leaves out
%   the ocean.
%
%   The bases (wavelets) come with priors weights. BAsically, the bases
%   cover all 816 1 deg squares that cover US. The grid-cells that have a
%   large prior weights are one where we have nightlights.
%
%  PROCEDURE
%
%    1. We consider the nighlights, and scale them so max nightlight = 100
%    2. We zero-pad the nightlights to be dyadic.
%    3. We wavelet trasnform it. The results is a dyadic matrix. We zero
%       out small wavelet coeffs (in constructNoisyField, with noise = 0.0)
%    4. we do an inverse wavelet transform and crop out the 27 x 64 matrix,
%       which is the reconstructed nighlights minus small wavelets
%    5. We zero out non-US areas (oceans)
%    6. We scale the nightlights to match EDGAR 2005 emissions
%    7. We use this as the "mask"  as input to creation of a masked wavelet 
%       coeff structure (MWCS). This identifies the non-ocean wavelets
%    8. We loop over each level and each wavlet in the empty MWCS, and set
%       wavelet to 1. We recreate the field due to it, and that makes it
%       the spatial base to stick into 'bases'. It also returns the weights
%       of the scaled up nightlight field that models emissions. These are
%       the prior weights.
%
%  INPUTS :
%   nlightsFile : a filename that has nightlights at 1deg X 1deg resoution
%                 in the domain (corners) [51N, 127W]  [24, 63W]. This is a
%                 27 x 64 matrix.
%   nlightsLatsFile : full filename, of file containing cell-centers 
%                     latitudes of nighlight cells
%   nlightsLongsFile : ditto, but longitudes
%   USBndryFile      : lat-long of US bndry, lower 48 states
%   edgarDir         : directory with EDGAR info
%   levelsToZero     : zero out the wavelets at the finest 'levelsToZero'
%                      levels.
%   visualize        : if true, show plots
%
%  OUTPUTS :       
%   bases   : a maxtrix, ncellsDomain X nRetainedWavelets in size.
%             ncellsDomain is the # of pixels in the domain of interest with 
%             cell-centers at [88.5N, -126.5W] and [24.5N, -63.5W] with a
%             cell-centered 64 x 64 mesh i.e., 4096 cells.
%             nRetainedWavelets is the size of the set of wavelets needed to model
%             emissions over US.
%
%   latlon  : ncellsDomain x 2 matrix, containing the latlon of the cells
%             with the emissions/fields given by bases * weights
%
%   weights : the wavelets weights, for the wavelets retained in 'bases'
%             that model the emissions using Em = c * nightlights model
%
% Jaideep Ray, 7/25/2012
%
% =========================================================================

%% Make the nightlight field
% ---- Read in nightlights 
[~, nlights    ] = hdrload(nlightsFile) ;
[~, nlightsLats]  = hdrload(nlightsLatsFile) ;
[~, nlightsLongs] = hdrload(nlightsLongsFile) ;
[~, borderCoords] = hdrload(USBndryFile);

% ---- Construct a version of nightlights, that has no nightlights in
% the ocean. noiseLevel = 0.0
noiseLevel = 0.0 ;
noisyNlights_27x64 = constructNoisyField(nlights, nlightsLats, nlightsLongs, ...
                                   borderCoords, noiseLevel) ;

% ---- Insert the noisyNightlights into a domain which is discrtized with 
% a square dyadic grid i.e., 64 x 64, with the same spatial resolution as 
% noisyNlights
targDomainTC = [87.5, -126.5] ; % top and bottom corners of the target domain
targDomainBC = [24.5, - 63.5] ;
[noisyNlights, lat2d, lon2d] = insertField(noisyNlights_27x64, nlightsLats, ...
                                nlightsLongs, targDomainTC, targDomainBC) ;
latlon = [lat2d(:), lon2d(:)] ;                            
                                 
% ---- Show the noisy nightlights
if (visualize == true)
    imagesc(noisyNlights) ; colorbar ; title('Noisy Nightlights'); 
    pause(2) ; % Pause 2 second
end

%% Scale up the nightlight field to match EDGAR 2005 total emissions
% ---- Get EDGAR2005 emissions 
topcorner = [ max(nlightsLats(:)), min(nlightsLongs(:)) ] ;
botcorner = [ min(nlightsLats(:)), max(nlightsLongs(:)) ] ;
[edgar2005, latlonEDGAR] = prepareEDGAR(edgarDir, topcorner, botcorner, ...
                                        false) ;
if  (visualize == true)
    figure(); 
    pcolor(latlonEDGAR(:, :, 2) - 360, latlonEDGAR(:, :, 1), edgar2005) ; 
    colorbar ; caxis([0, 5]); 
    xlabel('Longitude') ; ylabel('Latitude') ;
    title('EDGAR 2005 emissions; micromoles of CO2/m^2/sec') ;   
end

% ---- Scale up nlights to meet EDGAR emissions
totalEm = sum( edgar2005(:) ) ; totalNl = sum(noisyNlights(:)) ; 
c = totalEm / totalNl ;
emissionsFromNlights = c * noisyNlights ;

if (visualize == true)
    figure() ;
    pcolor(lon2d, lat2d, emissionsFromNlights) ; colorbar ;
    caxis([0, 0.8*max(emissionsFromNlights(:))]) ;
    title(' Emissions from linearly scaled nightlight radiance; micromoles of CO2/m^2/sec') ;
    xlabel('Longitude') ; ylabel('Latitude');
end

%% ---- Construct bases and priors
[bases, weights] = constructBasesAndPriors(emissionsFromNlights, levelsToZero) ;

% ---- Loop over bases and plot them
if (visualize == true)
    % Find the corners of the domain and add to points to be plotted
    westC  = min( lon2d(:) ) ;   eastC  = max( lon2d(:) ) ;
    northC = max( lat2d(:) ) ;   southC = min( lat2d(:) );
    Xpoints = [borderCoords(:, 1);  westC;  eastC] ;
    Ypoints = [borderCoords(:, 2); northC; southC] ;

    % Loop over bases
    figure() ;
    for ibase = 1 : size(bases, 2)
        base = bases(:, ibase) ;
        base2D = reshape(base, size(noisyNlights, 1), size(noisyNlights,2)) ; 
        pcolor(lon2d, lat2d, base2D) ; colorbar ; hold on
        scatter(Xpoints, Ypoints, 40, 'w') ;
        title('Basis for emissions');
        pause(1) ; % Pause a second
        hold off ;
    end
end
end % End of function
