function CUSUMPlots(CUSUMPost,CUSUMAverage,NumBurnin,CurrChain,PlotType,DSGEModel,CurrINI)
% CUSUMPlots: Presents CUSUM plots of certain statistics.
%
% USAGE:
%
%       CUSUMPlots(CUSUMPost,CUSUMAverage,NumBurnin,CurrChain,PlotType,DSGEModel,CurrINI)
%
% REQUIRED INPUT: CUSUMPost (matrix) of dimension NumIter-NumBurnin x NumParam. The CUSUMPost entries
%                           are calculated relative to the full sample mean, which is then added to
%                           the accumulated difference.
%
%                 CUSUMAverage (matrix) of dimension 0.9*(NumIter-NumBurnin)xNumParam.
%
%                 NumBurnin (integer) with the number of burn-in draws.
%
%                 CurrChain (integer) with the number of the current MCMC chain.
%
%                 PlotType (string vector) that takes on the values "logpost", "original", and "transformed".
%
%                 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.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: April 18, 2007.
%                        This version: December 19, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:
%
% * 9-11-2007:  Updated the documentation.
%
% * 10-1-2008:  Changed the last drawnow call to pause(0.02).
%
% * 10-3-2008:  Made sure that the "CurrChain" variable affects the filename for the figure.
%
% * 11-3-2008:  Took the changed dimension of CUSUMAverage into account.
%
% * 27-3-2008:  Made sure that the sequential grid has at least 2 values.
%
% * 3-4-2008:   Fixed an error in the size shown in the legend for the moving CUSUM window.
%
% * 18-4-2008:  Fixed a bug where the wrong x-axis values where used when plotting "CUSUMAverage".
%
% * 23-5-2008:  Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% first we sample from the draws
%
[NumIter NumParam] = size(CUSUMPost);
NumSum = size(CUSUMAverage,1);
%
% locate the controls structure
%
maingui = findobj('Type','figure','Tag','YADA');
controls = get(maingui,'UserData');
%
% sequential marginal likelihood start value and step length
%
SeqStartStr = get(controls.dsge.sequentialstart,'String');
Origin = str2num(strrep(StringTrim(SeqStartStr(DSGEModel.SequentialStartIteration,:)),',',''));
SeqStepStr = get(controls.dsge.sequentialstep,'String');
StepLength = str2num(strrep(StringTrim(SeqStepStr(DSGEModel.SequentialStepLength,:)),',',''));
%
% make sure that Origin is less than NumIter
%
if Origin==NumIter;
   Origin = floor(NumIter/10);
end;
%
% Create the sequential grid over which we will display the cusum plots
%
SeqGrid = (Origin:StepLength:NumIter);
if isempty(SeqGrid)==1;
   SeqGrid = NumIter;
elseif SeqGrid(length(SeqGrid))<NumIter;
   SeqGrid = [SeqGrid NumIter];
end;
SeqGrid2 = SeqGrid(SeqGrid<=NumSum);
%
% now we sample
%
CUSUMPost = CUSUMPost(SeqGrid,:);
CUSUMAverage = CUSUMAverage(SeqGrid2+1,:);
NumSum = 0.1*NumIter;
%
GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
%
if strcmp(lower(PlotType),'logpost')==0;
   %
   % 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;
   %
   ParameterNames = ModeData.ParameterNames;
   %
   % Compute the number of figures (max 16 plots per graph)
   %
   NumFigs = ceil(NumParam/16);
   %
   if NumFigs==1;
      [NumRows,NumCols] = VillaniConstructOptimalSubplot(NumParam);
   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',['CUSUM Plots ' PlotType ' Parameters - Chain ' num2str(CurrChain,'%0.0f')]);
      else;
         set(FigHandles{i},'Name',['CUSUM Plots ' PlotType ' Parameters - 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;
               k = k+1;
               subplot(NumRows,NumCols,k);
               %
               % first we plot the average
               %
               if isempty(CUSUMAverage)==0;
                  AverageHandle = plot(SeqGrid2+NumSum+NumBurnin,CUSUMAverage(:,NumPlots)','LineStyle','--','Color',[0.75 0 0],'LineWidth',2);
                  hold('on');
               end;
               %
               % plot the CUSUM plot
               %
               PostHandle = plot(SeqGrid+NumBurnin,CUSUMPost(:,NumPlots)','LineStyle','-','Color',[0 0 0.75],'LineWidth',2);
               hold('on');
               %
               SetAxesFonts(CurrINI,'Title',StringTrim(ParameterNames.estimated(NumPlots,:)));
               %
               % check the YLim
               %
               ylim = get(gca,'YLim');
               add_on = 0.5*(ylim(2)-ylim(1));
               ylim = [ylim(1)-add_on ylim(2)+add_on];
               set(gca,'XLim',[SeqGrid(1)+NumBurnin NumIter+NumBurnin],'YLim',ylim);
               %
               % check if we should add a zero line
               %
               if (ylim(1)<0)&(ylim(2)>0);
                  hold('on');
                  plot([NumBurnin NumIter+NumBurnin],[0 0],'-','Color',[0 0 0]);
               end;
               hold('off');
            end;
         end;
      end;
      %
      % add a legend
      %
      if isempty(CUSUMAverage)==0;
         legHandle = legend([PostHandle AverageHandle],strvcat('CUSUM',['Moving CUSUM window (size=' num2str(NumSum,'%0.0f') ')']));
      else;
         legHandle = legend(PostHandle,'CUSUM');
      end;
      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)==1;
         if NumFigs>1;
            file = [GraphDirectory '\CUSUM-' PlotType '-' num2str(NumIter,'%0.0f') '-' num2str(i,'%0.0f') '.fig'];
         else;
            file = [GraphDirectory '\CUSUM-' PlotType '-' num2str(NumIter,'%0.0f') '.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;
else;
   %
   % prepare a CUSUM plot for the log posterior values
   %
   FigHandle =figure('Visible','off', ...
          'Position',[(CurrINI.scrsz(3)-700)/2 (CurrINI.scrsz(4)-500)/2 700 500], ...
          'Tag','MargLikeCJ', ...
          'Name','CUSUM Plots - Log Posterior');
   %
   SetFigureProperties(FigHandle,CurrINI);
   %
   % first we plot the average
   %
   if isempty(CUSUMAverage)==0;
      AverageHandle = plot(SeqGrid2+NumSum,CUSUMAverage,'LineStyle','--','Color',[0.75 0 0],'LineWidth',2);
      hold('on');
   end;
   %
   % and now we plot the CUSUM plot
   %
   PostHandle = plot(SeqGrid+NumBurnin,CUSUMPost','LineStyle','-','Color',[0 0 0.75],'LineWidth',2);
   hold('on');
   %
   % check the YLim
   %
   ylim = get(gca,'YLim');
   add_on = 0.5*(ylim(2)-ylim(1));
   ylim = [ylim(1)-add_on ylim(2)+add_on];
   set(gca,'XLim',[SeqGrid(1)+NumBurnin NumIter+NumBurnin],'YLim',ylim);
   %
   % check if we should draw a horizontal zero line
   %
   if (ylim(1)<0)&(ylim(2));
      plot([SeqGrid(1)+NumBurnin NumIter+NumBurnin],[0 0],'-','Color',[0 0 0]);
      kids = get(gca,'Children');
      if isempty(CUSUMAverage)==0;
         set(gca,'Children',[PostHandle AverageHandle kids(1)]);
      else;
         set(gca,'Children',[PostHandle kids(1)]);
      end;
   end;
   %
   SetAxesFonts(CurrINI,'Title','CUSUM plot of log posterior','XLabel',['RWM draw number (with ' num2str(NumBurnin,'%0.0f') ' burn-in draws)'],'YLabel','Average log posterior');
   %
   % set the legend
   %
   if isempty(CUSUMAverage)==0;
      legHandle = legend([PostHandle AverageHandle],strvcat('CUSUM',['Moving CUSUM window (size=' num2str(NumSum-NumBurnin,'%0.0f') ')']),0);
   else;
      legHandle = legend(PostHandle,'CUSUM',0);
   end;
   SetAxesFonts(CurrINI,legHandle);
   hold('off');
   %
   % save the grasph
   %
   if DirExist(GraphDirectory);
      file = [GraphDirectory '\CUSUMLogPosterior-' num2str(NumIter,'%0.0f') '-' int2str(CurrChain) '.fig'];
      saveas(FigHandle,file);
   end;
   %
   % display the graph
   %
   SetToolbarImages(FigHandle,CurrINI.images);
   set(FigHandle,'Visible','on','CloseRequestFcn','delete(gcf); drawnow; pause(0.02);');
   drawnow;
   pause(0.02);
end;

%
% end of CUSUMPlots.m
%
