% MCMC_TD estimates paramters using a Monte Carlo Markov Chain.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: MCMC estimate the parameters of a given model given
%              inputs, outputs and noise using a Metropolis algorithm.
%              
%
% CALL:        b = mcmc(out,pl)
%
% INPUTS:      
%              out     - analysis objects with measured outputs
%              pl      - parameter list
%
% OUTPUTS:     b   - pest object contatining estimate information
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'mcmc_td')">Parameters Description</a>
%
% VERSION:    $Id: mcmc_td.m,v 1.6 2011/04/08 08:56:12 hewitson Exp $
%
% References:  "Catching supermassive black holes binaries without a net"
%              N.J. Cornish, E.K. Porter, Phys.Rev.D 75, 021301, 2007
%
% TODO: multiple chain option not implemented yet
%       metropolis/hastings not implemented
%       empty initial values not implemented
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = mcmc_td(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);
  
  % Method can not be used as a modifier
  if nargout == 0
    error('### metropolis2D cannot be used as a modifier. Please give an output variable.');
  end
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs smodels and plists
  [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);
  
  out = bs;
  
  
  % Combine plists
  pl = parse(pl, getDefaultPlist);
  
  
  % Get parameters
  N = find(pl, 'N');
  Tc = find(pl,'Tc');
  xi = find(pl,'heat');
  cvar  = find(pl, 'cov');
  rng = find(pl,'range');
  search = find(pl,'search');
  simplex = find(pl,'simplex');
  x0 = find(pl,'x0');
  mdl = find(pl,'model');
  %nse = find(pl,'noise');
  jumps = find(pl,'jumps');
  parplot = find(pl,'plot');
  debug = find(pl,'debug');
  
  in = find(pl,'Input');
  inNames = find(pl,'inNames');
  outNames = find(pl,'outNames');
  Noise = find(pl,'Noise');
  parnames = find(pl,'parnames');
  cutbefore = find(pl,'cutbefore');
  cutafter = find(pl,'cutafter');
  
  % Check input parameters
  if isempty(rng)
    error('### Please define a search range ''range''');
  end

  if isempty(parnames)
    error('### Please define the parameters ''parnames''');
  end
  nparam = numel(parnames);

  
  % Get range for parameters
  for i = 1:nparam
    range(:,i) = rng{i};
  end
  
  
  
  % do simplex
  if simplex
    if isempty(x0)
      error('### Simplex needs a starting guess. Please input a ''x0''.');
    else
      switch class(mdl)
        case 'smodel'
          %xo = fminsearch(@(x) utils.math.loglikehood(x,in,out,nse,mdl,param),xo,optimset('Display','iter'));
        case 'ssm'
          x0 = fminsearch(@(x)utils.math.loglikehood_ssm_td(x,in,out,parnames,mdl,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter),x0,optimset('Display','iter'));
        otherwise
          error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
      end
      
      for i = 1:numel(parnames)
        fprintf('###  Simplex estimate: %s = %d \n',parnames{i},x0(i))
        %save('parameters_simplex.txt','x0','-ASCII')
      end
    end
  end
  
  % sample distribution
  smpl = utils.math.mhsample_td(mdl,in,out,cvar,N,range,parnames,Tc,xi,x0,search,jumps,parplot,debug,inNames,outNames,Noise,cutbefore,cutafter);
  
  % statistics of the chain
  if isempty(Tc)
    initial =1;
  else
    initial = Tc(2)+1;
  end
  
  mn = mean(smpl(initial:end,:));
  cv = cov(smpl(initial:end,:));
  
  cr = utils.math.cov2corr(cv);
  
  % compute pdf
  for i = 1:nparam
    [count,bin] = hist(smpl(initial:end,i));
  end
  
  % create pest output
  p = pest(mn);
  p.setName('metropolis2D');
  p.setNames(parnames{:});
  % add statistical info
  p.setCov(cv);
  p.setCorr(cr);
  p.setDy(sqrt(diag(cv))); % need to change if non-gaussian
  p.setPdf([reshape(bin, 1, []), reshape(count, 1, [])]); % not yet able to store bins
  p.setChain(smpl);
  p.setModels(mdl);
  % set history
  p.addHistory(getInfo('None'), pl, [ao_invars(:)], [bs(:).hist mdl(:).hist]);
  
  % Set output
  varargout{1} = p;
  
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: mcmc_td.m,v 1.6 2011/04/08 08:56:12 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();
  
  % N
  p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
  pl.append(p);
  
  % Sigma
  p = param({'cov','covariance of the gaussian jumping distribution.'}, paramValue.DOUBLE_VALUE(1e-4));
  pl.append(p);
  
  % Param
  p = param({'parnames','A cell array of evaluated parameters.'}, paramValue.EMPTY_CELL);
  pl.append(p);
  
  
  % Model
  p = param({'model','The model for the system, at the moment only ssm is supported'}, paramValue.EMPTY_STRING);
  pl.append(p);
  
  % Search
  p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
  pl.append(p);
  
  % Simplex
  p = param({'simplex','Set to true to perform a simplex search to find the starting parameters of the MCMC chain.'}, paramValue.TRUE_FALSE);
  pl.append(p);
  
  % heat
  p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
  pl.append(p);
  
  % Tc
  p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_STRING);
  pl.append(p);
  
  % x0
  p = param({'x0','The proposed initial values.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % jumps
  p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
    'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
    'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % plot
  p = param({'plot','Select indexes of the parameters to be plotted.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % debug
  p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
  pl.append(p);
  
  % Input
  p = param({'Input','An array of analysis object containing the input signals'},paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % inNames
  p = param({'inNames','Cell array of input port names'}, paramValue.EMPTY_CELL);
  pl.append(p);
  
  % outNames
  p = param({'outNames','Cell array of output port names'}, paramValue.EMPTY_CELL);
  pl.append(p);
  
  
  % Noise
  p = param({'Noise','An array of analysis objects containg noise series'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % cutbefore
  p = param({'cutbefore','the data samples to cut at the starting of the series'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % cutafter
  p = param({'cutafter','the data samples to cut at the ending of the series'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  
end