function [ Gtowers, GCS, Phi64 ] = computeGainUsingH_v2( obsDir, hDir, info, ...
                                      periods, buffer, useNlights,  ...
                                      nlightsFile, nlightsLats, nlightsLongs, ...
                                      USBndryFile, edgarDir, levelsToZero, visualize)
%computeGainUsingH Compute gain and Phi using H matrices
%   [ Gtowers, GCS, Phi64 ] = computeGainUsingH( obsDir, hDir, info, ...
%                               periods, buffer, useNlights,  ...
%                                      nlightsFile, nlightsLats, nlightsLongs, ...
%                                      USBndryFile, edgarDir, levelsToZero, visualize)
% 
% This function returns the gain matrices for the H matrices. It can do so
% for a full wavelet basis set (Haars) or for those sparsified by
% nightlights.
%
% INPUT :
%  obsDir  : Directory with tower info
%  hDir    : Directory where the H matrices are kept
%  info    : inversion problem info, specifically
%            info.nUSCells      = 816 ;
%            info.nObsPerPeriod = 8*8 ;
%            info.nCompSamples  = 10 ;
%  periods : the 8-day periods of interest
%  buffer  : how much buffer to keep
%  useNlights  : Should I use nlights to sparsify the basis set?
%  nlightsFile : file with the radiances
%  nlightsLats, nlightsLongs: File with the lat/lons
%  USBndryFile : File with US boundary so that non-US lights can be zeroed
%  edgarDir    : edgar emissions so that we may also calculate a nightlight
%                emission model and calculate prior wavelet coeffs
%  levelToZero : keep at 0
%  visualize   : Do we plot the wavelets?
%
% OUTPUT :
%
%  Gtowers : Cell array, 35 long, contain arrays. Each array is
%            (nperiods * 64) x (nwavelets * nperiods) in size
%  GCS     : PsiM * PhiClower48, nCompSamples x (nwavelets * nperiods) in
%            size
%  Phi64   : The Haar bases over the dyadic grid.
%
% Jaideep Ray, 2/12/2013
% -------------------------------------------------------------------------

%% Get the tower names
[ ~, towerNames ] = getTowerObservations(obsDir, ...
                                [periods(1)-buffer, periods(2)+buffer], ...
                                info) ;
info.nTowers = length(towerNames) ;
info.buffer = buffer ; % buffer, in Periods, at the front and end of observations
info.nPeriods = periods(2) - periods(1) + 1 + 2*info.buffer;

%% Extract the Hmatrices for the same period
[Hmatrices, Hlatlon] = getHMatrices(hDir, ...
                                [periods(1)-buffer, periods(2)+buffer], ...
                                towerNames, info) ;

%% Construct a set of bases that provides emissions over US only, at
% full 1 deg x 1 deg resolution
tst = tic ;
if (useNlights == false )
    [Phi64, latlon] = constructBases() ;
else
    [Phi64, latlon, priorWeights] = ...
        constructBasesFromNlightEmissionModel(nlightsFile, nlightsLats, ...
               nlightsLongs, USBndryFile, edgarDir, levelsToZero, visualize);
   
end
fprintf('Constructing bases took %f secs\n', toc(tst));

% ---- Now split the basis set, once it's been scaled
[Phi64Lower48, latlonLower48, Phi64Clower48, latlonClower48] =  ...
                                    splitBasisSet(Phi64, latlon, Hlatlon) ;

%% Generate the random matrix for compressive sampling
R = genRandomMatrixForCompSampling(info.nCompSamples, size(Phi64Clower48, 1)) ;

%% Construct gain functions
disp('Computing gain matrices with H') ; tst = tic ;
[Gtowers, GCS] = computeAugmentedGainMatricesAfterNormalization(Hmatrices, ...
                        Phi64Lower48, R, Phi64Clower48, info) ;
fprintf('Gain matrices took %f secs to compute\n', toc(tst));

end

