diff --git a/gpflowopt/acquisition/__init__.py b/gpflowopt/acquisition/__init__.py index 4237d79..375a78f 100644 --- a/gpflowopt/acquisition/__init__.py +++ b/gpflowopt/acquisition/__init__.py @@ -13,7 +13,8 @@ # limitations under the License. # Framework components and interfaces -from .acquisition import Acquisition, AcquisitionAggregation, AcquisitionProduct, AcquisitionSum, MCMCAcquistion +from .acquisition import (Acquisition, ParallelBatchAcquisition, AcquisitionAggregation, AcquisitionProduct, \ + AcquisitionSum, MCMCAcquistion, IAcquisition, ParToSeqAcquisitionWrapper, setup_required) # Single objective from .ei import ExpectedImprovement @@ -24,5 +25,8 @@ # Multiobjective from .hvpoi import HVProbabilityOfImprovement +# Batch +from .qei import qExpectedImprovement + # Black-box constraint from .pof import ProbabilityOfFeasibility diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index c1b0fd4..d7f67fc 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -25,6 +25,8 @@ import copy from functools import wraps +from abc import ABCMeta, abstractmethod, abstractproperty +import six float_type = settings.dtypes.float_type @@ -36,7 +38,7 @@ def setup_required(method): """ @wraps(method) def runnable(instance, *args, **kwargs): - assert isinstance(instance, Acquisition) + assert isinstance(instance, ParallelBatchAcquisition) hp = instance.highest_parent if hp._needs_setup: # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition @@ -54,41 +56,25 @@ def runnable(instance, *args, **kwargs): return runnable -class Acquisition(Parameterized): +@six.add_metaclass(ABCMeta) +class IAcquisition(Parameterized): # pragma: no cover """ - An acquisition function maps the belief represented by a Bayesian model into a - score indicating how promising a point is for evaluation. + Interface for Acquisition functions mapping the belief represented by a Bayesian model into a score indicating how + promising a point is for evaluation. In Bayesian Optimization this function is typically optimized over the optimization domain - to determine the next point for evaluation. An object of this class holds a list of GPflow models. Subclasses - implement a build_acquisition function which computes the acquisition function (usually from the predictive - distribution) using TensorFlow. Optionally, a method setup can be implemented which computes some quantities which - are used to compute the acquisition, but do not depend on candidate points. + to determine the next point for evaluation. An objects satisfying this interface hold a list of GPflow models. + Implementing build_acquisition yields the score of the acquisition function (usually obtained by transforming the + predictive distribution) using TensorFlow. Optionally, a method _setup can be implemented which computes some + quantities which are used to compute the acquisition, do not vary with the candidate points. This method is not + automatically called, classes implementing this interface may implement their own mechanism on when _setup is + called. Acquisition functions can be combined through addition or multiplication to construct joint criteria. For instance, for constrained optimization. The objects then form a tree hierarchy. - - Acquisition models implement a lazy strategy to optimize models and run setup. This is implemented by a _needs_setup - attribute (similar to the _needs_recompile in GPflow). Calling :meth:`set_data` sets this flag to True. Calling methods - marked with the setup_require decorator (such as evaluate) optimize all models, then call setup if this flag is set. - In hierarchies, first acquisition objects handling constraint objectives are set up, then the objects handling - objectives. """ - def __init__(self, models=[], optimize_restarts=5): - """ - :param models: list of GPflow models representing our beliefs about the problem - :param optimize_restarts: number of optimization restarts to use when training the models - """ - super(Acquisition, self).__init__() - models = np.atleast_1d(models) - assert all(isinstance(model, (Model, ModelWrapper)) for model in models) - self._models = ParamList([DataScaler(m) for m in models]) - - assert (optimize_restarts >= 0) - self.optimize_restarts = optimize_restarts - self._needs_setup = True - + @abstractmethod def _optimize_models(self): """ Optimizes the hyperparameters of all models that the acquisition function is based on. @@ -103,81 +89,65 @@ def _optimize_models(self): As a special case, if optimize_restarts is set to zero, the hyperparameters of the models are not optimized. This is useful when the hyperparameters are sampled using MCMC. """ - if self.optimize_restarts == 0: - return - - for model in self.models: - runs = [] - for i in range(self.optimize_restarts): - if i > 0: - model.randomize() - try: - result = model.optimize() - runs.append(result) - except tf.errors.InvalidArgumentError: # pragma: no cover - print("Warning: optimization restart {0}/{1} failed".format(i + 1, self.optimize_restarts)) - if not runs: - raise RuntimeError("All model hyperparameter optimization restarts failed, exiting.") - best_idx = np.argmin([r.fun for r in runs]) - model.set_state(runs[best_idx].x) + pass - def build_acquisition(self, Xcand): - raise NotImplementedError + @abstractmethod + def _setup(self): + """ + Pre-calculation of quantities used later in the evaluation of the acquisition function for candidate points. - def enable_scaling(self, domain): + Subclasses can implement this method to compute quantities (such as fmin). The decision when to run this function + is governed by :class:`Acquisition`, based on the setup_required decorator on methods which require + setup to be run (e.g. set_data). """ - Enables and configures the :class:`.DataScaler` objects wrapping the GP models. - - Sets the _needs_setup attribute to True so the contained models are optimized and :meth:`setup` is run again - right before evaluating the :class:`Acquisition` function. Note that the models are modified directly and - references to them outside of the object will also point to scaled instances. + pass - :param domain: :class:`.Domain` object, the input transform of the data scalers is configured as a transform - from domain to the unit cube with the same dimensionality. + @abstractmethod + def _setup_constraints(self): """ - n_inputs = self.data[0].shape[1] - assert (domain.size == n_inputs) - for m in self.models: - m.input_transform = domain >> UnitCube(n_inputs) - m.normalize_output = True - self.highest_parent._needs_setup = True + Run only if some outputs handled by this acquisition are constraints. Used in aggregation. + """ + pass - def set_data(self, X, Y): + @abstractmethod + def _setup_objectives(self): """ - Update the training data of the contained models + Run only if all outputs handled by this acquisition are objectives. Used in aggregation. + """ + pass - Sets the _needs_setup attribute to True so the contained models are optimized and :meth:`setup` is run again - right before evaluating the :class:`Acquisition` function. + @abstractmethod + def get_suggestion(self, optimizer): + pass - Let Q be the the sum of the output dimensions of all contained models, Y should have a minimum of - Q columns. Only the first Q columns of Y are used while returning the scalar Q + @abstractmethod + def build_acquisition(self, *args): + pass - :param X: input data N x D - :param Y: output data N x R (R >= Q) - :return: Q (sum of output dimensions of contained models) + @abstractmethod + def enable_scaling(self, domain): """ - num_outputs_sum = 0 - for model in self.models: - num_outputs = model.Y.shape[1] - Ypart = Y[:, num_outputs_sum:num_outputs_sum + num_outputs] - num_outputs_sum += num_outputs + Enables and configures the :class:`.DataScaler` objects wrapping the GP models. - model.X = X - model.Y = Ypart + Sets the _needs_setup attribute to True so the contained models are optimized and :meth:`setup` is run again + right before evaluating the acquisition function. Note that the models are modified directly and + references to them outside of the object will also point to scaled instances. - self.highest_parent._needs_setup = True - return num_outputs_sum + :param domain: :class:`.Domain` object, the input transform of the data scalers is configured as a transform + from domain to the unit cube with the same dimensionality. + """ + pass - @property + @abstractproperty def models(self): """ The GPflow models representing our beliefs of the optimization problem. - - :return: list of GPflow models + + :return: list of GPflow models """ - return self._models.sorted_params + pass - @property + @abstractproperty def data(self): """ The training data of the models. @@ -187,85 +157,80 @@ def data(self): :return: tuple X, Y of tensors (if in tf_mode) or numpy arrays. """ - if self._tf_mode: - return self.models[0].X, tf.concat(list(map(lambda model: model.Y, self.models)), 1) - else: - return self.models[0].X.value, np.hstack(map(lambda model: model.Y.value, self.models)) + pass + @abstractproperty + def batch_size(self): + pass + + @abstractmethod + def set_data(self, X, Y): + """ + Update the training data of the contained models + + Sets the _needs_setup attribute to True so the contained models are optimized and :meth:`_setup` is run again + right before evaluating the acquisition function. + + Let Q be the the sum of the output dimensions of all contained models, Y should have a minimum of + Q columns. Only the first Q columns of Y are used while returning the scalar Q + + :param X: input data N x D + :param Y: output data N x R (R >= Q) + :return: Q (sum of output dimensions of contained models) + """ + pass + + @abstractmethod def constraint_indices(self): """ Method returning the indices of the model outputs which correspond to the (expensive) constraint functions. By default there are no constraint functions """ - return np.empty((0,), dtype=int) + pass + @abstractmethod def objective_indices(self): """ Method returning the indices of the model outputs which are objective functions. - + By default all outputs are objectives. - + :return: indices to the objectives, size R """ - return np.setdiff1d(np.arange(self.data[1].shape[1]), self.constraint_indices()) + pass + @abstractmethod def feasible_data_index(self): """ Returns a boolean array indicating which data points are considered feasible (according to the acquisition function(s) ) and which not. - + By default all data is considered feasible. - - :return: logical indices to the feasible data points, size N - """ - return np.ones(self.data[0].shape[0], dtype=bool) - def _setup(self): - """ - Pre-calculation of quantities used later in the evaluation of the acquisition function for candidate points. - - Subclasses can implement this method to compute quantities (such as fmin). The decision when to run this function - is governed by :class:`Acquisition`, based on the setup_required decorator on methods which require - setup to be run (e.g. set_data). + :return: logical indices to the feasible data points, size N """ pass - def _setup_constraints(self): + @abstractmethod + def evaluate(self, candidates): """ - Run only if some outputs handled by this acquisition are constraints. Used in aggregation. - """ - if self.constraint_indices().size > 0: - self._setup() + AutoFlow method to compute the acquisition scores for candidates, without returning the gradients. - def _setup_objectives(self): - """ - Run only if all outputs handled by this acquisition are objectives. Used in aggregation. + :return: acquisition scores, size N x 1 """ - if self.constraint_indices().size == 0: - self._setup() + pass - @setup_required - @AutoFlow((float_type, [None, None])) - def evaluate_with_gradients(self, Xcand): + @abstractmethod + def evaluate_with_gradients(self, candidates): """ AutoFlow method to compute the acquisition scores for candidates, also returns the gradients. - - :return: acquisition scores, size N x 1 - the gradients of the acquisition scores, size N x D - """ - acq = self.build_acquisition(Xcand) - return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] - @setup_required - @AutoFlow((float_type, [None, None])) - def evaluate(self, Xcand): - """ - AutoFlow method to compute the acquisition scores for candidates, without returning the gradients. - :return: acquisition scores, size N x 1 + the gradients of the acquisition scores, size N x D """ - return self.build_acquisition(Xcand) + pass + @abstractmethod def __add__(self, other): """ Operator for adding acquisition functions. Example: @@ -275,10 +240,9 @@ def __add__(self, other): >>> type(a1 + a2) """ - if isinstance(other, AcquisitionSum): - return AcquisitionSum([self] + other.operands.sorted_params) - return AcquisitionSum([self, other]) + pass + @abstractmethod def __mul__(self, other): """ Operator for multiplying acquisition functions. Example: @@ -288,56 +252,46 @@ def __mul__(self, other): >>> type(a1 * a2) """ - if isinstance(other, AcquisitionProduct): - return AcquisitionProduct([self] + other.operands.sorted_params) - return AcquisitionProduct([self, other]) - - def __setattr__(self, key, value): - super(Acquisition, self).__setattr__(key, value) - if key is '_parent': - self.highest_parent._needs_setup = True + pass -class AcquisitionAggregation(Acquisition): +@six.add_metaclass(ABCMeta) +class AcquisitionWrapper(IAcquisition): """ - Aggregates multiple acquisition functions, using a TensorFlow reduce operation. + Partial implementation of the Acquisition interface, useful for implementing acquisition functions which wrap + another acquisition function (and forward method calls). + + Typically used in combination with multiple inheritance: + + >>> class WrappedAcquisition(AcquisitionWrapper, Acquisition) + >>> def __init__(self, acquisition): + >>> super(WrappedAcquisition, self).__init__(acquisition, optimize_restarts=2) + >>> + >>> def build_acquisition(self, candidates): + >>> return tf.constant(1.0) """ - def __init__(self, operands, oper): + def __init__(self, wrap, *args, **kwargs): """ - :param operands: list of acquisition objects - :param oper: a tf.reduce operation (e.g., tf.reduce_sum) for aggregating the returned scores of each operand. + :param wrap: list of acquisition objects """ - super(AcquisitionAggregation, self).__init__() - assert (all([isinstance(x, Acquisition) for x in operands])) - self.operands = ParamList(operands) - self._oper = oper + self.wrapped = [] + super(AcquisitionWrapper, self).__init__(*args, **kwargs) + assert all([isinstance(x, IAcquisition) for x in wrap]) + assert all(wrap[0].batch_size == oper.batch_size for oper in wrap) + self.wrapped = ParamList(np.atleast_1d(wrap).tolist()) def _optimize_models(self): - for oper in self.operands: + for oper in self.wrapped: oper._optimize_models() - @Acquisition.models.getter - def models(self): - return [model for acq in self.operands for model in acq.models] - - def enable_scaling(self, domain): - for oper in self.operands: - oper.enable_scaling(domain) - - def set_data(self, X, Y): - offset = 0 - for operand in self.operands: - offset += operand.set_data(X, Y[:, offset:]) - return offset - def _setup_constraints(self): - for oper in self.operands: + for oper in self.wrapped: if oper.constraint_indices().size > 0: # Small optimization, skip subtrees with objectives only oper._setup_constraints() def _setup_objectives(self): - for oper in self.operands: + for oper in self.wrapped: oper._setup_objectives() def _setup(self): @@ -346,23 +300,297 @@ def _setup(self): # Then objectives as these might depend on the constraint acquisition self._setup_objectives() + @abstractmethod + def build_acquisition(self, *args): # pragma: no cover + pass + + @property + def models(self): + return [model for acq in self.wrapped for model in acq.models] + + def enable_scaling(self, domain): + for oper in self.wrapped: + oper.enable_scaling(domain) + + def set_data(self, X, Y): + offset = 0 + for operand in self.wrapped: + offset += operand.set_data(X, Y[:, offset:]) + return offset + def constraint_indices(self): offset = [0] idx = [] - for operand in self.operands: + for operand in self.wrapped: idx.append(operand.constraint_indices()) offset.append(operand.data[1].shape[1]) return np.hstack([i + o for i, o in zip(idx, offset[:-1])]) def feasible_data_index(self): - return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.operands)), axis=0) - - def build_acquisition(self, Xcand): - return self._oper(tf.concat(list(map(lambda operand: operand.build_acquisition(Xcand), self.operands)), 1), - axis=1, keep_dims=True, name=self.__class__.__name__) + return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.wrapped)), axis=0) def __getitem__(self, item): - return self.operands[item] + return self.wrapped[item] + + +@six.add_metaclass(ABCMeta) +class ParallelBatchAcquisition(IAcquisition): + """ + Abstract Acquisition class for batch acquisition functions, optimizing the batch in parallel (as opposed to + sequential or greedy approaches). + + Subclasses deriving this class should accept multiple candidates tensors for build_acquisition. Each tensor + represents points x1, x2, ... xn in the batch. + + ParallelBatchAcquisition objects implement a lazy strategy to optimize models and run setup. This is implemented by + a _needs_setup attribute (similar to the _needs_recompile in GPflow). Calling :meth:`set_data` sets this flag to + True. Calling methods marked with the setup_require decorator (such as evaluate) optimize all models, then call + setup if this flag is set. In hierarchies, first acquisition objects handling constraint objectives are set up, then + the objects handling objectives. + """ + + def __init__(self, models=[], optimize_restarts=5, batch_size=1): + """ + :param models: list of GPflow models representing our beliefs about the problem + :param optimize_restarts: number of optimization restarts to use when training the models + """ + super(ParallelBatchAcquisition, self).__init__() + models = np.atleast_1d(models) + assert all(isinstance(model, (Model, ModelWrapper)) for model in models) + self._models = ParamList([DataScaler(m) for m in models]) + + assert (optimize_restarts >= 0) + self.optimize_restarts = optimize_restarts + self._needs_setup = True + self._batch_size = batch_size + + def _optimize_models(self): + if self.optimize_restarts == 0: + return + + for model in self.models: + runs = [] + for i in range(self.optimize_restarts): + if i > 0: + model.randomize() + try: + result = model.optimize() + runs.append(result) + except tf.errors.InvalidArgumentError: # pragma: no cover + print("Warning: optimization restart {0}/{1} failed".format(i + 1, self.optimize_restarts)) + if not runs: + raise RuntimeError("All model hyperparameter optimization restarts failed, exiting.") + best_idx = np.argmin([r.fun for r in runs]) + model.set_state(runs[best_idx].x) + + def _inverse_acquisition(self, x): + return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) + + def get_suggestion(self, optimizer): + opt = copy.deepcopy(optimizer) + batch_domain = optimizer.domain.batch(self.batch_size) + opt.domain = batch_domain + result = opt.optimize(self._inverse_acquisition) + return np.vstack(np.split(result.x, self.batch_size, axis=1)) + + def enable_scaling(self, domain): + n_inputs = self.data[0].shape[1] + assert (domain.size == n_inputs) + for m in self.models: + m.input_transform = domain >> UnitCube(n_inputs) + m.normalize_output = True + self.highest_parent._needs_setup = True + + def set_data(self, X, Y): + num_outputs_sum = 0 + for model in self.models: + num_outputs = model.Y.shape[1] + Ypart = Y[:, num_outputs_sum:num_outputs_sum + num_outputs] + num_outputs_sum += num_outputs + + model.X = X + model.Y = Ypart + + self.highest_parent._needs_setup = True + return num_outputs_sum + + @property + def models(self): + return self._models.sorted_params + + @property + def data(self): + if self._tf_mode: + return self.models[0].X, tf.concat(list(map(lambda model: model.Y, self.models)), 1) + else: + return self.models[0].X.value, np.hstack(map(lambda model: model.Y.value, self.models)) + + @property + def batch_size(self): + return self._batch_size + + def constraint_indices(self): + return np.empty((0,), dtype=int) + + def objective_indices(self): + return np.setdiff1d(np.arange(self.data[1].shape[1]), self.constraint_indices()) + + def feasible_data_index(self): + return np.ones(self.data[0].shape[0], dtype=bool) + + def _setup(self): + pass + + def _setup_constraints(self): + if self.constraint_indices().size > 0: + self._setup() + + def _setup_objectives(self): + if self.constraint_indices().size == 0: + self._setup() + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate_with_gradients(self, Xcand): + acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate(self, Xcand): + return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + + @abstractmethod + def build_acquisition(self, *args): # pragma: no cover + pass + + def __add__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if isinstance(other, AcquisitionSum): + return NotImplemented + + if self.batch_size < other.batch_size: + return ParToSeqAcquisitionWrapper(self, other.batch_size) + other + + if self.batch_size > other.batch_size: + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + return AcquisitionSum([self, other]) + + def __mul__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if isinstance(other, AcquisitionProduct): + return NotImplemented + + if self.batch_size < other.batch_size: + return ParToSeqAcquisitionWrapper(self, other.batch_size) * other + + if self.batch_size > other.batch_size: + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + + return AcquisitionProduct([self, other]) + + def __radd__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if self.batch_size < other.batch_size: + return other + ParToSeqAcquisitionWrapper(self, other.batch_size) + + if self.batch_size > other.batch_size: + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + + return AcquisitionSum([other, self]) + + def __rmul__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if self.batch_size < other.batch_size: + return other * ParToSeqAcquisitionWrapper(self, other.batch_size) + + if self.batch_size > other.batch_size: + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + + return AcquisitionProduct([other, self]) + + +@six.add_metaclass(ABCMeta) +class Acquisition(ParallelBatchAcquisition): + """ + Class to be implemented for standard single-point acquisition functions like standard EI. + + This abstract class inherits the lazy setup mechanism from its ancestor. + """ + + def get_suggestion(self, optimizer): + result = optimizer.optimize(self._inverse_acquisition) + return result.x + + @abstractmethod + def build_acquisition(self, candidates): # pragma: no cover + pass + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate_with_gradients(self, Xcand): + acq = self.build_acquisition(Xcand) + return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate(self, Xcand): + return self.build_acquisition(Xcand) + + def __add__(self, other): + if not isinstance(other, Acquisition): + return NotImplemented + + return AcquisitionSum([self, other]) + + def __mul__(self, other): + if not isinstance(other, Acquisition): + return NotImplemented + + return AcquisitionProduct([self, other]) + + def __setattr__(self, key, value): + super(Acquisition, self).__setattr__(key, value) + if key is '_parent': + self.highest_parent._needs_setup = True + + +class AcquisitionAggregation(AcquisitionWrapper, ParallelBatchAcquisition): + """ + Aggregates multiple acquisition functions, using a TensorFlow reduce operation. + """ + + def __init__(self, operands, oper): + """ + :param operands: list of acquisition objects + :param oper: a tf.reduce operation (e.g., tf.reduce_sum) for aggregating the returned scores of each operand. + """ + super(AcquisitionAggregation, self).__init__(operands) + self._oper = oper + + @ParallelBatchAcquisition.batch_size.getter + def batch_size(self): + return self[0].batch_size + + @property + def operands(self): + return self.wrapped + + @operands.setter + def operands(self, value): + self.wrapped = value + + def build_acquisition(self, *args): + return self._oper(tf.concat(list(map(lambda operand: operand.build_acquisition(*args), self.operands)), 1), + axis=1, keep_dims=True, name=self.__class__.__name__) class AcquisitionSum(AcquisitionAggregation): @@ -373,10 +601,19 @@ def __init__(self, operands): super(AcquisitionSum, self).__init__(operands, tf.reduce_sum) def __add__(self, other): - if isinstance(other, AcquisitionSum): - return AcquisitionSum(self.operands.sorted_params + other.operands.sorted_params) + if self.batch_size == other.batch_size: + if isinstance(other, AcquisitionSum): + return AcquisitionSum(self.operands.sorted_params + other.operands.sorted_params) + else: + return AcquisitionSum(self.operands.sorted_params + [other]) + else: + return super(AcquisitionSum, self).__add__(other) + + def __radd__(self, other): + if self.batch_size == other.batch_size: + return AcquisitionSum([other] + self.operands.sorted_params) else: - return AcquisitionSum(self.operands.sorted_params + [other]) + return super(AcquisitionSum, self).__radd__(other) class AcquisitionProduct(AcquisitionAggregation): @@ -387,10 +624,19 @@ def __init__(self, operands): super(AcquisitionProduct, self).__init__(operands, tf.reduce_prod) def __mul__(self, other): - if isinstance(other, AcquisitionProduct): - return AcquisitionProduct(self.operands.sorted_params + other.operands.sorted_params) + if self.batch_size == other.batch_size: + if isinstance(other, AcquisitionProduct): + return AcquisitionProduct(self.operands.sorted_params + other.operands.sorted_params) + else: + return AcquisitionProduct(self.operands.sorted_params + [other]) + else: + return super(AcquisitionProduct, self).__mul__(other) + + def __rmul__(self, other): + if self.batch_size == other.batch_size: + return AcquisitionProduct([other] + self.operands.sorted_params) else: - return AcquisitionProduct(self.operands.sorted_params + [other]) + return super(AcquisitionProduct, self).__rmul__(other) class MCMCAcquistion(AcquisitionSum): @@ -431,7 +677,7 @@ def _optimize_models(self): for idx, draw in enumerate(self.operands): draw.set_state(hypers[idx, :]) - @Acquisition.models.getter + @ParallelBatchAcquisition.models.getter def models(self): # Only return the models of the first operand, the copies remain hidden. return self.operands[0].models @@ -443,9 +689,9 @@ def set_data(self, X, Y): offset = operand.set_data(X, Y) return offset - def build_acquisition(self, Xcand): + def build_acquisition(self, *args): # Average the predictions of the copies. - return 1. / len(self.operands) * super(MCMCAcquistion, self).build_acquisition(Xcand) + return 1. / len(self.operands) * super(MCMCAcquistion, self).build_acquisition(*args) def _kill_autoflow(self): """ @@ -456,3 +702,18 @@ def _kill_autoflow(self): """ super(MCMCAcquistion, self)._kill_autoflow() self._needs_new_copies = True + + +class ParToSeqAcquisitionWrapper(AcquisitionWrapper, ParallelBatchAcquisition): + + def __init__(self, acquisition, batch_size): + assert isinstance(acquisition, IAcquisition) + super(ParToSeqAcquisitionWrapper, self).__init__([acquisition], batch_size=batch_size) + + def build_acquisition(self, *args): + splits = self.batch_size // self.wrapped[0].batch_size + assert len(args) == splits + all_scores = self.wrapped[0].build_acquisition(tf.concat(args, 0)) + return tf.reduce_prod(tf.concat(tf.split(all_scores, axis=0, num_or_size_splits=splits), 1), axis=1) + + diff --git a/gpflowopt/acquisition/qei.py b/gpflowopt/acquisition/qei.py new file mode 100644 index 0000000..a6560cb --- /dev/null +++ b/gpflowopt/acquisition/qei.py @@ -0,0 +1,85 @@ +# Copyright 2017 Joachim van der Herten +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from .acquisition import ParallelBatchAcquisition + +from gpflow.model import Model +from gpflow.param import DataHolder +from gpflow import settings + +import numpy as np +import tensorflow as tf +import tensorflow.contrib.distributions as ds + +stability = settings.numerics.jitter_level +float_type = settings.dtypes.float_type + + +class qExpectedImprovement(ParallelBatchAcquisition): + def __init__(self, model, batch_size=4): + """ + :param model: GPflow model (single output) representing our belief of the objective + """ + super(qExpectedImprovement, self).__init__(model, batch_size=batch_size) + self.fmin = DataHolder(np.zeros(1)) + self._setup() + + def _setup(self): + super(qExpectedImprovement, self)._setup() + # Obtain the lowest posterior mean for the previous - feasible - evaluations + feasible_samples = self.data[0][self.highest_parent.feasible_data_index(), :] + samples_mean, _ = self.models[0].predict_f(feasible_samples) + self.fmin.set_data(np.min(samples_mean, axis=0)) + + def build_acquisition(self, *args): + # Obtain predictive distributions for candidates + N, D = tf.shape(args[0])[0], tf.shape(args[0])[1] + q = self.batch_size + + Xcand = tf.transpose(tf.stack(args, axis=0), perm=[1, 0, 2]) + m, sig = tf.map_fn(lambda x: self.models[0].build_predict(x, full_cov=True), Xcand, + dtype=(float_type, float_type)) # N x q x 1, N x q x q + + eye = tf.tile(tf.expand_dims(tf.eye(q, dtype=float_type), 0), [q, 1, 1]) + A = eye + A = A - tf.transpose(eye, perm=[1, 2, 0]) + A = A - tf.transpose(eye, perm=[2, 0, 1]) # q x q x q (k x q x q) + + mk = tf.tensordot(A, m, [[2], [1]]) # N x q(k) x q Mean of Zk (k x q) + sigk = tf.tensordot(A, sig, [[2], [1]]) # N x q x q x q + sigk = tf.reduce_sum(tf.expand_dims(sigk, 3) * tf.expand_dims(tf.expand_dims(A, 0), 2), axis=-1)# N x q(k) x q x q + + a = tf.tile(tf.expand_dims(tf.eye(self.batch_size, self.batch_size, dtype=float_type), 0), [self.batch_size, 1, 1]) + A1 = tf.gather_nd(a, [[[i, j + (j >= i)] for j in range(self.batch_size)] for i in range(self.batch_size)]) + # q(i) x (q-1) x q + + bk = -self.fmin * tf.eye(q, dtype=float_type) # q(k) x q + + Sigk_t = tf.expand_dims(tf.transpose(tf.tensordot(sigk, A1, [[-2], [2]]), [0, 1, 3, 4, 2]), -2) # N x q(k) x q(i) x (q-1) x 1 x q + Sigk = tf.reduce_sum(tf.expand_dims(tf.expand_dims(tf.expand_dims(A1, 0), 0), -3)* Sigk_t, axis=-1) + #Sigk = tf.einsum('ijklm,lnk->ijlmn', Sigk_t, A1) # N x q(k) x q(i) x q-1 x q-1 + c = tf.tensordot(tf.expand_dims(bk, 0) - mk, A1, [[2], [2]]) + c = tf.tensordot(tf.expand_dims(bk, 0) - mk, A1, [[2], [2]]) # N x q(k) x q(i) x q-1 + F = tf.expand_dims(tf.expand_dims(tf.expand_dims(bk, 0) - mk, -1) / tf.matrix_diag_part(sigk), + -1) # N x q(k) x q(i) x 1 + F *= tf.transpose(tf.squeeze(tf.matrix_diag_part(tf.transpose(Sigk_t, [0, 1, 3, 4, 5, 2])), 4), [0, 1, 3, 2]) + c -= F + + MVN = ds.MultivariateNormalFullCovariance(loc=mk, covariance_matrix=sigk) + MVN2 = ds.MultivariateNormalFullCovariance(loc=tf.zeros(tf.shape(c), dtype=float_type), covariance_matrix=Sigk) + UVN = ds.MultivariateNormalDiag(loc=mk, scale_diag=tf.sqrt(tf.matrix_diag_part(sigk))) + t1 = tf.reduce_sum((self.fmin - m) * MVN.cdf(bk), axis=1) + sigkk = tf.transpose(tf.matrix_diag_part(tf.transpose(sigk, perm=[0, 3, 1, 2])), perm=[0, 2, 1]) + t2 = tf.reduce_sum(sigkk * UVN.pdf(-tf.expand_dims(bk, 0)) * MVN2.cdf(c), axis=[1, 2]) + return tf.add(t1, t2, name=self.__class__.__name__) diff --git a/gpflowopt/bo.py b/gpflowopt/bo.py index 93c00f9..84f309c 100644 --- a/gpflowopt/bo.py +++ b/gpflowopt/bo.py @@ -19,7 +19,7 @@ import tensorflow as tf from gpflow.gpr import GPR -from .acquisition import Acquisition, MCMCAcquistion +from .acquisition import IAcquisition, MCMCAcquistion from .design import Design, EmptyDesign from .objective import ObjectiveWrapper from .optim import Optimizer, SciPyOptimizer @@ -89,20 +89,24 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr :class:`~.Acquisition` this allows several scenarios: do the optimization manually from the callback (optimize_restarts equals 0), or choose the starting point + some random restarts (optimize_restarts > 0). """ - assert isinstance(acquisition, Acquisition) + assert isinstance(acquisition, IAcquisition) assert hyper_draws is None or hyper_draws > 0 assert optimizer is None or isinstance(optimizer, Optimizer) assert initial is None or isinstance(initial, Design) super(BayesianOptimizer, self).__init__(domain, exclude_gradient=True) + # Configure MCMC self._scaling = scaling if self._scaling: acquisition.enable_scaling(domain) self.acquisition = acquisition if hyper_draws is None else MCMCAcquistion(acquisition, hyper_draws) + # Setup optimizer self.optimizer = optimizer or SciPyOptimizer(domain) self.optimizer.domain = domain + + # Setup initial evaluations initial = initial or EmptyDesign(domain) self.set_initial(initial.generate()) @@ -225,17 +229,14 @@ def _optimize(self, fx, n_iter): # Remove initial design for additional calls to optimize to proceed optimization self.set_initial(EmptyDesign(self.domain).generate()) - def inverse_acquisition(x): - return tuple(map(lambda r: -r, self.acquisition.evaluate_with_gradients(np.atleast_2d(x)))) - # Optimization loop for i in range(n_iter): # If a callback is specified, and acquisition has the setup flag enabled (indicating an upcoming - # compilation), run the callback. + # setup), run the callback. if self._model_callback and self.acquisition._needs_setup: self._model_callback([m.wrapped for m in self.acquisition.models]) - result = self.optimizer.optimize(inverse_acquisition) - self._update_model_data(result.x, fx(result.x)) + Xnew = self.acquisition.get_suggestion(self.optimizer) + self._update_model_data(Xnew, fx(Xnew)) return self._create_bo_result(True, "OK") diff --git a/gpflowopt/domain.py b/gpflowopt/domain.py index f7499be..e77a376 100644 --- a/gpflowopt/domain.py +++ b/gpflowopt/domain.py @@ -13,6 +13,7 @@ # limitations under the License. import numpy as np +import copy from itertools import chain from gpflow.param import Parentable @@ -53,6 +54,9 @@ def size(self): """ return sum(map(lambda param: param.size, self._parameters)) + def batch(self, size): + return np.sum([copy.deepcopy(self) for i in range(size)]) + def __setattr__(self, key, value): super(Domain, self).__setattr__(key, value) if key is not '_parent': diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 5994158..caf3120 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -7,25 +7,43 @@ domain = np.sum([gpflowopt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)]) - +####### TESTING CLASSES ####### class SimpleAcquisition(gpflowopt.acquisition.Acquisition): - def __init__(self, model): - super(SimpleAcquisition, self).__init__(model) + + def build_acquisition(self, Xcand): + return -self.models[0].build_predict(Xcand)[0] + + +class SimpleAcquisitionConstrained(gpflowopt.acquisition.Acquisition): + + def constraint_indices(self): + return np.arange(self.data[1].shape[1]) + + def build_acquisition(self, Xcand): + return -self.models[0].build_predict(Xcand)[0] + + +class SimpleParallelBatch(gpflowopt.acquisition.ParallelBatchAcquisition): + def __init__(self, model, batch_size=2): + super(SimpleParallelBatch, self).__init__(model, batch_size=batch_size) + self.n_args = 0 self.counter = 0 def _setup(self): - super(SimpleAcquisition, self)._setup() + super(SimpleParallelBatch, self)._setup() self.counter += 1 - def build_acquisition(self, Xcand): - return self.models[0].build_predict(Xcand)[0] + def build_acquisition(self, *args): + self.n_args = len(args) + return -tf.add(self.models[0].build_predict(args[0])[0], self.models[0].build_predict(args[1]+1)[0]) +############################## -class TestAcquisition(GPflowOptTestCase): +class TestParallelBatchAcquisition(GPflowOptTestCase): def setUp(self): self.model = create_parabola_model(domain) - self.acquisition = SimpleAcquisition(self.model) + self.acquisition = SimpleParallelBatch(self.model) def test_object_integrity(self): self.assertEqual(len(self.acquisition.models), 1, msg="Model list has incorrect length.") @@ -33,18 +51,20 @@ def test_object_integrity(self): def test_setup_trigger(self): with self.test_session(): + Xrand = np.hstack((gpflowopt.design.RandomDesign(10, domain).generate(), + gpflowopt.design.RandomDesign(10, domain).generate())) m = create_parabola_model(domain) self.assertTrue(np.allclose(m.get_free_state(), self.acquisition.models[0].get_free_state())) self.assertTrue(self.acquisition._needs_setup) self.assertEqual(self.acquisition.counter, 0) - self.acquisition.evaluate(gpflowopt.design.RandomDesign(10, domain).generate()) + self.acquisition.evaluate(Xrand) self.assertFalse(self.acquisition._needs_setup) self.assertEqual(self.acquisition.counter, 1) self.assertFalse(np.allclose(m.get_free_state(), self.acquisition.models[0].get_free_state())) self.acquisition._needs_setup = True self.acquisition.models[0].set_state(m.get_free_state()) - self.acquisition.evaluate_with_gradients(gpflowopt.design.RandomDesign(10, domain).generate()) + self.acquisition.evaluate_with_gradients(Xrand) self.assertFalse(self.acquisition._needs_setup) self.assertEqual(self.acquisition.counter, 2) @@ -88,25 +108,25 @@ def test_enable_scaling(self): def test_result_shape_tf(self): # Verify the returned shape of evaluate - design = gpflowopt.design.RandomDesign(50, domain) with self.test_session(graph=tf.Graph()): free_vars = tf.placeholder(tf.float64, [None]) l = self.acquisition.make_tf_array(free_vars) x_tf = tf.placeholder(tf.float64, shape=(50, 2)) with self.acquisition.tf_mode(): - tens = self.acquisition.build_acquisition(x_tf) + tens = self.acquisition.build_acquisition(x_tf, x_tf) self.assertTrue(isinstance(tens, tf.Tensor), msg="no Tensor was returned") def test_result_shape_np(self): with self.test_session(): design = gpflowopt.design.RandomDesign(50, domain) - res = self.acquisition.evaluate(design.generate()) + candidates = np.hstack((design.generate(), design.generate())) + res = self.acquisition.evaluate(candidates) self.assertTupleEqual(res.shape, (50, 1)) - res = self.acquisition.evaluate_with_gradients(design.generate()) + res = self.acquisition.evaluate_with_gradients(candidates) self.assertTrue(isinstance(res, tuple)) self.assertTrue(len(res), 2) self.assertTupleEqual(res[0].shape, (50, 1)) - self.assertTupleEqual(res[1].shape, (50, domain.size)) + self.assertTupleEqual(res[1].shape, (50, domain.size * 2)) def test_optimize(self): with self.test_session(): @@ -119,6 +139,36 @@ def test_optimize(self): self.acquisition._optimize_models() self.assertFalse(np.allclose(state, self.acquisition.get_free_state())) + def test_arg_splitting(self): + self.acquisition.evaluate(np.random.rand(10, 4)) + self.assertEqual(self.acquisition.n_args, 2) + + self.acquisition.n_args = 0 + X = np.random.rand(10, 4) + A = self.acquisition.evaluate_with_gradients(X) + self.assertEqual(self.acquisition.n_args, 2) + + def test_get_suggestion(self): + opt = gpflowopt.optim.StagedOptimizer([gpflowopt.optim.MCOptimizer(domain, 200), + gpflowopt.optim.SciPyOptimizer(domain)]) + Xnew = self.acquisition.get_suggestion(opt) + self.assertTupleEqual(Xnew.shape, (2, 2)) + self.assertTrue(np.allclose(Xnew[0, :], 0, atol=1e-4)) + self.assertTrue(np.allclose(Xnew[1, :], -1, atol=1e-4)) + + +class TestAcquisition(GPflowOptTestCase): + + def setUp(self): + self.model = create_parabola_model(domain) + self.acquisition = SimpleAcquisition(self.model) + + def test_get_suggestion(self): + opt = gpflowopt.optim.SciPyOptimizer(domain) + Xnew = self.acquisition.get_suggestion(opt) + self.assertTupleEqual(Xnew.shape, (1, 2)) + self.assertTrue(np.allclose(Xnew, 0)) + aggregations = list() aggregations.append(gpflowopt.acquisition.AcquisitionSum([ @@ -134,10 +184,11 @@ def test_optimize(self): ) -class TestAcquisitionAggregation(GPflowOptTestCase): +class TestAcquisitionAggregationImplementations(GPflowOptTestCase): @parameterized.expand(list(zip(aggregations))) def test_object_integrity(self, acquisition): + self.assertTrue(isinstance(acquisition, gpflowopt.acquisition.IAcquisition)) acquisition._kill_autoflow() with self.test_session(): for oper in acquisition.operands: @@ -196,15 +247,6 @@ def test_indices(self, acquisition): np.testing.assert_allclose(acquisition.objective_indices(), np.arange(2, dtype=int)) np.testing.assert_allclose(acquisition.constraint_indices(), np.arange(0, dtype=int)) - def test_generating_operators(self): - joint = gpflowopt.acquisition.ExpectedImprovement(create_parabola_model(domain)) + \ - gpflowopt.acquisition.ExpectedImprovement(create_parabola_model(domain)) - self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) - - joint = gpflowopt.acquisition.ExpectedImprovement(create_parabola_model(domain)) * \ - gpflowopt.acquisition.ExpectedImprovement(create_parabola_model(domain)) - self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) - @parameterized.expand(list(zip([aggregations[2]]))) def test_hyper_updates(self, acquisition): acquisition._kill_autoflow() @@ -247,32 +289,124 @@ def test_mcmc_acq(self): np.testing.assert_almost_equal(ei_mle, ei_mcmc, decimal=5) +class TestParallelToSequentialAdapter(GPflowOptTestCase): + + def test_wrapper(self): + acq = SimpleAcquisition(create_parabola_model(domain)) + adapter = gpflowopt.acquisition.ParToSeqAcquisitionWrapper(acq, 3) + self.assertEqual(acq.batch_size, 1) + self.assertEqual(adapter.batch_size, 3) + + a = np.array([[0.25, -0.25]]) + b = np.array([[0.67, 0]]) + c = np.array([[-0.1, 0.1]]) + + r1 = acq.evaluate(np.vstack((a, b, c))) + r2 = adapter.evaluate(np.hstack((a, b, c))) + self.assertTrue(np.allclose(np.prod(r1), r2)) + + class TestJointAcquisition(GPflowOptTestCase): - def test_constrained_ei(self): - with self.test_session(): - design = gpflowopt.design.LatinHyperCube(16, domain) - X = design.generate() - Yo = parabola2d(X) - Yc = -parabola2d(X) + 0.5 - m1 = gpflow.gpr.GPR(X, Yo, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) - m2 = gpflow.gpr.GPR(X, Yc, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) - ei = gpflowopt.acquisition.ExpectedImprovement(m1) - pof = gpflowopt.acquisition.ProbabilityOfFeasibility(m2) - joint = ei * pof - - # Test output indices - np.testing.assert_allclose(joint.objective_indices(), np.array([0], dtype=int)) - np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int)) + def test_simple_generating_operators(self): + joint = SimpleAcquisition(create_parabola_model(domain)) + SimpleAcquisition(create_parabola_model(domain)) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + + joint = SimpleAcquisition(create_parabola_model(domain)) * SimpleAcquisition(create_parabola_model(domain)) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + + def test_batch_generating_operators(self): + joint = SimpleParallelBatch(create_parabola_model(domain), 2) + \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + + joint = SimpleParallelBatch(create_parabola_model(domain), 2) * \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + + joint = SimpleParallelBatch(create_parabola_model(domain), 6) + \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + joint = SimpleParallelBatch(create_parabola_model(domain), 2) + \ + SimpleParallelBatch(create_parabola_model(domain), 6) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + joint = SimpleParallelBatch(create_parabola_model(domain), 6) * \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + joint = SimpleParallelBatch(create_parabola_model(domain), 2) * \ + SimpleParallelBatch(create_parabola_model(domain), 6) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + def test_mixed_generating_operators(self): + joint = SimpleAcquisition(create_parabola_model(domain)) + SimpleParallelBatch(create_parabola_model(domain), 3) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 3) + self.assertEqual(joint.operands[0].batch_size, 3) + self.assertEqual(joint.operands[1].batch_size, 3) + + joint = SimpleAcquisition(create_parabola_model(domain)) * SimpleParallelBatch(create_parabola_model(domain), 3) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 3) + self.assertEqual(joint.operands[0].batch_size, 3) + self.assertEqual(joint.operands[1].batch_size, 3) + + def test_incorrect_generating_operators(self): + with self.assertRaises(TypeError): + SimpleAcquisition(create_parabola_model(domain)) + 1 + + with self.assertRaises(TypeError): + 1 + SimpleAcquisition(create_parabola_model(domain)) + + with self.assertRaises(TypeError): + SimpleAcquisition(create_parabola_model(domain)) * 1 - # Test proper setup - joint._optimize_models() - joint._setup() - self.assertGreater(ei.fmin.value, np.min(ei.data[1]), msg="The best objective value is in an infeasible area") - self.assertTrue(np.allclose(ei.fmin.value, np.min(ei.data[1][pof.feasible_data_index(), :]), atol=1e-3), - msg="fmin computed incorrectly") + with self.assertRaises(TypeError): + 1 * SimpleAcquisition(create_parabola_model(domain)) - def test_hierarchy(self): + with self.assertRaises(TypeError): + SimpleParallelBatch(create_parabola_model(domain), 2) + 1 + + with self.assertRaises(TypeError): + 1 + SimpleParallelBatch(create_parabola_model(domain), 2) + + with self.assertRaises(TypeError): + SimpleParallelBatch(create_parabola_model(domain), 2) * 1 + + with self.assertRaises(TypeError): + 1 * SimpleParallelBatch(create_parabola_model(domain), 2) + + def test_simple_hierarchy(self): with self.test_session(): design = gpflowopt.design.LatinHyperCube(16, domain) X = design.generate() @@ -280,16 +414,16 @@ def test_hierarchy(self): m1 = create_parabola_model(domain, design) m2 = create_parabola_model(domain, design) m3 = gpflow.gpr.GPR(X, Yc, gpflow.kernels.RBF(2, ARD=True)) - joint = gpflowopt.acquisition.ExpectedImprovement(m1) * \ - (gpflowopt.acquisition.ProbabilityOfFeasibility(m3) - + gpflowopt.acquisition.ExpectedImprovement(m2)) + joint = SimpleAcquisition(m1) * (SimpleAcquisitionConstrained(m3) + SimpleAcquisition(m2)) np.testing.assert_allclose(joint.objective_indices(), np.array([0, 2], dtype=int)) np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int)) - def test_multi_aggr(self): + @parameterized.expand([([SimpleAcquisition(create_parabola_model(domain)) for i in range(4)],), + ([SimpleParallelBatch(create_parabola_model(domain)) for i in range(4)],)]) + def test_multi_aggr_simple(self, acq): + # In these tests, the batch sizes align with self.test_session(): - acq = [gpflowopt.acquisition.ExpectedImprovement(create_parabola_model(domain)) for i in range(4)] acq1, acq2, acq3, acq4 = acq joint = acq1 + acq2 + acq3 self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) @@ -325,6 +459,45 @@ def test_multi_aggr(self): self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3, acq4]) + def test_multi_aggr_mixed(self): + # In these tests, batch sizes do not align + with self.test_session(): + acq1 = SimpleAcquisition(create_parabola_model(domain)) + acq2 = SimpleAcquisition(create_parabola_model(domain)) + acq3 = SimpleParallelBatch(create_parabola_model(domain)) + + first = acq1 + acq2 + joint = first + acq3 + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[1], acq3) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[0].wrapped[0], first) + + joint = acq3 + first + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[0], acq3) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[1].wrapped[0], first) + + first = acq1 * acq2 + joint = first * acq3 + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[1], acq3) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[0].wrapped[0], first) + + joint = acq3 * first + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[0], acq3) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[1].wrapped[0], first) + + + class TestRecompile(GPflowOptTestCase): """ diff --git a/testing/test_domain.py b/testing/test_domain.py index 68b3e1b..a041173 100644 --- a/testing/test_domain.py +++ b/testing/test_domain.py @@ -68,6 +68,13 @@ def test_value(self): self.assertTupleEqual(p.value.shape, (1,), msg="Default value has incorrect shape.") self.assertTrue(np.allclose(p.value, 0.2), msg="Parameter has incorrect initialized value") + def test_batch(self): + p = gpflowopt.domain.ContinuousParameter("x1", 0, 1) + b = p.batch(3) + self.assertEqual(b.size, 3) + self.assertTrue(np.allclose(b.lower, 0)) + self.assertTrue(np.allclose(b.upper, 1)) + class TestHypercubeDomain(GPflowOptTestCase): diff --git a/testing/test_implementations.py b/testing/test_implementations.py index 77a5d9a..20f249e 100644 --- a/testing/test_implementations.py +++ b/testing/test_implementations.py @@ -1,4 +1,5 @@ import gpflowopt +import gpflow import numpy as np from parameterized import parameterized from .utility import create_parabola_model, create_plane_model, create_vlmop2_model, parabola2d, load_data, GPflowOptTestCase @@ -189,6 +190,32 @@ def test_HvPoI_validity(self): np.testing.assert_almost_equal(scores, self.data['scores'], decimal=2) +class TestConstrainedExpectedImprovement(GPflowOptTestCase): + + def test_constrained_ei(self): + with self.test_session(): + design = gpflowopt.design.LatinHyperCube(16, domain) + X = design.generate() + Yo = parabola2d(X) + Yc = -parabola2d(X) + 0.5 + m1 = gpflow.gpr.GPR(X, Yo, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m2 = gpflow.gpr.GPR(X, Yc, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + ei = gpflowopt.acquisition.ExpectedImprovement(m1) + pof = gpflowopt.acquisition.ProbabilityOfFeasibility(m2) + joint = ei * pof + + # Test output indices + np.testing.assert_allclose(joint.objective_indices(), np.array([0], dtype=int)) + np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int)) + + # Test proper setup + joint._optimize_models() + joint._setup() + self.assertGreater(ei.fmin.value, np.min(ei.data[1]), msg="The best objective value is in an infeasible area") + self.assertTrue(np.allclose(ei.fmin.value, np.min(ei.data[1][pof.feasible_data_index(), :]), atol=1e-3), + msg="fmin computed incorrectly") + + class TestMinValueEntropySearch(GPflowOptTestCase): def setUp(self): super(TestMinValueEntropySearch, self).setUp()