--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/spsdSubtraction.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,504 @@
+% SPSDSUBTRACTION makes a sPSD-weighted least-square iterative fit
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: SPSDSUBTRACTION makes a sPSD-weighted least-square iterative fit
+%
+% CALL: [MPest, plOut, aoResiduum, aoP, aoPini] = spsdSubtraction(ao_Y, [ao_U1, ao_U2, ao_U3 ...]);
+%       [MPest, plOut, aoResiduum, aoP, aoPini] = spsdSubtraction(ao_Y, [ao_U1, ao_U2, ao_U3 ...], pl);
+%
+%  The function finds the optimal M that minimizes the sum of the weighted sPSD of
+%  (ao_Y - M * [ao_U1 ao_U2 ao_U3 ...] )
+%  if ao_Y is a vector of aos, the use the matrix/spsdSubtraction is
+%  advised
+%
+%  OUTPUTS: - MPest: output PEST object with parameter estimates
+%           - aoResiduum: residuum times series
+%           - plOut: plist containing data like the parameter estimates
+%           - aoP: last weight used in the optimization (fater last
+%             Maximization/Expectation step)
+%           - aoPini: initial weight used in the optimization (before first
+%             Maximization/Expectation step)
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'spsdSubtraction')">Parameters Description</a>
+%
+% VERSION : $Id: spsdSubtraction.m,v 1.6 2011/08/03 19:21:10 adrien Exp $
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = spsdSubtraction(varargin)
+  
+  % use the caller is method flag
+  callerIsMethod = utils.helper.callerIsMethod;
+  
+  % Check if this is a call for parameters
+  if utils.helper.isinfocall(varargin{:})
+    varargout{1} = getInfo(varargin{3});
+    return
+  end
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  if ~nargin>1
+    error('optSubtraction requires at least the two input aos as first and second entries')
+  end
+  
+  %% retrieving the two input aos
+  [aos_in, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+  
+  aosY = varargin{1};
+  aosU = varargin{2};
+  if (~isa(aosY, 'ao')) || (~isa(aosU, 'ao'))
+    error('first two inputs should be two ao-arrays involved in the subtraction')
+  end
+  
+  % Collect plist
+  pl = utils.helper.collect_objects(varargin(:), 'plist');
+  
+  % Get default parameters
+  pl = combine(pl, getDefaultPlist);
+  
+  %% checking data sizes
+  NY = numel(aosY);
+  if NY==0
+    error('Nothing to subtract to!')
+  end
+  NU = size(aosU,2);
+  if NU==0
+    error('Nothing to subtract!')
+  end
+  if ~(size(aosY,2)==1)
+    error('The input ao Y array should be a column vector')
+  end
+  if ~(size(aosU,1)==NY)
+    error('The fields ''subtracted'' should be an array of aos with the height of numel(initial)')
+  end
+  
+  %% collecting history
+  if callerIsMethod
+    % we don't need the history of the input
+  else
+    inhist  = [aosY(:).hist aosU(:).hist];
+  end
+  
+  %% retrieving general quantities
+  ndata = numel(aosY(1).y);
+  Ts = 1/aosY(1).fs;
+  nFreqs = floor(ndata/2);
+  freqs = 1/(2*Ts) * linspace(0,1,nFreqs);
+  
+  %% produce window
+  Win = find(pl, 'Win');
+  if isa(Win, 'plist')
+    Win = ao( combine(plist( 'length', ndata), Win) );
+    W = Win.y;
+  elseif isa(Win, 'ao')
+    if ~isa(Win.data, 'tsdata')
+      error('An ao window should be a time series')
+    end
+    W = Win.y;
+    if ~length(W)==ndata
+      error('signals and windows don''t have the same length')
+    end
+  else
+    error('input option Win is not acceptable (not a plist nor an ao)!')
+  end
+  
+  %% get initial M coefficient matrix
+  M = pl.find('coefs');
+  if isempty(M)
+    M = zeros(1,NU);
+  end
+  
+  %% get criterion thinness
+  linCoef = pl.find('lincoef');
+  logCoef = pl.find('logcoef');
+  
+  %% getting the input data Y and taking FFT
+  Y = zeros(NY, nFreqs);
+  YLocNorm = zeros(NY,1);
+  
+  for ii=1:NY
+    if isempty(aosY(ii).data)
+      error('One ao for Y is empty!')
+    end
+    if ~(length(aosY(ii).y)==ndata)
+      error('various Y vectors do not have the same length')
+    end
+    yLoc = fft(aosY(ii).y .* W, ndata);
+    YLocNorm(ii) = norm(aosY(ii).y .* W)/sqrt(ndata);
+    Y(ii,:) = yLoc(1:nFreqs)/YLocNorm(ii);
+  end
+  
+  %% getting the data U norm
+  ULocNorm = zeros(NY,NU);
+  for iU=1:NU
+    for iY=1:NY
+      if ~isempty(aosU(iY,iU).data)
+        ULocNorm(iY,iU,:) = norm(aosU(iY,iU).y  .* W)/sqrt(ndata);
+      end
+    end
+  end
+  ULocNorm = max(ULocNorm,[],1);
+  
+  %% getting the input data U and taking FFT
+  U = zeros(NY,NU, nFreqs);
+  for iY=1:NY
+    for iU=1:NU
+      if ~isempty(aosU(iY,iU).data)
+        if ~(length(aosU(iY,iU).y)==ndata)
+          error('various U vectors do not have the same length as Y')
+        end
+        uLoc = fft(aosU(iY,iU).y .* W, ndata);
+        U(iY,iU,:) = uLoc(1:nFreqs)/ (YLocNorm(iY) * ULocNorm(iU));
+      end
+    end
+  end
+  
+  %% getting the weight powAvgWeight
+  weightingMethod =pl.find('weightingMethod');
+  switch lower(weightingMethod)
+    case 'pzmodel'
+      weightModel =pl.find('pzmodelWeight');
+      if numel(weightModel)~=NY
+        error('there should be as many pzmodels as weighted entries')
+      end
+      for ii=1:NY
+        weight = weightModel(ii).resp(freqs);
+        weight = abs(weight).^2;
+        pow = [0 ; weight.y(2:nFreqs)];
+        [freqsAvg, powAvgs, nFreqsAvg, nDofs, binningMatrix] = ltpda_spsd(freqs, pow, linCoef, logCoef);
+        powAvgWeight(ii,:) = powAvgs; %#ok<AGROW>
+      end
+    case 'ao'
+      weight = pl.find('aoWeight');
+      if numel(weight)~=NY
+        error('there should be as many AOs as weighted entries')
+      end
+      for ii=1:NY
+        if ~isa(weight(ii).data, 'fsdata')
+          error('if weight is an ao, it should be a FSdata')
+        elseif length(weight(ii).y)~=nFreqs
+          error(['length of FS weight is not length of the FFT vector : ' num2str(length(weight(ii).y)) ' instead of ' num2str(nFreqs)])
+        else
+          pow = weight(ii).y;
+          [freqsAvg, powAvgs, nFreqsAvg, nDofs, binningMatrix] = ltpda_spsd(freqs, pow, linCoef, logCoef);
+          powAvgWeight(ii,:) = powAvgs; %#ok<AGROW>
+          %% add unit check here!!
+        end
+      end
+    case 'residual'
+      [freqsAvg, powAvgWeight, nFreqsAvg, nDofs, binningMatrix] = computeWeight(Y, M, U, freqs, linCoef, logCoef);
+    otherwise
+      error('weighting method requested does not exist!')
+  end
+  powAvgInv = (powAvgWeight.*(nFreqsAvg.')./(nDofs.')).^-1;
+  
+  %% get ME iterations termination conditions
+  iterMax = pl.find('iterMax');
+  normCoefs = pl.find('normCoefs');
+  normCriterion = pl.find('normCriterion');
+  
+  %% Maximization Expectation iterations loop
+  for i_iter = 1:iterMax
+    utils.helper.msg(utils.const.msg.PROC3, ['starting iteration ', num2str(i_iter)]);
+    
+    %% initializing history
+    if i_iter==1 % storing intial weight
+      Pini = powAvgWeight;
+      MHist(1,:) = reshape(M, [1, numel(M)] );
+    end
+    fValIni = optimalCriterion(Y, M, U, powAvgInv, linCoef, logCoef);
+    
+    %% solving LSQ problem
+    [M, hessian] = solveProblem(M, Y, U, powAvgInv, nFreqsAvg, binningMatrix);
+    fval = optimalCriterion(Y, M, U, powAvgInv, linCoef, logCoef);
+    
+    %% store history
+    fValHist(i_iter) = fval/fValIni; %#ok<AGROW>
+    MHist(i_iter+1,:) = reshape(M, [1, numel(M)] ); %#ok<AGROW>
+    
+    %% updating weight, recomputing residuum power
+    [freqsAvg, powAvgWeight, nFreqsAvg, nDofs] = computeWeight(Y, M, U, freqs, linCoef, logCoef);
+    powAvgInv = (powAvgWeight.*(nFreqsAvg.')./(nDofs.')).^-1;
+    
+    %% deciding whether to pursue or not ME iterations
+    if strcmpi( weightingMethod, 'pzmodel')
+      display('One iteration for Pzmodel weighting only')
+      break
+    elseif strcmpi( weightingMethod, 'ao')
+      display('One iteration for ao weighting only')
+      break
+    elseif norm(fValHist(i_iter)-1) < normCriterion && norm(MHist(i_iter+1,:)-MHist(i_iter,:))<normCoefs
+      display(['Iterations stopped at iteration ' num2str(i_iter) ' because not enough progress was made (see parameter "normCriterion" and "normCoefs")'])
+      break
+    elseif i_iter == iterMax
+      display(['Iterations stopped at maximum number of iterations ' num2str(i_iter) ' (see parameter "iterMax")'])
+      break
+    end
+  end
+  
+  %% creating output pest
+  MVals = M * diag( ULocNorm.^-1 );
+  MStd = diag(diag(ULocNorm) * hessian * diag(ULocNorm)).^-0.5;
+  MCov = diag(ULocNorm)^-1 * hessian^-1 * diag(ULocNorm)^-1;
+  
+  % prepare model, units, names
+  model = [];
+  for jj = 1:NU
+    names{jj} = ['U' num2str(jj)]; %#ok<AGROW>
+    units{jj} = aosY(1).yunits / aosU(1,jj).yunits; %#ok<AGROW>
+    xunits{jj} = aosU(1,jj).yunits; %#ok<AGROW>
+    MNames{jj} = ['M' num2str(jj)]; %#ok<AGROW>
+    if jj == 1
+      model = ['M' num2str(jj) '*U' num2str(jj)];
+    else
+      model = [model ' + M' num2str(jj) '*U' num2str(jj)]; %#ok<AGROW>
+    end
+  end
+  
+  model = smodel(plist('expression', model, ...
+    'params', MNames, ...
+    'values', MVals.', ...
+    'xvar', names, ...
+    'xunits', xunits, ...
+    'yunits', aosY(1).yunits ...
+    ));
+  
+  % collect inputs names
+  argsname = aosY(1).name;
+  for jj = 1:numel(NU)
+    argsname = [argsname ',' aosU(jj).name];
+  end
+  
+  % Build the output pest object
+  MPest = pest;
+  MPest.setY( MVals.' );
+  MPest.setDy(MStd);
+  MPest.setCov(MCov);
+  MPest.setChi2(0);
+  MPest.setNames(names{:});
+  MPest.setYunits(units{:});
+  MPest.setModels(model);
+  MPest.name = sprintf('optSubtraction(%s)', argsname);
+  
+  % Set procinfo object
+  MPest.procinfo = plist('MPsdE', 0);
+  % Propagate 'plotinfo'
+  plotinfo = [aosY(:).plotinfo aosU(:).plotinfo];
+  if ~isempty(plotinfo)
+    MPest.plotinfo = combine(plotinfo);
+  end
+  
+  %% creating output plist
+  plOut = plist;
+  
+  p = param({ 'criterion' , 'last value of the criterion in the last optimization'}, fval );
+  plOut.append(p);
+  p = param({ 'M' , 'Best fitting value'}, MVals );
+  plOut.append(p);
+  p = param({ 'Mhist' , 'History of the best fit, through iteration'}, MHist * diag( ULocNorm.^-1 ) );
+  plOut.append(p);
+  p = param({ 'fValHist' , 'History of the criterion value, through iteration'}, fValHist );
+  plOut.append(p);
+  p = param({ 'hessian' , 'fitting hessian'},  diag(ULocNorm) * hessian * diag(ULocNorm) );
+  plOut.append(p);
+  %add history and use Mdata/Pest instead
+  
+  %% creating aos for the weights used
+  if nargout>2
+    aoP = ao.initObjectWithSize(NY, 1);
+    aoPini = ao.initObjectWithSize(NY, 1);
+    for ii=1:NY
+      aoP(ii).setData(fsdata( freqsAvg, YLocNorm(ii)^2 * powAvgWeight(ii,:) ));
+      aoP(ii).setName('final weight');
+      aoP(ii).setXunits('Hz');
+      aoP(ii).setYunits(aosY(ii).yunits^2 * unit('Hz^-1'));
+      aoP(ii).setDescription(['final weight in the channel "' aosY(ii).name '"']);
+      aoP(ii).setT0(aosY(ii).t0);
+      aoPini(ii).setData(fsdata( freqsAvg, YLocNorm(ii)^2 * Pini(ii,:) ));
+      aoPini(ii).setName('initial weight');
+      aoPini(ii).setXunits('Hz');
+      aoPini(ii).setYunits(aosY(ii).yunits^2 * unit('Hz^-1'));
+      aoPini(ii).setDescription(['initial weight in the channel "' aosY(ii).name '"']);
+      aoPini(ii).setT0(aosY(ii).t0);
+    end
+  end
+  
+  %% creating residuum time-series
+  if nargout>2
+    aoResiduum = ao.initObjectWithSize(NY,1);
+    for ii=1:NY
+      aoResiduumValue = aosY(ii).y;
+      for jj = 1:NU
+        aoResiduumValue = aoResiduumValue - MVals(jj)*aosU(ii,jj).y;
+      end
+      aoResiduum(ii).setData(tsdata( aoResiduumValue, aosY(ii).fs ));
+      aoResiduum(ii).setName(['residual in the channel "' aosY(ii).name '"' ]);
+      aoResiduum(ii).setXunits('s');
+      aoResiduum(ii).setYunits(aosY(ii).yunits);
+      aoResiduum(ii).setDescription(['residual corresponding to "' aosY(ii).description '"' ]);
+      aoResiduum(ii).setT0(aosY(ii).t0);
+    end
+  end
+  
+  %% adding history
+  if callerIsMethod
+    % we don't need to set the history
+  else
+    MPest.addHistory(getInfo('None'), pl, ao_invars, inhist);
+    if nargout>2
+      for ii=1:NY
+        aoP(ii).addHistory(getInfo('None'), pl, ao_invars, inhist);
+        aoPini(ii).addHistory(getInfo('None'), pl, ao_invars, inhist);
+        aoResiduum(ii).addHistory(getInfo('None'), pl, ao_invars, inhist);
+      end
+    end
+  end
+  
+  %% return coefficients and hessian and Jfinal and powAvgWeight
+  if nargout>2
+    varargout = {MPest, plOut, aoResiduum, aoP, aoPini};
+  else
+    varargout = {MPest, plOut};
+  end
+end
+
+%% weight for optimal criterion
+function [freqsAvg, powAvgWeight, nFreqsAvg, nDofs, binningMatrix] = computeWeight(Y, M, U, freqs, linCoef, logCoef)
+  errDft =  subtraction( Y, M, U);
+  errPow = real(errDft).^2 + imag(errDft).^2;
+  for ii=1:size(errDft,1)
+    [freqsAvg, powAvgs, nFreqsAvg, nDofs, binningMatrix] = ltpda_spsd(freqs, errPow, linCoef, logCoef);
+    powAvgWeight(ii,:) = powAvgs;  %#ok<AGROW>
+  end
+end
+
+%% optimal criterion
+function j = optimalCriterion(Y, M, U, powAvgInv, linCoef, logCoef)
+  errDft = subtraction(Y, M, U);
+  errPow = real(errDft).^2 + imag(errDft).^2;
+  j = 0;
+  for ii=1:size(errDft,1)
+    [freqsAvg, powAvgs, nFreqsAvg, nDofs] = ltpda_spsd([], errPow, linCoef, logCoef); %#ok<ASGLU>
+    powSum = powAvgs .* nDofs; % binning frequencies as in sPSD
+    j  = j + sum( powSum .* powAvgInv(:,ii) );
+    %     alpha = 4;
+    %     logProbaDensityFactor =  - nFreqsAvg * log(2) - gammaln(nFreqsAvg);
+    %     normlzChi2Sum = ((alpha*2)*powSum) .* powAvgInv(:,ii); % divide the sum by the expected average of each terms, so the chi2 is normalized
+    %     logProbaDensities = logProbaDensityFactor + (nFreqsAvg-1).*log(normlzChi2Sum) - normlzChi2Sum/2 ; % here computing log of probability
+    %     j = j - sum(logProbaDensities); % better than taking product of probabilities
+  end
+end
+
+%% time-series subtraction function
+function Y = subtraction( Y, M, U)
+  ndata = size(Y,2);
+  for ii=1:size(Y,1)
+    for j=1:numel(M)
+      Y(ii,:) = Y(ii,:) - reshape( M(j)*U(ii,j,:) , [1,ndata] );
+    end
+  end
+end
+
+%% Direct solver
+function [M, hessian] = solveProblem(M, Y, U, powAvgInv, nFreqsAvg, binningMatrix)
+  errDft = subtraction(Y, M, U);
+  NU = size(U,2);
+  NFreqs = size(binningMatrix,2);
+  ATB = zeros(NU,1);
+  ATA = zeros(NU,NU);
+  % matrix for frequency binning & Weighting & Summing :
+  matBSW = powAvgInv * binningMatrix;
+  for iiParam = 1:NU
+    Uii = reshape(U(1,iiParam,:), [1 NFreqs]);
+    ATB(iiParam) = 2 * ( matBSW * real( Uii .* conj(errDft) ).' );
+    for jjParam = 1:NU
+      Ujj = reshape(U(1,jjParam,:), [1 NFreqs]);
+      ATA(iiParam,jjParam) = 2 * ( matBSW * real( Uii .* conj(Ujj) ).' );
+    end
+  end
+  try
+    MUpdate = ATA^-1 * ATB;
+    M = M + MUpdate.';
+  catch 
+    warning('Numerical accuracy limited the number of iterations')
+  end
+  hessian = ATA;
+end
+
+
+%--------------------------------------------------------------------------
+% Get Info Object
+%--------------------------------------------------------------------------
+function ii = getInfo(varargin)
+  if nargin == 1 && strcmpi(varargin{1}, 'None')
+    sets = {};
+    pl   = [];
+  else
+    sets = {'Default'};
+    pl   = getDefaultPlist;
+  end
+  % Build info object
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: spsdSubtraction.m,v 1.6 2011/08/03 19:21:10 adrien Exp $', sets, pl);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function plout = getDefaultPlist()
+  persistent pl;
+  if exist('pl', 'var')==0 || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;
+end
+
+function pl = buildplist()
+  pl = plist;
+  
+  % initial coefficients for subtraction initialization
+  p = param({ 'coefs' , 'initial subtracted coefficients, must be a nY*nU double array. If not provided zeros are assumed'},  [] );
+  pl.append(p);
+  
+  % weighting scheme
+  p = param({ 'weightingMethod' , 'choose to define a frequency weighting scheme'},  {1, {'residual', 'ao', 'pzmodel'}, paramValue.SINGLE} );
+  pl.append(p);
+  
+  p = param({ 'aoWeight' , 'ao to define a frequency weighting scheme (if chosen in ''weightingMethod'')'},  ao.initObjectWithSize(0,0) );
+  pl.append(p);
+  
+  p = param({ 'pzmodelWeight' , 'pzmodel to define a frequency weighting scheme (if chosen in ''weightingMethod'')'},  pzmodel.initObjectWithSize(0,0) );
+  pl.append(p);
+  
+  p = param({ 'lincoef' , 'linear coefficient for scaling frequencies in chi2'},  5 );
+  pl.append(p);
+  
+  p = param({ 'logcoef' , 'logarithmic coefficient for scaling frequencies in chi2'},  0.3 );
+  pl.append(p);
+  
+  % iterations convergence stop criterion
+  p = param({ 'iterMax' , 'max number of Mex/Exp iterations'},  20 );
+  pl.append(p);
+  
+  p = param({ 'normCoefs' , 'tolerance on inf norm of coefficient update (used depending on ''CVCriterion'')'},  1e-15 );
+  pl.append(p);
+  
+  p = param({ 'normCriterion' , 'tolerance on norm of criterion variation (used depending on ''CVCriterion'')'},  1e-15 );
+  pl.append(p);
+  
+  % windowing options
+  p = param({ 'win' , 'window to operate FFT, may be a plist/ao'},  plist('win', 'levelledHanning', 'PSLL', 200, 'levelOrder', 4 ) );
+  pl.append(p);
+  
+  % display
+  p = param({ 'display' , 'choose how much to display of the optimizer output'},  {1, {'off', 'iter', 'final'}, paramValue.SINGLE} );
+  pl.append(p);
+  
+  % optimizer options
+  p = param({ 'maxcall' , 'maximum number of calls to the criterion function'},  5000 );
+  pl.append(p);
+  
+end
+
+