--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/bilinfit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,286 @@
+% BILINFIT is a linear fitting tool
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: BILINFIT linear fitting tool based on MATLAB's lscov
+% function. It solves an equation in the form
+%
+%     Y = X(1) * P(1) + X(2) * P(2) + ... + P(N+1)
+%
+% for the fit parameters P. It handles an arbitrary number of input vectors
+% and uncertainties on the dependent vector Y and input vectors X(1..N).
+% The output is a pest object where the fields are containing:
+% Quantity                              % Field
+% Fit coefficients                          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 = bilinfit(X1, X2, .., XN, Y, PL)
+%
+% INPUTS:     Y       - dependent variable
+%             X(1..N) - input variables
+%             PL      - parameter list
+%
+% OUTPUT:     P   - a pest object with the N+1 elements
+%
+%
+% PARAMETERS:
+%    'dy' - uncertainty on the dependent variable
+%    'dx' - uncertainties on the input variables
+%    'p0' - initial guess on the fit parameters to propagate uncertainities
+%           in the input variables X(1..N) to the dependent variable Y
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'bilinfit')">Parameters Description</a>
+%
+% VERSION:    $Id: bilinfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $
+%
+% EXAMPLES:
+%
+% % 1) 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', 'm'));
+%  x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
+%  n  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
+%  c = [ao(1,plist('yunits','m/m')) ao(2,plist('yunits','m/m'))];
+%  y = c(1)*x1 + c(2)*x2 + n;
+%  y.simplifyYunits;
+%
+% % Get a fit for the c coefficients and a constant term
+%  p = bilinfit(x1, x2, y)
+%
+% % Do linear combination: using eval
+%  pl_split = plist('times', [1 5]);
+%  yfit = eval(p, split(x1, pl_split), split(x2, pl_split));
+% 
+% % Plot (compare data with fit)
+% iplot(y, yfit, plist('Linestyles', {'-','--'}))
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = bilinfit(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);
+  pl               = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+  
+  if nargout == 0
+    error('### bilinfit can not be used as a modifier method. Please give at least one output');
+  end
+  
+  if numel(aos) < 2
+    error('### bilinfit needs at least 2 inputs AOs');
+  end
+  
+  % combine plists
+  pl = parse(pl, getDefaultPlist());
+  
+  % extract parameters
+  dy = find(pl, 'dy');
+  dx = find(pl, 'dx');
+  p0 = find(pl, 'p0');
+  
+  % collect inputs
+  Y = aos(end);
+  X = aos(1:end-1);
+  
+  % collect inputs names
+  argsname = aos(1).name;
+  for jj = 2:numel(aos)
+    argsname = [argsname ',' aos(jj).name];
+  end
+  
+  % get data from AOs
+  x = X(:).y;
+  y = Y.y;
+  
+  % vectors length
+  N = length(y);
+  
+  % uncertainty on Y
+  if isempty(dy)
+    dy = 1;
+  end
+  if isa(dy, 'ao')
+    % check units
+    if Y.yunits ~= dy.yunits
+      error('### Y and DY units are not compatible - %s %s', char(Y.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(N, 1) * dy;
+  end
+  
+  % squares
+  sigma2 = dy.^2;
+  sigma2y_rms = sqrt(sum(sigma2)/N);
+  
+  % extract values for initial guess
+  if (isa(p0, 'ao') || isa(p0, 'pest'))
+    p0 = p0.y;
+  end
+  
+  % uncertainty on X
+  if ~isempty(dx)
+    
+    for k = 1:length(dx)
+      dxi = dx(k);
+      
+      if ~isempty(dxi)
+        if isa(dxi, 'ao')
+          % check units
+          if X(k).yunits ~= dxi.yunits
+            error('### X and DX units are not compatible - %s %s', char(X.yunits), char(dxi.yunits));
+          end
+          % extract values from AO
+          dxi = dxi.y;
+        end
+        if isscalar(dxi)
+          % given a single value construct a vector
+          dxi = ones(N, 1) * dxi;
+        end
+        
+        % squares
+        sigma2xi = dxi.^2;
+        
+        % if A0 guess are not given
+        if isempty(p0(k))
+          % set it to obtain equal error contribution to the Y error
+          sigma2xi_rms = sqrt(sum(sigma2xi)/N);
+          p0(k) = sigma2y_rms/sigma2xi_rms;
+        end
+        
+        % add contribution to weights
+        sigma2 = sigma2 + sigma2xi * p0(k)^2;
+      end
+      
+    end
+  end
+  
+  % constant term
+  c = ones(N, 1);
+  
+  % build matrix
+  m = [x c];
+  
+  % solve
+  [p, stdx, mse, s] = lscov(m, y, 1./sigma2);
+  
+  % scale errors and covariance matrix
+  stdp = stdx ./ sqrt(mse);
+  s    = s ./ mse;
+  
+  % compute chi2
+  dof = N - length(p);
+  chi2 = sum((y - lincom(m, p)).^2 ./ sigma2) / dof;
+  
+  % prepare model, units, names
+  model = [];
+  for kk = 1:length(p)
+    switch kk
+      case 1
+        units(kk) = simplify(Y.yunits/X(kk).yunits);
+        model = ['P' num2str(kk) '*X' num2str(kk)];
+        xvar{kk} = ['X' num2str(kk)];
+        xunits{kk} = X(kk).yunits;
+      case length(p)
+        units(kk) = Y.yunits;
+        model = [model ' + P' num2str(kk)];
+      otherwise
+        units(kk) = simplify(Y.yunits/X(kk).yunits);
+        model = [model ' + P' num2str(kk) '*X' num2str(kk)];
+        xvar{kk} = ['X' num2str(kk)];
+        xunits{kk} = X(kk).yunits;
+    end
+    names{kk} = ['P' num2str(kk)];
+    
+  end
+  model = smodel(plist('expression', model, ...
+                       'params', names, ...
+                       'values', p, ...
+                       'xvar', xvar, ...
+                       'xunits', xunits, ...
+                       'yunits', Y.yunits));
+  
+  
+  % 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('bilinfit(%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
+
+% computes linear combination
+function out = lincom(x, p)
+  assert(size(x, 2) == length(p));
+  out = zeros(size(x, 1), 1);
+  for k = 1:length(p)
+    out = out + x(:,k) * p(k);
+  end
+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: bilinfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $', sets, pl);
+  ii.setModifier(false);
+  ii.setArgsmin(2);
+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.MULTILINEAR_FIT_PLIST;
+  
+end