function [thetaMargMode,UpdateMargMode] = GraphPosteriorDensity(thetaPostSample,NumBurnin,CurrChain,UpdatePostSample,UpdateMode,UpdateNames,UpdatePriorX,UpdatePriorY,UpdatePriorBounds,DSGEModel,CurrINI,controls)
% GraphPosteriorDensity: Creates plots of the posterior densities for all original parameters or estimate
%                        the posterior mode of the marginal distribution.
%
% USAGE:
%
%       [thetaMargMode,UpdateMargMode] = GraphPosteriorDensity(thetaPostSample,NumBurnin,CurrChain,UpdatePostSample,...
%                                                UpdateMode,UpdateNames,UpdatePriorX,UpdatePriorY,UpdatePriorBounds,...
%                                                DSGEModel,CurrINI,controls)
%
% REQUIRED INPUT:  thetaPostSample (matrix) of dimension NumIterxNumParam with all the posterior
%                                  draws of the original parameters.
%
%                  NumBurnin (integer) with the number of burn-in draws.
%
%                  CurrChain (integer) with the number for the used sampling chain.
%
%                  UpdatePostSample (matrix) of dimension (NumIter-NumBurnin x NumAdditional) with
%                                  the post burnin sample draws of the parameters that are functions
%                                  of the estimated parameters.
%
%                  UpdateMode (vector) of dimension NumAdditionalx1 with the posterior mode estimates
%                                  of the parameters that are functions of the estimated parameters.
%
%                  UpdateNames (string matrix) of dimension (NumAdditional x N1) with the names of the
%                                  additional parameters.
%
%                  UpdatePriorX (matrix) of dimension NumAdditional x NumValues with the values on the
%                                  X-axis for the prior densities of the NumAdditional additional
%                                  parameters.
%
%                  UpdatePriorY (matrix) of dimension NumAdditional x NumValues with the density values
%                                  for the Y-axis of the prior densities of the NumAdditional additional
%                                  parameters.
%
%                  UpdatePriorBounds (matrix) of dimension NumAdditional x 2 with suitable lower and
%                                  upper bounds for the additional parameters.
%
%                  DSGEModel (structure)  whose fields contain all required information
%                                  about the model, such as input files, variable and sample
%                                  data. Actual observations are not stored in this file.
%
%                  CurrINI (structure) with initialization information.
%
%                  controls (structure) whose fields provide handles to all the UI
%                                  controls in YADA's main dialog.
%
% OPTIONAL OUTPUT: thetaMargMode (vector) of dimension NumParamx1 with all estimated marginal posterior
%                                  mode values. The estimates are determined from the maximum of
%                                  the Kernel density values.
%
%                  UpdateMargMode (vector) of dimension NumAdditionalx1 with the estimated marginal
%                                  posterior mode values of the parameters that are a function of the
%                                  estimated parameters.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2009 European Central Bank.
%
%                       First version: December 15, 2006.
%                        This version: March 2, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 18-12-2006: Changed the wait dialog text such that the parameter information comes last. The function also
%               saves the figures to disk automatically.
%
% * 19-12-2006: Made sure that x and f vectors have the right dimension.
%
% * 08-01-2007: Placed a StringTrim command around parameter name is estimation wait dialog.
%
% * 08-02-2007: Added a second "drawnow" call at the end of the figure loop.
%
% * 18-04-2007: Changed the plot-style for the prior from dotted to dashed.
%
% * 09-11-2007: Updated the documentation.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 10-03-2008: The value of "CurrChain" is now used in the figure filenames.
%
% * 17-03-2008: Made sure that the y-axis height is meaningful relative to the posterior density height.
%
% * 04-04-2008: The function can now be used to estimate the posterior mode of the marginal distribution.
%
% * 23-05-2008: Updated the documentation.
%
% * 26-05-2008: The function now asks uses a subset of the posterior draws when "DSGEModel.PosteriorDensityValue"
%               is 3 or 4. The input variable controls is now also used.
%
% * 21-07-2008: Added the "LowerBound" input variable when calling "GraphPriorDensity".
%
% * 19-12-2008: Updated the documentation.
%
% * 25-02-2009: Added the input variables "UpdatePostSample", "UpdateMode" and "UpdateNames" so that the posterior
%               densities of additional parameters can also be prepared and plotted.
%
% * 26-02-2009: The function now has two optional output variables "thetaMargMode" and "UpdateMargMode". These
%               replace the old optional output variable "PostMode".
%
% * 02-03-2009: Added the three input variables "UpdatePriorX", "UpdatePriorY", and "UpdatePriorBounds" that
%               hold data for the implied prior densities of the parameters that are a function of the
%               estimated parameters. The function now checks if the user has "YADAksdensity" or "ksdensity"
%               on the path. If not, the normal kernel from the kernel density estimation toolbox is used
%               when the user has set "DSGEModel.PosteriorDensityValue" to unity, i.e., a Gaussian kernel.
%               Increased the number of points in the grid from 2^7 to 2^8 when calling the "YADAMHDensEst"
%               function.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

[NumIter,NumParam] = size(thetaPostSample);
%
%
% Try to load the posterior mode estimates
%
ModeFile = [DSGEModel.OutputDirectory '\mode\PosteriorMode-' DSGEModel.NameOfModel '.mat'];
ErrorStr = '';
try;
   ModeData = load(ModeFile);
catch;
   ErrorStr = ['Unable to load the file "' ModeFile '". Message caught is: ' lasterr];
end;
if isempty(ErrorStr)==0;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
theta = ModeData.theta;
thetaMode = ModeData.thetaMode;
ParameterNames = ModeData.ParameterNames;
thetaDist = ModeData.thetaDist;
PriorDist = ModeData.PriorDist;
LowerBound = ModeData.LowerBound;
%
% locate the controls structure
%
maingui = findobj('Type','figure','Tag','YADA');
controls = get(maingui,'UserData');
%
% directory for saving graphs
%
GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
%
% load the data for plotting the prior densities
%
[Xprior,Yprior,BoundsMatrix] = GraphPriorDensity(theta,thetaDist,LowerBound,PriorDist,ParameterNames,controls,DSGEModel,CurrINI);
%
% check if we should use a sub-sample
%
thetaPostSample = thetaPostSample(NumBurnin+1:NumIter,:);
if DSGEModel.PosteriorDensityValue>=3;
   %
   % select which posterior draws to use
   %
   UseDraws = KeepPosteriorDraws(DSGEModel,controls.posterior.usepostdraws,NumIter-NumBurnin);
   if length(UseDraws)<NumIter-NumBurnin;
      thetaPostSample = thetaPostSample(UseDraws,:);
      if isempty(UpdatePostSample)==0;
         UpdatePostSample = UpdatePostSample(UseDraws,:);
      end;
   end;
end;
%
ParameterNames = ParameterNames.estimated;
if isempty(UpdatePostSample)==1;
   NumAdditional = 0;
else;
   NumAdditional = size(UpdatePostSample,2);
   ParameterNames = strvcat(ParameterNames,UpdateNames);
   thetaPostSample = [thetaPostSample UpdatePostSample];
   thetaMode = [thetaMode;UpdateMode];
end;
%
% first we display a wait dialog
%
if DSGEModel.PosteriorDensityValue>2;
   WaittxtAdd = [' Estimation started at: ' StringTrim(datestr(now,14)) '. Using ' int2str(size(thetaPostSample,1)) ' draws of the ' int2str(NumIter-NumBurnin) ' post-burnin draws available. This may take a very long time...'];
   ShowCancel = 1;
else;
   WaittxtAdd = [' Estimation started at: ' StringTrim(datestr(now,14)) '.'];
   ShowCancel = 0;
end;
txt = ['Please wait while YADA computes kernel density estimates of the marginal posterior densities.' WaittxtAdd ' Current parameter: ' ParameterNames(1,:) '           - 1 out of ' num2str(NumParam+NumAdditional,'%0.0f') '.'];
WaitHandle = WaitDLG(txt,'information','Kernel Density Estimation - Posterior',500,150,CurrINI,ShowCancel);
WaitControls = get(WaitHandle,'UserData');
drawnow;
pause(0.02);
%
% now we try to construct the x-axis and y-axis data for the posterior densities
%
StopEstimation = 0;
for i=1:NumParam+NumAdditional;
   txt = ['Please wait while YADA computes kernel density estimates of the marginal posterior densities.' WaittxtAdd ' Current parameter: ' StringTrim(ParameterNames(i,:)) ' - ' num2str(i,'%0.0f') ' out of ' num2str(NumParam+NumAdditional,'%0.0f') '.'];
   set(WaitControls.text,'String',txt);
   drawnow;
   if DSGEModel.PosteriorDensityValue==1;
      if exist('YADAksdensity')==2;
         [f,x] = YADAksdensity(thetaPostSample(:,i),'npoints',256);
      else;
         if exist('ksdensity')==2;
            [f,x] = ksdensity(thetaPostSample(:,i),'npoints',256);
         else;
            [x,f] = KernelDensityPlot(thetaPostSample(:,i),'knorm');
         end;
      end;
   elseif DSGEModel.PosteriorDensityValue==2;
      [x,f] = YADAMHDensEst(thetaPostSample(:,i),2^8,0,'gaussian',[]);
   elseif DSGEModel.PosteriorDensityValue==3;
      [x,f] = YADAMHDensEst(thetaPostSample(:,i),2^8,-1,'gaussian',WaitControls.cancel);
   elseif DSGEModel.PosteriorDensityValue==4;
      [x,f] = YADAMHDensEst(thetaPostSample(:,i),2^8,-2,'gaussian',WaitControls.cancel);
   end;
   if isempty(x)==1;
      StopEstimation = 1;
      break;
   end;
   if i==1;
      Xposterior = zeros(NumParam+NumAdditional,length(x));
      Yposterior = zeros(NumParam+NumAdditional,length(x));
   end;
   if size(x,1)>size(x,2);
      Xposterior(i,:) = x';
   else;
      Xposterior(i,:) = x;
   end;
   if size(f,1)>size(f,2);
      Yposterior(i,:) = f';
   else;
      Yposterior(i,:) = f;
   end;
end;
%
% close the wai dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);
%
if StopEstimation==1;
   if nargout==2;
      thetaMargMode = [];
      UpdateMargMode = [];
   end;
   return;
end;
%
% determine the mode
%
ModeData = zeros(NumParam+NumAdditional,1);
for i=1:NumParam+NumAdditional;
   [MaxVal,MaxInd] = max(Yposterior(i,:));
   ModeData(i,1) = Xposterior(i,MaxInd);
end;
%
if nargout==2;
   thetaMargMode = ModeData(1:NumParam);
   if NumAdditional>0;
      UpdateMargMode = ModeData(NumParam+1:NumParam+NumAdditional);
   else;
      UpdateMargMode = [];
   end;
   return;
end;
%
% we are ready to do some real plotting!
%
%
% Compute the number of figures (max 16 plots per graph)
%
NumFigs = ceil((NumParam+NumAdditional)/16);
%
if NumFigs==1;
   [NumRows,NumCols] = VillaniConstructOptimalSubplot(NumParam+NumAdditional);
else;
   NumRows = 4;
   NumCols = 4;
end;
NumPlots = 0;
pos = [(CurrINI.scrsz(3)-1000)/2 (CurrINI.scrsz(4)-680)/2 1000 650];
%
%
% now we can start the loops
%
for i=1:NumFigs;
   if i>1;
      pos(1) = pos(1)+22;
      pos(2) = pos(2)-22;
      if pos(1)+1000>CurrINI.scrsz(3);
         pos(1) = 0;
      end;
      if pos(2)<=32;
         pos(2) = CurrINI.scrsz(4)-700;
      end;
   end;
   FigHandles{i} = figure('Visible','off','Position',pos);
   if NumFigs==1;
      set(FigHandles{i},'Name',['Marginal Posterior Densities - Chain ' num2str(CurrChain,'%0.0f')]);
   else;
      set(FigHandles{i},'Name',['Marginal Posterior Densities - Figure ' num2str(i,'%0.0f') ' Chain ' num2str(CurrChain,'%0.0f')]);
   end;
   SetFigureProperties(FigHandles{i},CurrINI);
   %
   k = 0;
   for r=1:NumRows;
      for c=1:NumCols;
         NumPlots = NumPlots+1;
         if NumPlots<=NumParam+NumAdditional;
            k = k+1;
            subplot(NumRows,NumCols,k);
            %
            % plot the prior density
            %
            if NumPlots<=NumParam;
               PriorHandle = plot(Xprior(NumPlots,:),Yprior(NumPlots,:),'--','Color',[0.75 0 0],'LineWidth',2);
            else;
               PriorHandle = plot(UpdatePriorX(NumPlots-NumParam,:),UpdatePriorY(NumPlots-NumParam,:),'--','Color',[0.75 0 0],'LineWidth',2);
            end;
            hold('on');
            %
            % plot the posterior density
            %
            PosteriorHandle = plot(Xposterior(NumPlots,:),Yposterior(NumPlots,:),'-','Color',[0 0 0.75],'LineWidth',2);
            hold('on');
            SetAxesFonts(CurrINI,'Title',StringTrim(ParameterNames(NumPlots,:)));
            axis('tight');
            xlim = [min(thetaPostSample(:,NumPlots)) max(thetaPostSample(:,NumPlots))];
            add_x = 0.1*(xlim(2)-xlim(1));
            x_lim = [xlim(1)-add_x xlim(2)+add_x];
            set(gca,'XLim',x_lim);
            %
            ylim = [0 1.2*max(Yposterior(NumPlots,:))];
            if NumPlots<=NumParam;
               if thetaDist(NumPlots)==1;
                  if x_lim(1)<0;
                     x_lim(1) = 0;
                  end;
               elseif thetaDist(NumPlots)==2;
                  if x_lim(1)<0;
                     x_lim(1) = 0;
                  end;
                  if x_lim(2)>1;
                     x_lim(2) = 1;
                  end;
               elseif thetaDist(NumPlots)==4;
                  if x_lim(1)<0;
                     x_lim(1) = 0;
                  end;
               elseif thetaDist(NumPlots)==5;
                  if x_lim(1)<BoundsMatrix(NumPlots,1);
                     x_lim(1) = BoundsMatrix(NumPlots,1);
                  end;
               elseif thetaDist(NumPlots)==6;
                  if x_lim(1)<BoundsMatrix(NumPlots,1);
                     x_lim(1) = BoundsMatrix(NumPlots,1);
                  end;
                  if x_lim(2)>BoundsMatrix(NumPlots,2);
                     x_lim(2) = BoundsMatrix(NumPlots,2);
                  end;
               end;
            else;
               if x_lim(1)<UpdatePriorBounds(NumPlots-NumParam,1);
                  x_lim(1) = UpdatePriorBounds(NumPlots-NumParam,1);
               end;
            end;
            %
            % plot the marginal mode estimate
            %
            MarginalModeHandle = plot(ModeData(NumPlots)*ones(1,2),ylim,'-','Color',[0 0 0],'LineWidth',1);
            hold('on');
            %
            % plot the joint posterior mode estimate
            %
            JointModeHandle = plot(thetaMode(NumPlots)*ones(1,2),ylim,':','Color',[0 0 0],'LineWidth',1);
            %
            set(gca,'XLim',x_lim,'YLim',ylim);
            hold('off');
         end;
      end;
   end;
   legHandle = legend([PriorHandle PosteriorHandle MarginalModeHandle JointModeHandle],strvcat('Prior','Posterior','Marginal Mode','Joint Mode'));
   %
   SetAxesFonts(CurrINI,legHandle);
   %
   % transform units temporarily into pixels
   %
   OldAxesUnits = get(gca,'Units');
   OldLegUnits = get(legHandle,'Units');
   set(gca,'Units','pixels');
   set(legHandle,'Units','pixels');
   %
   % check if we can move legend
   %
   LegPos = get(legHandle,'Position');
   AxPos = get(gca,'Position');
   if LegPos(1)<AxPos(1)+AxPos(3);
      if AxPos(1)+AxPos(3)+10+LegPos(3)<1000;
         LegPos(1) = AxPos(1)+AxPos(3)+10;
         set(legHandle,'Position',LegPos);
      end;
   end;
   set(gca,'Units',OldAxesUnits);
   set(legHandle,'Units',OldLegUnits);
   %
   % save it to disk before
   %
   if DirExist(GraphDirectory);
      if NumFigs>1;
         file = [GraphDirectory '\MarginalPosteriorDensity-' num2str(NumIter,'%0.0f') '-' int2str(CurrChain) '-' num2str(i,'%0.0f') '.fig'];
      else;
         file = [GraphDirectory '\MarginalPosteriorDensity-' num2str(NumIter,'%0.0f') '-' int2str(CurrChain) '.fig'];
      end;
      saveas(FigHandles{i},file);
   end;
   %
   % Fix the toolbar cdata
   %
   SetToolbarImages(FigHandles{i},CurrINI.images);
   %
   % Display the current figure
   %
   set(FigHandles{i},'Visible','on','CloseRequestFcn','delete(gcf); drawnow; pause(0.02);');
   drawnow;
   pause(0.02);
end;

%
% end of GraphPosteriorDensity.m
%
