% 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