Tutorials

What follows are 7 pySOT tutorials that are based on a short course given at CMWR 2016, Toronto. These notebooks are available at: https://people.cam.cornell.edu/~dme65/talks.html

Tutorial 1: Hello World!

First example to show how to use pySOT in serial and synchronous parallel for bound constrained optimization problems

Step 1: Import modules and create pySOT objects (1)-(4)

# Import the necessary modules
from pySOT import *
from poap.controller import SerialController, ThreadController, BasicWorkerThread
import numpy as np

# Decide how many evaluations we are allowed to use
maxeval = 500

# (1) Optimization problem
# Use the 10-dimensional Ackley function
data = Ackley(dim=10)
print(data.info)

# (2) Experimental design
# Use a symmetric Latin hypercube with 2d + 1 samples
exp_des = SymmetricLatinHypercube(dim=data.dim, npts=2*data.dim+1)

# (3) Surrogate model
# Use a cubic RBF interpolant with a linear tail
surrogate = RBFInterpolant(kernel=CubicKernel, tail=LinearTail, maxp=maxeval)

# (4) Adaptive sampling
# Use DYCORS with 100d candidate points
adapt_samp = CandidateDYCORS(data=data, numcand=100*data.dim)

Step 2: Launch a serial controller and use standard Surrogate Optimization strategy

# Use the serial controller (uses only one thread)
controller = SerialController(data.objfunction)

# (5) Use the sychronous strategy without non-bound constraints
strategy = SyncStrategyNoConstraints(
        worker_id=0, data=data, maxeval=maxeval, nsamples=1,
        exp_design=exp_des, response_surface=surrogate,
        sampling_method=adapt_samp)
controller.strategy = strategy

# Run the optimization strategy
result = controller.run()

# Print the final result
print('Best value found: {0}'.format(result.value))
print('Best solution found: {0}'.format(
    np.array_str(result.params[0], max_line_width=np.inf,
                precision=5, suppress_small=True)))

Possible output:

Best value found: 0.211036185111
Best solution found: [ 0.02193  0.00486  0.03323  0.03656 -0.00228 -0.00414  0.05239 -0.08511 -0.0002   0.00104]

Step 3: Make a progress plot

import matplotlib.pyplot as plt

# Extract function values from the controller
fvals = np.array([o.value for o in controller.fevals])

f, ax = plt.subplots()
ax.plot(np.arange(0,maxeval), fvals, 'bo')  # Points
ax.plot(np.arange(0,maxeval), np.minimum.accumulate(fvals), 'r-', linewidth=4.0)  # Best value found
plt.xlabel('Evaluations')
plt.ylabel('Function Value')
plt.title(data.info)
plt.show()

Possible output:

_images/tutorial1_pic1.png

Step 4: Launch a threaded controller with 4 workers and use standard Surrogate Optimization strategy allowing to do 4 simultaneous in parallel. Notice how similar the code in Step 3 is to the code in Step 2.

# Use the threaded controller
controller = ThreadController()

# (5) Use the sychronous strategy without non-bound constraints
# Use 4 threads and allow for 4 simultaneous evaluations
nthreads = 4
strategy = SyncStrategyNoConstraints(
        worker_id=0, data=data, maxeval=maxeval, nsamples=nthreads,
        exp_design=exp_des, response_surface=surrogate,
        sampling_method=adapt_samp)
controller.strategy = strategy

# Launch the threads and give them access to the objective function
for _ in range(nthreads):
    worker = BasicWorkerThread(controller, data.objfunction)
    controller.launch_worker(worker)

# Run the optimization strategy
result = controller.run()

# Print the final result
print('Best value found: {0}'.format(result.value))
print('Best solution found: {0}'.format(
    np.array_str(result.params[0], max_line_width=np.inf,
                precision=5, suppress_small=True)))

Possible output:

Best value found: 0.0143986694788
Best solution found: [-0.00869 -0.00018 -0.00311 -0.00211  0.00098  0.00161  0.00219 -0.00241  0.00285  0.00255]

Step 5 Make a progress plot

import matplotlib.pyplot as plt

# Extract function values from the controller
fvals = np.array([o.value for o in controller.fevals])

f, ax = plt.subplots()
ax.plot(np.arange(0,maxeval), fvals, 'bo')  # Points
ax.plot(np.arange(0,maxeval), np.minimum.accumulate(fvals), 'r-', linewidth=4.0)  # Best value found
plt.xlabel('Evaluations')
plt.ylabel('Function Value')
plt.title(data.info)
plt.show()

Possible output:

_images/tutorial1_pic2.png

Tutorial 2: Python objective function

This example shows how to define our own optimization problem in pySOT

Step 1: Define our own optimization problem

class SomeFun:
    def __init__(self, dim=10):
        self.xlow = -10 * np.ones(dim) # lower bounds
        self.xup = 10 * np.ones(dim) # upper bounds
        self.dim = dim # dimensionality
        self.info = "Our own " + str(dim)+"-dimensional function" # info
        self.integer = np.array([0]) # integer variables
        self.continuous = np.arange(1, dim) # continuous variables

    def objfunction(self, x):
        return np.sum(x) * np.cos(np.sum(x))

Step 2: Let’s make sure that our optimization problem follows the pySOT standard

import numpy as np
from pySOT import check_opt_prob
data = SomeFun(dim=10)
check_opt_prob(data)

Everything is fine as long as pySOT doesn’t complain!

Step 3: Import modules and create pySOT objects (1)-(4)

# Import the necessary modules
from pySOT import *
from poap.controller import SerialController, BasicWorkerThread

# Decide how many evaluations we are allowed to use
maxeval = 500

# (1) Optimization problem
# Use our 10-dimensional function
print(data.info)

# (2) Experimental design
# Use a symmetric Latin hypercube with 2d + 1 samples
exp_des = SymmetricLatinHypercube(dim=data.dim, npts=2*data.dim+1)

# (3) Surrogate model
# Use a cubic RBF interpolant with a linear tail
surrogate = RBFInterpolant(kernel=CubicKernel, tail=LinearTail, maxp=maxeval)

# (4) Adaptive sampling
# Use DYCORS with 100d candidate points
adapt_samp = CandidateDYCORS(data=data, numcand=100*data.dim)

Output:

Our own 10-dimensional function

Step 4: Run the optimization in serial

# Use the serial controller (uses only one thread)
controller = SerialController(data.objfunction)

# (5) Use the sychronous strategy without non-bound constraints
strategy = SyncStrategyNoConstraints(
        worker_id=0, data=data, maxeval=maxeval, nsamples=1,
        exp_design=exp_des, response_surface=surrogate,
        sampling_method=adapt_samp)
controller.strategy = strategy

# Run the optimization strategy
result = controller.run()

# Print the final result
print('Best value found: {0}'.format(result.value))
print('Best solution found: {0}'.format(
    np.array_str(result.params[0], max_line_width=np.inf,
                precision=5, suppress_small=True)))

Possible output:

Best value found: -72.2440613978
Best solution found: [ 9.       5.58049  9.34501  5.35848  9.26448  9.05695  5.45796  1.80559  8.16331  9.21498]

Step 5: Plot the progress:

import matplotlib.pyplot as plt

# Extract function values from the controller
fvals = np.array([o.value for o in controller.fevals])

f, ax = plt.subplots()
ax.plot(np.arange(0,maxeval), fvals, 'bo')  # Points
ax.plot(np.arange(0,maxeval), np.minimum.accumulate(fvals), 'r-', linewidth=4.0)  # Best value found
plt.xlabel('Evaluations')
plt.ylabel('Function Value')
plt.title(data.info)
plt.show()

Possible output:

_images/tutorial2_pic1.png

Tutorial 3: MATLAB objective function

This example shows how to use pySOT with an objective function that is written in MATLAB. You need the matlab_wrapper module or the MATLAB engine which is available for MATLAB R2014b or later. This example uses the matlab_wrapper module to work for older versions of MATLAB as well.

Step 1: The following shows an implementation of the Ackley function in a MATLAB script matlab_ackley.m that takes a variable x from the workspace and saves the value of the objective function as val:

dim = length(x);
val = -20*exp(-0.2*sqrt(sum(x.^2,2)/dim)) - ...
  exp(sum(cos(2*pi*x),2)/dim) + 20 + exp(1);

Step 2: This will create a MATLAB session. You may need to specify the root folder of your MATLAB installation. Type matlabroot in a MATLAB session to see what the root folder is.

import matlab_wrapper
matlab = matlab_wrapper.MatlabSession(matlab_root='/Applications/MATLAB_R2014a.app', options='-nojvm')

Step 3: Define Python optimization problem that uses our MATLAB objective function to do function evaluations:

# This is the path to the external MATLAB function, assuming it is in your current path
import os
mfile_location = os.getcwd()
matlab.workspace.addpath(mfile_location)

class AckleyExt:
    def __init__(self, dim=10):
        self.xlow = -15 * np.ones(dim)
        self.xup = 20 * np.ones(dim)
        self.dim = dim
        self.info = str(dim) + "-dimensional Ackley function \n" + \
            "Global optimum: f(0,0,...,0) = 0"
        self.min = 0
        self.integer = []
        self.continuous = np.arange(0, dim)

    def objfunction(self, x):
        matlab.put('x', x)
        matlab.eval('matlab_ackley')
        val = matlab.get('val')
        return val

** Step 4:** Optimize over our optimization problem

from pySOT import *
from poap.controller import SerialController, BasicWorkerThread
import numpy as np

maxeval = 500

data = AckleyExt(dim=10)
print(data.info)

# Use the serial controller for simplicity
# In order to run in parallel we need to maintain an array of MATLAB session
controller = SerialController(data.objfunction)
controller.strategy = \
    SyncStrategyNoConstraints(
        worker_id=0, data=data,
        maxeval=maxeval, nsamples=1,
        exp_design=LatinHypercube(dim=data.dim, npts=2*(data.dim+1)),
        response_surface=RBFInterpolant(kernel=CubicKernel, tail=LinearTail, maxp=maxeval),
        sampling_method=CandidateDYCORS(data=data, numcand=100*data.dim))

# Run the optimization strategy
result = controller.run()

# Print the final result
print('Best value found: {0}'.format(result.value))
print('Best solution found: {0}'.format(
    np.array_str(result.params[0], max_line_width=np.inf,
                 precision=5, suppress_small=True)))

Possible output:

10-dimensional Ackley function
Global optimum: f(0,0,...,0) = 0
Best value found: 0.00665167450159
Best solution found: [-0.00164  0.00162 -0.00122  0.0019  -0.00109  0.00197 -0.00102 -0.00124 -0.00194  0.00216]

Step 5: Plot the progress:

import matplotlib.pyplot as plt

# Extract function values from the controller
fvals = np.array([o.value for o in controller.fevals])

f, ax = plt.subplots()
ax.plot(np.arange(0,maxeval), fvals, 'bo')  # Points
ax.plot(np.arange(0,maxeval), np.minimum.accumulate(fvals), 'r-', linewidth=4.0)  # Best value found
plt.xlabel('Evaluations')
plt.ylabel('Function Value')
plt.title(data.info)
plt.show()

Possible output:

_images/tutorial3_pic1.png

Step 6: End the MATLAB session:

matlab.__del__()

Note: The example test_matlab_engine.py in pySOT.test shows how to use a MATLAB engine with more than 1 worker. The main idea is the give each worker its own MATLAB session that the worker can do for function evaluations.

Tutorial 4: Python objective function with inequality constraints

This example considers the Keane bump function which has two inequality constraints and takes the following form:

\[f(x_1,\ldots,x_d) = -\left| \frac{\sum_{j=1}^d \cos^4(x_j) - \ 2 \prod_{j=1}^d \cos^2(x_j)}{\sqrt{\sum_{j=1}^d jx_j^2}} \right|\]

subject to:

\[0 \leq x_i \leq 5\]
\[0.75 - \prod_{j=1}^d x_j < 0\]
\[\sum_{j=1}^d x_j - 7.5d < 0\]

The global optimum approaches -0.835 when d goes to infinity. We will use pySOT and the penalty method approach to optimize over the Keane bump function and we will use 4 workers in synchronous parallel. The code that achieves this is

from pySOT import *
from poap.controller import Threaded, BasicWorkerThread
import numpy as np

maxeval = 500

data = Keane(dim=10)
print(data.info)

controller = ThreadController()

# Use 4 threads and allow for 4 simultaneous evaluations
nthreads = 4
strategy = SyncStrategyPenalty(
        worker_id=0, data=data,
        maxeval=maxeval, nsamples=1,
        exp_design=LatinHypercube(dim=data.dim, npts=2*(data.dim+1)),
        response_surface=RBFInterpolant(kernel=CubicKernel, tail=LinearTail, maxp=maxeval),
        sampling_method=CandidateDYCORS(data=data, numcand=100*data.dim))
controller.strategy = strategy

# Launch the threads and give them access to the objective function
for _ in range(nthreads):
    worker = BasicWorkerThread(controller, data.objfunction)
    controller.launch_worker(worker)

# Returns f(x) is feasible, infinity otherwise
def feasible_merit(record):
    return record.value if record.feasible else np.inf

# Run the optimization strategy and ask the controller for the best FEASIBLE solution
result = controller.run(merit=feasible_merit)
best, xbest = result.value, result.params[0]

# Print the final result
print('Best value found: {0}'.format(result.value))
print('Best solution found: {0}'.format(
    np.array_str(result.params[0], max_line_width=np.inf,
                precision=5, suppress_small=True)))

# Check constraints
print('\nConstraint 1: 0.75 - prod(x) = {0}'.format(0.75 - np.prod(xbest)))
print('Constraint 2: sum(x) - 7.5*dim = {0}'.format(np.sum(xbest) - 7.5*data.dim))

Possible output:

Best value found: -0.683081148607
Best solution found: [ 3.11277  3.07498  2.91834  2.96004  2.84659  1.29008  0.17825  0.31923  0.19628  0.24831]

Constraint 1: 0.75 - prod(x) = -0.0921329885647
Constraint 2: sum(x) - 7.5*dim = -57.8551318917

A possible progress plot is:

_images/tutorial4_pic1.png

Tutorial 5: Equality constraints

The only way pySOT supports inequality constraints is via a projection method. That is, the user needs to supply a method that projects any infeasible point onto the feasible region. This is trivial in some cases such as when we have a normalization constraint of the form \(g(x) = 1 - \|x\| = 0\), in which case we can just rescale each infeasible point. The purpose of this example is to show how to use pySOT for such a constraint and we will modify the Ackley function by adding a constraint that the solution needs to have unit 2-norm. Our new objective function takes the form

class AckleyUnit:
    def __init__(self, dim=10):
        self.xlow = -1 * np.ones(dim)
        self.xup = 1 * np.ones(dim)
        self.dim = dim
        self.info = str(dim)+"-dimensional Ackley function on the unit sphere \n" +\
                             "Global optimum: f(1,0,...,0) = ... = f(0,0,...,1) = " +\
                             str(np.round(20*(1-np.exp(-0.2/np.sqrt(dim))), 3))
        self.min = 20*(1 - np.exp(-0.2/np.sqrt(dim)))
        self.integer = []
        self.continuous = np.arange(0, dim)
        check_opt_prob(self)

    def objfunction(self, x):
        n = float(len(x))
        return -20.0 * np.exp(-0.2*np.sqrt(np.sum(x**2)/n)) - \
            np.exp(np.sum(np.cos(2.0*np.pi*x))/n) + 20 + np.exp(1)

    def eval_eq_constraints(self, x):
        return np.linalg.norm(x) - 1

We next define a projection method as follows:

import numpy as np

def projection(x):
    return x / np.linalg.norm(x)

Optimizing over this function is done via

from pySOT import *
from poap.controller import Threaded, BasicWorkerThread
import numpy as np

maxeval = 500

data = AckleyUnit(dim=10)
print(data.info)

controller = ThreadController()

# Use 4 threads and allow for 4 simultaneous evaluations
nthreads = 4
strategy = SyncStrategyProjection(
        worker_id=0, data=data,
        maxeval=maxeval, nsamples=1,
        exp_design=LatinHypercube(dim=data.dim, npts=2*(data.dim+1)),
        response_surface=RBFInterpolant(kernel=CubicKernel, tail=LinearTail, maxp=maxeval),
        sampling_method=CandidateDYCORS(data=data, numcand=100*data.dim),
        proj_fun=projection)
controller.strategy = strategy

# Launch the threads and give them access to the objective function
for _ in range(nthreads):
    worker = BasicWorkerThread(controller, data.objfunction)
    controller.launch_worker(worker)

# Run the optimization strategy and ask the controller for the best FEASIBLE solution
result = controller.run()

# Print the final result
print('Best value found: {0}'.format(result.value))
print('Best solution found: {0}'.format(
    np.array_str(result.params[0], max_line_width=np.inf,
                precision=5, suppress_small=True)))

# Check constraint
print('\n ||x|| = {0}'.format(np.linalg.norm(result.params[0])))

Possible output:

Best value found: 1.22580826108
Best solution found: [-0.00017  0.00106  0.00172 -0.00126  0.0013  -0.00035  0.00133  0.99999 -0.00114  0.00138]
||x|| = 1.0

A possible progress plot if the following:

_images/tutorial5_pic1.png

Tutorial 6: C++ objective function

Stay patient!