function FittedStd = VillaniCheckOptimum(ObjFuncStr,CenterPoint,MinMaxGrid,InvHessian,NumberOfEval,SameVerticalScale,SubplotStructure,EstParTeXNames,FigTitle,CurrINI,GraphDirectory,ParamName,varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PURPOSE:     Computing and plotting perturbations of the log posterior around a given point (CenterPoint). The pertubations
%              are one parameter at the time keeping the other parameters fixed at CenterPoint. 
%              
%              NaN is allowed in both CenterPoint and MinMaxGrid, but not at the same time for a given row (parameter). 
%              If CenterPoint(i)=nan, then CenterPoint(i) becomes the midpoint of the interval specified by MinMaxGrid(i,:)
%              If any of the bounds in MinMaxGrid(i,:) is nan, then the bounds are +- 4 posterior standard deviations, using the InvHessian 
%              to compute approximate posterior standard deviations.
%
% INPUT:       ObjFuncStr      (string)            String with name of .m file containing the objective function
%              CenterPoint     (1-by-Npara)        Point where the normal approximation of the log posterior is centered.
%              MinMaxGrid      (Npara-by-2)        MinMaxGrid(i,:) contains the bounds for the i:th parameter. May contain nan's or MinMaxGrid=[]. See PURPOSE. 
%              InvHessian      (Npara-by-Npara)    InvHessian of the log posterior.
%              NumberOfEval    (scalar)            Number of point in the grid.
%              SameVerticalScale (scalar)          If SameVerticalScale==1, then the same scale is used on the vertical axis for each parameter.
%              EstParTeXNames  (Npara-by-Something)EstParTeXNames(i,:) contains the i:ths parameter's name in TeX code.
%              SubplotStructure (1-by-2)           Determines the structure of the Subplots (the number of rows and cols)
%              varargin                            List of additional arguments needed to compute the objective function in ObjFuncStr
%      
% OUTPUT:      FittedStd       (Npara-by-1)        Estimated standard deviation (local curvature) from fit of quadratic around optimum.
%
% AUTHOR:      Mattias Villani, Research Department, Sveriges Riksbank and
%              Department of Statistics, Stockholm University. 
%              E-mail: mattias.villani@riksbank.se
%
% REVISED:     2005-11-03
%
% ADAPTED for YADA by Anders Warne, August 11-29, September 22, November 10, 2006, February 8, 2007,
%              January 10, 2008, February 18-19, 2008, February 25-26, 2008, October 23, 2008, and
%              November 28, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
NumberOfParameters = length(CenterPoint);
if isempty(SubplotStructure);
   if NumberOfParameters<=16;
      SubplotStructure=[ceil(sqrt(NumberOfParameters)) ceil(sqrt(NumberOfParameters))];
   else;
      SubplotStructure = [4 4];
   end;
end
NplotsPerFigure=prod(SubplotStructure);
NFigures=ceil(NumberOfParameters/NplotsPerFigure);
PriorFigHandles=cell(NFigures,1);
%
if NFigures==1;
   [Nrows,Ncols] = VillaniConstructOptimalSubplot(NumberOfParameters);
else;
   Nrows=SubplotStructure(1);
   Ncols=SubplotStructure(2);
end;
%
FittedStd = zeros(NumberOfParameters,1);
%
% Computing approximate conditional standard deviations
%
ApproxStd = zeros(NumberOfParameters,1);
S = InvHessian;
FullSet = (1:NumberOfParameters);
for i=1:NumberOfParameters;
   iExcluded = setdiff(FullSet,i);
   ApproxStd(i) = sqrt(S(i,i)-S(i,iExcluded)*inv(S(iExcluded,iExcluded))*S(iExcluded,i));
end;
%
if isempty(MinMaxGrid);
   LowerBoundGrid=CenterPoint-4*ApproxStd;
   UpperBoundGrid=CenterPoint+4*ApproxStd;
   StepLength=(8*ApproxStd)/(NumberOfEval-1);
else;
   for i=1:NumberOfParameters;
      if any(isnan(MinMaxGrid(i,:)));
         LowerBoundGrid(i) = CenterPoint(i)-4*ApproxStd;
         UpperBoundGrid(i) = CenterPoint(i)+4*ApproxStd;
         StepLength(i) = (8*ApproxStd)/(NumberOfEval-1);
         if isnan(CenterPoint(i));
            return;
         end;
      else;
         LowerBoundGrid(i) = MinMaxGrid(i,1);
         UpperBoundGrid(i) = MinMaxGrid(i,2);
         StepLength(i) = (UpperBoundGrid(i)-LowerBoundGrid(i))/(NumberOfEval-1);
         if isnan(CenterPoint(i));
            CenterPoint(i) = (LowerBoundGrid(i)+UpperBoundGrid(i))/2;
         end;
         if (CenterPoint(i)<LowerBoundGrid(i))|(CenterPoint(i)>UpperBoundGrid(i));
            CenterPoint(i) = (LowerBoundGrid(i)+UpperBoundGrid(i))/2;
         end;
      end;
   end;
end;
%
if strcmp(ObjFuncStr,'logPosteriorPhiDSGE')==1;
   NameStr = 'Transformed';
else;
   NameStr = 'Original';
end;
%
% create a wait dialog
%
txt = ['Please wait while YADA computes the log posterior for "' StringTrim(EstParTeXNames(1,:)) '" with grid values (' num2str(MinMaxGrid(1,1),'%0.4f') ':' num2str(StepLength(1),'%0.4f') ':' num2str(MinMaxGrid(1,2),'%0.4f') ').'];
WaitHandle = WaitDLG(txt,'information',['Check Posterior Mode - ' NameStr ' Parameters'],500,150,CurrINI,0);
WaitControls = get(WaitHandle,'UserData');
%
Minimum = inf;
Maximum = -inf;
scrsz = get(0,'ScreenSize');
pos = [(scrsz(3)-1000)/2 (scrsz(4)-680)/2 1000 650];
if isempty(InvHessian)==0;
   LineHandles = zeros(NFigures,3);
else;
   LineHandles = zeros(NFigures,2);
end;
for j=1:NFigures;
   if j>1;
      pos(1) = pos(1)+22;
      pos(2) = pos(2)-22;
      if pos(1)+1000>scrsz(3);
         pos(1) = 0;
      end;
      if pos(2)<=32;
         pos(2) = scrsz(4)-700;
      end;
   end;
   PriorFigHandles{j} = figure('Visible','off','Position',pos,'CloseRequestFcn','delete(gcf); drawnow; pause(0.02);');
   SetFigureProperties(PriorFigHandles{j},CurrINI);
   k = 0;
   if j<NFigures;
      PlotCounter=1+(j-1)*NplotsPerFigure:j*NplotsPerFigure;
   else;
      PlotCounter=NplotsPerFigure*(NFigures-1)+1:NumberOfParameters;
   end;
   % 
   for i=PlotCounter;
      k = k+1;
      %
      % update the wait dialog
      %
      txt = ['Please wait while YADA computes the log posterior for "' StringTrim(EstParTeXNames(i,:)) '" with grid values (' num2str(MinMaxGrid(i,1),'%0.4f') ':' num2str(StepLength(i),'%0.4f') ':' num2str(MinMaxGrid(i,2),'%0.4f') ').'];
      set(WaitControls.text,'String',txt);
      drawnow;
      pause(0.02);
      %
      % make sure that PriorFigHandles{j} is the current figure
      %
      set(0,'CurrentFigure',PriorFigHandles{j});
      subplot(Nrows,Ncols,k);
      %
      % Computes the log posterior around the supposed optimum
      %
      GridAroundOpt=LowerBoundGrid(i):StepLength(i):UpperBoundGrid(i);
      if isempty(MinMaxGrid);
         if rem(NumberOfEval,2)==0;
            GridAroundOpt = [GridAroundOpt(1:NumberOfEval/2) CenterPoint(i) GridAroundOpt(NumberOfEval/2+1:NumberOfEval)];
            CenterPointLoc = NumberOfEval/2+1;
         else;
            GridAroundOpt = [GridAroundOpt(1:floor(NumberOfEval/2)) CenterPoint(i) GridAroundOpt(floor(NumberOfEval/2)+1:NumberOfEval)];
            CenterPointLoc = floor(NumberOfEval/2)+1;
         end;
      else;
         NLowerThanCenter = sum(GridAroundOpt<=CenterPoint(i));
         NHigherThanCenter = sum(GridAroundOpt>=CenterPoint(i));
         if NLowerThanCenter+NHigherThanCenter>length(GridAroundOpt); % Center is one of the grid points
            GridAroundOpt(GridAroundOpt~=CenterPoint(i)); % Exclude centerpoint
         end;
         GridAroundOpt = [GridAroundOpt(1:NLowerThanCenter) CenterPoint(i) GridAroundOpt(NLowerThanCenter+1:length(GridAroundOpt))];
         CenterPointLoc = NLowerThanCenter+1;
      end;
      LogPostOnGrid = zeros(1,length(GridAroundOpt));
      Solved = ones(1,length(GridAroundOpt));
      TempPoint = CenterPoint;
      u = 0;
      LogPostOnGridMode = -feval(ObjFuncStr,TempPoint,varargin{:});
      for q=GridAroundOpt;
         u = u+1;
         TempPoint(i)=q;
         LogPostDens = -feval(ObjFuncStr,TempPoint,varargin{:});
         LogPostOnGrid(u) = LogPostDens;
         if LogPostDens==-1000000;
            Solved(i) = 0;
         end;
      end;
      %
      LogPostOnGrid = LogPostOnGrid(Solved==1); % Discarding non-solvable grid points
      GridAroundOpt = GridAroundOpt(Solved==1);
      CloseToMode = (LogPostOnGridMode-LogPostOnGrid<10); % Plotting only grid point with log posterior
      LogPostOnGrid = LogPostOnGrid(CloseToMode==1);      % not more than 10 units smaller than Log posterior mode 
      GridAroundOpt = GridAroundOpt(CloseToMode==1);
      Minimum = min(Minimum,min(LogPostOnGrid));
      Maximum = max(Maximum,max(LogPostOnGrid));
      %
      % make sure the current figure is properly set
      %
      set(0,'CurrentFigure',PriorFigHandles{j});
      p1 = plot(GridAroundOpt,LogPostOnGrid,'.-','LineWidth',2,'Color',[0 0 0.75],'Tag','logposterior');
%      ApproxLogPDFMode = logNormalPDF(CenterPoint(i),CenterPoint(i),sqrt(InvHessian(i,i)));
      ApproxLogPDFMode = logNormalPDF(CenterPoint(i),CenterPoint(i),ApproxStd(i));
      hold('on');
      if isempty(InvHessian)==0;
%         ApproxLogPDF = logNormalPDF(GridAroundOpt,CenterPoint(i),sqrt(InvHessian(i,i)));
         ApproxLogPDF = logNormalPDF(GridAroundOpt,CenterPoint(i),ApproxStd(i));
         p2 = plot(GridAroundOpt,ApproxLogPDF-ApproxLogPDFMode+LogPostOnGridMode,':','LineWidth',2,'Color',[0.75 0 0],'Tag','hessian');
         hold('on');
      end;
      %LaplaceStuff=ApproxLogPDF-ApproxLogPDFMode+LogPostOnGridMode;
      %
      % Modified InvHessian by fitting quadratic to evaluated LogPost
      %
      y = LogPostOnGrid(:);
      x = GridAroundOpt(:);
      Z = [ones(length(x),1) x x.^2];
      Bols = inv(Z'*Z)*Z'*y;
      FittedStd(i) = 1/(sqrt(abs(2*Bols(3))));
      ApproxLogPDFMode = logNormalPDF(CenterPoint(i),CenterPoint(i),FittedStd(i));
      ApproxLogPDF = logNormalPDF(GridAroundOpt,CenterPoint(i),FittedStd(i));
      p3 = plot(GridAroundOpt,ApproxLogPDF-ApproxLogPDFMode+LogPostOnGridMode,'--','LineWidth',2,'Color',[0 0.75 0],'Tag','modhessian');
      hold('on');
   end;
   %
   % collect the last set of handles
   %
   if isempty(InvHessian)==0;
      LineHandles(j,:) = [p1 p2 p3];
   else;
      LineHandles(j,:) = [p1 p3];
   end;
end
%
% close the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);
%
ParamCount=1;
for j=1:NFigures;
   figure(PriorFigHandles{j});
   if NFigures>1;
      set(PriorFigHandles{j},'name',[FigTitle ' No. ' int2str(j)]);
   else;
      set(PriorFigHandles{j},'name',FigTitle);
   end;
   k = 0;
   for Row=1:Nrows;
      for Col=1:Ncols;
         if ParamCount<=NumberOfParameters;
            k = k+1;
            subplot(Nrows,Ncols,k);
            if SameVerticalScale;
               set(gca,'ylim',[Minimum Maximum]);
            else;
               YLimits = get(gca,'ylim');
               Minimum = YLimits(1);
               Maximum = YLimits(2);
            end;
            %
            % get the data for the log posterior and the modified hessian
            %
            LogPostHandle = findobj(gca,'Type','line','Tag','logposterior');
            ModHessianHandle = findobj(gca,'Type','line','Tag','modhessian');
            if (ishandle(LogPostHandle)==1)&(ModHessianHandle);
               YLogPostData = get(LogPostHandle,'YData');
               YModHessianData = get(ModHessianHandle,'YData');
               YLim = [min([YLogPostData YModHessianData]) max([YLogPostData YModHessianData])];
               Maximum = LogPostOnGridMode;
               AddYLim = 0.1*(YLim(2)-YLim(1));
               YLim = YLim+[0 AddYLim];
            else;
               YLim = [];
            end;
            %
            LineHandle = line([CenterPoint(ParamCount) CenterPoint(ParamCount)],[Minimum Maximum]);
            set(LineHandle,'color',[0 0 0],'LineWidth',1);
            SetAxesFonts(CurrINI,'Title',StringTrim(EstParTeXNames(ParamCount,:)));
            axis('tight');
            %
            % fix the y-axes
            %
            if isempty(YLim)==0;
               set(gca,'YLim',YLim);
            end;
         end;
         ParamCount = ParamCount+1;
      end;
   end;
   %
   % add a legend
   %
   if size(LineHandles,2)==3;
      legHandle = legend([LineHandle LineHandles(j,:)],strvcat('Mode','Log Posterior','Norm Approx - Hessian','Norm Approx - Modified Hessian'));
   else;
      legHandle = legend([LineHandle LineHandles(j,:)],strvcat('Mode','Log Posterior','Norm Approx - Modified Hessian'));
   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 NFigures>1;
         file = [GraphDirectory '\PosterorModeOptimum-' ParamName '-' num2str(j,'%0.0f') '.fig'];
      else;
         file = [GraphDirectory '\PosterorModeOptimum-' ParamName '.fig'];
      end;
      SetLogText(strvcat('   Saving plot to:',['   ' GraphDirectory '...']));
      saveas(PriorFigHandles{j},file);
   end;
   SetToolbarImages(PriorFigHandles{j},CurrINI.images);
   set(PriorFigHandles{j},'Visible','on');
   drawnow;
   pause(0.02);
end;

%
% end of VillaniCheckOptimum.m
%
