--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@matrix/mcmc.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,700 @@
+% 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)
\ No newline at end of file