% MCMC 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.
%              It handles (1 input / 1 output) systems, (2 input / 1 output) systems,
%              and (2 input / 2 output) systems.
%
% CALL:        b = mcmc(in,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')">Parameters Description</a>
%
% VERSION:    $Id: mcmc.m,v 1.24 2011/11/16 15:21:13 nikos 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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = mcmc(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('### mcmc 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
  out = copy(as, nargout);
  
  % Combine plists
  pl = parse(pl, getDefaultPlist);
  
  
  % Get parameters from plist
  mtxin = find(pl,'input');
  noi = find(pl,'noise');
  mdlin = find(pl,'model');
  
  % Decide on a deep copy or a modify
  in  = copy(mtxin, nargout);
  mdl = copy(mdlin, nargout);

  numexp = numel(in(1,:));
  numchannels_in = numel(in(:,1));
  numchannels_out = numel(out(:,1));

  if (numel(in) ~= numel(out))
      error('### Dimensions of input and output must be equal.');
  end

  % Redefining the ao inputs as matrix objects 
  for ii = 1:numexp
    matin(ii) = matrix(in(:,ii),plist('shape',[numchannels_in 1]));
    matout(ii) = matrix(out(:,ii),plist('shape',[numchannels_in 1]));
    matnoise(ii) = matrix(noi(:,ii),plist('shape',[numchannels_in 1]));
  end
    
  switch class(mdl)
    case 'smodel'
        matmodel = matrix(mdl);
    case 'ssm'
        matmodel = mdl;
    case 'matrix'
        matmodel = mdl;
  end

  pl = pset(pl,'model',matmodel,'input',matin,'noise',matnoise);
  mcmc_obj = mcmc(matout,pl);


  % Set output
  varargout{1} = mcmc_obj;
  
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.m,v 1.24 2011/11/16 15:21:13 nikos 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({'Fitparams','A cell array of evaluated parameters.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % Input
  p = param({'input','A matrix array of input signals.'}, paramValue.EMPTY_STRING);
  pl.append(p);
  
  % Frequencies
  p = param({'frequencies','Range of frequencies where the analysis is performed. If an array, only first and last are used'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % Noise
  p = param({'noise','An array of noise spectrum (PSD) used to compute the likelihood.<ul>', ...
    '<li>1 channel - Input one object</li>', ...
    '<li>2 channel - Input four objects defining the power spectrum matrix </li>',...
    '</ul>'}, paramValue.EMPTY_STRING);
  pl.append(p);
  
  % Model
  p = param({'model','An array of models.',...
    '<li>1 input / 1 output - Input one object</li>', ...
    '<li>2 input / 1 output - Input two objects</li>', ...
    '<li>2 channel (2 input/ 2 output) - Input four objects defining the transfer function matrix </li>',...
    '</ul>'}, 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_DOUBLE);
  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);
  
  p = param({'prior','Mean, sigma and normalization factor for priors. Still under test'}, paramValue.EMPTY_STRING);
  pl.append(p);

  p = param({'anneal',['Choose type of annealing during sampling. Default value is ',...
    'simulated annealing. Choose "thermo" for annealing with a thermostat.',...
    ' SNR is computed and if it is larger than a fixed value SNR0 (provided also in the plist), ',...
    'then the chains are heated by a factor of (SNR(1)/SNR0)^2. Choosing "simple" ',...
    'the deviation of the loglikelihood of every 10 points in the chains is stored. If this deviation ',...
    'is larger or smaller than two fixed values the chains are cooled or heated respectively.']}, {1, {'simul','thermo', 'simple'}, paramValue.SINGLE});
  pl.append(p);

  p = param({'SNR0','Fixed value for thermostated annealing.'},   {1, {200}, paramValue.OPTIONAL});
  pl.append(p);
          
  p = param({'DeltaL',['Deviation of Loglikelihood for 10 points of the chains. Used for',...
    'the "simple" choice of annealing with a thermostat.']},   {1, {[100 600 2 3]}, paramValue.OPTIONAL});
  pl.append(p);
  
end


% PARAMETERS:  J     - number of chains to compute.
%              sigma - dispersion of the jumping distribution.
%              range - range where the parameteters are sampled.
%              param - a cell array of evaluated parameters (from the
%                      smodel).
%              Tc    - an array with two values (default: 0, i.e. no annealing).
%              heat  - heat index flattening likelihood surface (default: 1)
%              x0    - proposed initial values, if empty selected randomly.
%                      (default: empty)