Benchmarking the Fractal Farmer algorithm¶
Sergio Hernández y Guillem Duran¶
Source code for the benchmarking process¶
In [2]:
from bug import *
from scipy import optimize
import numpy as np
import pandas as pd
In [3]:
%pylab inline
In [4]:
def calc_fun(function,X,Y,Z):
for i in np.arange(len(X[0])):
for j in np.arange(len(X[1])):
Z[i][j]=function(np.array([X[i][j],Y[i][j]]))
return Z.copy()
def mccornick(x):
return np.sin(x[0]+x[1])+(x[0]-x[1])**2-1.5*x[0]+2.5*x[1]+1
In [5]:
class MyBounds(object):
def __init__(self, xmax=[513,513], xmin=[-513,-513] ):
self.xmax = np.array(xmax)
self.xmin = np.array(xmin)
def __call__(self, **kwargs):
x = kwargs["x_new"]
tmax = bool(np.all(x <= self.xmax))
tmin = bool(np.all(x >= self.xmin))
return tmax and tmin
In [6]:
def run_benchmark(fun,x0,domain,nsa=100,nfa=100,nbh=100,maxticks=16,numbolts=1,seed=1):
xmin=domain[0][0][0]
xmax=domain[0][0][1]
ymin=domain[1][0][0]
ymax=domain[1][0][1]
myopts = {
'schedule' : 'boltzmann', # Non-default value.
'maxfev' : None, # Default, formerly `maxeval`.
'maxiter' : nsa, # Non-default value.
'maxaccept' : None, # Default value.
'ftol' : 1e-6, # Default, formerly `feps`.
'T0' : None, # Default value.
'Tf' : 1e-6, # Default value.
'boltzmann' : 1.0, # Default value.
'learn_rate' : 0.5, # Default value.
'quench' : 1.0, # Default value.
'm' : 1.0, # Default value.
'n' : 1.0, # Default value.
'lower' : -10, # Non-default value.
'upper' : +10, # Non-default value.
'dwell' : 250, # Non-default value.
'disp' : True # Default value.
}
cons = ({'type': 'ineq', 'fun': lambda x: -x[0] + xmin},
{'type': 'ineq', 'fun': lambda x: x[0] + xmax},
{'type': 'ineq', 'fun': lambda x: x[1] + ymax},
{'type': 'ineq', 'fun': lambda x: -x[1] + ymin})
np.random.seed(seed) # Seeded to allow replication.
anneal = optimize.minimize(fun, x0, method='Anneal',
options=myopts,constraints=cons)
np.random.seed(seed)
funcion=ObjFunc(fun,domain=domain)
fa=Farmer(funcion,maxticks=maxticks,
numbolts=numbolts,speed=1,
minmax=-1.0,maxfails=11000,
maxepoch=nfa,sens=8)
fa.optimize()
np.random.seed(seed)
mybounds = MyBounds(xmax=[xmax,ymax],xmin=[xmin,ymin])
minimizer_kwargs = {"method": "BFGS"}
bsh = optimize.basinhopping(fun, x0,
minimizer_kwargs=minimizer_kwargs,
niter=6000,niter_success=nbh,accept_test=mybounds)
return anneal,fa,bsh
In [7]:
def comparativa(sa,fa,bh,solution):
logs=fa.logger.data
fab=fa.best['val']
faev=logs['ticks'].dropna().cumsum().tail(1).values[0]*len(fa.bolts)
sab=sa.fun
saev=sa.nfev
bhb=bh.fun
bhev=bh.nfev
benchmark=pd.DataFrame(index=['Farmer','Simulated Annealing','Basin Hopping'],columns=['Evaluations','Minimum Found','Error'])
benchmark['Evaluations']=np.array([faev,saev,bhev])
benchmark['Minimum Found']=np.array([fab,sab,bhb])
benchmark['Error']=np.array([solution-fab,solution-sab,solution-bhb])
return benchmark
In [8]:
def create_index(samples,names):
algos=np.array(names)
tiled=np.tile(algos,samples)
iters=np.repeat(np.arange(samples),len(algos))
tuples=list(zip(iters,tiled))
multiindex = pd.MultiIndex.from_tuples(tuples, names=['iter', 'algorithm'])
return multiindex
In [9]:
names=['Farmer','Simulated Annealing','Basin Hopping']
def benchmark_df(samples,names):
return pd.DataFrame(index=create_index(samples,names),columns=['Evaluations','Minimum Found','Error'])
In [10]:
def results_df():
benchmark=pd.DataFrame(index=['Farmer','Simulated Annealing','Basin Hopping'],columns=['Evaluations','Minimum Found','Error'])
In [11]:
def record_benchmark(df,sa,fa,bh,solution,i):
logs=fa.logger.data
fab=fa.best['val']
faev=logs['ticks'].dropna().cumsum().tail(1).values[0]*len(fa.bolts)
sab=sa.fun
saev=sa.nfev
bhb=bh.fun
bhev=bh.nfev
df.loc[i,'Farmer']['Evaluations']=faev
df.loc[i,'Simulated Annealing']['Evaluations']=saev
df.loc[i,'Basin Hopping']['Evaluations']=bhev
df.loc[i,'Farmer']['Minimum Found']=fab
df.loc[i,'Simulated Annealing']['Minimum Found']=sab
df.loc[i,'Basin Hopping']['Minimum Found']=bhb
df.loc[i,'Farmer']['Error']=fab-solution
df.loc[i,'Simulated Annealing']['Error']=sab-solution
df.loc[i,'Basin Hopping']['Error']=bhb-solution
In [12]:
def comparativa(df):
return df.applymap(lambda x:float(x)).reset_index().drop('iter',axis=1).groupby('algorithm').describe()
In [13]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
def make_plot(fun,domain,points=100):
xmin=domain[0][0][0]
xmax=domain[0][0][1]
ymin=domain[1][0][0]
ymax=domain[1][0][1]
x = np.linspace(xmin, xmax, points)
y = np.linspace(ymin, ymax, points)
X,Y = np.meshgrid(x, y)
Z=X.copy()
Z=calc_fun(fun,X,Y,Z)
# `ax` is a 3D-aware axis instance because of the projection='3d' keyword argument to add_subplot
fig = plt.figure(figsize=(16,10))
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z, rstride=4, cstride=4, alpha=0.25,cmap=cm.gist_rainbow)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$y$', fontsize=18)
ax.set_xlim3d(xmin, xmax)
ax.set_ylim3d(ymin, ymax)
return ax,fig
In [14]:
def random_in_domain(domain):
xmin=domain[0][0][0]
xmax=domain[0][0][1]
ymin=domain[1][0][0]
ymax=domain[1][0][1]
return np.array([np.random.uniform(xmin,xmax),np.random.uniform(ymin,ymax)])
Benchmark¶
Schaffel function number 2¶
In [15]:
def schaffel2(X):
x=X[0]
y=X[1]
return 0.5+(np.sin((x**2-y**2))**2-0.5)/(1+0.001*(x**2+y**2))**2
In [19]:
%%capture
dsc4=[[(-100,100)],[(-100,100)]]
runs=100
bm=benchmark_df(runs,names)
for i in range(runs):
x0=random_in_domain(dsc4)
sa,fa,bh=run_benchmark(schaffel2,x0,dsc4,nsa=500,nfa=30,nbh=50,maxticks=8,numbolts=10,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=0.0,i=i)
In [20]:
ax,fig=make_plot(schaffel2,dsc4,500)
p=ax.set_title(r'$Schaffel2$',fontsize=18)
comparativa(bm)
Out[20]:
Schaffel function number 4¶
In [21]:
def schaffel4(X):
x=X[0]
y=X[1]
return 0.5+(np.cos(np.sin(np.abs(x**2-y**2)))**2-0.5)/(1+0.001*(x**2+y**2))**2
In [22]:
%%capture
dsc4=[[(-100,100)],[(-100,100)]]
runs=100
bm=benchmark_df(runs,names)
for i in range(runs):
x0=random_in_domain(dsc4)
sa,fa,bh=run_benchmark(schaffel4,x0,dsc4,nsa=500,nfa=60,nbh=110,maxticks=8,numbolts=10,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=0.292579,i=i)
In [23]:
ax,fig=make_plot(schaffel4,dsc4,500)
p=ax.set_title(r'$Schaffel 4$',fontsize=18)
comparativa(bm)
Out[23]:
Test function from scipy optimize¶
In [24]:
params = (2, 3, 7, 8, 9, 10, 44, -1, 2, 26, 1, -2, 0.5)
def f1(z,):
x, y = z
a, b, c, d, e, f, g, h, i, j, k, l, scale = params
return (a * x**2 + b * x * y + c * y**2 + d*x + e*y + f)
def f2(z,):
x, y = z
a, b, c, d, e, f, g, h, i, j, k, l, scale = params
return (-g*np.exp(-((x-h)**2 + (y-i)**2) / scale))
def f3(z,):
x, y = z
a, b, c, d, e, f, g, h, i, j, k, l, scale = params
return (-j*np.exp(-((x-k)**2 + (y-l)**2) / scale))
def f(z,):
x, y = z
a, b, c, d, e, f, g, h, i, j, k, l, scale = params
return f1(z) + f2(z) + f3(z)
In [25]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dsc4=[[(-100,100)],[(-100,100)]]
for i in range(runs):
x0=random_in_domain(dsc4)
sa,fa,bh=run_benchmark(f,x0,dsc4,nsa=500,nfa=100,nbh=150,maxticks=8,numbolts=40,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=-3.409,i=i)
In [26]:
ax,fig=make_plot(f,dsc4,500)
#p=ax.set_title(r'$Schaffel 4$',fontsize=18)
comparativa(bm)
Out[26]:
Easom Function¶
In [27]:
def easom(x):
return -np.cos(x[0])*np.cos(x[1])*np.exp(-((x[0]-np.pi)**2+(x[1]-np.pi)**2))
In [28]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dsc4=[[(-100,100)],[(-100,100)]]
for i in range(runs):
x0=random_in_domain(dsc4)
sa,fa,bh=run_benchmark(easom,x0,dsc4,nsa=40,nfa=50,nbh=500,maxticks=8,numbolts=100,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=-1,i=i)
In [29]:
dplot=[[(-40,40)],[(-40,40)]]
ax,fig=make_plot(easom,dplot,500)
ax.set_title(r'$Easom function$',fontsize=18)
#p=ax.set_title(r'$Schaffel 4$',fontsize=18)
comparativa(bm)
Out[29]:
Eggholder¶
In [30]:
def eggholder(x):
return -1*(x[1]+47)*np.sin(np.sqrt(np.abs(x[1]+x[0]/2.0+47.0)))-x[0]*np.sin(np.sqrt(np.abs(x[0]-(x[1]+47.0))))
In [31]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dom=[[(-512,512)],[(-512,512)]]
for i in range(runs):
x0=random_in_domain(dom)
sa,fa,bh=run_benchmark(eggholder,x0,dom,nsa=2500,nfa=1000,nbh=3500,maxticks=8,numbolts=120,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=-959.6407,i=i)
In [32]:
ax,fig=make_plot(eggholder,dom,1000)
miau=ax.set_title(r'$Eggholder$',fontsize=18)
comparativa(bm)
Out[32]:
Trefethen function¶
In [33]:
def trefethen(X):
x=X[0]
y=X[1]
return np.exp(np.sin(50.0*x))+np.sin(60.0*np.exp(y))+np.sin(70.0*np.sin(x))+np.sin(np.sin(80.0*y))-np.sin(10.0*(x+y))+0.25*(x**2+y**2)
In [34]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dom=dsc4
for i in range(runs):
x0=random_in_domain(dom)
sa,fa,bh=run_benchmark(trefethen,x0,dsc4,nsa=1500,nfa=500,nbh=1500,maxticks=8,numbolts=40,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=-3.3068686474,i=i)
In [35]:
ax,fig=make_plot(trefethen,dsc4,600)
miau=ax.set_title(r'$Trefethen$',fontsize=18)
comparativa(bm)
Out[35]:
Beale's Function¶
In [36]:
def beales(X):
x=X[0]
y=X[1]
return (1.5-x+x*y)**2+(2.25-x+x*y**2)**2+(2.625-x+x*y**3)**2
dbeale=[[(-4.5,4.5)],[(-4.5,4.5)]]
In [37]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dom=dbeale
for i in range(runs):
x0=random_in_domain(dom)
sa,fa,bh=run_benchmark(beales,x0,dbeale,nsa=500,nfa=100,nbh=1000,maxticks=8,numbolts=90,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=0,i=i)
In [38]:
ax,fig=make_plot(beales,dbeale,600)
miau=ax.set_title(r'$Beale$',fontsize=18)
comparativa(bm)
Out[38]:
Goldstein-Price Function¶
In [39]:
def gprice(X):
x=X[0]
y=X[1]
return (1+((x+y+1)**2)*(19-14*x+3*x**2-14*y+6*x*y+3*y**2))*(30+((2*x-3*y)**2)*(18-32*x+12*x**2+48*y-36*x*y+27*y**2))
dprice=[[(-2,2)],[(-2,2)]]
dpriceplot=[[(-2,1)],[(-2,1)]]
In [40]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dom=dprice
for i in range(runs):
x0=random_in_domain(dom)
sa,fa,bh=run_benchmark(gprice,x0,dprice,nsa=500,nfa=200,nbh=100,maxticks=8,numbolts=50,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=3,i=i)
In [41]:
ax,fig=make_plot(gprice,dpriceplot,700)
miau=ax.set_title(r'$Goldstein-Price$',fontsize=18)
comparativa(bm)
Out[41]:
Matyas Function¶
In [42]:
def matyas(x):
return 0.26*(x[0]**2+x[1]**2)-0.48*x[0]*x[1]
dmat=[[(-10,10)],[(-10,10)]]
In [43]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dom=dmat
for i in range(runs):
x0=random_in_domain(dom)
sa,fa,bh=run_benchmark(matyas,x0,dmat,nsa=500,nfa=20,nbh=100,maxticks=8,numbolts=50,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=0,i=i)
In [44]:
ax,fig=make_plot(matyas,dmat,700)
miau=ax.set_title(r'$Matyas$',fontsize=18)
comparativa(bm)
Out[44]:
Bukin Function number 6¶
In [45]:
def bukin6(X):
x=X[0]
y=X[1]
return 100.0*np.sqrt(np.abs(y-0.01*x**2))+0.01*np.abs(x+10)
dbuk=[[(-15,5)],[(-3,3)]]
In [46]:
%%capture
runs=100
bm=benchmark_df(runs,names)
dom=dbuk
for i in range(runs):
x0=random_in_domain(dom)
sa,fa,bh=run_benchmark(bukin6,x0,dbuk,nsa=2500,nfa=400,nbh=5000,maxticks=8,numbolts=100,seed=i*10+1)
record_benchmark(bm,sa,fa,bh,solution=0,i=i)
In [47]:
ax,fig=make_plot(bukin6,dbuk,700)
miau=ax.set_title(r'$Bukin6$',fontsize=18)
comparativa(bm)
Out[47]:
In [ ]: