Merge pull request #57 from GPflow/enhance_optimizers

Optimizer cleanup

Author: icouckuy

GPflowOpt/optim.py

@@ -12,12 +12,14 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import numpy as np
-from scipy.optimize import OptimizeResult, minimize
-from GPflow import settings
 import contextlib
-import sys
 import os
+import sys
+import warnings
+
+import numpy as np
+from GPflow import settings
+from scipy.optimize import OptimizeResult, minimize
 
 from .design import RandomDesign
 from .objective import ObjectiveWrapper
@@ -126,47 +128,31 @@ def silent(self):
         sys.stdout = save_stdout
 
 
-class CandidateOptimizer(Optimizer):
+class MCOptimizer(Optimizer):
     """
-    Optimization of an objective function by evaluating a set of pre-defined candidate points.
-
-    Returns the point with minimal objective value.
+    Optimization of an objective function by evaluating a set of random points.
 
-    For compatibility with the StagedOptimizer, the candidate points are concatenated with
-    the initial points and evaluated.
+    Note: each call to optimize, a different set of random points is evaluated.
     """
 
-    def __init__(self, domain, candidates, batch=False):
+    def __init__(self, domain, nsamples):
         """
-        :param domain: Optimization :class:`.domain.Domain`.
-        :param candidates: candidate points, should be within the optimization domain. 
-        :param batch: bool, if true evaluates the objective function on all points at once
+        :param domain: Optimization :class:`~.domain.Domain`.
+        :param nsamples: number of random points to use
         """
-        super(CandidateOptimizer, self).__init__(domain, exclude_gradient=True)
-        assert(candidates in domain)
-        self.candidates = candidates
-        self._batch_mode = batch
+        super(MCOptimizer, self).__init__(domain, exclude_gradient=True)
+        self._nsamples = nsamples
 
     @Optimizer.domain.setter
     def domain(self, dom):
-        # Attempt to transform candidates
-        t = self.domain >> dom
-        self.candidates = t.forward(self.candidates)
         self._domain = dom
-        self.set_initial(dom.value)
 
-    def get_initial(self):
-        return np.vstack((super(CandidateOptimizer, self).get_initial(), self.candidates))
-
-    def _evaluate_one_by_one(self, objective, X):
-        """
-        Evaluates each row of X individually.
-        """
-        return np.vstack(map(lambda x: objective(x), X))
+    def _get_eval_points(self):
+        return RandomDesign(self._nsamples, self.domain).generate()
 
     def _optimize(self, objective):
-        points = self.get_initial()
-        evaluations = objective(points) if self._batch_mode else self._evaluate_one_by_one(objective, points)
+        points = self._get_eval_points()
+        evaluations = objective(points)
         idx_best = np.argmin(evaluations, axis=0)
 
         return OptimizeResult(x=points[idx_best, :],
@@ -175,21 +161,39 @@ def _optimize(self, objective):
                               nfev=points.shape[0],
                               message="OK")
 
+    def set_initial(self, initial):
+        initial = np.atleast_2d(initial)
+        super(MCOptimizer, self).set_initial(initial)
+        if initial.size > 0:
+            warnings.warn("Initial points set in {0} are ignored.".format(self.__class__.__name__), UserWarning)
+
 
-class MCOptimizer(CandidateOptimizer):
+class CandidateOptimizer(MCOptimizer):
     """
-    Optimization of an objective function by evaluating a set of random points.
+    Optimization of an objective function by evaluating a set of pre-defined candidate points.
 
-    Note: each call to optimize, a different set of random points is evaluated.
+    Returns the point with minimal objective value.
     """
 
-    def __init__(self, domain, nsamples, batch=False):
-        super(MCOptimizer, self).__init__(domain, np.empty((0, domain.size)), batch=batch)
-        self._nsamples = nsamples
+    def __init__(self, domain, candidates):
+        """
+        :param domain: Optimization :class:`~.domain.Domain`.
+        :param candidates: candidate points, should be within the optimization domain.
+        """
+        super(CandidateOptimizer, self).__init__(domain, candidates.shape[0])
+        assert (candidates in domain)
+        self.candidates = candidates
+        # Clear the initial data points
+        self.set_initial(np.empty((0, self.domain.size)))
 
-    def _optimize(self, objective):
-        self.candidates = RandomDesign(self._nsamples, self.domain).generate()
-        return super(MCOptimizer, self)._optimize(objective)
+    def _get_eval_points(self):
+        return self.candidates
+
+    @MCOptimizer.domain.setter
+    def domain(self, dom):
+        t = self.domain >> dom
+        super(CandidateOptimizer, self.__class__).domain.fset(self, dom)
+        self.candidates = t.forward(self.candidates)
 
 
 class SciPyOptimizer(Optimizer):
@@ -231,6 +235,17 @@ def __init__(self, optimizers):
         no_gradient = any(map(lambda opt: not opt.gradient_enabled(), optimizers))
         super(StagedOptimizer, self).__init__(optimizers[0].domain, exclude_gradient=no_gradient)
         self.optimizers = optimizers
+        self.set_initial(np.empty((0, self.domain.size)))
+
+    @Optimizer.domain.setter
+    def domain(self, domain):
+        super(StagedOptimizer, self.__class__).domain.fset(self, domain)
+        for optimizer in self.optimizers:
+            optimizer.domain = domain
+
+    def _best_x(self, results):
+        best_idx = np.argmin([r.fun for r in results if r.success])
+        return results[best_idx].x, results[best_idx].fun
 
     def optimize(self, objectivefx):
         """
@@ -240,18 +255,21 @@ def optimize(self, objectivefx):
         """
 
         self.optimizers[0].set_initial(self.get_initial())
-        fun_evals = []
-        for current, next in zip(self.optimizers[:-1], self.optimizers[1:]):
+        results = []
+        for current, following in zip(self.optimizers[:-1], self.optimizers[1:]):
             result = current.optimize(objectivefx)
-            fun_evals.append(result.nfev)
+            results.append(result)
             if not result.success:
                 result.message += " StagedOptimizer interrupted after {0}.".format(current.__class__.__name__)
                 break
-            next.set_initial(result.x)
+            following.set_initial(self._best_x(results)[0])
 
         if result.success:
             result = self.optimizers[-1].optimize(objectivefx)
-            fun_evals.append(result.nfev)
-        result.nfev = sum(fun_evals)
-        result.nstages = len(fun_evals)
+            results.append(result)
+
+        result.nfev = sum(r.nfev for r in results)
+        result.nstages = len(results)
+        if any(r.success for r in results):
+            result.x, result.fun = self._best_x(results)
         return result

testing/test_optimizers.py

@@ -5,6 +5,7 @@
 import six
 import sys
 import os
+import warnings
 from contextlib import contextmanager
 from scipy.optimize import OptimizeResult
 
@@ -31,18 +32,20 @@ def __init__(self, iters_to_raise, f):
         self.iters_to_raise, self.f = iters_to_raise, f
         self.count = 0
 
-    def __call__(self, *a, **kw):
-        self.count += 1
+    def __call__(self, X):
         if self.count >= self.iters_to_raise:
             raise KeyboardInterrupt
-        return self.f(*a, **kw)
+        val = self.f(X)
+        self.count += X.shape[0]
+        return val
 
 
 class _TestOptimizer(object):
     _multiprocess_can_split_ = True
 
     def setUp(self):
         self.optimizer = None
+        warnings.simplefilter("once", category=UserWarning)
 
     @property
     def domain(self):
@@ -71,9 +74,21 @@ def setUp(self):
         design = GPflowOpt.design.FactorialDesign(4, self.domain)
         self.optimizer = GPflowOpt.optim.CandidateOptimizer(self.domain, design.generate())
 
+    def test_default_initial(self):
+        self.assertTupleEqual(self.optimizer._initial.shape, (0, 2), msg="Invalid shape of initial points array")
+
+    def test_set_initial(self):
+        # When run separately this test works, however when calling nose to run all tests on python 2.7 this records
+        # no warnings
+        with warnings.catch_warnings(record=True) as w:
+            super(TestCandidateOptimizer, self).test_set_initial()
+            assert len(w) == 1
+            assert issubclass(w[-1].category, UserWarning)
+
     def test_object_integrity(self):
         self.assertTupleEqual(self.optimizer.candidates.shape, (16, 2), msg="Invalid shape of candidate property.")
-        self.assertTupleEqual(self.optimizer.get_initial().shape, (17, 2), msg="Invalid shape of initial points")
+        self.assertTupleEqual(self.optimizer._get_eval_points().shape, (16, 2))
+        self.assertTupleEqual(self.optimizer.get_initial().shape, (0, 2), msg="Invalid shape of initial points")
         self.assertFalse(self.optimizer.gradient_enabled(), msg="CandidateOptimizer supports no gradients.")
 
     def test_set_domain(self):
@@ -83,17 +98,17 @@ def test_set_domain(self):
         self.assertNotEqual(self.optimizer.domain, self.domain)
         self.assertEqual(self.optimizer.domain, GPflowOpt.domain.UnitCube(2))
         rescaled_candidates = GPflowOpt.design.FactorialDesign(4, GPflowOpt.domain.UnitCube(2)).generate()
-        self.assertTrue(np.allclose(self.optimizer.get_initial(), np.vstack((0.5*np.ones((1,2)), rescaled_candidates))))
+        self.assertTrue(np.allclose(self.optimizer.candidates, rescaled_candidates))
 
     def test_optimize(self):
+        self.optimizer.candidates = np.vstack((self.optimizer.candidates, np.zeros((1,2))))
         result = self.optimizer.optimize(parabola2d)
         self.assertTrue(result.success, msg="Optimization should succeed.")
         self.assertTrue(np.allclose(result.x, 0), msg="Optimum should be identified")
         self.assertTrue(np.allclose(result.fun, 0), msg="Function value in optimum is 0")
         self.assertEqual(result.nfev, 17, msg="Number of function evaluations equals candidates + initial points")
 
     def test_optimize_second(self):
-        self.optimizer.set_initial([0.67, 0.67])
         result = self.optimizer.optimize(parabola2d)
         self.assertGreater(result.fun, 0, msg="Optimum is not amongst candidates and initial points")
         self.assertLess(result.fun, 2, msg="Function value not reachable within domain")
@@ -131,33 +146,51 @@ class TestStagedOptimizer(_TestOptimizer, unittest.TestCase):
     def setUp(self):
         super(TestStagedOptimizer, self).setUp()
         self.optimizer = GPflowOpt.optim.StagedOptimizer([GPflowOpt.optim.MCOptimizer(self.domain, 5),
+                                                          GPflowOpt.optim.MCOptimizer(self.domain, 5),
                                                           GPflowOpt.optim.SciPyOptimizer(self.domain, maxiter=10)])
 
+    def test_default_initial(self):
+        self.assertTupleEqual(self.optimizer._initial.shape, (0,2))
+
     def test_object_integrity(self):
-        self.assertEqual(len(self.optimizer.optimizers), 2, msg="Two optimizers expected in optimizerlist")
+        self.assertEqual(len(self.optimizer.optimizers), 3, msg="Two optimizers expected in optimizerlist")
         self.assertFalse(self.optimizer.gradient_enabled(), msg="MCOptimizer supports no gradients => neither "
                                                                 "does stagedoptimizer.")
 
     def test_optimize(self):
-        result = self.optimizer.optimize(parabola2d)
-        self.assertTrue(result.success)
-        self.assertLessEqual(result.nfev, 20, "Only 10 Iterations permitted")
-        self.assertTrue(np.allclose(result.x, 0), msg="Optimizer failed to find optimum")
-        self.assertTrue(np.allclose(result.fun, 0), msg="Incorrect function value returned")
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", category=UserWarning)
+            result = self.optimizer.optimize(parabola2d)
+            self.assertTrue(result.success)
+            self.assertLessEqual(result.nfev, 20, "Only 20 Iterations permitted")
+            self.assertTrue(np.allclose(result.x, 0), msg="Optimizer failed to find optimum")
+            self.assertTrue(np.allclose(result.fun, 0), msg="Incorrect function value returned")
 
     def test_optimizer_interrupt(self):
-        self.optimizer.set_initial([-1, -1])
-        result = self.optimizer.optimize(KeyboardRaiser(3, parabola2d))
-        self.assertFalse(result.success, msg="After two evaluations, a keyboard interrupt is raised, "
-                                             "non-succesfull result expected.")
-        self.assertFalse(np.allclose(result.x, 0.0), msg="After one iteration, the optimum will not be found")
-        self.assertEqual(result.nstages, 1, msg="Stage 1 should be in progress during interrupt")
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", category=UserWarning)
+            result = self.optimizer.optimize(KeyboardRaiser(0, parabola2d))
+            self.assertFalse(result.success, msg="non-succesfull result expected.")
+            self.assertEqual(result.nstages, 1, msg="Stage 2 should be in progress during interrupt")
+            self.assertEqual(result.nfev, 0)
+
+            result = self.optimizer.optimize(KeyboardRaiser(3, parabola2d))
+            self.assertFalse(result.success, msg="non-succesfull result expected.")
+            self.assertFalse(np.allclose(result.x, 0.0), msg="The optimum will not be found")
+            self.assertEqual(result.nstages, 2, msg="Stage 2 should be in progress during interrupt")
+            self.assertEqual(result.nfev, 5)
+
+            result = self.optimizer.optimize(KeyboardRaiser(12, parabola2d))
+            print(result)
+            self.assertFalse(result.success, msg="non-succesfull result expected.")
+            self.assertEqual(result.nfev, 12)
+            self.assertFalse(np.allclose(result.x[0, :], 0.0), msg="The optimum should not be found yet")
+            self.assertEqual(result.nstages, 3, msg="Stage 3 should be in progress during interrupt")
 
-        result = self.optimizer.optimize(KeyboardRaiser(8, parabola2d))
-        self.assertFalse(result.success, msg="After 7 evaluations, a keyboard interrupt is raised, "
-                                             "non-succesfull result expected.")
-        self.assertFalse(np.allclose(result.x[0, :], 0.0), msg="The optimum should not be found yet")
-        self.assertEqual(result.nstages, 2, msg="Stage 2 should be in progress during interrupt")
+    def test_set_domain(self):
+        super(TestStagedOptimizer, self).test_set_domain()
+        for opt in self.optimizer.optimizers:
+            self.assertEqual(opt.domain, GPflowOpt.domain.UnitCube(3))
 
 
 class TestBayesianOptimizer(_TestOptimizer, unittest.TestCase):