% DETREND detrends the input analysis object using a polynomial of degree N.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: DETREND detrends the input analysis object using a
%              polynomial of degree N.
%
% CALL:        b = detrend(a1,a2,a3,..., pl)
%
% INPUTS:      aN  - a list of analysis objects
%              pl  - a parameter list
%
% OUTPUTS:     b - array of analysis objects
%
%              If the last input argument is a parameter list (plist) it is used.
%              The following parameters are recognised.
%
% NOTE: detrend uses two possible algorithms. A fast C-code implementation
% is used for orders less than 11. If the order is higher than 10, then a
% MATLAB code is used which is typically much slower. You can also force
% the use of the MATLAB code using a plist option. When interpreting the
% resulting coefficients, you must be clear which algorithm was used. For
% the C-code algorithm, the coefficients are scaled from the original by
% 
%     z = 2.*ii/(n-1)-1;
% 
% such that the functional form (for a three-coefficient example) that is
% subtracted from the data can be recovered with:
% 
% fitted_c = c.y(3).*z.^2 + c.y(2).*z + c.y(1); 
% 
% 
% The procinfo field of the output AOs is filled with the following key/value
% pairs:
%
%    'COEFFS' - pest object with coefficients describing the subtracted trend function
%               Note that they correspond to physical trend coefficents
%               only if using the (slower) option 'M-FILE ONLY'. In the
%               case of using the C-code algorith, see the note above.
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'detrend')">Parameters Description</a>
%
% VERSION:     $Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = detrend(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
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  % Decide on a deep copy or a modify
  bs = copy(as, nargout);
  
  % Combine plists
  pl = parse(pl, getDefaultPlist);
  
  % Unpack parameter list
  % Leaving the "N" for backwards compatibility
  order = find(pl, 'N', find(pl, 'order'));
  
  % Loop over analysis objects
  for jj = 1:numel(bs)
    % check data
    if ~isa(bs(jj).data, 'tsdata')
      warning('!!! This method only works with time-series at the moment. Skipping %s', ao_invars{jj});
    else
      % detrend with polynomial
      if find(pl, 'M-FILE ONLY')
        if order >= 0
          [y,c] = polydetrend(bs(jj).x, bs(jj).y, order);
        else
          y = bs(jj).y;
          c = [];
        end
      else
        try
          [y,c] = ltpda_polyreg(bs(jj).y, order);
        catch
          warning('!!! failed to execture ltpda_polyreg.mex. Using m-file call.');
          if order >= 0
            [y,c] = polydetrend(bs(jj).x, bs(jj).y, order);
          else
            y = bs(jj).y;
            c = [];
          end
        end
      end
      % set the values in the y and procinfo fields
      bs(jj).data.setY(y);
      % build a pest object with the coefficients
      units = unit.initObjectWithSize(1, order + 1);
      names = [];
      mdl   = [];
      
      if order >= 0
        for kk = 1:order+1
          ll = order - kk + 1;
          units(kk) = bs(jj).yunits / simplify((unit(bs(jj).xunits)).^(ll));
          names{kk}= ['P' num2str(ll)];
          mdl = [mdl '+P' num2str(ll) '*t.^' num2str(ll) ' '];
        end
        
        coeff = pest(plist(...
          'names', names, ...
          'y', c, ...
          'dy', zeros(size(c)), ...
          'yunits', units ...
          ));
        bs(jj).procinfo = plist('coeffs', coeff.setModels(smodel(mdl)));
      end
      % add name
      bs(jj).name = sprintf('detrend(%s)', bs(jj).name);
      % add history
      bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
      % Clear the errors since they don't make sense anymore
      clearErrors(bs(jj));
    end
  end
  
  
  % Set output
  if nargout == numel(bs)
    % List of outputs
    for ii = 1:numel(bs)
      varargout{ii} = bs(ii);
    end
  else
    % Single output
    varargout{1} = bs;
  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.sigproc, '$Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $', sets, pl);
end

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

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

function pl = buildplist()
  pl = plist();
  
  % order
  p = param({'order',['The order of detrending:<ul>', ...
    '<li>-1 - no detrending</li>', ...
    '<li>0 - subtract mean</li>', ...
    '<li>1 - subtract linear fit</li>', ...
    '<li>N - subtract fit of polynomial, order N</li></ul>']}, {3, {-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9}, paramValue.SINGLE});
  pl.append(p);
  
  % m-file only
  p = param({'M-FILE ONLY','Using SLOWER m-file call'}, paramValue.FALSE_TRUE);
  pl.append(p);
  
end

function [y,p] = polydetrend(varargin)
  % POLYDETREND detrends the input data vector with a polynomial.
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % DESCRIPTION: POLYDETREND detrends the input data vector with a
  %              polynomial.
  %
  % CALL:        [y,p] = polydetrend(x,order);
  %
  % INPUTS:      x     - vector of x values
  %              order - order of polynomial to fit and subtract
  %
  % OUTPUTS:     y   - detrended data
  %              p   - coefficients fitted to the trend function
  %
  % VERSION:     $Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $
  %
  % HISTORY:     30-05-07 M Hewitson
  %                 Creation
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  if nargin == 2
    x     = varargin{1};
    order = varargin{2};
    t     = [1:length(x)].';
  elseif nargin == 3
    t     = varargin{1};
    x     = varargin{2};
    order = varargin{3};
  else
    error('### incorrect inputs.');
  end
  
  % fit polynomial
  p = polyfit(t, x, order);
  
  % make polynomial series
  py = polyval(p, t);
  
  % detrend
  y = x - py;
end
% END