Initialization from demo

.gitignore

@@ -87,3 +87,6 @@ ENV/
 
 # Rope project settings
 .ropeproject
+
+# Pycharm IDE directory
+.idea

.travis.yml

@@ -9,6 +9,6 @@ install:
   - pip install tensorflow==1.0.1 numpy scipy nose pandas codecov coverage
   - python setup.py install
 script:
-  - nosetests --nologcapture --with-coverage --cover-package=testing testing
+  - nosetests --nologcapture --with-coverage --cover-package=GPflowOpt testing
 after_success:
   - codecov

GPflowOpt/__init__.py

@@ -0,0 +1,19 @@
+# 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 . import acquisition
+from . import domain
+from .bo import BayesianOptimizer
+from . import optim
+from . import design

GPflowOpt/acquisition.py

@@ -0,0 +1,214 @@
+# 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 GPflow.param import Parameterized, AutoFlow, ParamList, DataHolder
+from GPflow.model import Model
+from GPflow._settings import settings
+
+import numpy as np
+
+import tensorflow as tf
+
+float_type = settings.dtypes.float_type
+stability = settings.numerics.jitter_level
+
+
+class Acquisition(Parameterized):
+    """
+    Class representing acquisition functions which map the predictive distribution of a Bayesian model into a score. In
+    Bayesian Optimization this function is typically optimized over the optimization domain to determine the next
+    point for evaluation.
+
+    The acquisition function object maintains a list of models. Subclasses should implement a build_acquisition
+    function which computes the acquisition function using tensorflow.
+
+    Acquisition functions can be added or multiplied to construct joint criteria (for instance for constrained
+    optimization)
+    """
+
+    def __init__(self, models=[], optimize_restarts=5):
+        super(Acquisition, self).__init__()
+        self.models = ParamList(np.atleast_1d(models).tolist())
+        self._default_params = map(lambda m: m.get_free_state(), self.models)
+
+        assert(optimize_restarts >= 0)
+        self._optimize_restarts = optimize_restarts
+        self._optimize_all()
+
+    def _optimize_all(self):
+        for model, hypers in zip(self.models, self._default_params):
+            runs = []
+            # Start from supplied hyperparameters
+            model.set_state(hypers)
+            for i in range(self._optimize_restarts):
+                if i > 0:
+                    model.randomize()
+                try:
+                    result = model.optimize()
+                    runs.append(result)
+                except tf.errors.InvalidArgumentError:
+                    print("Warning - optimization restart {0}/{1} failed".format(i + 1, self._optimize_restarts))
+            best_idx = np.argmin(map(lambda r: r.fun, runs))
+            model.set_state(runs[best_idx].x)
+
+    def _build_acquisition_wrapper(self, Xcand, gradients=True):
+        acq = self.build_acquisition(Xcand)
+        if gradients:
+            return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0]
+        else:
+            return acq
+
+    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._optimize_all()
+        self.setup()
+        return num_outputs_sum
+
+    @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))
+
+    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 setup(self):
+        """
+        Method triggered after calling set_data(). Override for pre-calculation of quantities used later in
+        the evaluation of the acquisition function for candidate points
+        """
+        pass
+
+    @AutoFlow((float_type, [None, None]))
+    def evaluate_with_gradients(self, Xcand):
+        return self._build_acquisition_wrapper(Xcand, gradients=True)
+
+    @AutoFlow((float_type, [None, None]))
+    def evaluate(self, Xcand):
+        return self._build_acquisition_wrapper(Xcand, gradients=False)
+
+    def __add__(self, other):
+        return AcquisitionSum(self, other)
+
+    def __mul__(self, other):
+        return AcquisitionProduct(self, other)
+
+
+class ExpectedImprovement(Acquisition):
+    """
+    Implementation of the Expected Improvement acquisition function for single-objective global optimization
+    (Mockus et al, 1975)
+    """
+
+    def __init__(self, model):
+        super(ExpectedImprovement, self).__init__(model)
+        assert (isinstance(model, Model))
+        self.fmin = DataHolder(np.zeros(1))
+        self.setup()
+
+    def setup(self):
+        super(ExpectedImprovement, self).setup()
+        # Obtain the lowest posterior mean for the previous evaluations
+        samples_mean, _ = self.models[0].predict_f(self.data[0])
+        self.fmin.set_data(np.min(samples_mean, axis=0))
+        # samples_mean, _ = self.models[0].build_predict(self.data[0])
+        # self.fmin = tf.reduce_min(samples_mean, axis=0)
+
+    def build_acquisition(self, Xcand):
+        # Obtain predictive distributions for candidates
+        candidate_mean, candidate_var = self.models[0].build_predict(Xcand)
+        candidate_var = tf.maximum(candidate_var, stability)
+
+        # Compute EI
+        normal = tf.contrib.distributions.Normal(candidate_mean, tf.sqrt(candidate_var))
+        t1 = (self.fmin - candidate_mean) * normal.cdf(self.fmin)
+        t2 = candidate_var * normal.pdf(self.fmin)
+        return tf.add(t1, t2, name=self.__class__.__name__)
+
+
+class ProbabilityOfFeasibility(Acquisition):
+    """
+    Probability of Feasibility acquisition function for learning  feasible regions
+    """
+
+    def __init__(self, model):
+        super(ProbabilityOfFeasibility, self).__init__(model)
+
+    def constraint_indices(self):
+        return np.arange(self.data[1].shape[1])
+
+    def build_acquisition(self, Xcand):
+        candidate_mean, candidate_var = self.models[0].build_predict(Xcand)
+        candidate_var = tf.maximum(candidate_var, stability)
+        normal = tf.contrib.distributions.Normal(candidate_mean, tf.sqrt(candidate_var))
+        return normal.cdf(tf.constant(0.0, dtype=float_type), name=self.__class__.__name__)
+
+
+class AcquisitionBinaryOperator(Acquisition):
+    """
+    Base class for implementing binary operators for acquisition functions, for defining joint acquisition functions
+    """
+
+    def __init__(self, lhs, rhs, oper):
+        super(AcquisitionBinaryOperator, self).__init__()
+        assert isinstance(lhs, Acquisition)
+        assert isinstance(rhs, Acquisition)
+        self.lhs = lhs
+        self.rhs = rhs
+        self._oper = oper
+
+    @Acquisition.data.getter
+    def data(self):
+        lhsX, lhsY = self.lhs.data
+        rhsX, rhsY = self.rhs.data
+        if self._tf_mode:
+            return lhsX, tf.concat((lhsY, rhsY), 1)
+        else:
+            return lhsX, np.hstack((lhsY, rhsY))
+
+    def set_data(self, X, Y):
+        offset_lhs = self.lhs.set_data(X, Y)
+        offset_rhs = self.rhs.set_data(X, Y[:, offset_lhs:])
+        return offset_lhs + offset_rhs
+
+    def constraint_indices(self):
+        offset = self.lhs.data[1].shape[1]
+        return np.hstack((self.lhs.constraint_indices(), offset + self.rhs.constraint_indices()))
+
+    def build_acquisition(self, Xcand):
+        return self._oper(self.lhs.build_acquisition(Xcand), self.rhs.build_acquisition(Xcand),
+                          name=self.__class__.__name__)
+
+
+class AcquisitionSum(AcquisitionBinaryOperator):
+    def __init__(self, lhs, rhs):
+        super(AcquisitionSum, self).__init__(lhs, rhs, tf.add)
+
+
+class AcquisitionProduct(AcquisitionBinaryOperator):
+    def __init__(self, lhs, rhs):
+        super(AcquisitionProduct, self).__init__(lhs, rhs, tf.multiply)

GPflowOpt/bo.py

@@ -0,0 +1,135 @@
+# 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.
+
+import numpy as np
+from scipy.optimize import OptimizeResult
+
+from .acquisition import Acquisition
+from .optim import Optimizer, SciPyOptimizer
+from .design import EmptyDesign
+
+
+class BayesianOptimizer(Optimizer):
+    """
+    Specific optimizer representing bayesian optimization. Takes a domain, an acquisition function and an optimization
+    strategy for the acquisition function.
+    """
+
+    def __init__(self, domain, acquisition, optimizer=None, initial=None):
+        assert isinstance(acquisition, Acquisition)
+        super(BayesianOptimizer, self).__init__(domain, exclude_gradient=True)
+        self.acquisition = acquisition
+        if initial is None:
+            initial = EmptyDesign(domain)
+        self.set_initial(initial.generate())
+        self.optimizer = SciPyOptimizer(domain) if optimizer is None else optimizer
+
+    def _update_model_data(self, newX, newY):
+        """
+        Update the underlying models of the acquisition function with new data
+        :param newX: samples (# new samples x indim)
+        :param newY: values obtained by evaluating the objective and constraint functions (# new samples x # targets)
+        """
+        assert self.acquisition.data[0].shape[1] == newX.shape[-1]
+        assert self.acquisition.data[1].shape[1] == newY.shape[-1]
+        assert newX.shape[0] == newY.shape[0]
+        X = np.vstack((self.acquisition.data[0], newX))
+        Y = np.vstack((self.acquisition.data[1], newY))
+        self.acquisition.set_data(X, Y)
+
+    def _acquisition_wrapper(self, x):
+        """
+        Negation of the acquisition score/gradient
+        :param x: candidate points (# candidates x indim)
+        :return: negative score and gradient (maximizing acquisition function vs minimization algorithm)
+        """
+        scores, grad = self.acquisition.evaluate_with_gradients(np.atleast_2d(x))
+        return -scores, -grad
+
+    def _evaluate(self, X, fxs):
+        fxs = np.atleast_1d(fxs)
+        if X.size > 0:
+            evaluations = np.hstack(map(lambda f: f(X), fxs))
+            assert evaluations.shape[1] == self.acquisition.data[1].shape[1]
+            return evaluations
+        else:
+            return np.empty((0, self.acquisition.data[1].shape[1]))
+
+    def _create_result(self, success, message):
+        """
+        Given all data evaluated after the optimization, analyse and return an OptimizeResult
+        :param success: Optimization successfull? (True/False)
+        :param message: return message
+        :return: OptimizeResult object
+        """
+        X, Y = self.acquisition.data
+
+        # Filter on constraints
+        # Extend with valid column - in case of no constrains
+        Yext = np.hstack((-np.ones((Y.shape[0], 1)), Y[:, self.acquisition.constraint_indices()]))
+        valid = np.all(Yext < 0, axis=1)
+
+        if not np.any(valid):
+            return OptimizeResult(success=False,
+                                  message="No evaluations satisfied the constraints")
+
+        valid_X = X[valid, :]
+        valid_Y = Y[valid, :]
+        valid_Yo = valid_Y[:, self.acquisition.objective_indices()]
+
+        # Here is the place to plug in pareto front if valid_Y.shape[1] > 1
+        # else
+        idx = np.argmin(valid_Yo)
+
+        return OptimizeResult(x=valid_X[idx, :],
+                              success=success,
+                              fun=valid_Yo[idx, :],
+                              message=message)
+
+    def optimize(self, objectivefx, n_iter=20):
+        """
+        Run Bayesian optimization for a number of iterations. Each iteration a new data point is selected by optimizing
+        an acquisition function. This point is evaluated on the objective and black-box constraints. This retrieved
+        information then updates the model
+        :param objectivefx: (list of) expensive black-box objective and constraint functions
+        :param n_iter: number of iterations to run
+        :return: OptimizeResult object
+        """
+
+        # Bayesian optimization is Gradient-free: provide wrapper. Also anticipate for lists and pass on a function
+        # which satisfies the optimizer.optimize interface
+        fx = lambda X: (self._evaluate(X, objectivefx), np.zeros((X.shape[0], 0)))
+        return super(BayesianOptimizer, self).optimize(fx, n_iter=n_iter)
+
+    def _optimize(self, fx, n_iter):
+        """
+        Internal optimization function. Receives and ObjectiveWrapper
+        :param fx: ObjectiveWrapper object wrapping expensive black-box objective and constraint functions
+        :param n_iter: number of iterations to run
+        :return: OptimizeResult object
+        """
+
+        # Evaluate and add the initial design (if any)
+        initial = self.get_initial()
+        values = fx(initial)
+        self._update_model_data(initial, values)
+        # Remove initial design for additional calls to optimize to proceed optimization
+        self.set_initial(EmptyDesign(self.domain).generate())
+
+        # Optimization loop
+        for i in range(n_iter):
+            result = self.optimizer.optimize(self._acquisition_wrapper)
+            self._update_model_data(result.x, fx(result.x))
+
+        return self._create_result(True, "OK")