"Brute force" optimization with Scipy

Official documentation:

Import required modules



In [ ]:

    
%matplotlib inline

import matplotlib
matplotlib.rcParams['figure.figsize'] = (8, 8)



In [ ]:

    
# Setup PyAI
import sys
sys.path.insert(0, '/Users/jdecock/git/pub/jdhp/pyai')



In [ ]:

    
import numpy as np
from scipy import optimize



In [ ]:

    
# Plot functions
from pyai.optimize.utils import plot_contour_2d_solution_space
from pyai.optimize.utils import plot_2d_solution_space

Define the objective function



In [ ]:

    
# Set the objective function
from pyai.optimize.functions import sphere1d
from pyai.optimize.functions import sphere2d

Minimize using the "Brute force" algorithm

Uses the "brute force" method, i.e. computes the function's value at each point of a multidimensional grid of points, to find the global minimum of the function.

See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brute.html#scipy.optimize.brute

First example: the 1D sphere function



In [ ]:

    
%%time

search_ranges = (slice(-3., 3.5, 0.5),)

res = optimize.brute(sphere1d,
                     search_ranges,
                     #args=params,
                     full_output=True,
                     finish=None)     # optimize.fmin)

print("x* =", res[0])
print("f(x*) =", res[1])



In [ ]:

    
print(res[2].shape)
print("tested x:", res[2])
print(res[3].shape)
print("tested f(x):", res[3])



In [ ]:

    
x_star = res[0]
y_star = res[1]

x = res[2]
y = res[3]

fig, ax = plt.subplots()

ax.set_title('Objective function')

ax.plot(x, y, 'k-', alpha=0.25, label="f")
ax.plot(x, y, 'g.', label="tested points")
ax.plot(x_star, y_star, 'ro', label="$x^*$")

ax.legend(fontsize=12);

Second example: the 2D sphere function



In [ ]:

    
%%time

search_ranges = (slice(-2., 2.5, 0.5), slice(-2., 2.5, 0.5))

res = optimize.brute(sphere2d,
                     search_ranges,
                     #args=params,
                     full_output=True,
                     finish=None)     # optimize.fmin)

print("x* =", res[0])
print("f(x*) =", res[1])



In [ ]:

    
print(res[2].shape)
print("tested x:", res[2])
print()
print(res[3].shape)
print("tested f(x):", res[3])



In [ ]:

    
# Setup data #########################

# Using the following 3 lines, pcolormesh won't display the last row and the last collumn...

#xx = res[2][0]
#yy = res[2][1]
#zz = res[3]

# Workaround to avoid pcolormesh ignoring the last row and last collumn...

x = res[2][0][:,0]
y = res[2][1][0,:]

x = np.append(x, x[-1] + x[-1] - x[-2])
y = np.append(y, y[-1] + y[-1] - y[-2])

# Make the meshgrid

xx, yy = np.meshgrid(x, y)   

# "Ideally the dimensions of X and Y should be one greater than those of C;
# if the dimensions are the same, then the last row and column of C will be ignored."
# https://stackoverflow.com/questions/44526052/can-someone-explain-this-matplotlib-pcolormesh-quirk
zz = res[3]

# Plot the image #####################

fig, ax = plt.subplots()

ax.set_title('Objective function')

# Shift to center pixels to data (workaround...)
# (https://stackoverflow.com/questions/43128106/pcolormesh-ticks-center-for-each-data-point-tile)

xx -= (x[-1] - x[-2])/2.
yy -= (y[-1] - y[-2])/2.

#im = ax.imshow(z, interpolation='bilinear', origin='lower')
im = ax.pcolormesh(xx, yy, zz, cmap='gnuplot2_r')
plt.colorbar(im)              # draw the colorbar

# Plot contours ######################

max_value = np.nanmax(zz)
levels = np.array([0.1*max_value, 0.3*max_value, 0.6*max_value])

# Shift back pixels for contours (workaround...)
# (https://stackoverflow.com/questions/43128106/pcolormesh-ticks-center-for-each-data-point-tile)

xx += (x[-1] - x[-2])/2.
yy += (y[-1] - y[-2])/2.

cs = plt.contour(xx[:-1,:-1], yy[:-1,:-1], zz, levels,
                 linewidths=(2, 2, 3), linestyles=('dotted', 'dashed', 'solid'),
                 alpha=0.5, colors='blue')
ax.clabel(cs, inline=False, fontsize=12)

# Plot x* ############################

ax.scatter(*res[0], c='red', label="$x^*$")

ax.legend(fontsize=12);