function [fh,xh,gh,H,itct,fcount,func_change,abort,OptResults,phiEstimates,retcodeh] = YADAcsminwel(fcn,x0,H0,grad,crit,nit,waitbar,ShowClock,waitbartext,CurrINI,varargin)
%[fhat,xhat,ghat,Hhat,itct,fcount,retcodeh] = YADAcsminwel(fcn,x0,H0,grad,crit,nit,varargin)
% fcn:   string naming the objective function to be minimized
% x0:    initial value of the parameter vector
% H0:    initial value for the inverse Hessian.  Must be positive definite.
% grad:  Either a string naming a function that calculates the gradient, or the null matrix.
%        If it's null, the program calculates a numerical gradient.  In this case fcn must
%        be written so that it can take a matrix argument and produce a row vector of values.
% crit:  Convergence criterion.  Iteration will cease when it proves impossible to improve the
%        function value by more than crit.
% nit:   Maximum number of iterations.
% varargin: A list of optional length of additional parameters that get handed off to fcn each
%        time it is called.
%        Note that if the program ends abnormally, it is possible to retrieve the current x,
%        f, and H from the files g1.mat and H.mat that are written at each iteration and at each
%        hessian update, respectively.  (When the routine hits certain kinds of difficulty, it
%        write g2.mat and g3.mat as well.  If all were written at about the same time, any of them
%        may be a decent starting point.  One can also start from the one with best function value.)
%
% NOTE:  YADAcsminwel is a minor rewrite of Christopher Sims' csminwel function. The changes
%        only reflect some needs of YADA. The actual computations have not been changed. All disp
%        calls have been commented out. The outputs func_change and abort was added with the absolute value
%        of the function change. The inputs waitbar, waitbartext, ShowClock, and CurrINI were also added to
%        allow for a progress dialog while minimizing.
%
%        First YADA version: August 9, 2006.
%         This YADA version: May 28, 2008.
%
%

[nx,no]=size(x0);
nx=max(nx,no);
Verbose=0;
NumGrad= isempty(grad);
done=0;
itct=0;
fcount=0;
snit=100;
func_change = NaN;
abort = '';
OptResults = [];
phiEstimates = [];
%
f0 = feval(fcn,x0,varargin{:});
if f0>1e50;
   disp('Bad initial parameter.');
   return;
end;
if NumGrad
   if length(grad)==0
      [g badg] = YADAnumgrad(fcn,x0, varargin{:});
   else
      badg=any(find(grad==0));
      g=grad;
   end
else
   [g badg] = feval(grad,x0,varargin{:});
end
retcode3=101;
x=x0;
f=f0;
H=H0;
cliff=0;
%
% setup a progress dialog is asked to
%
if waitbar==1;
   ProgressStructure.title = ['Progress for ' num2str(nit,'%0.0f') ' iterations'];
   ProgressStructure.facecolor = CurrINI.progress_facecolor;
   ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
   ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
   ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
   ProgressStructure.stop = 0;
   ProgressStructure.clock = ShowClock;
   ProgressStructure.label = 'Convergence Value:';
   %
   if strcmp(StringTrim(waitbartext),'')==0;
      ProgressStructure.name = waitbartext;
   else;
      ProgressStructure.name = 'Numerical Optimization with csminwel';
   end;
   ProgressStructure.CurrINI = CurrINI;
   mywait = ProgressDLG(0,ProgressStructure);
   set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
else;
   %
   % show a wait dialog
   %
   txt = ['Please wait while YADA attempts to estimate the posterior mode with csminwel. A maximum of ' int2str(nit) ' iterations are used.'];
   WaitHandle = WaitDLG(txt,'information','Numerical Optimization with csminwel',500,200,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
end;
%
while ~done
   g1=[]; g2=[]; g3=[];
   %addition fj. 7/6/94 for control
%   disp('-----------------')
%   disp('-----------------')
   %disp('f and x at the beginning of new iteration')
%   disp(sprintf('f at the beginning of new iteration, %20.10f',f))
   %-----------Comment out this line if the x vector is long----------------
%      disp([sprintf('x = ') sprintf('%15.8g %15.8g %15.8g %15.8g\n',x)]);
   %-------------------------
   itct=itct+1;
   %
   % update the progress dialog
   %
   if waitbar==1;
      abort = get(mywait,'UserData');
      if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
         break;
      else;
         if itct<=nit;
            ProgressDLG([itct/nit func_change]);
         end;
      end;
   end;
   [f1 x1 fc retcode1] = YADAcsminit(fcn,x,f,g,badg,H,varargin{:});
   %
   fcount = fcount+fc;
   % erased on 8/4/94
   % if (retcode == 1) | (abs(f1-f) < crit)
   %    done=1;
   % end
   % if itct > nit
   %    done = 1;
   %    retcode = -retcode;
   % end
   if retcode1 ~= 1
      if retcode1==2 | retcode1==4
         wall1=1; badg1=1;
      else
         if NumGrad
            [g1 badg1] = YADAnumgrad(fcn, x1,varargin{:});
         else
            [g1 badg1] = feval(grad,x1,varargin{:});
         end
         wall1=badg1;
         % g1
%         save([pwd '\tmp\g1.mat'],'g1','x1','f1');
      end
      if wall1 % & (~done) by Jinill
         % Bad gradient or back and forth on step length.  Possibly at
         % cliff edge.  Try perturbing search direction.
         %
         %fcliff=fh;xcliff=xh;
         Hcliff=H+diag(diag(H).*rand(nx,1));
%         disp('Cliff.  Perturbing search direction.')
         [f2 x2 fc retcode2] = YADAcsminit(fcn,x,f,g,badg,Hcliff,varargin{:});
         %
         fcount = fcount+fc; % put by Jinill
         if  f2 < f
            if retcode2==2 | retcode2==4
                  wall2=1; badg2=1;
            else
               if NumGrad
                  [g2 badg2] = YADAnumgrad(fcn, x2,varargin{:});
               else
                  [g2 badg2] = feval(grad,x2,varargin{:});
               end
               wall2=badg2;
               % g2
%               badg2
%               save g2 g2 x2 f2 varargin
            end
            if wall2
%               disp('Cliff again.  Try traversing')
               if norm(x2-x1) < 1e-13
                  f3=f; x3=x; badg3=1;retcode3=101;
               else
                  gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1);
                  if(size(x0,2)>1), gcliff=gcliff', end
                  [f3 x3 fc retcode3] = YADAcsminit(fcn,x,f,gcliff,0,eye(nx),varargin{:});
                  %
                  fcount = fcount+fc; % put by Jinill
                  if retcode3==2 | retcode3==4
                     wall3=1; badg3=1;
                  else
                     if NumGrad
                        [g3 badg3] = YADAnumgrad(fcn, x3,varargin{:});
                     else
                        [g3 badg3] = feval(grad,x3,varargin{:});
                     end
                     wall3=badg3;
                     % g3
%                     badg3
%                     save g3 g3 x3 f3 varargin;
                  end
               end
            else
               f3=f; x3=x; badg3=1; retcode3=101;
            end
         else
            f3=f; x3=x; badg3=1;retcode3=101;
         end
      else
         % normal iteration, no walls, or else we're finished here.
         f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101;
      end
   else 
      f2=f;f3=f;f1=f;retcode2=retcode1;retcode3=retcode1;
   end
   %how to pick gh and xh
   if f3 < f - crit & badg3==0
      ih=3;
      fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
   elseif f2 < f - crit & badg2==0
      ih=2;
      fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
   elseif f1 < f - crit & badg1==0
      ih=1;
      fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
   else
      [fh,ih] = min([f1,f2,f3]);
%      disp(sprintf('ih = %d',ih))
      %
      switch ih
         case 1
            xh=x1;
         case 2
            xh=x2;
         case 3
            xh=x3;
      end
      %
      retcodei=[retcode1,retcode2,retcode3];
      retcodeh=retcodei(ih);
      if exist('gh')
         nogh=isempty(gh);
      else
         nogh=1;
      end
      if nogh
         if NumGrad
            [gh badgh] = YADAnumgrad(fcn,xh,varargin{:});
         else
            [gh badgh] = feval(grad, xh,varargin{:});
         end
      end
      badgh=1;
   end
   %end of picking
   %ih
   %fh
   %xh
   %gh
   %badgh
   stuck = (abs(fh-f) < crit);
   if (~badg)&(~badgh)&(~stuck)
      H = YADAbfgsi(H,gh-g,xh-x);
   end
%   if Verbose
%      disp('----')
%      disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh))
%   end
   % if Verbose
      func_change = abs(fh-f);
      if itct > nit
%         disp('iteration count termination');
         done = 1;
      elseif stuck
%         disp('improvement < crit termination');
         done = 1;
      end
%      rc=retcodeh;
%      if rc == 1
%         disp('zero gradient');
%      elseif rc == 6
%         disp('smallest step still improving too slow, reversed gradient');
%      elseif rc == 5
%         disp('largest step still improving too fast');
%      elseif (rc == 4) | (rc==2)
%         disp('back and forth on step length never finished');
%      elseif rc == 3
%         disp('smallest step still improving too slow');
%      elseif rc == 7
%         disp('warning: possible inaccuracy in H matrix');
%      end
   % end
   f=fh;
   x=xh;
   g=gh;
   badg=badgh;
   %
   % Collect results
   %
   OptResults = [OptResults;[itct -fh func_change max(abs(gh))]];
   phiEstimates = [phiEstimates xh];
   drawnow;
end

if waitbar==1;
   if ishandle(mywait)==1;
      set(mywait,'UserData','done');
      close(mywait);
      drawnow;
      pause(0.02);
   end;
else;
   if ishandle(WaitHandle)==1;
      set(WaitControls.text,'UserData','done');
      delete(WaitHandle);
      drawnow;
      pause(0.02);
   end;
end;

% 
% end of YADAcsminwel.m
%
