"""Main optimization algorithm.
This module contains a class that implements the main optimization
algorithm. The class has a state so that optimization can be stopped
and resumed at any time.
Licensed under Revised BSD license, see LICENSE.
(C) Copyright Singapore University of Technology and Design 2016.
(C) Copyright International Business Machines Corporation 2016.
"""
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
import sys
import math
import time
import os
import pickle
import copy
import collections
import numpy as np
from multiprocessing import Pool
import rbfopt.rbfopt_utils as ru
import rbfopt.rbfopt_aux_problems as aux
import rbfopt.rbfopt_refinement as ref
from rbfopt.rbfopt_black_box import RbfoptBlackBox
from rbfopt.rbfopt_settings import RbfoptSettings
[docs]class RbfoptAlgorithm:
"""Optimization algorithm.
Implements the main optimization algorithm, and contains all its
state variables so that the algorithm can be warm-started from a
given state. Some of the logic of the algorithm is not inside this
class, because it can be run in parallel and therefore should not
modify the state.
Parameters
----------
settings : :class:`rbfopt_settings.RbfoptSettings`
Global and algorithmic settings.
black_box : :class:`rbfopt_black_box.BlackBox`
An object derived from class BlackBox, that describes the
problem.
init_node_pos : 2D numpy.ndarray[float] or None
Coordinates of points at which the function value is
known. These points will be added to the points generated by
the algorithm. If None, all the initial points will be
generated by the algorithm.
init_node_val : 1D numpy.ndarray[float] or None
Function values corresponding to the points given in
init_node_pos. Should be None if the previous argument is
None. If init_node_pos is not None but init_node_val is None,
the points in init_node_pos will be evaluated in the
initialization phase of the algorithm.
do_init_strategy : bool
Perform initialization strategy. If True, the algorithm will
generate an initial set of points on top of any user-provided
points. If False, we will rely on user-provided points
only. Notice that the algorithm may fail if not enough points
are provided, and do_init_strategy is False. Default True.
Attributes
----------
elapsed_time : float
Elapsed CPU time up to the point where state was last saved.
best_local_rbf : (string, float)
Best RBF type to construct model for local search, and the
corresponding shape parameter.
best_global_rbf : (string, float)
Best RBF type to construct model for global search, and the
corresponding shape parameter.
n : int
Dimension of the problem.
itercount : int
Iteration number.
evalcount : int
Total number of function evaluations in accurate mode.
noisy_evalcount : int
Total number of function evaluations in noisy mode.
current_step : int
Identifier of the current step within the cyclic optimization
strategy counter.
cyclecount : int
Number of cycles.
num_cons_refinement : int
Current number of consecutive refinement steps.
num_stalled_iter : int
Number of consecutive iterations without improvement.
num_noisy_restarts : int
Number of restarts in noisy mode.
discarded_iters : collections.deque
Rolling window to keep track of discarded iterations
do_init_strategy : bool
If True, perform initialization strategy on first cycle.
inf_step : int
Identifier of the InfStep.
local_search_step : int
Identifier of the LocalSearchStep.
refinement_step : int
Identifier of the Step.
cycle_length : int
Length of an optimization cycle.
restoration_step : int
Identifier of the RestorationStep.
first_step : int
Identifier of the first step of an optimization cycle.
two_phase_optimization : bool
Is the fast but noisy objective function is available?
eval_mode : string
Evaluation mode for the objective function at a given
stage. Can be either 'noisy' or 'accurate'.
node_pos : 2D numpy.ndarray[float]
Coordinates of the interpolation nodes (i.e. the points where
the objective function has already been evaluated). The
coordinates may be in the transformed space. This matrix only
includes points since the last restart.
node_val : 1D numpy.ndarray[float]
Objective function value at the points in node_pos. This array
only includes points since the last restart.
node_is_noisy : 1D numpy.ndarray[bool]
For each interpolation node in node_pos, was it evaluated in
'noisy' mode?
node_err_bounds : 2D numpy.ndarray[float]
The lower and upper variation of the function value for the
nodes in node_pos. The variation is assumed 0 for nodes
evaluated in accurate mode.
all_node_pos : 2D numpy.ndarray[float]
Coordinates of the interpolation nodes. This matrix contains all
evaluated points in the original space, and is persistent
across restarts.
all_node_val : 1D numpy.ndarray[float]
Objective function value at the points in all_node_pos.
all_node_is_noisy : 2D numpy.ndarray[bool]
For each interpolation node in all_node_pos, was it evaluated
in 'noisy' mode?
all_node_err_bounds : 2D numpy.ndarray[float]
The lower and upper variation of the function value for the
nodes in all_node_pos. The variation is assumed 0 for nodes
evaluated in accurate mode.
num_nodes_at_restart : int
Index of the first new node in all_node_pos after the latest
restart.
l_lower : 1D numpy.ndarray[float]
Variable lower bounds in the transformed space.
l_upper : 1D numpy.ndarray[float]
Variable upper bounds in the transformed space.
fmin_index : int
Index of the minimum value among the nodes since last restart.
fmin : float
Minimum value among the nodes since last restart.
fmax : float
Maximum value among the nodes since last restart.
fmin_stall_check : float
Best function value at the beginning of the most recent
optimization cycle.
fmin_last_refine : float
Best function value at the most recent refinement step.
iter_last_refine : int
Iteration number of the most recent refinement step.
fbest_index : int
Index in all_node_pos of the minimum value among all nodes.
fbest : float
Minimum value among all nodes.
best_gap_shown : float
Best gap shown on the log.
is_fmin_noisy : bool
Was the best known objective function since restart evaluated
in noisy mode?
is_fbest_noisy : bool
Was the best known objective function evaluated in noisy mode?
fixed_vars : List[(int, float)]
Indices and values of fixed variables. The indices are with
respect to the vector of all variables.
tr_model_set : 1D numpy.ndarray[int]
Indices of nodes (in node_pos) that are used to build the
linear model for the trust region refinement phase.
tr_radius : float
Radius of the trust region.
tr_iterate_index : int
Index of the node (in node_pos) that is the current iterate
for the trust region refinement method.
"""
def __init__(self, settings, black_box, init_node_pos=None,
init_node_val=None, do_init_strategy=True):
"""Constructor.
"""
assert(isinstance(black_box, RbfoptBlackBox))
assert(isinstance(settings, RbfoptSettings))
# Start timing
self.elapsed_time = 0.0
start_time = time.time()
# Save references of initial data
self.settings = settings
self.bb = black_box
dimension = black_box.get_dimension()
assert(dimension >= 1)
var_lower = black_box.get_var_lower()
var_upper = black_box.get_var_upper()
var_type = black_box.get_var_type()
assert(len(var_lower) == dimension)
assert(len(var_upper) == dimension)
assert(len(var_type) == dimension)
# We make sure these vectors are Numpy arrays
var_lower = np.array(var_lower, dtype=np.float_)
var_upper = np.array(var_upper, dtype=np.float_)
integer_vars = np.array([i for i in range(dimension)
if var_type[i].upper() == 'I'],
dtype=np.int_)
# Check for fixed variables
fixed_vars = np.isclose(var_lower, var_upper, 0,
settings.eps_zero).nonzero()[0]
self.fixed_vars = [(i, var_lower[i]) for i in fixed_vars]
# Adjust basic data if there are fixed vars
if (len(fixed_vars)):
unfixed_vars = np.setdiff1d(np.arange(dimension), fixed_vars,
assume_unique=True)
var_lower = var_lower[unfixed_vars]
var_upper = var_upper[unfixed_vars]
temp_array = np.array([False] * dimension)
temp_array[integer_vars] = True
integer_vars = np.where(temp_array[unfixed_vars] == True)[0]
dimension -= len(fixed_vars)
self.var_lower, self.var_upper = var_lower, var_upper
self.integer_vars = integer_vars
assert(len(var_lower) == dimension)
assert(len(var_upper) == dimension)
assert((len(integer_vars) == 0) or (max(integer_vars) < dimension))
assert(init_node_pos is None or init_node_val is None or
len(init_node_pos) == len(init_node_val))
assert(do_init_strategy or (init_node_pos is not None and
len(init_node_pos) >= dimension + 1))
# Set the value of 'auto' parameters if necessary
l_settings = settings.set_auto_parameters(dimension, var_lower,
var_upper, integer_vars)
# Save references
self.l_settings = l_settings
# Local and global RBF models are initially the same
self.best_local_rbf = (l_settings.rbf, l_settings.rbf_shape_parameter)
self.best_global_rbf = (l_settings.rbf, l_settings.rbf_shape_parameter)
self.best_local_rbf_list = list()
self.best_global_rbf_list = list()
# We use n to denote the dimension of the problem, same notation
# of the paper. This is redundant but it simplifies our life.
self.n = dimension
# Set random seed and external loggers. We always use numpy's
# random generator. This is taken care of by an appropriate
# function (useful for parallelization).
ru.init_environment(l_settings)
# Initialize counters
self.itercount = 0
self.evalcount = 0
self.noisy_evalcount = 0
self.current_step = 0
self.cyclecount = -1
self.num_cons_refinement = 0
self.num_stalled_iter = 0
self.num_noisy_restarts = 0
self.discarded_iters = collections.deque(
maxlen=l_settings.discarded_window_size*l_settings.num_cpus)
self.do_init_strategy = do_init_strategy
# Initialize identifiers of the search steps
self.inf_step = 0
self.local_search_step = (l_settings.num_global_searches + 1)
self.cycle_length = (l_settings.num_global_searches + 2)
self.refinement_step = (l_settings.num_global_searches + 2)
self.restoration_step = (l_settings.num_global_searches + 3)
# Determine which step is the first of each loop
self.first_step = (self.inf_step if l_settings.do_infstep
else self.inf_step + 1)
# Initialize settings for two-phase optimization.
if (self.bb.has_evaluate_noisy()):
self.two_phase_optimization = True
self.is_fmin_noisy = True
self.is_fbest_noisy = True
self.eval_mode = 'noisy'
else:
self.two_phase_optimization = False
self.is_fmin_noisy = False
self.is_fbest_noisy = False
self.eval_mode = 'accurate'
# Round variable bounds to integer if necessary
ru.round_integer_bounds(var_lower, var_upper, integer_vars)
# Save initial nodes if given
if (init_node_pos is None):
self.init_node_pos = np.array([])
else:
self.init_node_pos = np.array(init_node_pos)
for point in self.init_node_pos:
ru.round_integer_vars(point, self.integer_vars)
if (init_node_val is None):
self.init_node_val = np.array([])
else:
self.init_node_val = np.array(init_node_val)
# Initialize global lists
self.all_node_pos, self.all_node_val = np.array([]), np.array([])
self.node_pos, self.node_val = np.array([]), np.array([])
# Store if each function evaluation is noisy or accurate
self.node_is_noisy, self.all_node_is_noisy = np.array([]), np.array([])
# Error bounds
self.node_err_bounds = np.array([])
self.all_node_err_bounds = np.array([])
# We need to remember the index of the first node in all_node_pos
# after every restart
self.num_nodes_at_restart = 0
# Update domain bounds if necessary
(self.l_lower, self.l_upper) = ru.transform_domain_bounds(
l_settings, var_lower, var_upper)
# Current minimum value among the nodes, and its index
self.fmin_index = 0
self.fmin = float('+inf')
# Current maximum value among the nodes
self.fmax = float('+inf')
# Current best value found
self.fbest_index = 0
self.fbest = float('+inf')
# Best gap shown on the log so far
self.best_gap_shown = float('+inf')
# Best function value at the beginning of an optimization cycle
self.fmin_stall_check = self.fmin
# Best function value and itercount at the most recent refinement
self.fmin_last_refine = self.fmin
self.iter_last_refine = 0
# Trust region information
self.tr_model_set = np.array([])
self.tr_radius = float('inf')
self.tr_iterate_index = None
# Set default output stream
self.output_stream = sys.stdout
# Update timer
self.elapsed_time += time.time() - start_time
# -- end function
[docs] def set_output_stream(self, output_stream):
"""Set output stream for the log.
Parameters
----------
output_stream : file
Stream to be used for output. Must have a 'write' and a
'flush' method.
"""
assert('write' in dir(output_stream))
assert('flush' in dir(output_stream))
self.output_stream = output_stream
# -- end function
[docs] def print_init_line(self):
"""Print first line of the output, with headers.
"""
print(' {:>5s}'.format('Iter') +
' {:>6s}'.format('Cycle') +
' {:15s}'.format('Action') +
' {:>18s}'.format('Objective value') +
' {:>8s}'.format('Time') +
' {:>8s}'.format('Gap'),
file=self.output_stream)
print(' {:>5s}'.format('----') +
' {:>6s}'.format('-----') +
' {:15s}'.format('------') +
' {:>18s}'.format('---------------') +
' {:>8s}'.format('----') +
' {:>8s}'.format('---'),
file=self.output_stream)
self.output_stream.flush()
# -- end function
[docs] def update_log(self, tag, node_is_noisy=None, obj_value=None, gap=None):
"""Print a single line in the log.
Update the program's log, writing information about an
iteration of the optimization algorithm, or a special message.
Parameters
----------
tag : string
Iteration id tag, or unique message if at least one of the
other arguments are None.
node_is_noisy : bool or None
Is the objective function value to be printed associated
with a node evaluated in noisy mode?
obj_value : float or None
Objective function value to print.
gap : float or None
Relative distance from the optimum. This will be
multiplied by 100 before printing.
"""
if (node_is_noisy is None or obj_value is None or gap is None):
print(' {:5d}'.format(self.itercount) +
' {:6d}'.format(self.cyclecount) +
' {:35s}'.format(tag) +
' {:8.2f}'.format(time.time() - self.start_time),
file=self.output_stream)
else:
print(' {:5d}'.format(self.itercount) +
' {:6d}'.format(self.cyclecount) +
' {:15s}'.format(tag) +
' {:18.6f}'.format(obj_value) +
'{:s}'.format('~' if node_is_noisy else ' ') +
' {:8.2f}'.format(time.time() - self.start_time) +
' {:8.2f}'.format(gap*100) +
' {:s}'.format('*' if (gap < self.best_gap_shown and
obj_value == self.fbest) else ' '),
file=self.output_stream)
if (gap < self.best_gap_shown and obj_value == self.fbest):
self.best_gap_shown = gap
self.output_stream.flush()
# -- end function
[docs] def print_summary_line(self, node_is_noisy, gap):
"""Print summary line of the algorithm.
Parameters
----------
node_is_noisy : bool
Is the objective function value to be printed associated
with a node evaluated in noisy mode?
gap : float
Relative distance from the optimum. This will be
multiplied by 100 before printing.
"""
print('Summary: iters {:3d}'.format(self.itercount) +
' evals {:3d}'.format(self.evalcount) +
' noisy_evals {:3d}'.format(self.noisy_evalcount) +
' cycles {:3d}'.format(self.cyclecount) +
' opt_time {:7.2f}'.format(time.time() - self.start_time) +
' tot_time {:7.2f}'.format(self.elapsed_time) +
' obj{:s}'.format('~' if node_is_noisy else ' ') +
' {:15.6f}'.format(self.fbest) +
' gap {:6.2f}'.format(100*gap),
file=self.output_stream)
self.output_stream.flush()
# -- end function
[docs] def remove_node(self, index, all_node_shift=0):
"""Remove a node from the lists of interpolation nodes.
Given the index of a node, remove its references from all
relevant places.
Parameters
----------
index : int
Index of the node to be removed, in the list self.node_pos
all_node_shift : int
A shift that has to be applied to the index to find the
corresponding node in self.all_node_pos. Typically, this
is the size of self.all_node_pos at the latest restart.
"""
assert(0 <= index <= len(self.node_pos))
assert(0 <= index + all_node_shift <= len(self.all_node_pos))
self.node_pos = np.delete(np.atleast_2d(self.node_pos), index, axis=0)
self.node_val = np.delete(self.node_val, index)
self.node_is_noisy = np.delete(self.node_is_noisy, index)
self.all_node_pos = np.delete(np.atleast_2d(self.all_node_pos),
all_node_shift + index, axis=0)
self.all_node_val = np.delete(self.all_node_val,
all_node_shift + index)
self.all_node_is_noisy = np.delete(self.all_node_is_noisy,
all_node_shift + index)
# -- end function
[docs] def add_node(self, point, orig_point, value):
"""Add a node to the all relevant data structures.
Given the data corresponding to a node, add it to all relevant
places: the list of current nodes, the list of all
nodes. Also, update function minimum and maximum.
Parameters
----------
point : 1D numpy.ndarray[float]
Coordinates of the node.
orig_point : 1D numpy.ndarray[float]
The point coordinates in the original space.
value : float
Objective function value of the node
"""
if not self.node_pos.size:
self.node_pos = point.copy()
self.node_err_bounds = np.array([[0, 0]])
else:
self.node_pos = np.vstack((self.node_pos, point))
self.node_err_bounds = np.vstack((self.node_err_bounds,
np.array([[0, 0]])))
self.node_val = np.append(self.node_val, value)
self.node_is_noisy = np.append(self.node_is_noisy, False)
if not self.all_node_pos.size:
self.all_node_pos = orig_point.copy()
self.all_node_err_bounds = np.array([[0, 0]])
else:
self.all_node_pos = np.vstack((self.all_node_pos, orig_point))
self.all_node_err_bounds = np.vstack((self.all_node_err_bounds,
np.array([[0, 0]])))
self.all_node_val = np.append(self.all_node_val, value)
self.all_node_is_noisy = np.append(self.all_node_is_noisy, False)
# Update fmin and fmax
self.fmax = max(self.fmax, value)
if (value < self.fmin):
self.fmin_index = len(self.node_pos) - 1
self.fmin = value
self.is_fmin_noisy = False
if (value < self.fbest):
self.fbest_index = len(self.all_node_pos) - 1
self.fbest = value
self.is_fbest_noisy = False
# -- end function
[docs] def add_noisy_node(self, point, orig_point, value, err_l, err_u):
"""Add a noisy node to the all relevant data structures.
Given the data corresponding to a node, add it to all relevant
places: the list of current nodes, the list of all
nodes. Also, update function minimum and maximum.
Parameters
----------
point : 1D numpy.ndarray[float]
Coordinates of the node.
orig_point : 1D numpy.ndarray[float]
The point coordinates in the original space.
value : float
Objective function value of the node
err_l : float
Lower variation for the error interval of the node.
err_u : float
Upper variation for the error interval of the node
"""
if not self.node_pos.size:
self.node_pos = point.copy()
self.node_err_bounds = np.array([[err_l, err_u]])
else:
self.node_pos = np.vstack((self.node_pos, point))
self.node_err_bounds = np.vstack((self.node_err_bounds,
np.array([[err_l, err_u]])))
self.node_val = np.append(self.node_val, value)
self.node_is_noisy = np.append(self.node_is_noisy, True)
if not self.all_node_pos.size:
self.all_node_pos = orig_point.copy()
self.all_node_err_bounds = np.array([[err_l, err_u]])
else:
self.all_node_pos = np.vstack((self.all_node_pos, orig_point))
self.all_node_err_bounds = np.vstack((self.all_node_err_bounds,
np.array([err_l, err_u])))
self.all_node_val = np.append(self.all_node_val, value)
self.all_node_is_noisy = np.append(self.all_node_is_noisy, True)
# Update fmin and fmax
self.fmax = max(self.fmax, value)
if (value < self.fmin):
self.fmin_index = len(self.node_pos) - 1
self.fmin = value
self.is_fmin_noisy = True
if (value < self.fbest):
self.fbest_index = len(self.all_node_pos) - 1
self.fbest = value
self.is_fbest_noisy = True
# -- end function
[docs] def save_to_file(self, filename):
"""Save object on file, with its state.
Saves the current state of the algorithm on file. The
optimization can be subsequently resumed reading the state
from file. This function will also attempt to save the state
of the random number generators so that if resumed on the same
machine, the optimization process is identical to an
uninterrupted process.
Parameters
----------
filename : string
Name of the file that the state will be saved to.
"""
# Save state of RNG
self.random_state = np.random.get_state()
# We cannot pickle this attribute. Erase it.
output_stream = self.output_stream
self.output_stream = None
# Dump to file
with open(filename, 'wb') as pickle_file:
pickle.dump(self, pickle_file, pickle.HIGHEST_PROTOCOL)
# Restore erased attribute
self.output_stream = output_stream
# -- end function
[docs] @classmethod
def load_from_file(cls, filename):
"""Load object from file, with its state.
Read the current state from file, and return an object of this
class. The optimization can be resumed immediately. This
function will attempt to set the random number generators to
the state they were in. Note that the output stream is set to
stdout, regardless of the output stream when the state was
saved, so the caller may have to set the desired output
stream.
Parameters
----------
filename : string
Name of the file from which the state will be read.
Returns
-------
RbfoptAlgorithm
An object of this class.
"""
assert(os.path.isfile(filename))
with open(filename, 'rb') as pickle_file:
alg = pickle.load(pickle_file)
np.random.set_state(alg.random_state)
# Set default output stream
alg.output_stream = sys.stdout
return alg
# -- end function
[docs] def optimize(self, pause_after_iters=sys.maxsize):
"""Optimize a black-box function.
Optimize an unknown function over a box using an RBF-based
algorithm. This function will select the serial or parallel
version of the optimizer, depending on settings.
Parameters
----------
pause_after_iters : int
Number of iterations after which the optimization process
should pause. This allows the user to do other activities
and resume optimization at a later time. Default
sys.maxsize, which is larger than any practical integer.
Returns
-------
(float, 1D numpy.ndarray[float], int, int, int)
A quintuple (value, point, itercount, evalcount,
noisy_evalcount) containing the objective function value of
the best solution found, the corresponding value of the
decision variables, the number of iterations of the
algorithm, the total number of function evaluations, and
the number of these evaluations that were performed in
'noisy' mode.
"""
self.start_time = time.time()
self.print_init_line()
if (self.n == 0 and len(self.fixed_vars[0])):
# There is nothing to do in this case
self.update_log('All variables fixed, nothing to do!')
self.add_node(self.var_lower, self.var_lower,
objfun([self.bb, self.var_lower, self.fixed_vars]))
self.evalcount += 1
elif (self.l_settings.num_cpus == 1):
self.optimize_serial(pause_after_iters)
else:
self.optimize_parallel(pause_after_iters)
start_time_retrieve_min = time.time()
# Find best point and return it
i = self.all_node_val.argmin()
gap = ru.compute_gap(self.l_settings, self.fbest +
self.all_node_err_bounds[i, 1])
# Update timer
self.elapsed_time += time.time() - start_time_retrieve_min
# Print summary and return
self.print_summary_line(self.all_node_is_noisy[i], gap)
return (self.all_node_val[i], self.all_node_pos[i],
self.itercount, self.evalcount, self.noisy_evalcount)
[docs] def optimize_serial(self, pause_after_iters=sys.maxsize):
"""Optimize a black-box function. Serial engine.
Optimize an unknown function over a box using an RBF-based
algorithm. This is the serial version of the optimization
routine.
Parameters
----------
pause_after_iters : int
Number of iterations after which the optimization process
should pause. Default sys.maxsize.
"""
start_time = time.time()
# Localize variables that will be used often and that will
# never be erased through the algorithm. These are all arrays
# or objects, so the reference points to the original.
var_lower, var_upper = self.var_lower, self.var_upper
l_lower, l_upper = self.l_lower, self.l_upper
integer_vars = self.integer_vars
settings, l_settings = self.settings, self.l_settings
# The dimension will not be changed so it is safe to localize
n = self.n
if (l_settings.algorithm == 'MSRSM'):
Amatinv = None
# Save number of iterations at start
itercount_at_start = self.itercount
# If this is the first iteration, initialize the algorithm
if (self.itercount == 0):
self.restart()
# We need to update the gap
gap = ru.compute_gap(l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
# Main loop
while (self.itercount - itercount_at_start < pause_after_iters and
self.itercount < l_settings.max_iterations and
self.evalcount < l_settings.max_evaluations and
self.cyclecount < l_settings.max_cycles and
time.time() - start_time < l_settings.max_clock_time and
gap > l_settings.eps_opt):
# If the user wants to skip inf_step, we proceed to the
# next iteration.
if (self.current_step == self.inf_step and
not l_settings.do_infstep):
self.advance_step_counter()
continue
# Check if we should restart.
if (self.current_step <= self.first_step and
((len(self.discarded_iters) >=
l_settings.discarded_window_size/2 and
self.discarded_iters.count(True) >=
l_settings.max_fraction_discarded *
len(self.discarded_iters)) or
(self.num_stalled_iter >= l_settings.max_stalled_iterations
and self.evalcount + 3*(n + 1 + self.cycle_length) <
l_settings.max_evaluations and
self.cyclecount < l_settings.max_cycles - 1))):
self.update_log('Restart')
self.restart()
# Number of nodes at current iteration
k = len(self.node_pos)
# If function scaling is automatic, determine which one to use
if (settings.function_scaling == 'auto' and
self.current_step <= self.first_step):
sorted_node_val = np.sort(self.node_val)
if (sorted_node_val[len(sorted_node_val)//2] -
sorted_node_val[0] > l_settings.log_scaling_threshold):
l_settings.function_scaling = 'log'
else:
l_settings.function_scaling = 'off'
# Rescale nodes if necessary
tfv = ru.transform_function_values(l_settings, self.node_val,
self.fmin, self.fmax,
self.node_err_bounds)
(scaled_node_val, scaled_fmin, scaled_fmax,
scaled_err_bounds, rescale_function) = tfv
# If RBF selection is automatic, at the beginning of each
# cycle check if a different RBF yields a better model
if (settings.rbf == 'auto' and k > n+2 and
self.current_step <= self.first_step and
len(self.best_local_rbf_list) <
l_settings.max_cross_validations):
loc_iter = int(math.ceil(k*0.1))
glob_iter = int(math.ceil(k*0.7))
best_local_rbf = ru.get_best_rbf_model(
l_settings, n, k, self.node_pos,
scaled_node_val, loc_iter)
best_global_rbf = ru.get_best_rbf_model(
l_settings, n, k, self.node_pos,
scaled_node_val, glob_iter)
self.best_local_rbf = best_local_rbf
self.best_global_rbf = best_global_rbf
self.best_local_rbf_list.append(best_local_rbf)
self.best_global_rbf_list.append(best_global_rbf)
# If we reached the maximum number of cross
# validations, find the best model by looking at
# frequency
if (len(self.best_local_rbf_list) ==
l_settings.max_cross_validations):
best_local_rbf = ru.get_most_common_element(
self.best_local_rbf_list, ['failed'])
self.best_local_rbf = (best_local_rbf
if best_local_rbf is not None
else ('cubic', 0.1))
best_global_rbf = ru.get_most_common_element(
self.best_global_rbf_list, ['failed'])
self.best_global_rbf = (best_global_rbf
if best_global_rbf is not None
else ('cubic', 0.1))
# If we are in local search or just before local search, use a
# local model.
if (self.current_step >= (self.local_search_step - 1)):
l_settings.rbf = self.best_local_rbf[0]
l_settings.rbf_shape_parameter = self.best_local_rbf[1]
# Otherwise, global.
else:
l_settings.rbf = self.best_global_rbf[0]
l_settings.rbf_shape_parameter = self.best_global_rbf[1]
if (self.current_step <= self.local_search_step):
try:
# Compute the matrices necessary for the algorithm
Amat = ru.get_rbf_matrix(l_settings, n, k, self.node_pos)
if (l_settings.algorithm == 'Gutmann'):
Amatinv = ru.get_matrix_inverse(l_settings, Amat)
# Compute RBF interpolant at current stage
if (np.any(self.node_is_noisy)):
# Get coefficients for the exact RBF
rbf_l, rbf_h = ru.get_rbf_coefficients(
l_settings, n, k, Amat, scaled_node_val)
# RBF with some noisy function evaluations
rbf_l, rbf_h = aux.get_noisy_rbf_coefficients(
l_settings, n, k, Amat[:k, :k], Amat[:k, k:],
scaled_node_val, scaled_err_bounds, rbf_l, rbf_h)
else:
# Fully accurate RBF
rbf_l, rbf_h = ru.get_rbf_coefficients(
l_settings, n, k, Amat, scaled_node_val)
except np.linalg.LinAlgError:
# Record statistics on failed model
if (self.current_step >= (self.local_search_step - 1)):
self.best_local_rbf_list[-1] = 'failed'
else:
self.best_global_rbf_list[-1] = 'failed'
# Error in the solution of the linear system. We must
# switch to a restoration phase.
self.current_step = self.restoration_step
elif (self.current_step == self.refinement_step):
if (self.num_cons_refinement == 0):
# For the first refinement step iteration, some data
# has to be computed
self.tr_model_set, self.tr_radius = ref.init_trust_region(
l_settings, n, k, self.node_pos,
self.node_pos[self.fmin_index])
self.tr_iterate_index = self.fmin_index
# Compute linear model for trust region method
try:
h, b, rank_deficient = ref.get_linear_model(
l_settings, n, k, self.node_pos,
scaled_node_val, self.tr_model_set)
except np.linalg.LinAlgError:
# Error in the solution of the linear system. We
# go back to the usual search.
self.current_step = self.first_step
self.cyclecount += 1
continue
# For displaying purposes, record what type of iteration
# we are performing
iteration_id = ''
# Initialize the new point to None
next_p = None
if (self.current_step == self.inf_step):
# Infstep: explore the parameter space
next_p = pure_global_step(l_settings, n, k, l_lower,
l_upper, integer_vars,
self.node_pos, Amatinv)
iteration_id = 'InfStep'
elif (self.current_step == self.restoration_step):
# Restoration
next_p = self.restoration_search()
if (next_p is None ):
self.update_log('Restoration phase failed. Restart.')
self.update_log('Restart')
self.restart()
continue
k = len(self.node_pos)
iteration_id = 'Restoration'
elif (self.current_step == self.local_search_step):
# Local search
(adj, next_p, ind) = local_step(
l_settings, n, k, l_lower, l_upper, integer_vars,
self.node_pos, rbf_l, rbf_h, tfv[:4],
Amat, Amatinv, self.fmin_index,
self.two_phase_optimization, self.eval_mode,
self.node_is_noisy)
# Re-evaluate point if necessary
if (ind is not None):
self.remove_node(ind, self.num_nodes_at_restart)
# We must update k here to make sure it is consistent
# until the start of the next iteration.
k = len(self.node_pos)
if (adj):
iteration_id = 'AdjLocalStep'
else:
iteration_id = 'LocalStep'
elif (self.current_step == self.refinement_step):
# Refinement
next_p, model_impr, tr_to_replace = refinement_step(
l_settings, n, k, l_lower, l_upper, integer_vars,
self.node_pos, tfv[:4], h, b, rank_deficient,
self.tr_model_set, self.tr_iterate_index, self.tr_radius)
iteration_id = 'RefinementStep'
else:
# Global search
next_p = global_step(l_settings, n, k, l_lower, l_upper,
integer_vars, self.node_pos, rbf_l,
rbf_h, tfv[:4], Amatinv,
self.fmin_index, self.current_step)
iteration_id = 'GlobalStep'
# -- end if
# If the optimization failed or the point is too close to
# current nodes, discard it. Otherwise, add it to the list.
if ((next_p is None) or
(ru.get_min_distance(next_p, self.node_pos) <=
l_settings.min_dist)):
self.advance_step_counter()
self.discarded_iters.append(True)
self.num_cons_refinement = 0
self.update_log('Discarded')
else:
# Transform back to original space if necessary
next_p_orig = ru.transform_domain(l_settings, var_lower,
var_upper, next_p, True)
# Evaluate the new point, in accurate mode or noisy
# mode. If we performed a restart, we also check if the
# same node was evaluated before the restart happened,
# to make sure we do not evaluate it twice.
next_val = None
if (self.eval_mode == 'accurate'):
if (self.num_nodes_at_restart):
min_dist, i = ru.get_min_distance_and_index(
next_p_orig,
self.all_node_pos[:self.num_nodes_at_restart])
if (min_dist < l_settings.eps_zero and
not self.all_node_is_noisy[i]):
next_val = self.all_node_val[i]
if (next_val is None):
next_val = objfun([self.bb, next_p_orig,
self.fixed_vars])
self.evalcount += 1
curr_is_noisy = False
else:
if (self.num_nodes_at_restart):
min_dist, i = ru.get_min_distance_and_index(
next_p_orig,
self.all_node_pos[:self.num_nodes_at_restart])
if (min_dist < l_settings.eps_zero):
next_val = self.all_node_val[i]
curr_is_noisy = self.all_node_is_noisy[i]
if (next_val is None):
next_val, err_l, err_u = objfun_noisy(
[self.bb, next_p_orig, self.fixed_vars])
self.noisy_evalcount += 1
curr_is_noisy = True
if (curr_is_noisy and
self.require_accurate_evaluation(next_val + err_l)):
self.update_log(iteration_id, True, next_val, gap)
next_val = objfun([self.bb, next_p_orig,
self.fixed_vars])
self.evalcount += 1
curr_is_noisy = False
# Add to the lists
if (curr_is_noisy):
self.add_noisy_node(next_p, next_p_orig, next_val,
err_l, err_u)
else:
self.add_node(next_p, next_p_orig, next_val)
gap = ru.compute_gap(
l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
self.update_log(iteration_id, self.node_is_noisy[-1],
next_val, gap)
if (self.current_step == self.refinement_step):
real_impr = (scaled_node_val[self.tr_iterate_index] -
rescale_function(next_val))
self.refinement_update(model_impr, real_impr,
tr_to_replace)
self.iter_last_refine = self.itercount
elif ((self.current_step == self.local_search_step) and
(self.itercount >= self.iter_last_refine +
self.l_settings.refinement_frequency *
(self.cycle_length -
(1 if self.l_settings.do_infstep else 2)) or
self.unlimited_refinement_active())):
self.current_step = self.refinement_step
self.iter_last_refine = self.itercount
else:
self.advance_step_counter()
self.discarded_iters.append(False)
# Update iteration number
self.itercount += 1
# At the beginning of each loop of the cyclic optimization
# strategy, check if the main loop is stalling
self.stalling_update()
# Check if we should switch to the second phase of
# two-phase optimization.
self.phase_update()
# Check if we should save the state.
if (self.itercount % l_settings.save_state_interval == 0):
self.save_to_file(l_settings.save_state_file)
# -- end while
# Update timer
self.elapsed_time += time.time() - start_time
# -- end function
[docs] def optimize_parallel(self, pause_after_iters=sys.maxsize):
"""Optimize a black-box function using parallel evaluations.
Optimize an unknown function over a box using an RBF-based
algorithm, using as many CPUs as requested.
Parameters
----------
pause_after_iters : int
Number of iterations after which the optimization process
should pause. This allows the user to do other activities
and resume optimization at a later time. Default
sys.maxsize, which is larger than any practical integer.
"""
start_time = time.time()
# Localize variables that will be used often and that will
# never be erased through the algorithm. These are all arrays
# or objects, so the reference points to the original.
var_lower, var_upper = self.var_lower, self.var_upper
l_lower, l_upper = self.l_lower, self.l_upper
integer_vars = self.integer_vars
settings, l_settings = self.settings, self.l_settings
# The dimension will not be changed so it is safe to localize
n = self.n
# Create pool of workers
pool = Pool(l_settings.num_cpus, ru.init_environment,
(l_settings, ))
# List of new point evaluations. A new point evaluation has
# the format: [result, point, is_node_noisy, iteration_id],
# where result is an object of class AsyncResult, and point is
# in the transformed space.
res_eval = list()
# List of new points to be explored. A new search request has
# the format: [result, is_node_noisy, iteration_id], where
# result is an object of class AsyncResult.
res_search = list()
# List of point evaluations removed from the computation
# because numerically unstable. Same format as res_eval above.
res_removed = list()
# Save number of iterations at start
itercount_at_start = self.itercount
# We will keep here temporary nodes submitted for
# evaluation. A position will be none if unfilled.
temp_node_pos = np.array([])
temp_node_val = np.array([])
temp_node_is_noisy = np.array([])
if (l_settings.algorithm == 'MSRSM'):
Amatinv = None
# Number of iterations spent in refinement. We used this
# strange initialization statement to make sure that the
# number is consistent (if not accurate) in case the algorithm
# is stopped and restarted.
refinement_itercount = self.itercount/l_settings.refinement_frequency
# Is a refinement iteration in progress?
refinement_in_progress = False
# If this is the first iteration, initialize the algorithm
if (self.itercount == 0):
self.restart(pool=pool)
# We need to update the gap
gap = ru.compute_gap(l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
# Main loop
while (self.itercount - itercount_at_start < pause_after_iters and
self.itercount < l_settings.max_iterations and
self.evalcount < l_settings.max_evaluations and
self.cyclecount < l_settings.max_cycles and
time.time() - start_time < l_settings.max_clock_time and
gap > l_settings.eps_opt):
itercount_before_rescheck = self.itercount
# If some evaluations are complete, update data structures
while (ru.results_ready(res_eval)):
# Obtain one point evaluation that is ready
j = ru.get_one_ready_index(res_eval)
res, next_p, node_is_noisy, iid = res_eval.pop(j)
if (node_is_noisy):
next_val, err_l, err_u = res.get()
else:
next_val = res.get()
# If distance is too small for whatever reason, better
# to skip this point
if (ru.get_min_distance(next_p, self.node_pos) <=
l_settings.min_dist):
continue
# Transform back to original space if necessary
next_p_orig = ru.transform_domain(l_settings, var_lower,
var_upper, next_p, True)
if (node_is_noisy and
self.require_accurate_evaluation(next_val + err_l)):
# Evaluate again
node_is_noisy = False
new_res = pool.apply_async(
objfun,
([self.bb, next_p_orig, self.fixed_vars], ))
self.evalcount += 1
res_eval.append([new_res, next_p, node_is_noisy, iid])
# Update temporary node values. First, remove old node.
temp_node_pos = np.delete(np.atleast_2d(temp_node_pos),
j, axis=0)
temp_node_val = np.delete(temp_node_val, j)
temp_node_is_noisy = np.delete(temp_node_is_noisy, j)
# Then add again at the end of the array
if (temp_node_pos.size):
temp_node_pos = np.vstack((temp_node_pos, next_p))
else:
temp_node_pos = np.atleast_2d(next_p)
temp_node_val = np.append(
temp_node_val, np.clip(next_val, self.fmin, self.fmax))
temp_node_is_noisy = np.append(temp_node_is_noisy,
node_is_noisy)
else:
# Add to data structures.
if (node_is_noisy):
self.add_noisy_node(next_p, next_p_orig, next_val,
err_l, err_u)
else:
self.add_node(next_p, next_p_orig, next_val)
gap = ru.compute_gap(
l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
self.update_log(iid, node_is_noisy, next_val, gap)
# Perform refinement updates if necessary
if (iid == 'RefinementStep'):
real_impr = (ref_rescale_function(
node_val[self.tr_iterate_index]) -
ref_rescale_function(next_val))
self.refinement_update_parallel(model_impr, real_impr,
tr_to_replace)
refinement_itercount += 1
refinement_in_progress = False
# Remove from list of temporary points
temp_node_pos = np.delete(np.atleast_2d(temp_node_pos),
j, axis=0)
temp_node_val = np.delete(temp_node_val, j)
temp_node_is_noisy = np.delete(temp_node_is_noisy, j)
# Perform main updates
self.itercount += 1
self.stalling_update()
self.phase_update()
# Check if we should save the state.
if (self.itercount % l_settings.save_state_interval == 0):
self.save_to_file(l_settings.save_state_file)
# -- end while
if (itercount_before_rescheck < self.itercount):
# If some points have been evaluated, force another
# loop to check termination conditions
continue
# At this point, the model could be updated. Are there
# points that should be submitted for evaluation?
while (ru.results_ready(res_search)):
# Obtain one search point that is ready
j = ru.get_one_ready_index(res_search)
(res, node_is_noisy, iteration_id) = res_search.pop(j)
# Local search is treated differently, because it
# may require re-evaluation of previous points
if (iteration_id == 'LocalStep'):
adj, next_p, ind = res.get()
# Re-evaluate point if necessary
if (ind is not None):
# Because of parallelism, there is a chance
# that the node position changed. Make sure
# we have the right one.
if (not np.array_equal(self.node_pos[ind], next_p)):
ind = np.where(np.all(
self.node_pos == next_p, axis=1))[0][0]
self.remove_node(ind, self.num_nodes_at_restart)
# Since node_pos is a copy of self.node_pos
# and temp_node_pos, there is no need to
# adjust k
if (adj):
iteration_id = 'AdjLocalStep'
else:
iteration_id = 'LocalStep'
elif (iteration_id == 'RefinementStep'):
next_p, model_impr, tr_to_replace = res.get()
else:
next_p = res.get()
# Verify that we have an actual point available
if ((next_p is None) or
(ru.get_min_distance(next_p, self.node_pos) <=
l_settings.min_dist) or
(temp_node_pos.size and
(ru.get_min_distance(next_p, temp_node_pos) <=
l_settings.min_dist))):
self.discarded_iters.append(True)
if (iteration_id == 'RefinementStep'):
refinement_itercount += 1
refinement_in_progress = False
self.num_cons_refinement = 0
# Update iteration number
self.itercount += 1
self.update_log('Discarded')
else:
# Transform back to original space if necessary
next_p_orig = ru.transform_domain(
l_settings, var_lower, var_upper, next_p, True)
# If we performed a restart, we also check if
# the same node was evaluated before the
# restart happened, to make sure we do not
# evaluate it twice.
next_val = None
if (self.num_nodes_at_restart):
min_dist, i = ru.get_min_distance_and_index(
next_p_orig,
self.all_node_pos[:self.num_nodes_at_restart])
if (min_dist < l_settings.eps_zero and
(node_is_noisy or not self.all_node_is_noisy[i])):
next_val = self.all_node_val[i]
node_is_noisy = self.all_node_is_noisy[i]
# Add to the data structures.
if (node_is_noisy):
self.add_noisy_node(
next_p, next_p_orig, next_val,
err_l, err_u)
else:
self.add_node(next_p, next_p_orig, next_val)
gap = ru.compute_gap(
l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
self.update_log(iid, node_is_noisy, next_val, gap)
# Update iteration number
self.itercount += 1
# Check if we should save the state.
if (self.itercount %
l_settings.save_state_interval == 0):
self.save_to_file(l_settings.save_state_file)
# Perform main updates
self.stalling_update()
self.phase_update()
if (next_val is None):
if (node_is_noisy):
new_res = pool.apply_async(
objfun_noisy,
([self.bb, next_p_orig, self.fixed_vars], ))
self.noisy_evalcount += 1
else:
new_res = pool.apply_async(
objfun,
([self.bb, next_p_orig, self.fixed_vars], ))
self.evalcount += 1
res_eval.append([new_res, next_p, node_is_noisy,
iteration_id])
self.discarded_iters.append(False)
# Append point to the list of temporary nodes.
# When evaluating the surrogate model, use the
# RbfoptSettings object that was used when
# rbf_l and rbf_h were last computed.
# Furthermore, ensure we only pick
# nodes used to compute coefficients.
val = ru.evaluate_rbf(
rbf_coeff_settings, next_p, n, len(rbf_l),
node_pos[:len(rbf_l)], rbf_l, rbf_h)
if (temp_node_pos.size):
temp_node_pos = np.vstack((temp_node_pos, next_p))
else:
temp_node_pos = np.atleast_2d(next_p)
temp_node_val = np.append(
temp_node_val, np.clip(val, self.fmin, self.fmax))
temp_node_is_noisy = np.append(temp_node_is_noisy,
node_is_noisy)
# -- end while
# If no CPUs are available, wait a bit and try again.
if (len(res_eval) + len(res_search) >= l_settings.num_cpus):
time.sleep(l_settings.parallel_wakeup_time)
continue
# At this point, we know that there are free workers, so
# we should try to generate new search points.
# If the user wants to skip inf_step, we proceed to the
# next iteration.
if (self.current_step == self.inf_step and
not l_settings.do_infstep):
self.advance_step_counter()
continue
# Check if we should restart.
if (self.current_step <= self.first_step and
((len(self.discarded_iters) >=
l_settings.discarded_window_size*l_settings.num_cpus/2
and self.discarded_iters.count(True) >=
l_settings.max_fraction_discarded *
len(self.discarded_iters)) or
(self.num_stalled_iter >= l_settings.max_stalled_iterations
and self.evalcount + 3*(n + 1 + self.cycle_length) <
l_settings.max_evaluations and
self.cyclecount < l_settings.max_cycles - 1))):
# If there are still results in the pipeline, wait for
# them, then try again.
if (res_eval or res_search):
[result[0].wait() for result in res_eval]
[result[0].wait() for result in res_search]
continue
# If there were no results in the pipeline, we should
# restart.
self.update_log('Restart')
self.restart(pool=pool)
temp_node_pos = np.array([])
temp_node_val = np.array([])
temp_node_is_noisy = np.array([])
# Nodes at current iteration, including temporary ones
if (temp_node_pos.size):
node_pos = np.vstack((self.node_pos, temp_node_pos))
node_err_bounds = np.vstack(
(self.node_err_bounds, np.zeros((len(temp_node_pos), 2))))
else:
node_pos = np.array(self.node_pos)
node_err_bounds = np.array(self.node_err_bounds)
node_val = np.append(self.node_val, temp_node_val)
node_is_noisy = np.append(self.node_is_noisy, temp_node_is_noisy)
k = len(node_pos)
# Make a copy as there might be parallel access
l_settings = copy.deepcopy(l_settings)
# If function scaling is automatic, determine which one to use
if (settings.function_scaling == 'auto' and
self.current_step <= self.first_step):
sorted_node_val = np.sort(node_val)
if (sorted_node_val[len(sorted_node_val)//2] -
sorted_node_val[0] > l_settings.log_scaling_threshold):
l_settings.function_scaling = 'log'
else:
l_settings.function_scaling = 'off'
# Rescale nodes if necessary
tfv = ru.transform_function_values(l_settings, node_val,
self.fmin, self.fmax,
node_err_bounds)
(scaled_node_val, scaled_fmin, scaled_fmax,
scaled_err_bounds, rescale_function) = tfv
# If RBF selection is automatic, at the beginning of each
# cycle check if a different RBF yields a better model
if (settings.rbf == 'auto' and k > n+2 and
self.current_step <= self.first_step and
len(self.best_local_rbf_list) <
l_settings.max_cross_validations):
loc_iter = int(math.ceil(k*0.1))
glob_iter = int(math.ceil(k*0.7))
best_local_rbf = ru.get_best_rbf_model(
l_settings, n, k, node_pos, scaled_node_val, loc_iter)
best_global_rbf = ru.get_best_rbf_model(
l_settings, n, k, node_pos, scaled_node_val, glob_iter)
self.best_local_rbf = best_local_rbf
self.best_global_rbf = best_global_rbf
self.best_local_rbf_list.append(best_local_rbf)
self.best_global_rbf_list.append(best_global_rbf)
# If we reached the maximum number of cross
# validations, find the best model by looking at
# frequency
if (len(self.best_local_rbf_list) ==
l_settings.max_cross_validations):
best_local_rbf = ru.get_most_common_element(
self.best_local_rbf_list, ['failed'])
self.best_local_rbf = (best_local_rbf
if best_local_rbf is not None
else ('cubic', 0.1))
best_global_rbf = ru.get_most_common_element(
self.best_global_rbf_list, ['failed'])
self.best_global_rbf = (best_global_rbf
if best_global_rbf is not None
else ('cubic', 0.1))
# If we are in local search or just before local search, use a
# local model.
if (self.current_step >= (self.local_search_step - 1)):
l_settings.rbf = self.best_local_rbf[0]
l_settings.rbf_shape_parameter = self.best_local_rbf[1]
# Otherwise, global.
else:
l_settings.rbf = self.best_global_rbf[0]
l_settings.rbf_shape_parameter = self.best_global_rbf[1]
try:
# Compute the matrices necessary for the algorithm
Amat = ru.get_rbf_matrix(l_settings, n, k, node_pos)
if (l_settings.algorithm == 'Gutmann'):
Amatinv = ru.get_matrix_inverse(l_settings, Amat)
# Compute RBF interpolant at current stage
if (np.any(self.node_is_noisy)):
# Get coefficients for the exact RBF
rbf_l, rbf_h = ru.get_rbf_coefficients(
l_settings, n, k, Amat, scaled_node_val)
# RBF with some noisy function evaluations
rbf_l, rbf_h = aux.get_noisy_rbf_coefficients(
l_settings, n, k, Amat[:k, :k], Amat[:k, k:],
scaled_node_val, scaled_err_bounds, rbf_l, rbf_h)
else:
# Fully accurate RBF
rbf_l, rbf_h = ru.get_rbf_coefficients(
l_settings, n, k, Amat, scaled_node_val)
# Copy in case something goes wrong
rbf_coeff_settings = copy.deepcopy(l_settings)
except np.linalg.LinAlgError:
# Record statistics on failed model
if (self.current_step >= (self.local_search_step - 1)):
self.best_local_rbf_list[-1] = 'failed'
else:
self.best_global_rbf_list[-1] = 'failed'
# Error in the solution of the linear system. We must
# switch to a restoration phase.
self.current_step = self.restoration_step
# For displaying purposes, record what type of iteration
# we are performing
iteration_id = ''
# Initialize the new point to None
next_p = None
curr_is_noisy = (self.eval_mode == 'noisy')
if (self.current_step == self.restoration_step):
# If there are still results in the pipeline, wait for
# them, then try again.
if (res_eval or res_search):
[result[0].wait() for result in res_eval]
[result[0].wait() for result in res_search]
continue
# If we failed, we restart because in parallel mode
# many things can go wrong when deleting points
self.update_log('Restoration phase failed. Restart.')
self.update_log('Restart')
self.restart(pool=pool)
temp_node_pos = np.array([])
temp_node_val = np.array([])
temp_node_is_noisy = np.array([])
continue
# If no refinement is in progress, and the frequency is
# not excessive, start it or keep it going
elif (not refinement_in_progress and
refinement_itercount <= self.itercount /
l_settings.refinement_frequency):
if (self.num_cons_refinement == 0 or
(self.fmin <= node_val[self.tr_iterate_index] -
l_settings.eps_impr *
max(1.0, abs(node_val[self.tr_iterate_index])))):
# For the first refinement step iteration, some data
# has to be computed
self.tr_model_set, self.tr_radius = ref.init_trust_region(
l_settings, n, k, node_pos, node_pos[self.fmin_index])
self.tr_iterate_index = self.fmin_index
try:
# Compute linear model for trust region method
h, b, rank_deficient = ref.get_linear_model(
l_settings, n, k, node_pos,
scaled_node_val, self.tr_model_set)
new_res = pool.apply_async(
refinement_step,
(l_settings, n, k, l_lower, l_upper, integer_vars,
node_pos, tfv[:4], h, b, rank_deficient,
self.tr_model_set, self.tr_iterate_index,
self.tr_radius))
iteration_id = 'RefinementStep'
ref_rescale_function = rescale_function
res_search.append([new_res, curr_is_noisy, iteration_id])
refinement_in_progress = True
except np.linalg.LinAlgError:
# Error in the solution of the linear system. We
# ignore and continue the search.
pass
elif (self.current_step == self.inf_step):
# Infstep: explore the parameter space
new_res = pool.apply_async(
pure_global_step,
(l_settings, n, k, l_lower, l_upper, integer_vars,
node_pos, Amatinv))
iteration_id = 'InfStep'
res_search.append([new_res, curr_is_noisy, iteration_id])
elif (self.current_step == self.local_search_step):
# Local search
new_res = pool.apply_async(
local_step,
(l_settings, n, k, l_lower, l_upper, integer_vars,
node_pos, rbf_l, rbf_h, tfv[:4], Amat, Amatinv,
self.fmin_index, self.two_phase_optimization,
self.eval_mode, node_is_noisy))
iteration_id = 'LocalStep'
res_search.append([new_res, curr_is_noisy, iteration_id])
else:
# Global search
new_res = pool.apply_async(
global_step,
(l_settings, n, k, l_lower, l_upper, integer_vars,
node_pos, rbf_l, rbf_h, tfv[:4], Amatinv,
self.fmin_index, self.current_step))
iteration_id = 'GlobalStep'
res_search.append([new_res, curr_is_noisy, iteration_id])
# -- end if
# Move forward, without waiting for results
if (iteration_id != 'RefinementStep'):
self.advance_step_counter()
# -- end while
# Wait for all results to complete
pool.close()
pool.join()
# Obtain all point evaluations that are ready
for (res, next_p, node_is_noisy, iid) in res_eval + res_removed:
next_val = res.get()
min_dist = ru.get_min_distance(next_p, self.node_pos)
# Transform back to original space if necessary
next_p_orig = ru.transform_domain(l_settings, var_lower,
var_upper, next_p, True)
# Add to the lists.
if (node_is_noisy):
self.add_noisy_node(next_p, next_p_orig, next_val,
err_l, err_u)
else:
self.add_node(next_p, next_p_orig, next_val)
gap = ru.compute_gap(
l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
self.update_log(iid, self.node_is_noisy[-1], next_val, gap)
# Update iteration number
self.itercount += 1
# -- end for
# Update timer
self.elapsed_time += time.time() - start_time
# -- end function
[docs] def restart(self, pool=None):
"""Perform a complete restart of the optimization.
Restart the optimization algorithm, i.e. discard the current
RBF model, and select new sample points to start the algorithm
from scratch. Previous point evaluations are ignored, but they
are still recorded in the appropriate arrays.
Parameters
----------
pool : multiprocessing.Pool()
A pool of workers to evaluate the initialization points in
parallel. If None, parallel evaluation will not be
performed.
Raises
------
RuntimeError
If the routine does not find enough good initialization points.
"""
# We update the number of noisy restarts here, so that if
# we hit the limit on noisy restarts, we can evaluate
# points in accurate mode after restarting (even if
# eval_mode is updated in a subsequent block of code)
self.num_noisy_restarts += (1 if self.eval_mode == 'noisy'
else 0)
# Store the current number of nodes
self.num_nodes_at_restart = len(self.all_node_pos)
if (self.itercount > 0 or self.do_init_strategy):
# Compute a new set of starting points
node_pos = ru.initialize_nodes(self.l_settings, self.var_lower,
self.var_upper, self.integer_vars)
# Check if some points were previously evaluated, i.e. before
# the restart. If so, remove them from the list of nodes that
# will be evaluated now.
if (self.num_nodes_at_restart):
dist = ru.bulk_get_min_distance(
node_pos, self.all_node_pos[:self.num_nodes_at_restart])
indices = np.where(dist <= self.l_settings.eps_zero)[0]
prev_node_pos = list()
prev_node_val = list()
prev_node_is_noisy = list()
prev_node_err_bounds = list()
for i in reversed(sorted(indices)):
val, j = ru.get_min_distance_and_index(
node_pos[i],
self.all_node_pos[:self.num_nodes_at_restart])
if (self.eval_mode == 'noisy' or
not self.all_node_is_noisy[j]):
prev_node_pos.append(self.all_node_pos[j])
prev_node_val.append(self.all_node_val[j])
prev_node_is_noisy.append(self.all_node_is_noisy[j])
prev_node_err_bounds.append(
self.all_node_err_bounds[j])
node_pos = np.delete(np.atleast_2d(node_pos), i,
axis=0)
# Evaluate new points
if (self.eval_mode == 'accurate' or
self.num_noisy_restarts > self.l_settings.max_noisy_restarts
or (self.noisy_evalcount + self.n + 1 >=
self.l_settings.max_noisy_evaluations)):
if (pool is None):
node_val = np.array([objfun([self.bb, point,
self.fixed_vars])
for point in node_pos])
else:
map_arg = [[self.bb, point, self.fixed_vars]
for point in node_pos]
node_val = np.array(pool.map(objfun, map_arg))
node_err_bounds = np.zeros((len(node_val), 2))
self.evalcount += len(node_val)
else:
if (pool is None):
res = np.array([objfun_noisy([self.bb, point,
self.fixed_vars])
for point in node_pos])
else:
map_arg = [[self.bb, point, self.fixed_vars]
for point in node_pos]
res = np.array(pool.map(objfun_noisy, map_arg))
node_val = res[:, 0]
node_err_bounds = res[:, 1:]
self.noisy_evalcount += len(node_val)
node_is_noisy = np.array([self.eval_mode == 'noisy'
for val in node_val])
# Add previously evaluated points, if any
if (self.num_nodes_at_restart and prev_node_pos):
node_pos = np.vstack((node_pos, np.array(prev_node_pos)))
node_val = np.append(node_val, np.array(prev_node_val))
node_is_noisy = np.append(
node_is_noisy, np.array(prev_node_is_noisy, dtype=bool))
node_err_bounds = np.vstack((node_err_bounds,
np.array(prev_node_err_bounds)))
else:
# We did not perform the initialization strategy, so use
# empty arrays
node_pos = np.zeros(shape=(0, self.n))
node_val = np.zeros(shape=(0, 1))
node_is_noisy = np.zeros(shape=(0, 1), dtype=bool)
node_err_bounds = np.zeros(shape=(0, 2))
# Add user-provided points if this is the first iteration
if (self.init_node_pos.size and self.itercount == 0):
# Determine which points can be added, based on distance
# from other points
self_dist = ([float('+inf')] +
[ru.get_min_distance(self.init_node_pos[i],
self.init_node_pos[:i])
for i in range(1, len(self.init_node_pos))])
if (node_pos.size > 0):
dist = ru.bulk_get_min_distance(self.init_node_pos,
node_pos)
selected = [i for i in range(len(self.init_node_pos))
if (dist[i] > self.l_settings.min_dist and
self_dist[i] > self.l_settings.min_dist)]
else:
selected = [i for i in range(len(self.init_node_pos))
if (self_dist[i] > self.l_settings.min_dist)]
init_node_pos = np.array(self.init_node_pos[selected])
if (self.init_node_val.size < len(self.init_node_pos)):
# Not enough values are provided, so we evaluate each node.
if (pool is None):
init_node_val = np.array([objfun([self.bb, point,
self.fixed_vars])
for point in init_node_pos])
else:
map_arg = [[self.bb, point, self.fixed_vars]
for point in init_node_pos]
init_node_val = np.array(pool.map(objfun, map_arg))
self.evalcount += len(init_node_val)
else:
init_node_val = self.init_node_val[selected]
if (init_node_pos.size):
node_pos = np.vstack((node_pos, init_node_pos))
node_val = np.append(node_val, init_node_val)
node_is_noisy = np.append(node_is_noisy,
np.zeros(len(init_node_val), dtype=bool))
node_err_bounds = np.vstack((node_err_bounds,
np.zeros((len(init_node_val), 2))))
if (self.all_node_pos.size == 0):
self.all_node_pos = node_pos.copy()
self.all_node_err_bounds = node_err_bounds.copy()
else:
self.all_node_pos = np.vstack((self.all_node_pos, node_pos))
self.all_node_err_bounds = np.vstack((self.all_node_err_bounds,
node_err_bounds))
self.all_node_val = np.append(self.all_node_val, node_val)
self.all_node_is_noisy = np.append(self.all_node_is_noisy,
node_is_noisy)
# Rescale the domain of the function
node_pos = ru.bulk_transform_domain(self.l_settings, self.var_lower,
self.var_upper, node_pos)
# Update references
self.node_pos, self.node_val = node_pos, node_val
self.node_is_noisy = node_is_noisy
self.node_err_bounds = node_err_bounds
# Update all counters and values to restart properly
self.fmin_index = np.argmin(node_val)
self.fmin = node_val[self.fmin_index]
self.fmax = np.max(node_val)
self.fmin_stall_check = self.fmin
self.fmin_last_refine = self.fmin
self.iter_last_refine = self.itercount
self.num_stalled_iter = 0
self.discarded_iters.clear()
self.num_cons_refinement = 0
self.is_fmin_noisy = self.node_is_noisy[self.fmin_index]
if (self.fmin < self.fbest):
self.fbest = self.fmin
self.fbest_index = self.fmin_index + self.num_nodes_at_restart
self.is_fbest_noisy = self.node_is_noisy[self.fmin_index]
gap = ru.compute_gap(self.l_settings, self.fbest +
self.all_node_err_bounds[self.fbest_index, 1])
self.current_step = self.first_step
self.cyclecount += 1
# Print the initialization points
for (i, val) in enumerate(self.node_val):
min_dist = ru.get_min_distance(self.node_pos[i],
np.vstack((self.node_pos[:i],
self.node_pos[(i+1):])))
self.update_log('Initialization', self.node_is_noisy[i], val, gap)
if (len(self.node_pos) < self.n + 1):
raise RuntimeError('Not enough initialization points; try ' +
'again with do_init_strategy=True.')
# -- end function
[docs] def restoration_search(self):
"""Perform restoration step to repair RBF matrix.
Try to repair an ill-conditioned RBF matrix by selecting
points far enough from current interpolation nodes, until
numerical stability is restored.
Returns
-------
1D numpy.ndarray[float] or None
The next point to be evaluated, or None if it cannot be
found.
"""
restoration_done = False
cons_restoration = 0
self.remove_node(len(self.node_pos) - 1, self.num_nodes_at_restart)
while (not restoration_done and cons_restoration <
self.l_settings.max_consecutive_restoration):
if (cons_restoration == 0):
# First, try to get the next point through something
# similar to a global search, using the MSRSM
# algorithm.
temp_settings = RbfoptSettings(
rbf='linear', algorithm='MSRSM',
global_search_method='genetic')
next_p = pure_global_step(temp_settings,
self.n, len(self.node_pos),
self.l_lower, self.l_upper,
self.integer_vars,
self.node_pos, None)
else:
# If that does not work (unlikely), generate a random point
next_p = (np.random.rand(self.n) *
(self.var_upper - self.var_lower) + self.var_lower)
node_pos = np.vstack((self.node_pos, next_p))
ru.round_integer_vars(next_p, self.integer_vars)
if (ru.get_min_distance(next_p, self.node_pos) >
self.l_settings.min_dist):
# Try inverting the RBF matrix to see if
# nonsingularity is restored
try:
Amat = ru.get_rbf_matrix(self.l_settings, self.n,
len(node_pos), node_pos)
Amatinv = ru.get_matrix_inverse(self.l_settings, Amat)
restoration_done = True
except np.linalg.LinAlgError:
cons_restoration += 1
else:
cons_restoration += 1
return (next_p if restoration_done else None)
# -- end function
[docs] def phase_update(self):
"""Check if we should switch phase in two-phase optimization.
Check if we should switch to the second phase of two-phase
optimization. The conditions for switching are:
1) Optimization in noisy mode restarted too many times.
2) We reached the limit of noisy mode iterations.
If both are met, the switch is performed.
"""
if ((self.two_phase_optimization == True) and
(self.eval_mode == 'noisy') and
((self.num_noisy_restarts > self.l_settings.max_noisy_restarts) or
(self.itercount >= self.l_settings.max_noisy_iterations) or
(self.noisy_evalcount >= self.l_settings.max_noisy_evaluations))):
self.update_log('Switching to accurate mode')
self.eval_mode = 'accurate'
# -- end function
[docs] def unlimited_refinement_active(self):
"""Should the usual limits of refinement phase be ignored?
Verify if, based on the unlimited_refinement_thresh parameter,
the usual limitations on the refinement phase should be
ignored.
Returns
-------
bool
True if unlimited refinement is active, False otherwise.
"""
return ((self.itercount >= self.l_settings.max_iterations *
self.l_settings.thresh_unlimited_refinement) or
(self.evalcount >= self.l_settings.max_evaluations *
self.l_settings.thresh_unlimited_refinement) or
(self.elapsed_time >= self.l_settings.max_clock_time *
self.l_settings.thresh_unlimited_refinement) or
(self.num_stalled_iter >=
self.l_settings.max_stalled_iterations *
self.l_settings.thresh_unlimited_refinement_stalled))
# -- end function
[docs] def refinement_update(self, model_impr, real_impr, to_replace):
"""Perform updates to refinement step and decide if continue.
Update the radius of the trust region and the iterate for the
refinement step. Also, decide if the next step should be
another refinement step, or if we should go back to global
search.
Parameters
----------
model_impr : float
Improvement in quadratic model value.
real_impr : float
Improvement in the real function value.
to_replace : int
Index in tr_model_set of the point to replace.
"""
self.num_cons_refinement += 1
# Update trust region information
self.tr_radius, tr_move = ref.update_trust_region_radius(
self.l_settings, self.tr_radius, model_impr, real_impr)
if (tr_move or
self.tr_model_set[to_replace] == self.tr_iterate_index):
# Accept new iterate
self.tr_iterate_index = len(self.node_pos) - 1
if (self.tr_radius < self.l_settings.tr_min_radius or
(self.num_cons_refinement >=
self.l_settings.max_consecutive_refinement and
not self.unlimited_refinement_active())):
# Do no continue refinement
self.num_cons_refinement = 0
self.current_step = self.first_step
self.cyclecount += 1
else:
# Continue refinement. Update set of points
# used to build model.
self.tr_model_set[to_replace] = len(self.node_pos) - 1
if (self.num_cons_refinement <
self.l_settings.max_consecutive_refinement):
# Update value at last refinement, unless we stopped
# refinement because of limit of iterations reached
self.fmin_last_refine = self.fmin
# -- end function
[docs] def refinement_update_parallel(self, model_impr, real_impr, to_replace):
"""Perform updates to refinement step and decide if continue.
Update the radius of the trust region and the iterate for the
refinement step. This is the version that should be used for
parallel computation, because the update phase is different.
Parameters
----------
model_impr : float
Improvement in quadratic model value.
real_impr : float
Improvement in the real function value.
to_replace : int
Index in the tr_model_set of the point to replace.
"""
self.num_cons_refinement += 1
# Update trust region information
self.tr_radius, tr_move = ref.update_trust_region_radius(
self.l_settings, self.tr_radius, model_impr, real_impr)
if (tr_move or
self.tr_model_set[to_replace] == self.tr_iterate_index):
# Accept new iterate
self.tr_iterate_index = len(self.node_pos) - 1
if (self.tr_radius < self.l_settings.tr_min_radius):
# Do no continue refinement
self.num_cons_refinement = 0
else:
# Continue refinement. Update set of points
# used to build model.
self.tr_model_set[to_replace] = len(self.node_pos) - 1
if (self.num_cons_refinement <
self.l_settings.max_consecutive_refinement):
# Update value at last refinement, unless we stopped
# refinement because of limit of iterations reached
self.fmin_last_refine = self.fmin
# -- end function
[docs] def advance_step_counter(self):
"""Advance the step counter of the optimization algorithm.
Advance the step counter of the optimization algorithm, and
update cycle number.
"""
# Advance to next step.
if (self.current_step > (self.current_step + 1) % self.cycle_length):
self.cyclecount += 1
self.current_step = (self.current_step + 1) % self.cycle_length
# -- end function
[docs] def stalling_update(self):
"""Check if the algorithm is stalling.
Check if the algorithm is stalling, and perform the
corresponding updates.
"""
if (self.fmin <= (self.fmin_stall_check -
self.l_settings.eps_impr *
max(1.0, abs(self.fmin_stall_check)))):
self.num_stalled_iter = 0
self.fmin_stall_check = self.fmin
else:
self.num_stalled_iter += 1
# -- end function
[docs] def require_accurate_evaluation(self, noisy_val):
"""Check if a given noisy value qualifies for accurate evaluation.
Verify if a point with the given objective function value in
noisy mode qualifies for an immediate accurate re-evaluation.
Parameters
----------
noisy_val : float
Value of the point to be tested, in noisy mode.
Returns
-------
bool
True if the point should be re-evaluated in accurate mode
immediately.
"""
# Check if the point improves over existing points,
# or if it could be optimal according to tolerances.
# In this case, perform a double evaluation.
best_possible = self.fmin
if ((noisy_val <= best_possible -
self.l_settings.eps_impr*max(1.0, abs(best_possible))) or
(noisy_val <= self.l_settings.target_objval +
self.l_settings.eps_opt*abs(self.l_settings.target_objval))):
return True
else:
return False
# -- end function
# -- end class
[docs]def pure_global_step(settings, n, k, var_lower, var_upper, integer_vars,
node_pos, mat):
"""Perform the pure global search step.
Parameters
----------
mat : numpy.matrix
The matrix necessary for the computation. This is the inverse
of the matrix [Phi P; P^T 0]. Must be a square numpy.matrix of
appropriate dimension.
settings : :class:`rbfopt_settings.RbfoptSettings`
Global and algorithmic settings.
n : int
The dimension of the problem, i.e. size of the space.
k : int
Number of nodes, i.e. interpolation points.
var_lower : 1D numpy.ndarray[float]
Vector of variable lower bounds.
var_upper : 1D numpy.ndarray[float]
Vector of variable upper bounds.
integer_vars : 1D numpy.ndarray[int]
A list containing the indices of the integrality constrained
variables. If empty list, all variables are assumed to be
continuous.
node_pos : List[List[float]]
List of coordinates of the nodes
mat : numpy.matrix
The matrix necessary for the computation. This is the inverse
of the matrix [Phi P; P^T 0], see paper as cited above. Must
be a square numpy.matrix of appropriate dimension. Can be None
when using the MSRSM algorithm.
Returns
-------
List[float] or None
The point to be evaluated next, or None if errors occurred.
"""
assert(len(var_lower)==n)
assert(len(var_upper)==n)
assert(len(node_pos)==k)
assert((mat is None and settings.algorithm == 'MSRSM') or
isinstance(mat, np.matrix))
assert(isinstance(settings, RbfoptSettings))
assert(isinstance(var_lower, np.ndarray))
assert(isinstance(var_upper, np.ndarray))
assert(isinstance(integer_vars, np.ndarray))
assert(isinstance(node_pos, np.ndarray))
# Infstep: explore the parameter space
return aux.pure_global_search(settings, n, k, var_lower, var_upper,
integer_vars, node_pos, mat)
# -- end function
[docs]def local_step(settings, n, k, var_lower, var_upper, integer_vars,
node_pos, rbf_lambda, rbf_h, tfv, Amat, Amatinv,
fmin_index, two_phase_optimization, eval_mode,
node_is_noisy):
"""Perform local search step, possibly adjusted.
Perform a local search step. This typically accepts the
minimum of the RBF model as the next point if it is a viable
option; if this is not viable, it will perform an adjusted
local search and try to generate a different candidate. It
also verifies if it is better to evaluate a brand new point,
or re-evaluate a previously known point. The test is based on
bumpiness of the resulting interpolant.
Parameters
----------
settings : :class:`rbfopt_settings.RbfoptSettings`
Global and algorithmic settings.
n : int
The dimension of the problem, i.e. size of the space.
k : int
Number of nodes, i.e. interpolation points.
var_lower : 1D numpy.ndarray[float]
Vector of variable lower bounds.
var_upper : 1D numpy.ndarray[float]
Vector of variable upper bounds.
integer_vars: 1D numpy.ndarray[int]
A list containing the indices of the integrality constrained
variables. If empty list, all variables are assumed to be
continuous.
node_pos : 2D numpy.ndarray[float]
List of coordinates of the nodes.
rbf_lambda : 1D numpy.ndarray[float]
The lambda coefficients of the RBF interpolant, corresponding
to the radial basis functions. List of dimension k.
rbf_h : 1D numpy.ndarray[float]
The h coefficients of the RBF interpolant, corresponding to
the polynomial. List of dimension n+1.
tfv : (1D numpy.ndarray[float], float, float, 2D numpy.ndarray[float])
Transformed function values: scaled node values, scaled
minimum, scaled maximum, scaled node error bounds.
Amat : numpy.matrix
RBF matrix, i.e. [Phi P; P^T 0].
Amatinv : numpy.matrix or None
Inverse of the RBF matrix, i.e. [Phi P; P^T 0]^{-1}. Can be
None if the algorithm is MSRSM.
fmin_index : int
Index of the minimum value among the nodes.
two_phase_optimization : bool
Is the noisy but fast objective function is available?
eval_mode : string
Evaluation mode for the objective function at a given
stage. Can be either 'noisy' or 'accurate'.
node_is_noisy : List[bool]
For each interpolation node in node_pos, was it evaluated in
'noisy' mode?
Returns
-------
(bool, 1D numpy.ndarray[float] or None, int or None)
A triple (adjusted, point, index) where adjusted is True if
the local search was adjusted rather than a pure local search,
point is the point to be evaluated next (or None if errors
occurred), and index is the index that this point should
replace, or None if the point should be appended.
See also
--------
rbfopt_utils.transform_function_values()
"""
assert(isinstance(var_lower, np.ndarray))
assert(isinstance(var_upper, np.ndarray))
assert(isinstance(integer_vars, np.ndarray))
assert(isinstance(node_pos, np.ndarray))
assert(isinstance(rbf_lambda, np.ndarray))
assert(isinstance(rbf_h, np.ndarray))
assert(len(var_lower)==n)
assert(len(var_upper)==n)
assert(len(rbf_lambda)==k)
assert(len(node_pos)==k)
assert(len(node_is_noisy)==k)
assert(0 <= fmin_index < k)
assert((eval_mode=='noisy') or (eval_mode=='accurate'))
assert(isinstance(settings, RbfoptSettings))
assert(isinstance(Amat, np.matrix))
assert((Amatinv is None and settings.algorithm == 'MSRSM') or
isinstance(Amatinv, np.matrix) )
scaled_node_val, scaled_fmin, scaled_fmax, scaled_err_bounds = tfv
# Local search: compute the minimum of the RBF.
min_rbf = aux.minimize_rbf(settings, n, k, var_lower, var_upper,
integer_vars, node_pos, rbf_lambda, rbf_h,
node_pos[fmin_index])
if (min_rbf is not None):
min_rbf_val = ru.evaluate_rbf(settings, min_rbf, n, k,
node_pos, rbf_lambda, rbf_h)
# If the RBF cannot me minimized, or if the minimum is
# larger than the node with smallest value, just take the
# node with the smallest value.
if (min_rbf is None or
(min_rbf_val >= scaled_fmin + settings.eps_zero)):
min_rbf = node_pos[fmin_index]
min_rbf_val = scaled_fmin
# Check if point can be accepted: is there an improvement?
if (min_rbf_val <= (scaled_fmin - settings.eps_impr *
max(1.0, abs(scaled_fmin)))):
target_val = min_rbf_val
next_p = min_rbf
adjusted = False
else:
# If the point is not improving, we solve a global
# search problem, rescaling the search box to enforce
# some sort of local search
target_val = scaled_fmin - 0.01*max(1.0, abs(scaled_fmin))
dist_weight = 0.05
local_varl = np.maximum(min_rbf - settings.local_search_box_scaling *
0.33 * (var_upper - var_lower),
var_lower)
local_varu = np.minimum(min_rbf + settings.local_search_box_scaling *
0.33 * (var_upper - var_lower),
var_upper)
ru.round_integer_bounds(local_varl, local_varu, integer_vars)
next_p = aux.global_search(settings, n, k, local_varl,
local_varu, integer_vars,
node_pos, rbf_lambda, rbf_h,
Amatinv, target_val, dist_weight,
scaled_fmin, scaled_fmax)
adjusted = True
# If previous points were evaluated in low quality and we are
# now in high-quality local search mode, then we should verify
# if it is better to evaluate a brand new point or re-evaluate
# a previously known point.
if ((two_phase_optimization == True) and (eval_mode == 'accurate')):
ind, bump = aux.get_min_bump_node(settings, n, k, Amat,
scaled_node_val,
scaled_err_bounds, target_val)
if (ind is not None and next_p is not None):
# Check if the newly proposed point is very close to an
# existing one.
dist, ind = ru.get_min_distance_and_index(next_p, node_pos)
if (dist > settings.min_dist):
# If not, compute bumpiness of the newly proposed point.
n_bump = aux.get_bump_new_node(settings, n, k, node_pos,
scaled_node_val, next_p,
scaled_err_bounds,
target_val)
else:
# If yes, we will simply reevaluate the existing point
# (if it can be reevaluated).
n_bump = (float('inf') if node_is_noisy[ind] else
float('-inf'))
if (n_bump > bump):
# In this case we want to put the new point at the
# same location as one of the old points.
return (True, node_pos[ind], ind)
# -- end if
# -- end if
return (adjusted, next_p, None)
# -- end function
[docs]def refinement_step(settings, n, k, var_lower, var_upper, integer_vars,
node_pos, tfv, h, b, rank_deficient, model_set,
start_point_index, tr_radius):
"""Perform a refinement step.
Perform a refinement step to locally improve the solution.
Parameters
----------
settings : :class:`rbfopt_settings.RbfoptSettings`
Global and algorithmic settings.
n : int
The dimension of the problem, i.e. size of the space.
k : int
Number of nodes, i.e. interpolation points.
var_lower : 1D numpy.ndarray[float]
Vector of variable lower bounds.
var_upper : 1D numpy.ndarray[float]
Vector of variable upper bounds.
integer_vars: 1D numpy.ndarray[int]
A list containing the indices of the integrality constrained
variables. If empty list, all variables are assumed to be
continuous.
node_pos : 2D numpy.ndarray[float]
List of coordinates of the nodes.
h : 1D numpy.ndarray[float]
Linear coefficients of the quadratic model.
b : float
Constant term of the quadratic model.
rank_deficient : bool
True if the points used to build the linear model do not span
the space.
model_set : 1D numpy.ndarray[int]
Indices of points in node_pos to be used to compute model.
start_point_index : int
Index in node_pos of the starting point for the descent.
tr_radius : float
Radius of the trust region.
Returns
-------
(1D numpy.ndarray[float], float, int)
Next candidate point for the search, the corresponding model
value difference, and the index in model_set of the point to
replace.
"""
assert(isinstance(var_lower, np.ndarray))
assert(isinstance(var_upper, np.ndarray))
assert(isinstance(h, np.ndarray))
assert(len(var_lower)==n)
assert(len(var_upper)==n)
assert(start_point_index < k)
assert(len(h)==n)
assert(tr_radius>=0)
assert(isinstance(settings, RbfoptSettings))
scaled_node_val, scaled_fmin, scaled_fmax, scaled_err_bounds = tfv
found_model_impr = False
if (rank_deficient):
point, found_model_impr, to_replace = ref.get_model_improving_point(
settings, n, k, var_lower, var_upper, node_pos, model_set,
start_point_index, tr_radius, integer_vars)
model_impr = np.dot(h, point) + b
if (not rank_deficient or not found_model_impr):
start_point = node_pos[start_point_index]
point, model_impr, grad_norm = ref.get_candidate_point(
settings, n, k, var_lower, var_upper, h, start_point, tr_radius)
to_replace = np.argmax(scaled_node_val[model_set])
if (len(integer_vars)):
# Get a rounding of the point
point, model_impr_adj = ref.get_integer_candidate(
settings, n, k, h, start_point, tr_radius, point, integer_vars)
model_impr += model_impr_adj
if (grad_norm <= settings.tr_min_grad_norm):
return None, model_impr, to_replace
return point, model_impr, to_replace
# -- end function
[docs]def global_step(settings, n, k, var_lower, var_upper, integer_vars,
node_pos, rbf_lambda, rbf_h, tfv, Amatinv,
fmin_index, current_step):
"""Perform global search step.
Perform a global search step, with a different methodology
depending on the algorithm chosen.
Parameters
----------
settings : :class:`rbfopt_settings.RbfoptSettings`
Global and algorithmic settings.
n : int
The dimension of the problem, i.e. size of the space.
k : int
Number of nodes, i.e. interpolation points.
var_lower : 1D numpy.ndarray[float]
Vector of variable lower bounds.
var_upper : 1D numpy.ndarray[float]
Vector of variable upper bounds.
integer_vars: 1D numpy.ndarray[int]
A list containing the indices of the integrality constrained
variables. If empty list, all variables are assumed to be
continuous.
node_pos : 2D numpy.ndarray[float]
List of coordinates of the nodes.
rbf_lambda : 1D numpy.ndarray[float]
The lambda coefficients of the RBF interpolant, corresponding
to the radial basis functions. List of dimension k.
rbf_h : 1D numpy.ndarray[float]
The h coefficients of the RBF interpolant, corresponding to
the polynomial. List of dimension n+1.
tfv : (List[float], float, float, List[(float, float)])
Transformed function values: scaled node values, scaled
minimum, scaled maximum, and node error bounds.
Amatinv : numpy.matrix or None
The matrix necessary for the computation. This is the inverse
of the matrix [Phi P; P^T 0], see paper as cited above. Must
be a square numpy.matrix of appropriate dimension. Can be None
if algorithm is MSRSM.
fmin_index : int
Index of the minimum value among the nodes.
current_step : int
Identifier of the current step within the cyclic optimization
strategy counter.
Returns
-------
1D numpy.ndarray[float] or None
The point to be evaluated next, or None if errors occurred.
"""
assert(isinstance(var_lower, np.ndarray))
assert(isinstance(var_upper, np.ndarray))
assert(isinstance(integer_vars, np.ndarray))
assert(isinstance(rbf_lambda, np.ndarray))
assert(isinstance(rbf_h, np.ndarray))
assert(isinstance(node_pos, np.ndarray))
assert(len(var_lower)==n)
assert(len(var_upper)==n)
assert(len(rbf_lambda)==k)
assert(len(node_pos)==k)
assert(0 <= fmin_index < k)
assert((Amatinv is None and settings.algorithm == 'MSRSM') or
isinstance(Amatinv, np.matrix))
assert(isinstance(settings, RbfoptSettings))
assert(0 <= current_step <= settings.num_global_searches)
scaled_node_val, scaled_fmin, scaled_fmax, scaled_err_bounds = tfv
assert (isinstance(scaled_node_val, np.ndarray))
# Global search: compromise between finding a good value of the
# objective function, and improving the model.
if (settings.algorithm == 'Gutmann'):
# If we use Gutmann's algorithm, we need the minimum of the
# RBF interpolant to choose the target value.
min_rbf = aux.minimize_rbf(settings, n, k, var_lower, var_upper,
integer_vars, node_pos, rbf_lambda,
rbf_h, node_pos[fmin_index])
if (min_rbf is not None):
min_rbf_val = ru.evaluate_rbf(settings, min_rbf, n, k,
node_pos, rbf_lambda, rbf_h)
# If the RBF cannot me minimized, or if the minimum is larger
# than the node with smallest value, just take the node with
# the smallest value.
if (min_rbf is None or
min_rbf_val >= scaled_fmin + settings.eps_zero):
min_rbf = node_pos[fmin_index]
min_rbf_val = scaled_fmin
else:
# For the Metric SRSM method, pick the node with the smallest
# value as minimum of the RBF.
min_rbf = node_pos[fmin_index]
min_rbf_val = scaled_fmin
# Compute the function value used to determine the target
# value. This is given by the sorted value in position sigma_n,
# where sigma_n is a function described in the paper by Gutmann
# (2001). If clipping is disabled, we simply take the largest
# function value.
if (settings.targetval_clipping):
local_fmax = ru.get_fmax_current_iter(settings, n, k,
current_step, scaled_node_val)
else:
local_fmax = scaled_fmax
# For Gutmann's RBF method, the scaling factor is 1 - h/kappa,
# where h goes from 0 to kappa-1 over the course of one global
# search cycle, and kappa is the number of global searches.
scaling = (1 - ((current_step - 1) / settings.num_global_searches))**2
target_val = (min_rbf_val - scaling * (local_fmax - min_rbf_val))
# For Metric SRSM, The weighting factor for the distance is is
# (h+1)/kappa, where h goes from 0 to kappa-1 over the course of
# one global search cycle, and kappa is the number of global
# searches. If dist_weight would be 0, we set it to 0.05.
dist_weight = (1 - (current_step / settings.num_global_searches)
if (current_step < settings.num_global_searches)
else 0.05)
# If the global search is almost a local search, we restrict the
# search to a box following Regis and Shoemaker (2007)
if (scaling <= settings.local_search_threshold):
local_varl = np.maximum(min_rbf - settings.local_search_box_scaling *
0.33 * (var_upper - var_lower),
var_lower)
local_varu = np.minimum(min_rbf + settings.local_search_box_scaling *
0.33 * (var_upper - var_lower),
var_upper)
ru.round_integer_bounds(local_varl, local_varu, integer_vars)
else:
# Otherwise, use original bounds
local_varl = var_lower
local_varu = var_upper
next_p = aux.global_search(settings, n, k, local_varl, local_varu,
integer_vars, node_pos, rbf_lambda,
rbf_h, Amatinv, target_val, dist_weight,
scaled_fmin, scaled_fmax)
return next_p
# -- end function
[docs]def objfun(data):
"""Call the evaluate() method of a RbfoptBlackBox object.
Apply the evaluate() method of the given RbfoptBlackBox object to
the given point. This way of calling the method indirectly is
necessary for parallelization.
Parameters
----------
data : (rbfopt_black_box.RbfoptBlackBox, 1D numpy.ndarray[float],
List[(int, float)])
A triple or list with three elements (black_box, point,
fixed_vars) containing an object derived from class
RbfoptBlackBox, that describes the problem, the point at which
we want to apply the evaluate() method, and a list of fixed
variables given as pairs (index, value).
Returns
-------
float
The value of the function evaluate() at the point.
"""
assert(len(data)==3)
assert(isinstance(data[0], RbfoptBlackBox))
if (data[2]):
point = data[1].tolist()
for (i, val) in data[2]:
point.insert(i, val)
point = np.array(point)
else:
point = data[1]
return data[0].evaluate(point)
# -- end function
[docs]def objfun_noisy(data):
"""Call the evaluate_noisy() method of a RbfoptBlackBox object.
Apply the evaluate_noisy() method of the given RbfoptBlackBox object to
the given point. This way of calling the method indirectly is
necessary for parallelization.
Parameters
----------
data : (rbfopt_black_box.RbfoptBlackBox, 1D numpy.array[float],
List[(int, float)])
A triple or list with three elements (black_box, point,
fixed_vars) containing an object derived from class
RbfoptBlackBox, that describes the problem, the point at which
we want to apply the evaluate() method, and a list of fixed
variables given as pairs (index, value).
Returns
-------
(float, float, float)
The value of the function evaluate_noisy() at the point, and
the possible variation given as lower and upper variation.
"""
assert(len(data)==3)
assert(isinstance(data[0], RbfoptBlackBox))
if (data[2]):
point = data[1].tolist()
for (i, val) in data[2]:
point.insert(i, val)
point = np.array(point)
else:
point = data[1]
return data[0].evaluate_noisy(point)
# -- end function