--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/linfit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,294 @@
+% 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
+