%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Metropolis algorithm
%
% M Nofrarias 15-06-09
%
% $Id: mhsample_td.m,v 1.1 2010/11/26 13:48:48 luigi Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function smpl = mhsample_td(model,in,out,cov,number,limit,parnames,Tc,xi,xo,search,jumps,parplot,dbg_info,inNames,outNames,Noise,cutbefore,cutafter)
  
  import utils.const.*
  
  % compute mean and range
  mn = mean(limit);
  range = diff(limit);
  % initialize
  acc = [];
  nacc = 1;
  nrej = 1;
  loop = 0;
  oldacc = 0;
  oldrej = 0;
  chgsc = 0;
  nparnames = length(parnames);
  
  
  % switch depending on model class (initial likelihood)
  switch class(model)
    case 'matrix'
%       loglk1 = utils.math.loglikehood_matrix(xo,in,out,nse,model,parnames);
%       loglk1 = utils.math.loglikehood_matrix(xo,in,out,nse,model,parnames,inModel,outModel);
    case 'smodel'
%       loglk1 = utils.math.loglikehood(xo,in,out,nse,model,parnames);
    case 'ssm'
      loglk1 = utils.math.loglikehood_ssm_td(xo,in,out,parnames,model,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter);
    otherwise
      error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
  end
  
  % accept the first sample if no search active
  if ~search
    smpl(1,:) = xo;
    nacc = 1;
  else
    smpl = [];
    nacc = 0;
  end
  nrej = 0;
  
  utils.helper.msg(msg.IMPORTANT, 'Starting Monte Carlo sampling', mfilename('class'), mfilename);
  
  
  hjump = 0;
  % main loop
  while(nacc<number)
    %    jumping criterion during search phase
    %    - 2-sigma  by default
    %    - 10-sigma  jumps at mod(10) samples
    %    - 100-sigma jumps at mod(25) samples
    %    - 1000-sigma jumps at mod(100) samples
    if search
      if nacc <= Tc(1)
        if(mod(nacc,10) == 0 && mod(nacc,25) ~= 0 && mod(nacc,100) ~= 0 && hjump  ~= 1)
          hjump = 1;
          xn = mvnrnd(xo,jumps(2)^2*cov);
        elseif(mod(nacc,25) == 0 && mod(nacc,100) ~= 0 && hjump  ~= 1)
          hjump = 1;
          xn = mvnrnd(xo,jumps(3)^2*cov);
        elseif(mod(nacc,100) == 0 && hjump  ~= 1)
          hjump = 1;
          xn = mvnrnd(xo,jumps(4)^2*cov);
          %           xn = xo + range/2+rand(size(xo));
        else
          hjump = 0;
          xn = mvnrnd(xo,jumps(1)^2*cov);
        end
      else
        xn = mvnrnd(xo,cov);
      end
    else
      xn = mvnrnd(xo,cov);
    end
    
    % check if out of limits
    if( any(xn < limit(1,:)) || any(xn > limit(2,:)))
      logprior =  inf;
      loglk2 = inf;
      betta = inf;
    else
      % compute annealing
      if ~isempty(Tc)
        if nacc <= Tc(1)
          betta = 1/2 * 10^(-xi*(1-Tc(1)/Tc(2)));
        elseif Tc(1) < nacc  && nacc <= Tc(2)
          betta = 1/2 * 10^(-xi*(1-nacc/Tc(2)));
        else
          betta = 1/2;
        end
      else
        betta = 1/2;
      end
      % switch depending on model class (new likelihood)
      switch class(model)
        case 'matrix'
%           loglk2 = utils.math.loglikehood_matrix(xn,in,out,nse,model,parnames,inModel,outModel);
        case 'smodel'
%           loglk2 = utils.math.loglikehood(xn,in,out,nse,model,parnames);
        case 'ssm'
          loglk2 = utils.math.loglikehood_ssm_td(xn,in,out,parnames,model,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter);
        otherwise
          error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
      end
    end
    
    % likelihood ratio
    logr = betta*(loglk2 - loglk1) ;
    % decide if sample is accepted or not
    if logr < 0
      xo = xn;
      nacc = nacc+1;
      smpl(nacc,:) = xn;
      loglk1 =loglk2;
      if dbg_info
        utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d  betta: %d  ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
      end
    elseif rand(1) > (1 - exp(-logr))
      xo = xn;
      nacc = nacc+1;
      smpl(nacc,:) = xn;
      loglk1 =loglk2;
      if dbg_info
        utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d  betta: %d  ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
      end
    else
      nrej = nrej+1;
      if dbg_info
        utils.helper.msg(msg.IMPORTANT, sprintf('rej.\t loglik: %d -> %d  betta: %d  ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
      end
    end
    
    % display and save
    if(mod(nacc,100) == 0 && nacc ~= oldacc)
      updacc = nacc-oldacc;
      updrej = nrej-oldrej;
      ratio(nacc/10,:) = updacc/(updacc+updrej);
      utils.helper.msg(msg.IMPORTANT, sprintf('accepted: %d   rejected : %d   acc. rate: %4.2f',nacc,updrej,ratio(end)), mfilename('class'), mfilename);
      for i = 1:numel(parnames)
        fprintf('###  Parameters: %s = %d \n',parnames{i},xn(i))
      end
      oldacc = nacc;
      oldrej = nrej;
      save('chain.txt','smpl','-ASCII')
      save('acceptance.txt','ratio','-ASCII')
    end
    
    % plot 
    if ((mod(nacc,100) == 0) && (nacc ~= 0) && ~isempty(parplot))
      for i = 1:numel(parplot)
        figure(parplot(i))
        plot(smpl(:,parplot(i)))
      end
    end
    
  end
end