Week 8 - Search

In [355]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
from sympy import *
import numpy as np

Plot Rosenbrock function

In [356]:
def rosenbrock(x,y):
    return (1 - x)**2 + 100*(y - x**2)**2
In [357]:
def plotRosen(end, res):
    x = np.arange(-end, end, res)
    y = np.arange(-end, end, res)
    X, Y = np.meshgrid(x, y) 
    z = rosenbrock(X, Y)
    
    fig = plt.figure(figsize = (12, 10))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(X, Y, z, cmap = 'jet', alpha = .6, 
                    edgecolor = 'none' )
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.view_init(50, 35)
    plt.show()
In [358]:
plotRosen(1, 0.01)
In [359]:
plotRosen(0.1, 0.01) # for fun

b)

Pick stopping criterion and use downhill simplex (Nelder-Mead) method to search for its minimum starting x = y = -1, plotting the search path, and counting the number of algorithm iterations and function evaluations used.

This post helped in understanding the simplex method

To do's:

  • create a simplex
  • calc center of defined simplex = centroid
  • reflect the point across the face
  • check if new point is better than 2nd worst

Ordering scheme of points: worst = xh, second worst = xs, best/last = xl

In [360]:
# in this case use just one point as input - 
# might be a little dirty but works for now
def rosen(point):
    x = point[0]; y = point[1]
    return (1 - x)**2 + 100*(y - x**2)**2
In [361]:
# Downhill simplex / Nelder-Mead method

# possible update moves: reflect, reflect+grow,
# reflect+shrink, shrink, shrinkToBest 

# factor along the axes to define the initial simplex this 
# also influences the number of evaluations
u = 0.12 

# starting points for x & y
xStart = -1.0 
yStart = -1.0


# run the downhill simplex method
def simplexDownhill():
    initialSimplex = np.array([[xStart,yStart],
                               [xStart+u,yStart],
                               [xStart,yStart+u]])
    S = []
    newSimplex = stepSimplex(initialSimplex,rosen)
    S.append(initialSimplex)
    S.append(newSimplex)
    for i in range(100):
        if newSimplex is None:
            print("Number of evaluations = " + str(i))
            ##print("Minimum at = " + str(rosen(S[-2])))
            return S[0:-2]
        newSimplex = stepSimplex(newSimplex,rosen)
        S.append(newSimplex)
    return S
 
    
# Do all the steps in the method for each simplex you feed it    
def stepSimplex(simplex, f):
        
    # calculate the simplex points to be sorted 
    z = np.array([f(simplex[0]),f(simplex[1]),f(simplex[2])])
    zMax = np.argmax(z)
    zMin = np.argmin(z)

    # find the mean of the lowest points
    xMean = np.average(np.delete(z,zMax)) 
    
    # assign xh as worst, xs as second worst and xl as "best" points based on z value
    if zMax == 0: 
        xl = simplex[0]; xh = simplex[1]; xs = simplex[2]
    elif zMax == 1:
        xl = simplex[1]; xh = simplex[0]; xs = simplex[2]
    else:
        xl = simplex[2]; xh = simplex[0]; xs = simplex[1]

    # find vector from xl to xh        
    v1 = xh - xl 
    # find vector from xl to xs
    v2 = xs - xl 

    # calc the defined steps: reflect, reflect+grow, reflect+shrink, shrinkAll
    xRefl = xl + v1 + v2
    xGrow = xl + 2*v1 + 2*v2
    xShrink = xl + 0.75*v1 + 0.75*v2
    xShrinkToMin = xl + 0.5*v1 + 0.5*v2

    # assign those to new array and second array for Rosenbrock values
    X = np.array([xRefl,xGrow,xShrink,xShrinkToMin])
    Xf = np.array([f(xRefl),f(xGrow),f(xShrink),f(xShrinkToMin)])
    
    # now check if the conditions hold and assign values if they work
    if Xf[0] < xMean: #if reflect worked keep going
        if Xf[1] < Xf[0]: # if grow worked assign it
            return np.array([xh,xs,xGrow])
        elif Xf[2] < Xf[0]:  # now shrink
            return np.array([xh,xs,xShrink])
        elif Xf[3] < xMean: # and now check shrink to best point
            return np.array([xh,xs,xShrinkToMin])

    # and when all of those are done give me the points I want
    else: 
        if zMin == 0: # assign xh, xs and xl to the points
            low = simplex[0]; xl1 = simplex[1]; xl2 = simplex[2]
        elif zMin == 1:
            low = simplex[1]; xl1 = simplex[0]; xl2 = simplex[2]
        else:
            low = simplex[2]; xl1 = simplex[0]; xl2 = simplex[1]
        u1 = xl1 - low
        u2 = xl2 - low
        return np.array([low,low+0.75*u1,low+0.75*u2])

S = simplexDownhill()

Number of evaluations = 98
In [362]:
res = 0.01
end = 1;
x = np.arange(-end, end, res)
y = np.arange(-end, end, res)
X, Y = np.meshgrid(x, y) 
z = rosenbrock(X, Y)
    
fig = plt.figure(figsize = (12, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, z, cmap = 'jet', alpha = .6, 
                    edgecolor = 'none' )

for i in range(len(S)):
    x1 = S[i][0][0]; y1 = S[i][0][1];
    x2 = S[i][1][0]; y2 = S[i][1][1];
    x3 = S[i][2][0]; y3 = S[i][2][1];
    ax.plot([x1,x2],[y1,y2],[rosen(S[i][0]),rosen(S[i][1])],color="black")
    ax.plot([x2,x3],[y2,y3],[rosen(S[i][1]),rosen(S[i][2])],color="black")
    ax.plot([x3,x1],[y3,y1],[rosen(S[i][2]),rosen(S[i][0])],color="black")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.view_init(50, 35)
        
plt.show()

Try to show this in 2D

In [363]:
# 2D contour plot just for fun 
x = np.arange(-1, 1, res)
y = np.arange(-1, 1, res)
X, Y = np.meshgrid(x, y) 
z = rosenbrock(X, Y)

fig = plt.figure(figsize = (8, 8))

for i in range(len(S)):
    x1 = S[i][0][0]; y1 = S[i][0][1];
    x2 = S[i][1][0]; y2 = S[i][1][1];
    x3 = S[i][2][0]; y3 = S[i][2][1];
    plt.plot([x1,x2], [y1,y2], color="black")
    plt.plot([x2,x3], [y2,y3], color="black")
    plt.plot([x3,x1], [y3,y1], color="black")

plt.contour(X,Y,z, 150, cmap = 'jet')
Out[363]:
<matplotlib.contour.QuadContourSet at 0x12648fcf8>

2D doesn't help so much in viewing the simplex march