csminwel1

[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,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.
 method: integer scalar, 2, 3 or 5 points formula.
 epsilon: scalar double, numerical differentiation increment
 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.)

0001 function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin)
0002 %[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin)
0003 % fcn:   string naming the objective function to be minimized
0004 % x0:    initial value of the parameter vector
0005 % H0:    initial value for the inverse Hessian.  Must be positive definite.
0006 % grad:  Either a string naming a function that calculates the gradient, or the null matrix.
0007 %        If it's null, the program calculates a numerical gradient.  In this case fcn must
0008 %        be written so that it can take a matrix argument and produce a row vector of values.
0009 % crit:  Convergence criterion.  Iteration will cease when it proves impossible to improve the
0010 %        function value by more than crit.
0011 % nit:   Maximum number of iterations.
0012 % method: integer scalar, 2, 3 or 5 points formula.
0013 % epsilon: scalar double, numerical differentiation increment
0014 % varargin: A list of optional length of additional parameters that get handed off to fcn each
0015 %        time it is called.
0016 %        Note that if the program ends abnormally, it is possible to retrieve the current x,
0017 %        f, and H from the files g1.mat and H.mat that are written at each iteration and at each
0018 %        hessian update, respectively.  (When the routine hits certain kinds of difficulty, it
0019 %        write g2.mat and g3.mat as well.  If all were written at about the same time, any of them
0020 %        may be a decent starting point.  One can also start from the one with best function value.)
0021 
0022 % Original file downloaded from:
0023 % http://sims.princeton.edu/yftp/optimize/mfiles/csminwel.m
0024 
0025 % Copyright (C) 1993-2007 Christopher Sims
0026 % Copyright (C) 2006-2011 Dynare Team
0027 %
0028 % This file is part of Dynare.
0029 %
0030 % Dynare is free software: you can redistribute it and/or modify
0031 % it under the terms of the GNU General Public License as published by
0032 % the Free Software Foundation, either version 3 of the License, or
0033 % (at your option) any later version.
0034 %
0035 % Dynare is distributed in the hope that it will be useful,
0036 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0038 % GNU General Public License for more details.
0039 %
0040 % You should have received a copy of the GNU General Public License
0041 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0042 
0043 global bayestopt_
0044 
0045 fh = [];
0046 xh = [];
0047 [nx,no]=size(x0);
0048 nx=max(nx,no);
0049 Verbose=1;
0050 NumGrad= isempty(grad);
0051 done=0;
0052 itct=0;
0053 fcount=0;
0054 snit=100;
0055 %tailstr = ')';
0056 %stailstr = [];
0057 % Lines below make the number of Pi's optional.  This is inefficient, though, and precludes
0058 % use of the matlab compiler.  Without them, we use feval and the number of Pi's must be
0059 % changed with the editor for each application.  Places where this is required are marked
0060 % with ARGLIST comments
0061 %for i=nargin-6:-1:1
0062 %   tailstr=[ ',P' num2str(i)  tailstr];
0063 %   stailstr=[' P' num2str(i) stailstr];
0064 %end
0065 
0066 [f0,junk1,junk2,cost_flag] = feval(fcn,x0,varargin{:});
0067 
0068 if ~cost_flag
0069     disp('Bad initial parameter.') 
0070     return 
0071 end
0072 
0073 if NumGrad
0074     switch method
0075       case 2
0076         [g,badg] = numgrad2(fcn, f0, x0, epsilon, varargin{:});
0077       case 3
0078         [g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:});
0079       case 5
0080         [g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:});    
0081     end
0082 elseif ischar(grad)
0083     [g,badg] = feval(grad,x0,varargin{:});
0084 else
0085     g=junk1;
0086     badg=0;
0087 end
0088 retcode3=101;
0089 x=x0;
0090 f=f0;
0091 H=H0;
0092 cliff=0;
0093 while ~done
0094     bayestopt_.penalty = f;
0095     g1=[]; g2=[]; g3=[];
0096     %addition fj. 7/6/94 for control
0097     disp('-----------------')
0098     disp('-----------------')
0099     %disp('f and x at the beginning of new iteration')
0100     disp(sprintf('f at the beginning of new iteration, %20.10f',f))
0101     %-----------Comment out this line if the x vector is long----------------
0102     %   disp([sprintf('x = ') sprintf('%15.8g %15.8g %15.8g %15.8g\n',x)]);
0103     %-------------------------
0104     itct=itct+1;
0105     [f1 x1 fc retcode1] = csminit1(fcn,x,f,g,badg,H,varargin{:});
0106     %ARGLIST
0107     %[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,...
0108     %           P8,P9,P10,P11,P12,P13);
0109     % itct=itct+1;
0110     fcount = fcount+fc;
0111     % erased on 8/4/94
0112     % if (retcode == 1) || (abs(f1-f) < crit)
0113     %    done=1;
0114     % end
0115     % if itct > nit
0116     %    done = 1;
0117     %    retcode = -retcode;
0118     % end
0119     if retcode1 ~= 1
0120         if retcode1==2 || retcode1==4
0121             wall1=1; badg1=1;
0122         else
0123             if NumGrad
0124                 switch method 
0125                   case 2
0126                     [g1 badg1] = numgrad2(fcn, f1, x1, epsilon, varargin{:});
0127                   case 3
0128                     [g1 badg1] = numgrad3(fcn, f1, x1, epsilon, varargin{:});
0129                   case 5
0130                     [g1,badg1] = numgrad5(fcn, f1, x1, epsilon, varargin{:});             
0131                 end
0132             elseif ischar(grad),
0133                 [g1 badg1] = feval(grad,x1,varargin{:});
0134             else
0135                 [junk1,g1,junk2, cost_flag] = feval(fcn,x1,varargin{:});
0136                 badg1 = ~cost_flag;
0137             end
0138             wall1=badg1;
0139             % g1
0140             save g1.mat g1 x1 f1 varargin;
0141             %ARGLIST
0142             %save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
0143         end
0144         if wall1 % && (~done) by Jinill
0145                  % Bad gradient or back and forth on step length.  Possibly at
0146                  % cliff edge.  Try perturbing search direction.
0147                  %
0148                  %fcliff=fh;xcliff=xh;
0149             Hcliff=H+diag(diag(H).*rand(nx,1));
0150             disp('Cliff.  Perturbing search direction.')
0151             [f2 x2 fc retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,varargin{:});
0152             %ARGLIST
0153             %[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,...
0154             %     P5,P6,P7,P8,P9,P10,P11,P12,P13);
0155             fcount = fcount+fc; % put by Jinill
0156             if  f2 < f
0157                 if retcode2==2 || retcode2==4
0158                     wall2=1; badg2=1;
0159                 else
0160                     if NumGrad
0161                         switch method
0162                           case 2
0163                             [g2 badg2] = numgrad2(fcn, f2, x2, epsilon, varargin{:});
0164                           case 3
0165                             [g2 badg2] = numgrad3(fcn, f2, x2, epsilon, varargin{:});
0166                           case 5
0167                             [g2,badg2] = numgrad5(fcn, f2, x2, epsilon, varargin{:});                   
0168                         end
0169                     elseif ischar(grad),
0170                         [g2 badg2] = feval(grad,x2,varargin{:});
0171                     else
0172                         [junk1,g2,junk2, cost_flag] = feval(fcn,x1,varargin{:});
0173                         badg2 = ~cost_flag;
0174                     end
0175                     wall2=badg2;
0176                     % g2
0177                     badg2
0178                     save g2.mat g2 x2 f2 varargin
0179                     %ARGLIST
0180                     %save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
0181                 end
0182                 if wall2
0183                     disp('Cliff again.  Try traversing')
0184                     if norm(x2-x1) < 1e-13
0185                         f3=f; x3=x; badg3=1;retcode3=101;
0186                     else
0187                         gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1);
0188                         if(size(x0,2)>1), gcliff=gcliff', end
0189                         [f3 x3 fc retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),varargin{:});
0190                         %ARGLIST
0191                         %[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,...
0192                         %         P4,P5,P6,P7,P8,...
0193                         %      P9,P10,P11,P12,P13);
0194                         fcount = fcount+fc; % put by Jinill
0195                         if retcode3==2 || retcode3==4
0196                             wall3=1; badg3=1;
0197                         else
0198                             if NumGrad
0199                                 switch method
0200                                   case 2
0201                                     [g3 badg3] = numgrad2(fcn, f3, x3, epsilon, varargin{:});
0202                                   case 3
0203                                     [g3 badg3] = numgrad3(fcn, f3, x3, epsilon, varargin{:});
0204                                   case 5
0205                                     [g3,badg3] = numgrad5(fcn, f3, x3, epsilon, varargin{:});                         
0206                                 end
0207                             elseif ischar(grad),
0208                                 [g3 badg3] = feval(grad,x3,varargin{:});
0209                             else
0210                                 [junk1,g3,junk2, cost_flag] = feval(fcn,x1,varargin{:});
0211                                 badg3 = ~cost_flag;
0212                             end
0213                             wall3=badg3;
0214                             % g3
0215                             save g3.mat g3 x3 f3 varargin;
0216                             %ARGLIST
0217                             %save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
0218                         end
0219                     end
0220                 else
0221                     f3=f; x3=x; badg3=1; retcode3=101;
0222                 end
0223             else
0224                 f3=f; x3=x; badg3=1;retcode3=101;
0225             end
0226         else
0227             % normal iteration, no walls, or else we're finished here.
0228             f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101;
0229         end
0230     else 
0231         f2=f;f3=f;f1=f;retcode2=retcode1;retcode3=retcode1;
0232     end
0233     %how to pick gh and xh
0234     if f3 < f - crit && badg3==0 && f3 < f2 && f3 < f1
0235         ih=3;
0236         fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
0237     elseif f2 < f - crit && badg2==0 && f2 < f1
0238         ih=2;
0239         fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
0240     elseif f1 < f - crit && badg1==0
0241         ih=1;
0242         fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
0243     else
0244         [fh,ih] = min([f1,f2,f3]);
0245         %disp(sprintf('ih = %d',ih))
0246         %eval(['xh=x' num2str(ih) ';'])
0247         switch ih
0248           case 1
0249             xh=x1;
0250           case 2
0251             xh=x2;
0252           case 3
0253             xh=x3;
0254         end %case
0255             %eval(['gh=g' num2str(ih) ';'])
0256             %eval(['retcodeh=retcode' num2str(ih) ';'])
0257         retcodei=[retcode1,retcode2,retcode3];
0258         retcodeh=retcodei(ih);
0259         if exist('gh')
0260             nogh=isempty(gh);
0261         else
0262             nogh=1;
0263         end
0264         if nogh
0265             if NumGrad
0266                 switch method
0267                   case 2
0268                     [gh,badgh] = numgrad2(fcn, fh, xh, epsilon, varargin{:});
0269                   case 3
0270                     [gh,badgh] = numgrad3(fcn, fh, xh, epsilon, varargin{:});
0271                   case 5
0272                     [gh,badgh] = numgrad5(fcn, fh, xh, epsilon, varargin{:});
0273                 end
0274             elseif ischar(grad),
0275                 [gh badgh] = feval(grad, xh,varargin{:});
0276             else
0277                 [junk1,gh,junk2, cost_flag] = feval(fcn,x1,varargin{:});
0278                 badgh = ~cost_flag;
0279             end
0280         end
0281         badgh=1;
0282     end
0283     %end of picking
0284     %ih
0285     %fh
0286     %xh
0287     %gh
0288     %badgh
0289     stuck = (abs(fh-f) < crit);
0290     if (~badg) && (~badgh) && (~stuck)
0291         H = bfgsi(H,gh-g,xh-x);
0292     end
0293     if Verbose
0294         disp('----')
0295         disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh))
0296     end
0297     % if Verbose
0298     if itct > nit
0299         disp('iteration count termination')
0300         done = 1;
0301     elseif stuck
0302         disp('improvement < crit termination')
0303         done = 1;
0304     end
0305     rc=retcodeh;
0306     if rc == 1
0307         disp('zero gradient')
0308     elseif rc == 6
0309         disp('smallest step still improving too slow, reversed gradient')
0310     elseif rc == 5
0311         disp('largest step still improving too fast')
0312     elseif (rc == 4) || (rc==2)
0313         disp('back and forth on step length never finished')
0314     elseif rc == 3
0315         disp('smallest step still improving too slow')
0316     elseif rc == 7
0317         disp('warning: possible inaccuracy in H matrix')
0318     end
0319     % end
0320     f=fh;
0321     x=xh;
0322     g=gh;
0323     badg=badgh;
0324 end
0325 % what about making an m-file of 10 lines including numgrad.m
0326 % since it appears three times in csminwel.m