Code covered by the BSD License

# OpenOpt

### Dmitrey Kroshko (view profile)

25 Nov 2006 (Updated )

nonSmoothSolve (similar to fsolve), non-smooth & noisy local + some global solvers; works in Octave

buscarnd(prob)
```function r = buscarnd(prob)
% line 186 "if norm(mem_x(:,k)-x(:,j)) < metric" EATS VERY MUCH CPU!!!//Dmitrey
% function [xo,Ot,nS]=buscarnd(S,x0,ip,nOt,samples,Lb,Ub,problem,tol,mxit,R,red,mem)
%   Unconstrained global optimization using adaptive random search.
%
%   [xo,Ot,nS]=buscarnd(S,x0,ip,nOt,samples,Lb,Ub,problem,tol,mxit,R,red,mem)
%
%   S: objective function
%   x0: initial point
%   ip: (0): no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
%   nOt: maximum number of optimal points (default = 1)
%   samples: number of samples per stage (default = 3*size(x0(:),1))
%   Lb, Ub: lower and upper bound vectors to plot (default = x0*(1+/-2))
%   problem: (-1): minimum (default), (1): maximum
%   tol: tolerance (default = 1e-4)//WHAT TOL?
%   mxit: maximum number of stages (default = 50*(1+4*~(ip>0)))
%   R: axis vector of the hyperellipse centered in x0 (default = max(0.1*abs(x0+~x0),1))
%   red: factor to reduce the size of the axis vector (0,1) (default = 0.2)
%   mem: number of stored points during the exploration phase [0, min(18,samples-1)]
%   xo: matrix of optimal points (column vectors)
%   Ot: vector of optimal values of S
%   nS: number of objective function evaluations
%   Copyright (c) 2001 by LASIM-DEQUI-UFRGS
%   \$Revision: 1.0 \$  \$Date: 2001/07/10 19:40:05 \$
%   Argimiro R. Secchi (arge@enq.ufrgs.br)
%
%   Based on the algorithm of the same authors written in C
%   for constrained global optimization in 1989/09/29 and published as:
%   A.R. Secchi and C.A. Perlingeiro, "Busca Aleatoria Adaptativa",
%   in Proc. of XII Congresso Nacional de Matematica Aplicada e
%   Computacional, Sao Jose do Rio Preto, SP, pp. 49-52 (1989).
%
% Modified by Giovani Tonel(giovani.tonel@ufrgs.br) on September 2006
% Modified by Dmitrey Kroshko(openopt@ukr.net) on November 2006

r = oor('authors: Argimiro R. Secchi arge@enq.ufrgs.br, Giovani Tonel giovani.tonel@ufrgs.br', ...
'alg:A.R. Secchi and C.A. Perlingeiro, "Busca Aleatoria Adaptativa"',...
'lastChanger: Dmitrey openopt@ukr.net');

rand('state',sum(100*clock)); 	%rand('state',sum(100*clock)) resets it to
%		         a different state each time.

x0 = prob.x0;
% S = prob.f;<-don't use such names pls!//Dmitrey

prob.assert(all(isinf(prob.lb)) && all(isinf(prob.ub)), 'buscarnd can''t handle lb-ub bounds')
samples=3*size(x0,1);

tol=prob.TolX;

r.ff = inf;
r.xf = nan*ones(prob.n,1);

if ~isfield(prob,'r0') || isempty(prob.r0);
prob.err('buscarnd: You must supply prob.r0 such that norm(x0-x_opt)<= r0')
else
R=prob.r0;
end
red=0.2;
mem=min(18,samples-1);
if mem >= samples,
mem=samples-1;
end
% initialization

R=abs(R(:));
n=size(x0,1);
distrib=1;          % distr =  1  --> exploration phase
% distr >  1  --> progression phase
asym1=-1*ones(n,1); % asym1 =  1  --> asymmetric distribution to reduce r.xf
% asym1 = -1  --> asymmetric distribution to increase r.xf
asym2=2*ones(n,1);  % asym2 =  2  --> symmetric distribution
% asym2 = 1.5 --> asymmetric distribution
metric=0.5*norm(R);
a0=0.5*metric;
a2=0.1*metric;
a3=10*tol;
n3=1;
holdm=floor(1/red+0.5)*n;
nhold=0;            % counter to keep distribution constant during holdm times

top=0;
if mem > 0,
n4=10*mem;
idx=1:n4+1;
last=n4;
first=0;
idx(last+1)=first;
mem_x=zeros(n,n4);
mem_S=-ones(1,n4)*inf;
end

lim_distrib=floor(0.5*(1+log10(metric/tol)/log10(2.0)));

xOt={};

r.xf=x0;
yo=-prob.f0;
nS=1;
it=0;
opt=0;
x=zeros(n,samples);
y=zeros(1,samples);

while 1
l3=0;

for j=1:samples,   % sampling
a=min(max(rand(n,1),0.05),0.95);
x(:,j)=r.xf+(R.*asym1.*(asym2.*a-1).^distrib)/distrib;
y(j)=-feval(prob.f,x(:,j),prob);
nS=nS+1;
end
if ~it;x(:,1) = r.xf; y(1) = yo; end

it=it+1;
[r.f m_ind] = min(-y);
r.x = x(:, m_ind);
if r.f<r.ff
r.xf = r.x;
r.ff = r.f;
end
l4=mem & distrib == 1;  % sorting the samples
if samples > 1
if ~l4
[y(1),i]=max(y);
x(:,1)=x(:,i);
else
[y,i]=sort(-y);
y=-y;
x=x(:,i);
end
end
if l4 % memorize samples
n2=top;
l1=0;
for j=1:samples
n1=0;
if opt
for i=1:opt
if norm(xOt{i}-x(:,j)) < a2,
n1=1;
if j == 1,
l1=1;
end
break
end
end
end
if ~n1
for i=1:j-1
if norm(x(:,j)-x(:,i)) < metric,
n1=1;
break
end
end
if ~n1
k=idx(1);
%CHANGED BY Dmitrey
if norm(mem_x(:,k)-x(:,j)) < metric; n1=1; else
for i=1:top
k=idx(k+1);
end
end
if ~n1
first=idx(first+1);
if ~first              % expand memory
mem_x(end,end+n4)=0;
mem_S(end+1:end+n4)=-inf;
idx(end+1:end+n4)=n2+1:n2+n4;
first=n2;
idx(last+1)=first;
last=n2+n4;
idx(last+1)=0;
end
k=first;
n2=n2+1;
end
if ~n1 || y(j) > mem_S(k)
mem_x(:,k)=x(:,j);
mem_S(:,k)=y(j);
end
if l1
l1=0;
x(:,1)=x(:,j);
y(1)=y(j);
end
end
end
if n2-top == mem
break
end
end
top=n2;
end
% analysis of best sampled point
a1=norm(x(:,1)-r.xf)/(0.1+norm(r.xf))+abs(y(1)-yo)/(0.1+abs(yo));
l1=y(1) > yo;
if l1
if norm(x(:,1)-r.xf) > a0
nhold=0;
n3=1;
end

z=r.xf; r.xf=x(:,1); x(:,1)=z;
a4=yo; yo=y(1); y(1)=a4;
end
if a1 < tol && distrib >= lim_distrib
l2=1;
if opt
for i=1:opt
if norm(xOt{i}-r.xf) < a3
l2=0;
break
end
end
end

if l2
opt=opt+1;
%             ff(end+1)=yo;
xOt{end+1}=r.xf;
end

% rescue another point from memory
if ~top
if opt
r.xf=xOt{end};
end
break
end

if l2 && ~l1
n2=0;
for j=1:top,
k=idx(n2+1);
if norm(mem_x(:,k)-r.xf) < a2,
top=top-1;
idx(last+1)=k;
last=k;
idx(n2+1)=idx(k+1);
idx(k+1)=0;
if k == first,
first=n2;
end
else
n2=k;
end
end
end

if ~top
if opt,
r.xf=xOt{end};
end
break
end

k=idx(1);
r.xf=mem_x(:,k);
yo=mem_S(k);
idx(last+1)=k;
last=k;
idx(1)=idx(k+1);
idx(k+1)=1;
top=top-1;

if k == first,
first=1;
end

l3=1;

elseif opt && (~l1 || (l1 && ~l4))
n1=0;
for i=1:opt,
if norm(xOt{i}-r.xf) < a2,
n1=1;
break
end
end

if n1
if ~top,
r.xf=xOt{end};
break
end

k=idx(1);
r.xf=mem_x(:,k);
yo=mem_S(k);
idx(last+1)=k;
last=k;
idx(1)=idx(k+1);
idx(k+1)=1;
top=top-1;

if k == first,
first=1;
end

l3=1;
end
end
% adjust search direction and criterion
if l3,
n3=1;
distrib=1;
nhold=0;
elseif (a1 < a2 || ~l1) && (nhold+n3 > holdm)
n3 = n3+~mod(distrib,3);
distrib=distrib+2;
nhold=0;
else
nhold=nhold+n3;
end

for i=1:n
if l3 || (~l3 && x(i,1) == r.xf(i)),
asym2(i)=2;
else
asym2(i)=1.5;
if x(i,1) < r.xf(i),
asym1(i)=1;
else
asym1(i)=-1;
end
end
end

r = prob.iterfcn(prob, r);
if r.istop; return; end
end

% if it == mxit
%     %disp('Warning Buscarnd: reached maximum number of stages!');
%     if opt
%         yo=r.ff;
%         r.xf=xOt;
%     end
% end
%
% if opt == nOt
%     yo=r.ff;
%     r.xf=xOt;
% end
%
% if opt > 1
%     [yo,i]=sort(-yo); yo=-yo;
%     r.xf=r.xf(:,i);
% end
%
% r.ff=-yo;
```