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