numgrad5

 Computes the gradient of the objective function fcn using a five points
 formula if possible.

 Adapted from Sims' numgrad.m routine.

 See section 25.3.6 Abramovitz and Stegun (1972, Tenth Printing, December) Handbook of Mathematical Functions.
 http://www.math.sfu.ca/~cbm/aands/ 

 TODO Try Four points formula when cost_flag3=0 or cost_flag4=0.

0001 function [g, badg, f0, f1, f2, f3, f4] = numgrad5(fcn,f0,x,epsilon,varargin)
0002 % Computes the gradient of the objective function fcn using a five points
0003 % formula if possible.
0004 %
0005 % Adapted from Sims' numgrad.m routine.
0006 %
0007 % See section 25.3.6 Abramovitz and Stegun (1972, Tenth Printing, December) Handbook of Mathematical Functions.
0008 % http://www.math.sfu.ca/~cbm/aands/
0009 %
0010 % TODO Try Four points formula when cost_flag3=0 or cost_flag4=0.
0011 
0012 % Original file downloaded from:
0013 % http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m
0014 
0015 % Copyright (C) 1993-2007 Christopher Sims
0016 % Copyright (C) 2008-2010 Dynare Team
0017 %
0018 % This file is part of Dynare.
0019 %
0020 % Dynare is free software: you can redistribute it and/or modify
0021 % it under the terms of the GNU General Public License as published by
0022 % the Free Software Foundation, either version 3 of the License, or
0023 % (at your option) any later version.
0024 %
0025 % Dynare is distributed in the hope that it will be useful,
0026 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0027 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0028 % GNU General Public License for more details.
0029 %
0030 % You should have received a copy of the GNU General Public License
0031 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0032 
0033 f1 = NaN;
0034 f2 = NaN;
0035 f3 = NaN;
0036 f4 = NaN;
0037 
0038 delta = epsilon;
0039 n=length(x);
0040 tvec=delta*eye(n);
0041 g=zeros(n,1);
0042 
0043 badg=0;
0044 goog=1;
0045 scale=1;
0046 
0047 for i=1:n
0048     if size(x,1)>size(x,2)
0049         tvecv=tvec(i,:);
0050     else
0051         tvecv=tvec(:,i);
0052     end
0053     [f1,junk1,junk2,cost_flag1] = feval(fcn, x+scale*transpose(tvecv), varargin{:});
0054     [f2,junk1,junk2,cost_flag2] = feval(fcn, x-scale*transpose(tvecv), varargin{:});
0055     if cost_flag1==0 || cost_flag2==0
0056         cost_flag3 = 0;
0057         cost_flag4 = 0;
0058         disp('numgrad:: I cannot use the five points formula!!')
0059     else
0060         [f3,junk1,junk2,cost_flag3] = feval(fcn, x+2*scale*transpose(tvecv), varargin{:});
0061         [f4,junk1,junk2,cost_flag4] = feval(fcn, x-2*scale*transpose(tvecv), varargin{:});
0062     end
0063     if cost_flag1 && cost_flag2 && cost_flag3 && cost_flag4% Five Points formula
0064         g0 = (8*(f1 - f2)+ f4-f3) / (12*scale*delta);
0065     elseif cost_flag3==0 || cost_flag4==0
0066         if cost_flag1 && cost_flag2% Three points formula
0067             g0 = (f1-f2)/(2*scale*delta);
0068         else
0069             if cost_flag1% Two points formula
0070                 g0 = (f1-f0) / (scale*delta);
0071             elseif cost_flag2% Two points formula
0072                 g0 = (f0-f2) / (scale*delta);
0073             else% Bad gradient!
0074                 goog=0;
0075             end
0076         end       
0077     end
0078     if goog && abs(g0)< 1e15 
0079         g(i)=g0;
0080     else
0081         disp('bad gradient ------------------------')
0082         g(i)=0;
0083         badg=1;
0084     end
0085 end