function [ Phi64 ] = scaleByPriorWeights(weights, priorWeights, Phi64) 
%scaleByPriorWeights Scale Basis sets by prior weights
%   scaleByPriorWeights(weights, priorWeights, Phi64) scales the basis set
%   Phi64 by prior weights. The reason is as follows
%
%   StOMP solves the problem min L1(x), such that L2(y - Ax) < epsilon. Let
%   us say that we have a prior belief that x ~ x_pr. Then we can write the
%   same as min L1(z) such that L2(y - x_pr ** A z) < eps, z = x / x_pr,
%   were x_pr ** A means each element of x_pr multiplies a column of A. In
%   such a case, if x_pr(i) is small and z(i) is important, then z(i) has
%   to be very large indeed; L1 will try to remove it very quickly. x_pr
%   becomes like a normalization for x
%
%   x_pr should be wavelet weights of flux. Unfortunately, we have wavelet
%   weights of nighlights, which have different units. Somehow we have to
%   do the unit conversion.
%
%  INPUTS :
%    weights : wavelet weights, corresponding to basis set Phi64, for an
%              emission field. This is a vector size(Phi64, 2) long
%    priorWeights : Weights from a wavelet transform of nightlights
%    Phi64        : Basis set
%
%  OUTPUTS
%    Phi64 : the scaled Phi64 i.e. x_pr ** Phi64
%    
%  Jaideep Ray, 5/31/2012
% =========================================================================

% % %tmp = diag( (abs(priorWeights)) ) * transpose(Phi64) ;
tmp = diag( 0.1*(abs(priorWeights)).^(0.75) ) * transpose(Phi64) ;
Phi64 = transpose(tmp);

% % Propose a model w_e = (w_n/C)^a, where w_e are absolute values of
% % emission wavelets, w_n are the counterparts from nightlights and a is an
% % exponent. we need to fit a model log(w_e) = a log(w_n) -a log(C)
% % or y = ax + b
% 
% % Find those wavelets which matter
% aw = abs(weights) ;
% an = abs(priorWeights) ;
% ratio = an ./ aw ; 
% ii = find( (ratio > 1) & (ratio < 1.0e6) ) ;
% 
% % Set up the fitting problem y = ax + b
% y = log(aw(ii)) ; x = log(an(ii)) ;
% 
% N = length(x) ;
% a = (N * sum( x .* y) - sum(x)*sum(y)) / (N*sum(x.*x) - sum(x)*sum(x)) ;
% b = sum(y - a*x) / N;
% C =  exp( b / (-a) ) ;
% sc = ( abs(priorWeights) / C ).^a ;
% 
% % Scale the basis set
% tmp = diag(sc) * transpose(Phi64) ;
% Phi64 = transpose(tmp) ;
end

