--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@matrix/modelselect.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,482 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%   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
+
+