%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   modelselect.m - method to compute the Bayes Factor using
%                   reverse jump mcmc algorithm
%
%            call - Bxy = modelselect(out,pl)
%
%          inputs - out: matrix objects with measured outputs
%                   pl:  parameter list
%
%         outputs - Bxy: an array of AOs containing the evolution
%                   of the Bayes factors. (comparing each model 
%                   with each other)
%
%      
% <a href="matlab:utils.helper.displayMethodInfo('matrix','modelselect')">Parameters Description</a>
%
%   N. Karnesis 27/09/2011 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = modelselect(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');
param = find(pl,'FitParams');
cvar  = find(pl, 'cov');
rng = find(pl,'range');
flim = find(pl,'frequencies');
fsout = find(pl,'fsout');
search = find(pl,'search');
xo = find(pl,'x0');
mtxin = find(pl,'input');
mdlin = find(pl,'models');
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,'inModel1');
outNames = find(pl,'outNames');
inNames = find(pl,'inNames');
anneal = find(pl,'anneal');
DeltaL = find(pl,'DeltaL');
SNR0 = find(pl,('SNR0'));

nummodels = numel(mdlin);

% Decide on a deep copy or a modify
in  = copy(mtxin, nargout);
out = copy(mtxs, nargout);
for k = 1:nummodels
  mdl(k) = copy(mdlin{1,k}, nargout);
  nparam(k) = numel(param{1,k});
  % lighten the model
  mdl(k).clearHistory;
  if isa(mdl(k), 'ssm')
  mdl(k).clearNumParams;
  mdl(k).clearAllUnits;
  mdl(k).params.simplify;
  end
end



% Check input parameters
if isempty(rng)
  error('### Please define a search range ''range''');
end
if isempty(param)
  error('### Please define the parameters ''param''');
end
if size(mdl(1)) ~= 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 jj = 1:nummodels
 for i = 1:numel(param{1,jj})   
  rang{1,jj}(:,i) = rng{1,jj}{i}';
 end
end
% Get number of experiments
nexp = numel(in);

switch class(mdl)
  case 'matrix'
    % Check model sizes
    for k = 1:nummodels
    if (isempty(outModel) && isempty(inModel))
      if (~(numel(in(1).objs) == mdl(k).ncols) || ~(numel(out(1).objs) == mdl(k).nrows))
        error('Check model or input/output sizes');
      end
    elseif isempty(inModel)
      if (~(numel(in(1).objs) == mdl(k).ncols) || ~(size(outModel,2) == mdl(k).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(k).ncols) || ~(numel(out(1).objs) == mdl(k).nrows))
        error('Check model or input/output sizes');
      end
    else
      if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl(k).ncols) || ~(size(outModel,2) == mdl(k).nrows) || ~(numel(out(1).objs) == size(outModel,1)))
        error('Check model or input/output sizes');
      end
    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]));
    
  elseif((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))  
    % 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 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);
    
  else
      error('### Sorry, implemented only for 2 inputs - 2 outputs and 1 input - 1 output only...')
  end

    switch class(mdl)
        
      case 'matrix'
       for jj = 1:nummodels
        for i = 1:numel(mdl(jj).objs)
          if (mdlFreqDependent)
            % set Xvals
            mdl(jj).objs(i).setXvals(freqs);
            % set alias
            mdl(jj).objs(i).assignalias(mdl(jj).objs(i),plist('xvals',freqs));
          else
            mdl(jj).objs(i).setXvals(1);
          end
        end
       end
      case 'ssm'
       if k < 2
        for jj = 1:nummodels
         mdl(jj).clearNumParams;
         spl = plist('set', 'for bode', ...
           'outputs', outNames, ...
           'inputs', inNames);
         % first optimise our model for the case in hand
         mdl(jj).reorganize(spl);
         % make it lighter
         mdl(jj).optimiseForFitting();
        end 
       end
    end
    % add lines to the matrix of models. Each line corresponds to one
    % experiment
    mmdl(k,:) = mdl;
end

  [Bxy LogLambda] = utils.math.rjsample(mmdl,fin,fout,mnse,cvar,N,rang,param,Tc,xi,xo,search,jumps,parplot,debug,inNames,outNames,anneal,SNR0,DeltaL,inModel,outModel); 

% Set output
numfactors = 0;
for gg = 1:(nummodels)
  for dd = gg+1:(nummodels)
    numfactors = numfactors + 1;
    BayesFactor(numfactors) = ao(Bxy{gg,dd});
    BayesFactor(numfactors).setName(sprintf('B%d%d',gg,dd));  
    BayesFactor(numfactors).procinfo = plist('Likelihoods',LogLambda);
  end
end

output = {BayesFactor};
varargout = output(1:nargout);

  
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.33 2011/07/31 17:38:34 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','An array input of models.'}, paramValue.EMPTY_STRING);
pl.append(p);

% Param
p = param({'FitParams','A cell array of evaluated parameters for each model.'}, 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','Cell array containing the covariances of the gaussian jumping distribution for each model.'}, 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);

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

% initial values
p = param({'x0','The proposed initial values (cell array again).'}, paramValue.EMPTY_DOUBLE);
pl.append(p);

% ranges
p = param({'range','ranges'}, 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','If set equal to one, the evolution of the Bayes factor is plotted every 500 steps.'}, 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);

% outModel
p = param({'outModel','Output model. 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