% LINFIT is a linear fitting tool
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
% DESCRIPTION: LINFIT is a linear fitting tool based on MATLAB's
% lscov function. It solves an equation in the form
%
%     Y = P(1) + X * P(2)
%
% for the fit parameters P. 
% The output is a pest object where the fields are containing: 
% Quantity                              % Field
% Fit parameters                            y
% Uncertainties on the fit parameters 
% (given as standard deviations)            dy
% The reduced CHI2 of the fit              chi2
% The covariance matrix                    cov
% The degrees of freedom of the fit        dof
%  
% CALL:       P = linfit(X, Y, PL)
%             P = linfit(A, PL)
%  
% INPUTS:     Y   - dependent variable
%             X   - input variable
%             A   - data ao whose x and y fields are used in the fit
%             PL  - parameter list
% 
% OUTPUT:     P   - a pest object with the fitting coefficients
%
%
% PARAMETERS:
%    'dy' - uncertainty on the dependent variable
%    'dx' - uncertainties on the input variable
%    'p0' - initial guess on the fit parameters used ONLY to propagate
%           uncertainities in the input variable X to the dependent variable Y
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'linfit')">Parameters Description</a>
%
% VERSION:     $Id: linfit.m,v 1.23 2011/05/15 22:52:57 mauro Exp $
%
% EXAMPLES:
%
%  %% Make fake AO from polyval
%   nsecs = 100;
%   fs    = 10;
%   
%   u1 = unit('fm s^-2');
%   u2 = unit('nT');
%   
%   pl1 = plist('nsecs', nsecs, 'fs', fs, ...
%     'tsfcn', 'polyval([10 1], t) + randn(size(t))', ...
%     'xunits', 's', 'yunits', u1);
% 
%   pl2 = plist('nsecs', nsecs, 'fs', fs, ...
%     'tsfcn', 'polyval([-5 0.2], t) + randn(size(t))', ...
%     'xunits', 's', 'yunits', u2);
%
%   a1 = ao(pl1);
%   a2 = ao(pl2); 
% 
%  %% 1) Determine dependance from time of a time-series
%  %% Fit a stright line the a1 dependance from time
%   p1 = linfit(a1, plist());
%   p2 = linfit(a1, plist('dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.y)), 'P0', ao([0 0])));
%   p3 = linfit(a1, plist('dx', ao(0.1, plist('yunits', a1.xunits)), 'dy', ao(0.1, plist('yunits', a1.yunits)), 'P0', p1));
%
%  %% Compute fit: evaluating pest
%   
%   b1 = p1.eval(plist('type', 'tsdata', 'XData', a1, 'xfield', 'x'));
%   b2 = p2.eval(plist('type', 'tsdata', 'XData', a1.x));
%   b3 = p3.eval(plist('type', 'tsdata', 'XData', a1.x));
%   
%  %% Plot fit
%   iplot(a1, b1, b2, b3, plist('LineStyles', {'', '--', ':', '-.'}));
%
%  %% Remove linear trend
%   c = a1 - b1;
%   iplot(c)
%
%  %% 2) Determine dependance of a time-series from another time-series
%  %% Fit with a straight line the a1 dependance from a2
%   
%   p1 = linfit(a1, a2, plist());
%   p2 = linfit(a1, a2, plist('dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.x)), 'P0', ao([0 0])));
%   p3 = linfit(a1, a2, plist('dx', ao(0.1, plist('yunits', a1.yunits)), 'dy', ao(0.1, plist('yunits', a2.yunits)), 'P0', p1));
%   
%  %% Compute fit: evaluating pest
%   
%   b1 = p1.eval(plist('type', 'xydata', 'XData', a1.y, 'xunits', a1.yunits));
%   b2 = p2.eval(plist('type', 'xydata', 'XData', a1));
%   b3 = p3.eval(plist('type', 'xydata', 'XData', a1.y, 'xunits', a1.yunits));
%
%  %% Build reference object
%   a12 = ao(plist('xvals', a1.y, 'yvals', a2.y, ...
%     'xunits', a1.yunits, 'yunits', a2.yunits));
%   
%  %% Plot fit
%   iplot(a12, b1, b2, b3, plist('LineStyles', {'', '--', ':', '-.'}));
%   
%   %% Remove linear trend
%   c = a12 - b3;
%   iplot(c)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = linfit(varargin)
  
  % check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
 
  % tell the system we are runing
  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
  [aos, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pli              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  if nargout == 0
    error('### linfit can not be used as a modifier method. Please give at least one output');
  end
  
  % combine plists, making sure the user input is not empty
  pli = combine(pli, plist());
  pl = parse(pli, getDefaultPlist());
  
  % extract arguments
  if (length(aos) == 1)
    % we are using x and y fields of the single ao we have
    x = aos(1).x;
    dx = aos(1).dx;
    y = aos(1).y;
    dy = aos(1).dy;
    xunits = aos(1).xunits;
    yunits = aos(1).yunits;
    argsname = aos(1).name;
  elseif (length(aos) == 2)
    % we are using y fields of the two aos we have
    x = aos(1).y;
    dx = aos(1).dy;
    y = aos(2).y;
    dy = aos(2).dy;
    xunits = aos(1).yunits;
    yunits = aos(2).yunits;
    argsname = [aos(1).name ',' aos(2).name];
  else
    error('### linfit needs one or two input AOs');
  end

  % extract plist parameters. For dx and dy we check the user input plist before
  dx = find(pli, 'dx', dx);
  dy = find(pli, 'dy', dy);
  p0 = find(pl, 'p0');
  
  % vectors length
  len = length(y);
  
  % uncertainty on Y  
  if isempty(dy)
    dy = 1;    
  end
  if isa(dy, 'ao')
    % check units
    if yunits ~= dy.yunits
      error('### Y and DY units are not compatible - %s %s', char(yunits), char(dy.yunits));
    end
    % extract values from AO
    dy = dy.y;
  end
  if isscalar(dy)
    % given a single value construct a vector
    dy = ones(len, 1) * dy;
  end
  
  % weights
  sigma2 = dy.^2;
  
  % extract values for initial guess
  if (isa(p0, 'ao') || isa(p0, 'pest'))
    p0 = p0.y;
  end

  % uncertainty on X
  if ~isempty(dx)       
    if length(p0) ~= 2
      error('### initial parameters guess p0 is mandatory for proper handling of X uncertainties');
    end
    
    if isa(dx, 'ao')
      % check units
      if xunits ~= dx.yunits
        error('### X and DX units are not compatible - %s %s', char(xunits), char(dx.yunits));
      end
      % extract values from AO
      dx = dx.y;
    end
    if isscalar(dx)
      % given a single value construct a vector
      dx = ones(len, 1) * dx;
    end
    
    % add contribution to weights
    sigma2 = sigma2 + p0(2)^2 .* dx.^2;    
  end
      
  % construct matrix
  m = [ ones(len, 1) x ];

  % solve
  [p, stdp, mse, s] = lscov(m, y, 1./sigma2);

  % scale errors and covariance matrix  
  stdp = stdp ./ sqrt(mse);
  s = s ./ mse;
  
  % compute chi2
  dof = len - 2;
  chi2 = sum((y - p(1)-p(2)*x).^2 ./ sigma2) / dof;
  
  % prepare model, units, names
  names = {'P1', 'P2'};  
  model = 'P1 + P2*X';
  model = smodel(plist('expression', model, ...
    'params', names, ...
    'values', p, ...
    'xvar', 'X', ...
    'xunits', xunits, ...
    'yunits', yunits ...
    ));
  units = [yunits simplify(yunits/xunits)];  

  % Build the output pest object
  out = pest;
  out.setY(p);
  out.setDy(stdp);
  out.setCov(s); 
  out.setChi2(chi2);
  out.setDof(dof);
  out.setNames(names{:}); 
  out.setYunits(units);
  out.setModels(model); 
  out.name = sprintf('linfit(%s)', argsname);
  out.addHistory(getInfo('None'), pl,  ao_invars, [aos(:).hist]);
  % Set procinfo object
  out.procinfo = plist('MSE', mse);
  
  % set outputs
  varargout{1} = out;
  
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: linfit.m,v 1.23 2011/05/15 22:52:57 mauro Exp $', sets, pl);
  ii.setModifier(false);
  ii.setArgsmin(1);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------

function plout = getDefaultPlist()
  persistent pl;
  if ~exist('pl', 'var') || isempty(pl)
    pl = buildplist();
  end
  plout = pl;
end

function pl = buildplist()
  
  % default plist for linear fitting
  pl = plist.LINEAR_FIT_PLIST;
  
end