function [ y ] = projectionsUsingAugmentedGain(mode, m, n, x, I, dim, UserData)
%projectionUsingAugmentedGain Computes y = A(:,I)*x if mode = 1, or y = A(:,I)'*x if
%mode = 2
%   [ y ] = projectionUsingGain(mode, m, n, x, I, dim, UserData) computes
%   the projections y = A(:,I)*x or y = A(:,I)'*x. Here x are the wavelet
%   coeffs (for entire period) and I denotes the coeffs that have been
%   selected by the StOMP in this round of the iteration. 
%   
%   INPUTS :
%   mode : 1 or 2, depending on whether we'll use the gain or gain' matrix
%   m    : m by dim is the implicit A matrix. Also, length(y) = m in mode 1
%   n    : The number of wavelet coeffs selected in this iter
%   I    : Indices of the wavelet coeffs selected. length(I) = n
%   x    : Selected wavelet coeffs, length(x) = n, if mode 1. Observation
%          residuals, length(x) = m, if mode = 2.
%   dim  : dimensionality of the problem. Total number of wavelet
%          coefficients, which need to be sparsified.
%   UserData : Struct, made by the user, and passed into this function
%              untouched by StOMP.
%  
%   OUTPUT :
%   y    : if mode = 1, co2 concetration predictions corresponding to the
%          guessed, sparse wavelet coeffs sent in. length(y) = m
%          if mode = 2, projections of the residuals onto the basis set
%          represented by the gain matrix. length(y) = n
%
% -------------------------------------------------------------------------

%% Check whether UserData has been passed in
global PROB_GLOBAL ;
if (nargin < 7)              % If UserData not passed in
    UserData = PROB_GLOBAL ; % Use the global version
end

%% Checks

% max value of I should be <= ncoeffsPerPeriod * nPeriods
assert( max(I(:)) <= (size(UserData.Phi, 2)*UserData.info.nPeriods) ) ;

% dimensionality should be equal to number of wavelet coeffs
assert( dim == (size(UserData.Phi, 2)*UserData.info.nPeriods) ) ;

%% Work

if ( mode == 1 ) % Given wavelet coeffs, return co2 concentrations

    % ---- m should be equal to the number of observations
    assert(m == ...
        UserData.info.nObsPerPeriod*UserData.info.nPeriods*UserData.info.nTowers + ...
        UserData.info.nCompSamples) ;

    % ---- number of coeffs should be less than those we started with
    assert( length(x) <= (size(UserData.Phi, 2)*UserData.info.nPeriods) ) ;

    % ---- Allocate output
    y = zeros(m, 1) ; % concetrations over nTowers and nPeriod

    % ---- Loop over towers
    for itower = 1 : UserData.info.nTowers
        G = UserData.Gtowers{itower} ;
        Gpart = G(:, I) ;
        tmp = Gpart * x ; % concentrations at this tower due to indicated coeffs for nPeriods
        len = length(tmp) ;
        iend = itower * len ;
        y(iend - len + 1 : iend) = tmp ;
    end
    
    % --- Now for the compressively sensed samples of fluxes in non-Lower48
    Gpart = UserData.GCS(:, I) ;
    y((end - UserData.info.nCompSamples + 1) : end) = Gpart * x ;
    return ;
end

if ( mode == 2 ) % Given residuals at each tower, return wavelet coeffs
    y = zeros(n, 1) ; % These will contain wavelet coeffs
    
    % ---- number of residuals == number of observations 
    assert(length(x) <= ....
        (UserData.info.nObsPerPeriod*UserData.info.nPeriods*UserData.info.nTowers + ...
         UserData.info.nCompSamples)) ;

    nobsPerTower = UserData.info.nObsPerPeriod * UserData.info.nPeriods ;
    for i = 1 : UserData.info.nTowers
        towerObs = x(i*nobsPerTower - nobsPerTower + 1 : i*nobsPerTower) ;
        G = UserData.Gtowers{i} ;
        Gpart = G(:, I) ;
        tmp = transpose(Gpart) * towerObs ; % tmp contains contribs to wavelet coeffs from each tower
        y = y + tmp ;
    end
    
    % ---- Finally, additions to wavelets coeffs due to compressive
    % samples. 
    obs     = x( (end - UserData.info.nCompSamples + 1) : end );
    Gpart   = UserData.GCS(:, I) ;
    tmp     = transpose(Gpart) * obs ;
    y       = y + tmp ;
    return ;
end

end


