--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@pest/eval.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,391 @@
+% EVAL evaluate a pest object
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: EVAL evaluate a pest model.
+%
+% CALL:        b = eval(p, pl)
+%              b = eval(p, x, pl)
+%              b = eval(p, x1, ... , xN, pl)
+%              b = eval(p, [x1 ... xN], pl)
+%
+% INPUTS:      p  - input pest(s) containing parameter values.
+%              xi - input ao(s) containing x values (as x or y fields, depending on the 'xfield' parameter)
+%              pl - parameter list (see below)
+%
+% OUTPUTs:     b  - an AO containing the model evaluated at the given X
+%                   values, with the given parameter values.
+%
+% <a href="matlab:utils.helper.displayMethodInfo('pest', 'eval')">Parameters Description</a>
+%
+% VERSION:     $Id: eval.m,v 1.28 2011/05/15 22:47:14 mauro Exp $
+%
+% EXAMPLES:
+%
+% % 1)
+% % Prepare the symbolic model
+% mdl = smodel(plist('expression', 'a1.*x + a2.*x.^2 + a0', 'xvar', 'x', 'yunits', 'V'));
+%
+% % Prepare the pest object
+%
+% p = pest(plist('paramnames', {'a0','a1','a2'}, 'y', [1 2 3], 'models', mdl));
+%
+% % Evaluate the object
+% a1 = eval(p, plist('xdata', ao([1:10])))
+% a2 = eval(p, ao([1:10]))
+%
+% % 2)
+% % Prepare the symbolic model
+% mdl = smodel(plist('expression', 'a1.*x1 + a2.*x2 + a0', 'xvar', {'x1', 'x2'}, 'yunits', 'm', 'xunits', {'T', 'K'}));
+%
+% % Prepare the pest object
+%
+% p = pest(plist('paramnames', {'a0','a1','a2'}, 'y', [1 2 3], 'yunits', {'m', 'T/m', 'K/m'}, 'models', mdl));
+%
+% % Evaluate the object
+% x1 = ao(plist('yvals', [1:10], 'fs', 1, 'yunits', 'T'));
+% x2 = ao(plist('yvals', [1:10], 'fs', 1, 'yunits', 'K'));
+% a1 = eval(p, plist('xdata', [x1 x2]))
+% a2 = eval(p, [x1 x2])
+% a3 = eval(p, x1, x2)
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = eval(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
+  [psts, pst_invars, rest] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
+  [pl, pl_invars, rest]   = utils.helper.collect_objects(rest, 'plist', in_names);
+  [as, as_invars, rest]   = utils.helper.collect_objects(rest, 'ao', in_names);
+  [c, c_invars]           = utils.helper.collect_objects(rest, 'cell', in_names);
+  
+  if nargout == 0
+    error('### eval can not be used as a modifier method. Please give at least one output');
+  end
+  
+  % combine plists
+  pl = parse(pl, getDefaultPlist());
+  
+  % Extract necessary parameters
+  index     = find(pl, 'index');
+  x = find(pl, 'Xdata');
+  if ~isempty(as)
+    x = as;
+  end
+  % I don't know how to deal with the history of a cell array of aos
+  if ~isempty(c)
+    error('Please pass the arguments in a vector or a list')
+  end
+  
+  % Extract the information about the x field, if necessary
+  xfield = pl.find('xfield');
+  switch xfield
+    case 'x'
+      op = str2func('x');
+    case 'y'
+      op = str2func('y');
+    otherwise
+      error('oops')
+  end
+  
+  % Extract the xvals for the smodel, and the output x for the ao
+  switch class(x)
+    case 'cell'
+      switch class(x{1})
+        case 'ao'
+          data_type = find(pl, 'type');
+          if isempty(data_type)
+            % Nothing to do, output data type will be inherited from the first AO
+            if ~isa(x{1}.data, 'cdata')
+              out_x      = x{1};
+            else
+              out_x      = [];
+            end
+            out_xunits = '';
+          else
+            % In this case I have to extract the vector with the output x and xunits
+            out_x = feval(op, x{1});
+            out_xunits = feval(str2func([xfield 'units']), x{1});
+          end
+        case 'double'
+          data_type = find(pl, 'type');
+          xvals = x;
+          switch data_type
+            case {'tsdata', 'fsdata', 'xydata'}
+              out_x      = x;
+              out_xunits = find(pl, 'xunits');
+            case {'cdata', ''}
+              out_x      = [];
+              out_xunits = '';
+            otherwise
+              error('LTPDA error')
+          end
+        otherwise
+      end
+    case 'ao'
+      data_type = find(pl, 'type');
+      if isempty(data_type)
+        % Nothing to do, output data type will be inherited from the first AO
+        if ~isa(x(1).data, 'cdata')
+          out_x      = x(1);
+        else
+          out_x      = [];
+        end
+        out_xunits = '';
+      else
+        % In this case I have to extract the vector with the output x and xunits
+        out_x = feval(op, x(1));
+        out_xunits = feval(str2func([xfield 'units']), x(1));
+      end
+    case 'double'
+      data_type = find(pl, 'type');
+      xvals = x;
+      switch data_type
+        case {'tsdata', 'fsdata', 'xydata'}
+          out_x      = x;
+        case {'cdata', ''}
+          out_x      = [];
+        otherwise
+          error('LTPDA error')
+      end
+      out_xunits = find(pl, 'xunits');
+    otherwise
+  end
+  
+  % If the user wants to override the pest/smodel yunits, let's get them
+  % This works only if the user sets them to something not empty
+  yunits = find(pl, 'yunits');
+  
+  % If we have AOs in a cell....
+  if iscell(x) && all(cellfun(@(x)isa(x, 'ao'), x))
+    % if we have multiple x, we need to convert the y values into a cell array
+    xvals = cellfun(op, x, 'UniformOutput', false);
+  elseif isa(x, 'ao')
+    % we put the y values in to a cell array
+    if numel(x) > 1
+      xvals = cellfun(op, num2cell(x), 'UniformOutput', false);
+    else
+      % ... we take the x values, as per the help
+      xvals = feval(op, x);
+    end
+  end
+  
+  % Loop over input objects
+  for jj = 1:numel(psts)
+    
+    pst = psts(jj);
+    % evaluate models
+    m = copy(pst.models, true);
+    switch class(m)
+      case 'smodel'
+        % Make sure the smodel parameters are named the same as the pest
+        m(index).setParams(pst.names, pst.y);
+        % If the user provided the x vector(s), override the smodel x with these
+        if ~isempty(xvals)
+          m(index).setXvals(xvals);
+        end
+        % Go for the model evaluation
+        out(jj) = eval(m(index), plist(...
+          'output type', data_type, 'output x', out_x, 'output xunits', out_xunits));
+
+        % Setting the units of the evaluated model
+        if ~isempty(yunits)
+          out(jj).setYunits(yunits);
+        end
+        
+      case 'ao'
+        % do linear combination: using lincom
+        out(jj) = lincom(m, pst);
+        out(jj).simplifyYunits;
+        
+  
+      case 'matrix'
+        % check objects of the matrix and switch
+        switch class(m(1).objs)
+          case 'smodel'
+            % Make sure the smodel parameters are named the same as the pest
+            for ii = 1:numel(m.objs)
+              m.objs(ii).setParams(pst.names, pst.y);
+            end
+            % If the user provided the x vector(s), override the smodel x with these
+            if ~isempty(x)
+              for ii = 1:numel(m.objs)
+                m.objs(ii).setXvals(x);
+              end
+            end
+            % Go for the model evaluation
+            tout = ao.initObjectWithSize(size(m.objs,1),size(m.objs,2));
+            for ii=1:size(m.objs,1)
+              for kk=1:size(m.objs,2)
+                tout(ii,kk) = eval(m.objs(ii,kk), plist(...
+                  'output type', data_type, 'output x', out_x));
+                % Setting the units of the evaluated model
+                if ~isempty(yunits)
+                  tout(ii,kk).setYunits(yunits);
+                end
+              end
+            end
+            out(jj) = matrix(tout);
+          case 'ao'
+            % get params from the pest object
+            prms = pst.y;
+            % build cdata aos
+            prmsao = ao.initObjectWithSize(numel(prms),1);
+            for ii = 1:numel(prms)
+              prmsao(ii) = ao(cdata(prms(ii)));
+              prmsao(ii).setYunits(pst.yunits(ii));
+            end
+            % build matrix for parameters
+            prm = matrix(prmsao);
+            % build matrix for the model
+            mm = ao.initObjectWithSize(numel(m(1).objs),numel(prms));
+            for ii = 1:numel(m(1).objs)
+              for kk = 1:numel(prms)
+                mm(ii,kk) = m(kk).getObjectAtIndex(ii);
+              end
+            end
+            mmat = matrix(mm);
+            % eval model
+            tout = mmat*prm;
+            out(jj) = tout;
+        end
+      otherwise
+        error('### current version of pest/eval needs the ''models'' field to be a smodel')
+    end
+    
+    
+    
+    % uncertainties for the evaluated model: calculate them from covariance matrix
+    if ~isempty(pst.cov) && utils.prog.yes2true(pl.find('errors'));
+      switch class(m)
+        case 'smodel'
+          C = pst.cov;
+          p = pst.names;
+          % here we need a matrix of "functions" which are the derivatives wrt parameters,
+          % evaluated at each point x:
+
+          F = [];
+          for kk = 1:length(p)
+            md = eval(diff(m(index), plist('var', p{kk})));
+            F = [F md.y];
+          end
+          % The formula is:
+          % D = F * C * F';
+          % and then we need to take
+          % dy = sqrt(diag(D))
+          if size(md.y, 1) > 1
+            % Make sure we work with columns
+            out(jj).setDy(sqrt(sum((F * C)' .* F'))');
+          else
+            out(jj).setDy(sqrt(sum((F' * C)' .* F)));
+          end
+          
+        otherwise
+          warning('Propagation of the errors on the model not yet implemented')
+      end
+    end
+    
+    % Set output AO name
+    name = sprintf('eval(%s,', pst.name);
+    for kk = 2:numel(pst)
+      name = [name pst(kk).name ','];
+    end
+    name = [name(1:end-1) ')'];
+    out(jj).name = name;
+    % Add history
+    if isempty(as)
+      out(jj).addHistory(getInfo('None'), pl, pst_invars, pst(:).hist);
+    else
+      out(jj).addHistory(getInfo('None'), pl, {pst_invars as_invars}, [pst(:).hist as(:).hist]);
+    end
+  end
+  
+  % Set output
+  varargout = utils.helper.setoutputs(nargout, 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, 'pest', 'ltpda', utils.const.categories.sigproc, '$Id: eval.m,v 1.28 2011/05/15 22:47:14 mauro Exp $', sets, pl);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function plout = getDefaultPlist()
+  persistent pl;  
+  if ~exist('pl', 'var') || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;  
+end
+
+function pl = buildplist()
+  
+  pl = plist();
+  
+  % INDEX
+  p = param({'index', 'Select which model must be evaluated if more than one.'}, paramValue.DOUBLE_VALUE(1));
+  pl.append(p);
+  
+  % XDATA
+  p = param({'Xdata', ['The X values to evaluate the model at. This can be:<ul>'...
+    '<li>a double vector </li>' ...
+    '<li>a cell array of double vectors</li>' ...
+    '<li>a single AO (from which the Y data will be extracted)</li>' ...
+    '<li>a cell array of AOs (from which the Y data will be extracted)</li></ul>' ...
+    ]}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % XFIELD
+  p = param({'xfield', 'Choose the field to extract the x values from when inputting AOs for parameter ''xdata''.'},  {2, {'x', 'y'}, paramValue.SINGLE});
+  pl.append(p);
+  
+  % TYPE
+  pv = paramValue.DATA_TYPES;
+  % Add an 'empty' on top of the list
+  pv{2} = [{''} pv{2}];
+  p = param({'type', ['Choose the data type for the output ao.<br>'...
+    'If empty, and if the user input AOs as ''XDATA'', the type will be inherited.']},  pv);
+  p.val.setValIndex(1);
+  pl.append(p);
+  
+  % YUNITS
+  p = param({'yunits','Unit on Y axis.'},  paramValue.STRING_VALUE(''));
+  pl.append(p);
+  
+  % XUNITS
+  p = param({'xunits','Unit on X axis.'},  paramValue.STRING_VALUE(''));
+  pl.append(p);
+  
+  % ERRORS
+  p = param({'errors', ['Estimate the uncertainty of the output values based <br>' ...
+    'on the parameters covariance matrix']}, paramValue.TRUE_FALSE);
+  pl.append(p);
+  
+end
+% END