--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@pest/mcmcPlot.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,223 @@
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION:  simple tool that plots mcmc pest objects 
+%
+% CALL: mcmcPlot(pest_obj,pl)
+%
+% Parameters: - pest_obj: pest object
+%             - pl: plist
+%
+%    example: - mcmcPlot(p,plist('plotmatrix',true,'burnin',5000,'pdfs',true,'chain',[1 2 3 4 5 6]))
+%
+%<a href="matlab:utils.helper.displayMethodInfo('pest', 'mcmcPlot')">ParametersDescription</a>
+%
+% Nikos Oct 2011
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = mcmcPlot(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);
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  % Collect all AOs and plists
+  [pests, pest_invars] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
+  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+  
+  % Decide on a deep copy or a modify
+  p = copy(pests, nargout);
+  
+  % combine plists
+  pl = parse(pl, getDefaultPlist());
+  BurnIn = find(pl, 'burnin');
+  nbins = find(pl, 'nbins');
+  paramarray = find(pl, 'chain');
+  %colorm = find(pl, 'colormap');
+    
+  if ~all(isa(pests, 'pest'))
+    error('### mcmcPlot must be only applied to pest objects.');
+  end
+  
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ 
+ outfigs = [];
+ N = numel(p); 
+ 
+ if (BurnIn == 1 && ((find(pl, 'plotmatrix'))))
+     utils.helper.msg(msg.IMPORTANT, sprintf(['The burn-in field is left empty or equal to one. For '...
+         'better and more accurate display the burn-in section of the chains should be discarded.']));
+ elseif (BurnIn == 1 && ((find(pl, 'pdfs'))))
+     utils.helper.msg(msg.IMPORTANT, sprintf(['The burn-in field is left empty or equal to one. For '...
+         'better and more accurate display the burn-in section of the chains should be discarded.']));
+ end
+ 
+ for numpest=1:N     
+     
+    % compute PDF
+    chain=p(numpest).chain(:,2:size(p(numpest).chain,2));
+    p(numpest).computePdf(plist('BurnIn',BurnIn,'nbins',nbins));
+     
+    if isempty(paramarray)
+      % plot chain field (skip 1st column where the Loglikelihood is stored)
+      outfigs = [outfigs ; figure]; 
+      data = plot(chain);
+    else
+      ch = p(numpest).chain(:,paramarray); 
+      outfigs = [outfigs ; figure];  
+      data = plot(ch);
+    end 
+    
+    if (find(pl, 'plotmatrix'));
+      chn = p(numpest).chain(BurnIn:size(chain,1),2:size(p(numpest).chain,2));    
+      outfigs = [outfigs ; figure];     
+      plotmatrix(chn);
+    end
+    
+    if (find(pl, 'results'));
+      chainn = chain(BurnIn:size(chain,1),:);    
+      utils.helper.msg(msg.IMPORTANT, sprintf('Results:'));
+      for ii = 1:(size(chainn,2))
+        mu = mean(chainn(:,ii)); 
+        sigma = std(chainn(:,ii));
+        res = [mu sigma];     
+        utils.helper.msg(msg.IMPORTANT, sprintf(' %d \t',res));
+      end
+    end
+    
+    if (find(pl, 'pdfs'));
+      outfigs = [outfigs ; figure];   
+      
+      if ~(find(pl, 'plotmatrix'))  
+        chn = p.chain(BurnIn:size(p.chain(:,:),1),:);
+      end
+      
+      a=p(numpest).pdf;
+      a(:,1) = [];
+      a(:,1) = [];
+      
+      for kk =1:size(chn,2)
+      subplot(2,4,kk)
+      x = linspace(min(a(:,2*kk-1)),max(a(:,2*kk-1)),10);
+      h = bar(a(:,2*kk-1),a(:,2*kk));
+      hold on;
+      y=normpdf(x,mean(chn(:,kk)),std(chn(:,kk)));       
+      s=sum(y);
+      y=y/s;
+      plot(x,y,'r-','LineWidth',2);
+      hold off;
+
+      shading interp          % Needed to graduate colors
+        
+      ch = get(h,'Children');
+      fvd = get(ch,'Faces');
+      fvcd = get(ch,'FaceVertexCData');
+      n = 10;
+      [zs, izs] = sortrows(a(:,2*kk),1);
+      k = 128;                % Number of colors in color table
+      colormap(summer(k));    % Expand the previous colormap
+      shading interp          % Needed to graduate colors
+      for i = 1:n
+          color = floor(k*i/n);       % Interpolate a color index
+          row = izs(i);               % Look up actual row # in data
+          fvcd(fvd(row,1)) = 1;       % Color base vertices 1st index
+          fvcd(fvd(row,4)) = 1;    
+          fvcd(fvd(row,2)) = color;   % Assign top vertices color
+          fvcd(fvd(row,3)) = color;
+      end
+      set(ch,'FaceVertexCData', fvcd);  % Apply the vertex coloring
+      set(ch,'EdgeColor','k')           % Give bars black borders
+      end
+        
+        
+        
+    end
+    
+    
+ end
+ 
+ 
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ 
+  if nargout == 0
+      out = outfigs;
+    else
+      error('### mcmcPlot cannot be used as a modifier!');   
+  end
+       
+  % Set outputs
+  if nargout > 0
+    varargout{1} = out;
+  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, 'pest', 'ltpda', utils.const.categories.helper, '$Id: computePdf.m,v 1.2 2011/06/06 14:02:12 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();
+
+p = param({'chain',['Insert an array containing the parameters to plot. If left empty,'...
+    'then by default will plot the chains of every parameter. (note: The loglikelihood is stored '...
+    'in the first column)']}, paramValue.DOUBLE_VALUE([]));
+pl.append(p);
+
+p = param({'BurnIn',['Number of samples (of the chains) to be discarded for the computation of the PDFs of the parameters. Also used'...
+                    'for producing the plotmatrix figure.']}, paramValue.DOUBLE_VALUE(1));
+pl.append(p);
+
+p = param({'nbins','Number of bins of the pdf histogram computed for every parameter (used again for the computation of the PDFs of the parameters)'}, paramValue.DOUBLE_VALUE(10));
+pl.append(p);
+
+p = param({'plotmatrix','Boolean to determine if a plotmatrix is desired'}, {1, {false,true}, paramValue.OPTIONAL});
+pl.append(p);
+
+p = param({'pdfs','Boolean to determine if a plot of the PDFs of each parameter is desired'}, {1, {false,true}, paramValue.OPTIONAL});
+pl.append(p);
+
+%p = param({'colormap','Choose a default matlab colormap for the parameter histogarms.'}, paramValue.DOUBLE_VALUE(summer));
+%pl.append(p);
+
+p = param({'results',['Set to "true" if a table of the results of the estimated parameters is desired.'...
+    'The results are printed on screen in 2 columns: the 1st contains the mean value'....
+    'and the second the sigma. Burn-in field is requiered.']}, {1, {false,true}, paramValue.OPTIONAL});
+pl.append(p);
+
+end
+