function [Gtowers, GCS] = ...
    computeAugmentedGainMatricesAfterNormalization(Hmatrices, ...
                                      Phi64Lower48, R, Phi64Clower48, info) 
%computeAugmentedGainMatricesAfterNormalization() Computes the gain matrix for towers and
%compressive sampling
%   [Gtowers, GCS] = computeAugmentedGainMatricesAfterNormalization(Hmatrices, ...
%                                      Phi64Lower48, R, Phi64Clower48, info) 
%  Computes the gain matrices for all towers and compressive sampled
%  emissions. A gain matrix is such that y^(I) = G^(I) w, where y^(I) is a
%  vector containing the [CO2] concentrations measured at tower I over a
%  whole year, G^(I) = H^(I) * Phi64Lower48 is the gain matrix for tower I
%  and w are wavelet weights. For the "fake" tower, which is really a
%  compressive sampling of emissions overs canada, mexico and oceans,
%  G^(II) = R * Phi64Clower48 and 0 = R * Phi64Clower48 * w.
%
%  INPUTS :
%   Hmatrices : cell array, 35 long, for the H matrices for 35 towers
%   Phi64Lower48 : 816 x 1054 martrix, mapping wavelet coeffs to emissions
%   in lower48 states
%   R         : random matrix to act as H matrix to samples emissions
%               outside lower48 states
%   Phi64Clower48 : bases vector for gridcells outside lower48 states
%   info : struct with some info about # of cells inside US, # of
%                 observations per periods etc
%
%  In this code, we make sure we normalize the H and Phi64Lower48 and
%  Phi64Clower48 matrices to compute incoherence.
%
% Jaideep Ray, 2/11/2013
% =========================================================================

%% Perform normalization
%
% H matrices first
% The H matrix, for a given tower, for, say 12 8-day periodsm has the form
% H, for 12 8-day periods is
% H = [ H_11  H_12 .... H_1,12
%       H_21  H_22      H_2,12
%       H_31            H_3,12
%
%       H_12,1          H_12,12]
%
% where  H_i,j is  nObsPerPeriod x nUSCells in size.
%
% Our job is makse sure that each row of H has magnitude 1
%
for itower = 1 : length(Hmatrices)
   tmp = Hmatrices{itower} ;
   norms = sqrt( sum(tmp.^2, 2) ) ; % Normalize each row
   invnorms = 1 ./ norms ;
   diagInvNorms = diag(invnorms) ;
   Hmatrices{itower} = diagInvNorms * tmp ;
end

% Phi64 next. We'll have to concatenate Phi64Lower48 and Phi64Clower64
Phi64         = [Phi64Lower48 ; Phi64Clower48] ; % Concatenate the Phi64 matrix
norms         = sqrt( sum(Phi64.^2, 1) ) ;  % Normalize each column
invnorms      = 1 ./ norms ;
diagInvNorms  = diag(invnorms) ;
Phi64         = diagInvNorms * transpose(Phi64) ;
Phi64         = transpose(Phi64) ;
nrowsLower48  = size(Phi64Lower48, 1) ;
Phi64Lower48  = Phi64(1 : nrowsLower48, :) ; % Disaggregate the Phi64 matrix
Phi64Clower48 = Phi64( (nrowsLower48+1) : end, :) ;

% Don't bother with normalizing R; it's already normalized

%% Compute the gain matrix for the towers
prob.H    = Hmatrices ;
prob.Phi  = Phi64Lower48 ;
prob.info = info ;
Gtowers = computeGain(prob) ;

% Compute the gain matrix for 1 period
tmp = R * Phi64Clower48 ;

% Replicate it for nperiods
nperiods = size(Hmatrices{1}, 2) / info.nUSCells ;
GCS     = repmat(tmp, 1, nperiods) ;
GCS = sparse(GCS) ;
end

