% LSCOV is a wrapper for MATLAB's lscov function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LSCOV is a wrapper for MATLAB's lscov function. It solves a
% set of linear equations by performing a linear least-squares fit. It
% solves the problem
%
%        Y = HX
%
%   where X are the parameters, Y the measurements, and H the linear
%   equations relating the two.
%
% CALL:        X = lscov([C1 C2 ... CN], Y, pl)
%              X = lscov(C1,C2,C3,...,CN, Y, pl)
%
% INPUTS:      C1...CN - AOs which represent the columns of H.
%              Y       - AO which represents the measurement set
%
% Note: the length of the vectors in Ci and Y must be the same.
% Note: the last input AO is taken as Y.
%
%              pl - parameter list (see below)
%
% OUTPUTs:     X  - A pest object with fields: 
%                   y   - the N fitting coefficients to y_i 
%                   dy  - the parameters' standard deviations (lscov 'STDX' vector)
%                   cov - the parameters' covariance matrix (lscov 'COV' vector)
%                   
% The procinfo field of the output PEST object is filled with the following key/value
% pairs:
%     'MSE' - the mean-squared errors
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'lscov')">Parameters Description</a>
%
% VERSION:     $Id: lscov.m,v 1.42 2011/04/08 08:56:17 hewitson Exp $
%
% EXAMPLES:
%
% % 1) Determine the coefficients of a linear combination of noises:
%
% % Make some data
%  fs    = 10;
%  nsecs = 10;
%  B1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
%  B2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
%  B3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
%  n  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
%  c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T')) ao(3,plist('yunits','m T^-1'))];
%  y = c(1)*B1 + c(2)*B2 + c(3)*B3 + n;
%  y.simplifyYunits;
%  % Get a fit for c
%  p_s = lscov(B1, B2, B3, y);
%  % do linear combination: using lincom
%  yfit1 = lincom(B1, B2, B3, p_s);
%  yfit1.simplifyYunits;
%  % do linear combination: using eval
%  yfit2 = p_s.eval(B1, B2, B3);
% 
%  % Plot (compare data with fit)
%  iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
%
% % 2) Determine the coefficients of a linear combination of noises:
%
% % Make some data
%  fs    = 10;
%  nsecs = 10;
%  x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
%  x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
%  x3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'C'));
%  n  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
%  c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/m')) ao(3,plist('yunits','m C^-1'))];
%  y = c(1)*x1 + c(2)*x2 + c(3)*x3 + n;
%  y.simplifyYunits;
%  % Get a fit for c
%  p_m = lscov(x1, x2, x3, y);
%  % do linear combination: using lincom
%  yfit1 = lincom(x1, x2, x3, p_m);
%  % do linear combination: using eval
%  pl_split = plist('times', [1 5]);
%  yfit2 = p_m.eval(plist('Xdata', {split(x1, pl_split), split(x2, pl_split), split(x3, pl_split)}));
%  % Plot (compare data with fit)
%  iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = lscov(varargin)
  
  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
  
  import utils.const.*
  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs and plists
  [A, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl             = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  if numel(A) < 2
    error('### lscov needs at least 2 inputs AOs');
  end
  
  if nargout == 0
    error('### lscov can not be used as a modifier method. Please give at least one output');
  end
  
  % combine plists
  pl = parse(pl, getDefaultPlist());
  
  % collect inputs names
  argsname = A(1).name;
  for jj = 2:numel(A)
    argsname = [argsname ',' A(jj).name];
  end
  
  % Extract parameters
  W = find(pl, 'weights');
  V = find(pl, 'cov');
  
  % Build matrices for lscov
  C = A(1:end-1);
  Y = A(end);
  
  H = C(:).y;
  y = Y.y;
  if isa(W,'ao'), W = W.y; end;
  if isa(V,'ao'), V = V.y; end;
    
  if isempty(V)
    [p,stdx,mse,s] = lscov(H, y, W);
  else
    [p,stdx,mse,s] = lscov(H, y, W, V);
  end
      
  % prepare model, units, names
  model = [];
  for jj = 1:numel(p)
    names{jj} = ['C' num2str(jj)];
    units{jj} = Y.yunits / C(jj).yunits;
    xunits{jj} = C(jj).yunits;
    xvar{jj} = ['X' num2str(jj)];
    if jj == 1
      model = ['C' num2str(jj) '*X' num2str(jj) ' '];
    else
      model = [model ' + C' num2str(jj) '*X' num2str(jj)];
    end
  end

  model = smodel(plist('expression', model, ...
    'params', names, ...
    'values', p, ...
    'xvar', xvar, ...
    'xunits', xunits, ...
    'yunits', Y.yunits ...
    ));  
  
  % Build the output pest object
  X = pest;
  X.setY(p);
  X.setDy(stdx);
  X.setCov(s);
  X.setChi2(mse);
  X.setNames(names{:});
  X.setYunits(units{:});
  X.setModels(model);
  X.name = sprintf('lscov(%s)', argsname);
  X.addHistory(getInfo('None'), pl, ao_invars, [A(:).hist]);
  % Set procinfo object
  X.procinfo = plist('MSE', mse);
  % Propagate 'plotinfo'
  plotinfo = [A(:).plotinfo];
  if ~isempty(plotinfo)
    X.plotinfo = combine(plotinfo);
  end
  
  % Set outputs
  varargout{1} = X;
    
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: lscov.m,v 1.42 2011/04/08 08:56:17 hewitson Exp $', sets, pl);
  ii.setModifier(false);
  ii.setArgsmin(2);
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();
  
  % Weights
  p = param({'weights','An ao containing weights for the fit.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % Cov
  p = param({'cov','An ao containing a covariance matrix for the fit.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
end
% END