function [x, histout] = hooke(x0, f, budget, scales, tol, v) % % Hooke-Jeeves optimization % % C. T. Kelley, July 10, 1998 % % % This code comes with no guarantee or warranty of any kind. % % function [x, histout] = hooke(x0, f, budget, scales, tol, v) % % inputs: % % x0= initial iterate % f = objective function % budget = max f evals (default=50*number of variables) % The iteration will terminate after the iteration that % exhausts the budget % scales = the decreasing sequence of stencil sizes % This is an optional argument and the default is % 1, 1/2, ... 1/128 % tol = termination tolerance % difference is function values at successive iterations % (default = 1.d-6) % v = matrix of search directions (default $v = I$) % target = value of f at which the iteration should be terminated % This is an optional argument, which you SHOULD set to % something reasonable for your problem. The default is % as close to no limit as we can get, -1.d8 % % % outputs: % % final results = x % iteration histor = histout itout x 5 % histout = iteration history, updated after each nonlinear iteration % = lhist x 2 array, the rows are % [fcount, fval, sval] % fcount = cumulative function evals % fval = current best function value % sval = current scale % % set debug = 1 to print iteration stats % debug=0; % % n=length(x0); fcount=0; bud=50*n; scal=-(0:8)'; scal=2.^scal; dir=eye(n); tolf=1.d-6; targ=-1.d8; % % cache control paramters % global cache_size cache cache_ptr cache_fvals cache_size=4*n; cache=zeros(n,cache_size); cache_ptr=0; cache_fvals=zeros(cache_size,1); % % count the incoming arguments % if nargin > 2 bud=budget; end if nargin > 3 scal=scales; end if nargin > 4 tolf=tol; end if nargin > 5 targ=target; end if nargin > 6 targ=target; end if nargin > 7 dir = v; end % % main loop to sweep the scales % fcount=1; h=scal(1); fv=geval(f,x0); histout=[fcount,fv,h]; fbest=fv; x=x0; nscal=length(scal); ns=0; while (ns < nscal & fcount <= bud) ns=ns+1; [x, sf, lhist, fcount] = hjsearch(x, f, scal(ns),dir, fcount,bud); histout=[histout',lhist']'; % % end main loop % end % % Search with a single scale % function [x, sf, lhist, finc] = hjsearch(xb, f, h,dir, fcount,bud) x=xb; xc=x; sf=0; finc=fcount; lhist=[]; [x, fv, sf, numf] = hjexplore(xb, xc, f, h, dir); finc=finc+numf; if sf==1 thist=[finc, fv, h]; lhist=[lhist',thist']'; end while sf==1 & finc < bud % % pattern move % d=x-xb; xb=x; xc=x+d; fb=fv; [x, fv, sf, numf] = hjexplore(xb, xc, f, h, dir,fb); finc=finc+numf; if sf == 0 % pattern move fails! [x, fv, sf, numf] = hjexplore(xb, xb, f, h, dir, fb); finc=finc+numf; end if sf==1 thist=[finc, fv, h]; lhist=[lhist',thist']'; end end % % exploratory move % function [x, fv, sf, numf] = hjexplore(xb, xc, f, h, dir, fbold); global cache_size cache cache_ptr cache_fvals n=length(xb); x=xb; numf=0; if nargin == 5 [fb,ctr]=geval(f,x); numf=numf+ctr; else fb=fbold; end fv=fb; xt=xc; sf=0; dirh=h*dir; fvold=fv; for k=1:n p=xt+dirh(:,k); ft=feval(f,p); numf=numf+1; if(ft >= fb) p=xt-dirh(:,k); [ft,ctr]=geval(f,p); numf=numf+ctr; end if(ft < fb) sf=1; xt=p; fb=ft; end end if sf==1 x=xt; fv=fb; end % % function [fs,ctr]=geval(fh,xh) global cache_size cache cache_ptr cache_fvals for i=1:cache_size nz(i)=norm(xh-cache(:,i)); end [vz,iz]=min(nz); if vz == 0 & cache_ptr ~=0 fs=cache_fvals(iz); ctr=0; else fs=feval(fh,xh); ctr=1; cache_ptr=mod(cache_ptr,cache_size)+1; cache(:,cache_ptr)=xh; cache_fvals(cache_ptr)=fs; end
function [x, histout] = hooke(x0, f, budget, scales, tol, v)
ReplyDelete%
% Hooke-Jeeves optimization
%
% C. T. Kelley, July 10, 1998
%
%
% This code comes with no guarantee or warranty of any kind.
%
% function [x, histout] = hooke(x0, f, budget, scales, tol, v)
%
% inputs:
%
% x0= initial iterate
% f = objective function
% budget = max f evals (default=50*number of variables)
% The iteration will terminate after the iteration that
% exhausts the budget
% scales = the decreasing sequence of stencil sizes
% This is an optional argument and the default is
% 1, 1/2, ... 1/128
% tol = termination tolerance
% difference is function values at successive iterations
% (default = 1.d-6)
% v = matrix of search directions (default $v = I$)
% target = value of f at which the iteration should be terminated
% This is an optional argument, which you SHOULD set to
% something reasonable for your problem. The default is
% as close to no limit as we can get, -1.d8
%
%
% outputs:
%
% final results = x
% iteration histor = histout itout x 5
% histout = iteration history, updated after each nonlinear iteration
% = lhist x 2 array, the rows are
% [fcount, fval, sval]
% fcount = cumulative function evals
% fval = current best function value
% sval = current scale
%
% set debug = 1 to print iteration stats
%
debug=0;
%
%
n=length(x0); fcount=0; bud=50*n; scal=-(0:8)'; scal=2.^scal; dir=eye(n);
tolf=1.d-6; targ=-1.d8;
%
% cache control paramters
%
global cache_size cache cache_ptr cache_fvals
cache_size=4*n; cache=zeros(n,cache_size); cache_ptr=0;
cache_fvals=zeros(cache_size,1);
%
% count the incoming arguments
%
if nargin > 2 bud=budget; end
if nargin > 3 scal=scales; end
if nargin > 4 tolf=tol; end
if nargin > 5 targ=target; end
if nargin > 6 targ=target; end
if nargin > 7 dir = v; end
%
% main loop to sweep the scales
%
fcount=1; h=scal(1); fv=geval(f,x0);
histout=[fcount,fv,h];
fbest=fv;
x=x0; nscal=length(scal); ns=0;
while (ns < nscal & fcount <= bud)
ns=ns+1;
[x, sf, lhist, fcount] = hjsearch(x, f, scal(ns),dir, fcount,bud);
histout=[histout',lhist']';
%
% end main loop
%
end
%
% Search with a single scale
%
function [x, sf, lhist, finc] = hjsearch(xb, f, h,dir, fcount,bud)
x=xb; xc=x; sf=0; finc=fcount; lhist=[];
[x, fv, sf, numf] = hjexplore(xb, xc, f, h, dir);
finc=finc+numf;
if sf==1 thist=[finc, fv, h]; lhist=[lhist',thist']'; end
while sf==1 & finc < bud
%
% pattern move
%
d=x-xb; xb=x; xc=x+d; fb=fv;
[x, fv, sf, numf] = hjexplore(xb, xc, f, h, dir,fb);
finc=finc+numf;
if sf == 0 % pattern move fails!
[x, fv, sf, numf] = hjexplore(xb, xb, f, h, dir, fb);
finc=finc+numf;
end
if sf==1 thist=[finc, fv, h]; lhist=[lhist',thist']'; end
end
%
% exploratory move
%
function [x, fv, sf, numf] = hjexplore(xb, xc, f, h, dir, fbold);
global cache_size cache cache_ptr cache_fvals
n=length(xb); x=xb; numf=0;
if nargin == 5
[fb,ctr]=geval(f,x);
numf=numf+ctr;
else
fb=fbold;
end
fv=fb;
xt=xc; sf=0; dirh=h*dir;
fvold=fv;
for k=1:n
p=xt+dirh(:,k); ft=feval(f,p); numf=numf+1;
if(ft >= fb) p=xt-dirh(:,k); [ft,ctr]=geval(f,p); numf=numf+ctr; end
if(ft < fb) sf=1; xt=p; fb=ft; end
end
if sf==1 x=xt; fv=fb; end
%
%
function [fs,ctr]=geval(fh,xh)
global cache_size cache cache_ptr cache_fvals
for i=1:cache_size
nz(i)=norm(xh-cache(:,i));
end
[vz,iz]=min(nz);
if vz == 0 & cache_ptr ~=0
fs=cache_fvals(iz);
ctr=0;
else
fs=feval(fh,xh);
ctr=1;
cache_ptr=mod(cache_ptr,cache_size)+1;
cache(:,cache_ptr)=xh;
cache_fvals(cache_ptr)=fs;
end
need the input values..pls help
ReplyDelete