% 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.
%
% CALL:        [b smplr] = mcmc(out,pl)
%
% INPUTS:      out     -  matrix objects with measured outputs
%              pl      -  parameter list
%
% OUTPUTS:     b   - pest object contatining estimate information
%              smplr - matrix containing info about the rejected points
%
% <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mcmc')">Parameters Description</a>
%
% VERSION:    $Id: mcmc.m,v 1.35 2011/11/16 15:21:13 nikos Exp $
%
% References:  M Nofrarias et al. Phys. Rev. D 82, 122002 (2010)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TODO: multiple chain option not implemented yet
%       metropolis/hastings not implemented
%       empty initial values 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
[mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);

% 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');
flim = find(pl,'frequencies');
fsout = find(pl,'fsout');
search = find(pl,'search');
simplex = find(pl,'simplex');
mhsampleTrue = find(pl,'mhsample');
xo = find(pl,'x0');
mtxin = find(pl,'input');
mdlin = find(pl,'model');
mdlFreqDependent = find(pl,'modelFreqDependent');
nse = find(pl,'noise');
jumps = find(pl,'jumps');
parplot = find(pl,'plot');
debug = find(pl,'debug');
outModel = find(pl,'outModel');
inModel = find(pl,'inModel');
outNames = find(pl,'outNames');
inNames = find(pl,'inNames');
fpars = find(pl,'prior');
anneal = find(pl,'anneal');
DeltaL = find(pl,'DeltaL');
SNR0 = find(pl,('SNR0'));

% Decide on a deep copy or a modify
in  = copy(mtxin, nargout);
out = copy(mtxs, nargout);
mdl = copy(mdlin, nargout);

% checking if input or noise are aos
if isa(in, 'ao')
    in = matrix(in,plist('shape',[size(in,1) size(in,2)]));
end
if isa(nse, 'ao')
    nse = matrix(nse,plist('shape',[size(nse,1) size(nse,2)]));
end

% lighten the model
mdl.clearHistory;
if isa(mdl, 'ssm')
  mdl.clearNumParams;
  mdl.clearAllUnits;
  mdl.params.simplify;
end

% Check input parameters
if isempty(rng)
  error('### Please define a search range ''range''');
end
param = find(pl,'FitParams');
if isempty(param)
  error('### Please define the parameters ''param''');
end
nparam = numel(param);
if size(mdl) ~= size(nse)
  error('### Parameters ''model'' and ''noise'' must be the same dimension');
end
if size(in) ~= size(out)
  error('### Number of input and output experiments must be the same');
end


% Get range for parameters
for i = 1:nparam
  rang(:,i) = rng{i};
end

% Get number of experiments
nexp = numel(in);

switch class(mdl)
  case 'matrix'
    % Check model sizes
    if (isempty(outModel) && isempty(inModel))
      if (~(numel(in(1).objs) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
        error('Check model or input/output sizes');
      end
    elseif isempty(inModel)
      if (~(numel(in(1).objs) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(size(outModel,1) == numel(out(1).objs)))
        error('Check model or input/output sizes');
      end
    elseif isempty(outModel)
      if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
        error('Check model or input/output sizes');
      end
    else
      if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(numel(out(1).objs) == size(outModel,1)))
        error('Check model or input/output sizes');
      end
    end
  case 'ssm'
    inNames = find(pl,'inNames');
    outNames = find(pl,'outNames');
    if( (numel(inNames) ~= numel(in(1).objs)) || numel(outNames) ~= numel(out(1).objs))
      error('Check model inNames and outNames, they do not match with the input objects')
    end
  otherwise
    error('### Model must be either from the ''matrix'' or the ''ssm'' class. Please check the inputs')
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Frequency domain pre-process  %%%%%%%%%%%%%%%%%%%%%%%%

utils.helper.msg(msg.IMPORTANT, 'Computing frequency domain objects', mfilename('class'), mfilename);
for k = 1:nexp
  
  if((numel(out(1).objs) == 2) && (numel(in(1).objs) == 2))  % previous code done by miquel
    % optional resampling (not matrix/resample method yet)
    if ~isempty(fsout)
      for i = 1:numel(in(k).objs(:))
        in(k).objs(i).resample(plist('fsout',fsout));
        out(k).objs(i).resample(plist('fsout',fsout));
      end
    end
    
    % fft
    fin(k) = fft(in(k));
    fout(k) = fft(out(k));
    
    Nsamples = length(fin(k).getObjectAtIndex(1).x);
    fs = fin(k).getObjectAtIndex(1).fs;
    
    
    if (isempty(outModel))
      % Get rid of fft f =0, reduce frequency range if needed
      if ~isempty(flim)
        fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
        fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
      end
    elseif (~isempty(outModel))
      % Get rid of fft f =0, reduce frequency range if needed
      if ~isempty(flim)
        fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
        fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
        if(~isempty(outModel))
          for lll=1:size(outModel,1)
            for kkk=1:size(outModel,2)
              outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
            end
          end
        end
        if(~isempty(inModel))
          inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
        end
      end
    end
    
    % use signal fft to get frequency vector. Take into account signal
    % could be empty or set to zero
    
    % 1st channel
    if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
      i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
      % rebuild input with new zeroed signal
      fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
    else
      freqs = fin(k).getObjectAtIndex(1,1).x;
    end
    % 2nd channel
    if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
      i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
      % rebuild input with new zeroed signal
      fin(k) = matrix(fin(k).getObjectAtIndex(1,1),i2,plist('shape',[2 1]));
    else
      freqs = fin(k).getObjectAtIndex(2,1).x;
    end
    
    % Compute psd
    n1  = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
    n2  = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(2,1), pl);
    n12 = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(1,1),nse(k).getObjectAtIndex(2,1), pl);
    
    % interpolate noise to fft frequencies
    S11 = interp(n1,plist('vertices',freqs));
    S12 = interp(n12,plist('vertices',freqs));
    S22 = interp(n2,plist('vertices',freqs));
    S21 = conj(S12);
    % build cross-spectrum matrix
    mnse(k) = matrix(S11,S21,S12,S22,plist('shape',[2 2]));
    
    switch class(mdl)
      case 'matrix'
        for i = 1:numel(mdl.objs)
          if (mdlFreqDependent)
            % set Xvals
            mdl.objs(i).setXvals(freqs);
            % set alias
            mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
          else
            mdl.objs(i).setXvals(1);
          end
        end
      case 'ssm'
        if k<2  
        mdl.clearNumParams;
        spl = plist('set', 'for bode', ...
          'outputs', outNames, ...
          'inputs', inNames);
        % first optimise our model for the case in hand
        mdl.reorganize(spl);
        % make it lighter
        
        %for k = 1:nexp
            
          mdl(k).optimiseForFitting();
        %end
        end
    end
    % set Xvals for model (frequency vector could change for each experiment)
    
    mmdl(k) = mdl;
    
  elseif ((numel(out(1).objs) == 3) && (numel(in(1).objs) == 4))
    
    % here we are implementing only the magnetic case
    % We have 4 inputs (the 4 conformator waveforms of the magnetic
    % analysis and
    % 3 outputs (that correspond to the IFO.x12 and IFO.ETA1 and
    % IFO.PHI1
    
    
    % fft
    fin(k) = fft(in(k));
    fout(k) = fft(out(k));
    Nsamples = length(fin(k).getObjectAtIndex(1,1).x);
    fs = fin(k).getObjectAtIndex(1).fs;
    
    % Get rid of fft f =0, reduce frequency range if needed
    if ~isempty(flim)
      fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
      fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
      
      if(~isempty(outModel))
        for lll=1:size(outModel,1)
          for kkk=1:size(outModel,2)
            outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
          end
        end
      end
      if(~isempty(inModel))
        inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
      end
    end
    
    freqs = fin(k).getObjectAtIndex(1,1).x;
    
    % Build noise model
    for ii = 1:numel(out.objs)
      for jj = ii:numel(out.objs)
        % Compute psd
        if (ii==jj)
          n(ii,jj)  = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(ii), pl);
          S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
        else
          n(ii,jj) = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(ii),nse(k).getObjectAtIndex(jj),pl);
          S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
          S(jj,ii) = conj(S(ii,jj));
        end
      end
    end
    % build cross-spectrum matrix
    mnse(k) = matrix(S,plist('shape',[numel(out.objs) numel(out.objs)]));
    
    % Check model sizes
    
    switch class(mdl)
      case 'matrix'
        % set Xvals for model (frequency vector could change for each experiment)
        for i = 1:numel(mdl.objs)
          % set Xvals
          if (mdlFreqDependent)
            mdl.objs(i).setXvals(freqs);
          else
            mdl.objs(i).setXvals(1);
          end
        end
      case 'ssm'
        spl = plist('set', 'for bode', ...
          'outputs', outNames, ...
          'inputs', inNames);
        % first optimise our model for the case in hand
        mdl.reorganize(spl);
        % make it lighter
        for k = 1:nexp
          mmdl(k).optimiseForFitting();
        end
    end
    mmdl(k) = mdl;
    
  elseif ((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))
      
      % case 1 input, 1 output used mostly for passing
      % ao objects from ao.mcmc.
      
    % optional resampling (not matrix/resample method yet)
    if ~isempty(fsout)
        in(k).objs(1).resample(plist('fsout',fsout));
        out(k).objs(1).resample(plist('fsout',fsout));
    end
    
    % fft
    fin(k) = fft(in(k));
    fout(k) = fft(out(k));
    
    Nsamples = length(fin(k).getObjectAtIndex(1).x);
    fs = fin(k).getObjectAtIndex(1).fs;
    
    
    if (isempty(outModel))
      % Get rid of fft f =0, reduce frequency range if needed
      if ~isempty(flim)
        fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
        fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
      end
    elseif (~isempty(outModel))
      % Get rid of fft f =0, reduce frequency range if needed
      if ~isempty(flim)
        fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
        fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
        if(~isempty(outModel))
          for lll=1:size(outModel,1)
            for kkk=1:size(outModel,2)
              outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
            end
          end
        end
        if(~isempty(inModel))
          inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
        end
      end
    end
    
    % use signal fft to get frequency vector. Take into account signal
    % could be empty or set to zero
    
    % 1st and only channel
    if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
      i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
      % rebuild input with new zeroed signal
      fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
    else
      freqs = fin(k).getObjectAtIndex(1,1).x;
    end
    
    % Compute psd
    n1  = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
    
    % interpolate noise to fft frequencies
    S11 = interp(n1,plist('vertices',freqs));

    % build cross-spectrum matrix
    mnse(k) = matrix(S11);
    
    switch class(mdl)
      case 'matrix'
        for i = 1:numel(mdl.objs)
          if (mdlFreqDependent)
            % set Xvals
            mdl.objs(i).setXvals(freqs);
            % set alias
            mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
          else
            mdl.objs(i).setXvals(1);
          end
        end
      case 'ssm'
       if k<2   
        mdl.clearNumParams;
        spl = plist('set', 'for bode', ...
          'outputs', outNames, ...
          'inputs', inNames);
        % first optimise our model for the case in hand
        mdl.reorganize(spl);
        % make it lighter
          mdl(k).optimiseForFitting();
       end
    end
    % set Xvals for model (frequency vector could change for each experiment)
    
    mmdl(k) = mdl;

    
  else
    error('Implemented cases: 1 inputs / 1output, 2 inputs / 2outputs (TN3045 analysis), and 4 inputs / 3 outpus (magnetic complete analysis model. Other cases have not been implemented yet. Sorry for the inconvenience)');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% do simplex
if simplex
  if isempty(xo)
    error('### Simplex needs a starting guess. Please input a ''x0''.');
  else
    switch class(mmdl)
      case 'matrix'
        xo = fminsearch(@(x) utils.math.loglikelihood_matrix(x,fin,fout,mnse,mmdl,param,inModel,outModel),xo,optimset('Display','iter'));
        %         case 'smodel'
        %           xo = fminsearch(@(x) utils.math.loglikelihood(x,in,out,nse,mmdl,param),xo,optimset('Display','iter'));
      case 'ssm'
        Amats = mdl.amats;
        Bmats = mdl.bmats;
        Cmats = mdl.cmats;
        Dmats = mdl.dmats;
        xo = fminsearch(@(x) utils.math.loglikelihood_ssm(x,fin,fout,mnse,mmdl,param,inNames,outNames,spl,Amats,Bmats,Cmats,Dmats),xo,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(param)
      fprintf('###  Simplex estimate: %s = %d \n',param{i},xo(i))
      save('parameters_simplex.txt','xo','-ASCII')
    end
  end
  
  p = pest(xo);
  p.setName('SimplexEstimate');
  p.setNames(param{:});
  p.setModels(mmdl);
  
  
end

if mhsampleTrue
  
  % Sampling and saving rejected points. smplr is a matrix of
  % the size of (# of samples)x(# of params + 2) wich contains the whole
  % chain for each parameter, a column of zeros and ones (rejected or 
  % accepted point) and the value of the loglikelihood.
  [smpl smplr] = utils.math.mhsample(mmdl,fin,fout,mnse,cvar,N,rang,param,Tc,xi,xo,search,jumps,parplot,debug,inNames,outNames,fpars,anneal,SNR0,DeltaL,inModel,outModel);

  % 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,:));
  
  for i = 1:nparam
    for j = 1:nparam
      cr(i,j) = cv(i,j)/sqrt(cv(i,i)*cv(j,j));
    end
  end
  
  % compute dof
  ndata = 0;
  for i = 1:length(mtxs)
    for j = 1:length(mtxs(1).objs)
      ndata = ndata + mtxs(i).objs(j).len;
    end
  end
  dof = ndata-nparam;
  
  % create pest output
  p = pest(mn);
  p.setName('mcmc');
  p.setNames(param{:});
  % add statistical info
  p.setCov(cv);
  p.setCorr(cr);
  p.setDy(sqrt(diag(cv)));
  p.setChain(smpl);
  %p.computePdf;
  p.setDof(dof);
  p.setModels(mdl);
  % set history
  %p.addHistory(getInfo('None'), pl, mtxs_invars(:), [out(:).hist mdl(:).hist]);
end
% Set output
% Implemented it this way in order to leave the pest class as it is. 
output = {p,smplr};
%varargout{1} = p;
varargout = output(1:nargout);

end

% %--------------------------------------------------------------------------
% % Compute priors means and sigmas
% %--------------------------------------------------------------------------
% function fit = fitPrior(prior,nparam) 
% fit = [];
% for ii=1:nparam
% [muhat,sigmahat] = normfit(prior(:,2*ii-1),prior(:,2*ii));
% fit = [fit ,[muhat,sigmahat]'];
% 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: mcmc.m,v 1.35 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();

% % inNames
% p = param({'inNames','Input names. Used for ssm models'}, paramValue.EMPTY_STRING);
% pl.append(p);
% 
% % outNames
% p = param({'outNames','Output names. Usde for ssm models'}, paramValue.EMPTY_STRING);
% pl.append(p);
%
% Model
% p = param({'model','A matrix array of models.'}, paramValue.EMPTY_STRING);
% 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);

pl = plist.MCH_FIT_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);

% Noise
p = param({'noise','A matrix array of noise spectrum (PSD) used to compute the likelihood.'}, paramValue.EMPTY_STRING);
pl.append(p);

% Search
p = param({'modelFreqDependent','Set to true to use frequency dependent s models, set to false when using constant models'}, paramValue.TRUE_FALSE);
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);

% 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);

% Resample
p = param({'fsout','Desired sampling frequency to resample the input time series'}, paramValue.EMPTY_DOUBLE);
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);

% mhsampleTrue
p = param({'mhsample','Set to true to perform a mhsample search. This is set to true by default. Only to be set to false by the user if we does not want to perform the mcmc search'}, 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);

% heat
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);

% heat
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);

% inModel
p = param({'inModel','Input model. Still under test'}, paramValue.EMPTY_STRING);
pl.append(p);

% outModel
p = param({'outModel','Output model. Still under test'}, paramValue.EMPTY_STRING);
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)