In [1]:
from matplotlib import pyplot as plt
from mpl_toolkits import mplot3d
import numpy as np
import copy as cp
%matplotlib inline
15.1a¶
In [4]:
x = np.linspace(-1,1,1000)
y = np.linspace(-1,1,1000)
def f(x,y):
return (1-x)**2 +100*(y-x**2)**2
X, Y = np.meshgrid(x, y)
Z = f(X,Y)
In [5]:
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_wireframe(X, Y, Z, color='grey')
ax.set_title('wireframe');
ax.view_init(60, 35)
15.1b¶
In [234]:
def reflect(point,mean,factor):
return(np.add(mean,np.multiply(factor,np.subtract(mean,point))))
def shrink(points,best):
return(.5*np.add(points,best))
In [235]:
def simplex2D(x,f,evals):
#X is an array of points
#f is the function to be evaluated
tempx, tempy = np.hsplit(x, 2)
tempx = tempx.flatten()
tempy = tempy.flatten()
#Evaluate the function
f_of_x = f(tempx,tempy).transpose()
f_max = np.max(f_of_x)
f_max_loc = np.where(f_of_x==f_max)
#Get the second worst
secondmax = sorted(f_of_x)[-2]
#Get the worst point coordinates
worst = np.array([tempx[f_max_loc].item(), tempy[f_max_loc].item()])
#Get the mean of the other point
others = np.concatenate([np.delete(tempx,f_max_loc).reshape((2,1)), np.delete(tempy,f_max_loc).reshape((2,1))],axis=1)
means = others.mean(axis=0)
#Reflect
new = reflect(worst,means,1)
tempx, tempy = np.split(new, 2)
new = new.reshape((1,2))
evals += 1
if(f(tempx,tempy).item()<secondmax):
new2 = reflect(worst,means,2)
tempx, tempy = np.split(new2,2)
new2 = new2.reshape((1,2))
evals +=1
if(f(tempx,tempy).item()<secondmax):
return(np.append(others,new2,axis=0),evals)
else:
return(np.append(others,new,axis=0),evals)
#Reflect and shrink
new = reflect(worst,means,.5)
tempx, tempy = np.split(new, 2)
new = new.reshape((1,2))
evals +=1
if(f(tempx,tempy).item()<secondmax):
return(np.append(others,new,axis=0),evals)
#Shrink
f_min = np.min(f_of_x)
f_min_loc = np.where(f_of_x==f_min)
tempx, tempy = np.hsplit(x, 2)
best = np.array([tempx[f_min_loc].item(), tempy[f_min_loc].item()])
new = shrink(x,best)
evals += 1
return(new,evals)
In [237]:
tolerance = 1e-8
diff = 100
evals = 0
iters = 0
orig = np.array([[-1,-1],[-.9,-1],[-1,-.9]])
tracking = orig
while(diff>tolerance):
new, evals = simplex2D(orig,f,evals)
iters += 1
temp = np.sum(abs(new-orig),axis=1).max()
#print(temp)
tracking = np.append(tracking,orig,axis=0)
orig = new
if(temp != 0):
diff = temp
else:
diff = 1
print(f"""Final: {orig[0,]},
Iterations: {iters},
Evals: {evals}""")
In [238]:
plotting = np.append(tracking,f(tracking[:,0],tracking[:,1]).reshape(351,1),axis=1)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_wireframe(X, Y, Z, color='grey')
ax.view_init(60, 35)
ax.plot3D(plotting[:,0],plotting[:,1],plotting[:,2],color='blue')
Out[238]:
15.1c¶
In [14]:
def f10(x,d):
sum = 0
for i in range(1,d):
sum +=(1-x[i-1])**2+100*(x[i]-x[i-1]**2)**2
return sum
In [246]:
def simplexinit(x,startsize):
xs = x
xs = xs.reshape((1,10))
for i in range(0,10):
xtemp = cp.deepcopy(x)
xtemp[i] = xtemp[i] + startsize
xs = np.concatenate((xs,xtemp.reshape((1,10))),axis=0)
return(xs)
def simplex10d(x,f,evals):
#X is an 10x10 of points
#f is the function to be evaluated
#Evaluate the function
f_of_x = f10(x.transpose(),10)
f_max = np.max(f_of_x)
f_max_loc = np.where(f_of_x==f_max)[0][0]
#Get the second worst
secondmax = sorted(f_of_x)[-2]
#Get the worst point coordinates
worst = x[f_max_loc,]
#Get the best point
f_min = np.min(f_of_x)
f_min_loc = np.where(f_of_x==f_min)[0][0]
best = x[f_min_loc,]
#Get the mean of the other point
others = np.delete(x,f_max_loc,axis=0)
means = others.mean(axis=1)
#Reflect
new = reflect(worst,means,1)
new = new.reshape((10,))
evals += 1
if(f10(new,10)<secondmax):
new2 = reflect(worst,means,2)
new2 = new2.reshape((10,))
evals +=1
if(f10(new2,10)<secondmax):
new2 = new2.reshape((1,10))
#print("reflect2")
return(np.concatenate((others,new2),axis=0),evals,f_max)
else:
new = new.reshape((1,10))
#print("reflect")
return(np.concatenate((others,new),axis=0),evals,f_max)
#Reflect and shrink
new = reflect(worst,means,.5)
new = new.reshape((10,))
evals +=1
if(f10(new,10)<secondmax):
new = new.reshape((1,10))
#print("reflectshrink")
return(np.concatenate((others,new),axis=0),evals,f_max)
#Shrink
new = reflect(worst,means,(-.5))
new = new.reshape((10,))
evals +=1
if(f10(new,10)<secondmax):
new = new.reshape((1,10))
#print("shrink")
return(np.concatenate((others,new),axis=0),evals,f_max)
#Shrink all
others = np.delete(x,f_min_loc,axis=0)
new = shrink(others,best)
evals += 1
best = best.reshape((1,10))
#print("shrinkall")
return(np.concatenate((best,new),axis=0),evals,f_max)
In [249]:
x0 = -1*np.ones(10)
orig = simplexinit(x0,.05)
tolerance = 1e-8
diff = 100
evals = 0
iters = 0
tracking = cp.deepcopy(x)
best = []
for iters in range(2500):
#while diff>tolerance:
new, evals,besttemp = simplex10d(orig,f10,evals)
iters += 1
best.append(besttemp)
temp = np.sum(abs(new[10:]-orig[10:]),axis=1).max()
#tracking = np.append(tracking,orig,axis=0)
orig = new
# if iters < 3:
# diff = 100
# else:
# diff = abs(best[iters-1]-best[iters-2])
if(temp != 0):
diff = temp
else:
diff = 1
print(f"""Final: {orig[10,]},
Iterations: {iters},
Evals: {evals},
Final Function eval: {best[iters-1]}""")
Clearly, something isn't working, but I'm not sure what.
15.1 d¶
In [250]:
#Redefine Rosenbrock function for simplicity.
def f(x):
return (1-x[0])**2 +100*(x[1]-x[0]**2)**2
def dfdy(x):
return (2*(1-x[0])+200*(x[1]-x**2)*-2*x[0],200*(x[1]-x[0]**2)**2)
In [ ]: