function [bases, prior] = constructBasesAndPriors(field, levelsToZero)
%constructBasesAndPriors() Construct bases and prior beliefs about their
% weights from a given field
%
%    1. We use this field as the "mask" input to a masked wavelet coeff 
%       structure (MWCS). This identifies the non-ocean wavelets
%    2. We loop over each level and each wavlet in the empty MWCS, and set
%       wavelet to 1. We recreate the emission due to it, and that makes it
%       the spatial base to stick into 'bases'
%    3. At the end of it all, we do a masked wavelet transform of the field
%       and extract the weights into a vector. These are the priors.
%    4. We normalize bases to be unit magnitude. We multiply priors
%       accordingly
%
% INPUTS:
%  field : a 2D field. does not have to be dyadic
%  levelsToZero : How many levels, from the finest, where we can zero out
%  bases.
%
% OUTPUTS:
%   bases : a matrix ncells X nbases in size. ncells = numel(field). Each
%           column is a base, is unit magnitude, and orthogonal to every 
%           one else. 
%   prior : a vector nbases in size.
%
% Jaideep Ray, 5/18/2012
% =========================================================================

% ---- Make its masked wavelet transform
small = 1.0e-6 ;
mwcs = createEmptyMWCS(field, 'Haar', 2,  field, small) ;
mwcs = computeMaskedWaveletTransform(field, mwcs) ;

% ---- Zero out all wavelets. Also, count them
numActiveLevels = mwcs.nlevels - levelsToZero ;
prior = mwcs.baseVal ; % first entry to priors
mwcs.baseVal = 0.0 ;
for l = 1 : mwcs.nlevels 
    mwcv = extractMaskedWaveletCoeffs(mwcs, l) ;
    if ( l <= numActiveLevels )
        prior = [prior; mwcv] ;
        fprintf('constructBasesAndPriors() : Level = %d  num active wavelets = %d\n', ...
                 l, numel(mwcv)) ;
    end
    mwcv = zeros(numel(mwcv), 1) ;
    mwcs = insertMaskedWaveletCoeffs(mwcs, mwcv, l) ; 
end
numParams = numel(prior) ;
fprintf('constructBasesAndPriors() : Number of free wavelet coefficients = %d\n', ...
         numParams) ;

% ---- Construct the matrix to hold the bases. It'll have as many rows as
% there are grid-cells - 27x64, and as many colms as there are bases (+1
% since we have to account for baseVal too)
bases = zeros( numel(field), numParams) ;

% ---- Loop over all weights in the the MWCS. Set each to one, reconstruct
% the emission field and save the base.
% ---- Turn on baseVal
mwcs.baseVal = 1.0;
E = inverseWaveletTransform(mwcs, mwcs.nlevels) ;
mwcs.baseVal = 0.0;     % Reset the baseVal ;

% Save the first base and increment paramNo
paramNo = 0 ;  bases(:, paramNo+1) = E(:) ;  % save field in W matrix
paramNo = paramNo + 1 ;

for l = 1 : numActiveLevels
    
    % Extract all wavelet coeffs on this level, and zero all out
    mwcv = extractMaskedWaveletCoeffs(mwcs, l) ;
    mwcv(:) = 0.0 ;      
    
    % For each of the wavelet coeffs on this level, set to 1 and generate
    % field
    for p = 1 : numel(mwcv)
        mwcv(p) = 1.0 ;       % Set the particular one to 1
        mwcs = insertMaskedWaveletCoeffs(mwcs, mwcv, l) ; % Stick it in.
        E = inverseWaveletTransform(mwcs, mwcs.nlevels) ; % Reconstruct emissions
        
        % Stick into W matrix, after normalizing
        E = E/norm(E(:)) ;
        bases(:, paramNo+1) = E(:) ;
        prior(paramNo+1) = norm(E(:)) * prior(paramNo+1) ; % adjust prior

        % Increment parameter count, reset mwcs
        mwcv(p) = 0.0 ;
        mwcs = insertMaskedWaveletCoeffs(mwcs, mwcv, l) ;
        paramNo = paramNo + 1 ;
        
    end
end


