added COBYLA (no REST interface yet)

moe/optimal_learning/python/cpp_wrappers/domain.py

@@ -75,6 +75,14 @@ def get_bounding_box(self):
         """Return a list of ClosedIntervals representing a bounding box for this domain."""
         return copy.copy(self._domain_bounds)
 
+    def get_constraint_list(self, start_index=0):
+        """Return a list of lambda functions expressing the domain bounds as linear constraints. Used by COBYLA."""
+        constraints = []
+        for i, interval in enumerate(self._domain_bounds):
+            constraints.append((lambda x: x[i + start_index] - interval.min))
+            constraints.append((lambda x: interval.max - x[i + start_index]))
+        return constraints
+
     def generate_random_point_in_domain(self, random_source=None):
         """Generate ``point`` uniformly at random such that ``self.check_point_inside(point)`` is True.
 

moe/optimal_learning/python/interfaces/domain_interface.py

@@ -31,6 +31,11 @@ def get_bounding_box(self):
         """Return a list of ClosedIntervals representing a bounding box for this domain."""
         pass
 
+    @abstractmethod
+    def get_constraint_list(self, start_index=0):
+        """Return a list of lambda functions expressing the domain bounds as linear constraints. Used by COBYLA."""
+        pass
+
     @abstractmethod
     def generate_random_point_in_domain(self, random_source=None):
         """Generate ``point`` uniformly at random such that ``self.check_point_inside(point)`` is True.

moe/optimal_learning/python/python_version/domain.py

@@ -85,6 +85,25 @@ def get_bounding_box(self):
         """Return a list of ClosedIntervals representing a bounding box for this domain."""
         return copy.copy(self._domain_bounds)
 
+    def get_constraint_list(self, start_index=0):
+        """Return a list of lambda functions expressing the domain bounds as linear constraints. Used by COBYLA.
+
+        Since COBYLA in scipy only optimizes arrays, we flatten out our points while doing multipoint EI optimization.
+        But in order for the constraints to access the correct index, the RepeatedDomain class has to signal which index
+        the TensorProductDomain should start from, using the start_index optional parameter.
+
+        :param start_index: the dimension this tensor product domain should start indexing from
+        :type start_index: integer (>= 0)
+        :return: a list of lambda functions corresponding to constraints
+        :rtype: array of lambda functions with shape (dim * 2)
+
+        """
+        constraints = []
+        for i, interval in enumerate(self._domain_bounds):
+            constraints.append((lambda x: x[i + start_index] - interval.min))
+            constraints.append((lambda x: interval.max - x[i + start_index]))
+        return constraints
+
     def generate_random_point_in_domain(self, random_source=None):
         """Generate ``point`` uniformly at random such that ``self.check_point_inside(point)`` is True.
 

moe/optimal_learning/python/python_version/optimization.py

@@ -324,6 +324,32 @@ class LBFGSBParameters(_BaseLBFGSBParameters):
     __slots__ = ()
 
 
+# See ConstrainedDFOParameters (below) for docstring.
+_BaseConstrainedDFOParameters = collections.namedtuple('_BaseConstrainedDFOParameters', [
+    'rhobeg',
+    'rhoend',
+    'maxfun',
+    'catol',
+])
+
+
+class ConstrainedDFOParameters(_BaseConstrainedDFOParameters):
+
+    r"""Container to hold parameters that specify the behavior of COBYLA.
+
+    Suggested values come from scipy documentation for scipy.optimize.fmin_cobyla:
+    http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.optimize.fmin_cobyla.html
+
+    :ivar rhobeg: (*float64 > 0.0*) reasonable initial changes to the variables (suggest: 1.0)
+    :ivar rhoend: (*float64 > 0.0*) final accuracy in the optimization (not precisely guaranteed), which is a lower bound on the size of the trust region (suggest: 1.0e-4)
+    :ivar maxfun: (*int > 0*) maximum number of objective function calls to make (suggest: 1000)
+    :ivar catol: (*float64 > 0.0*) absolute tolerance for constraint violations (suggest: 2.0e-4)
+
+    """
+
+    __slots__ = ()
+
+
 class NullOptimizer(OptimizerInterface):
 
     """A "null" or identity optimizer: this does nothing. It is used to perform "dumb" search with MultistartOptimizer."""
@@ -599,15 +625,16 @@ def __init__(self, domain, optimizable, optimization_parameters, num_random_samp
     def _scipy_decorator(self, func, **kwargs):
         """Wrapper function for expected improvement calculation to feed into BFGS.
 
-        func should be of the form compute_* in interfaces.optimization_interface.OptimizableInterface.
+        func should be of the form `compute_*` in :class:`moe.optimal_learning.python.interfaces.optimization_interface.OptimizableInterface`.
         """
         def decorated(point):
             """Decorator for compute_* functions in interfaces.optimization_interface.OptimizableInterface.
 
-            Converts the point to proper format and sets the current point before calling the compute function.
+            Converts the point to proper format (array with dim (self._num_points, self.domain.dim) instead of flat array)
+            and sets the current point before calling the compute function.
 
             :param point: the point on which to do the calculation
-            :type point: array of float64 with shape (self._num_points * self.domain.dim)
+            :type point: array of float64 with shape (self._num_points * self.domain.dim, )
             """
             shaped_point = point.reshape(self._num_points, self.domain.dim)
             self.objective_function.current_point = shaped_point
@@ -646,3 +673,81 @@ def optimize(self, **kwargs):
         else:
             shaped_point = unshaped_point.reshape(self._num_points, self.domain.dim)
         self.objective_function.current_point = shaped_point
+
+
+class ConstrainedDFOOptimizer(OptimizerInterface):
+
+    r"""Optimizes an objective function over the specified contraints with the COBYLA method.
+
+    .. Note:: See :func:`~moe.optimal_learning.python.python_version.optimization.ConstrainedDFOOptimizer.optimize()` docstring for more details.
+
+    """
+
+    def __init__(self, domain, optimizable, optimization_parameters, num_random_samples=None):
+        """Construct a ConstrainedDFOOptimizer.
+
+        :param domain: the domain that this optimizer operates over
+        :type domain: interfaces.domain_interface.DomainInterface subclass. Only supports TensorProductDomain for now.
+        :param optimizable: object representing the objective function being optimized
+        :type optimizable: interfaces.optimization_interface.OptimizableInterface subclass
+        :param optimization_parameters: parameters describing how to perform optimization (tolerances, iterations, etc.)
+        :type optimization_parameters: python_version.optimization.ConstrainedDFOParameters object
+
+        """
+        self.domain = domain
+        self.objective_function = optimizable
+        self.optimization_parameters = optimization_parameters
+        self._num_points = 1
+        # Check if this is a repeated domain, and if so set points equal to number of repeats.
+        if hasattr(self.domain, 'num_repeats'):
+            self._num_points = self.domain.num_repeats
+
+    def _scipy_decorator(self, func, **kwargs):
+        """Wrapper function for expected improvement calculation to feed into COBYLA.
+
+        func should be of the form compute_* in :class:`moe.optimal_learning.python.interfaces.optimization_interface.OptimizableInterface`.
+        """
+        def decorated(point):
+            """Decorator for compute_* functions in interfaces.optimization_interface.OptimizableInterface.
+
+            Converts the point to proper format (array with dim (self._num_points, self.domain.dim) instead of flat array)
+            and sets the current point before calling the compute function.
+
+            :param point: the point on which to do the calculation
+            :type point: array of float64 with shape (self._num_points * self.domain.dim)
+            """
+            shaped_point = point.reshape(self._num_points, self.domain.dim)
+            self.objective_function.current_point = shaped_point
+            value = -func(**kwargs)
+            if isinstance(value, (numpy.ndarray)):
+                return value.flatten()
+            else:
+                return value
+
+        return decorated
+
+    def optimize(self, **kwargs):
+        """Perform a COBYLA optimization given the parameters in optimization_parameters.
+
+        For more information, visit the scipy docs page and the original paper by Powell:
+        http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.optimize.fmin_cobyla.html
+        http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2007_03.pdf
+
+        objective_function.current_point will be set to the optimal point found.
+        """
+        # Parameters defined above in :class:`~moe.optimal_learning.python.python_version.optimization.LBFGSBParameters` class.
+        unshaped_point = scipy.optimize.fmin_cobyla(
+            func=self._scipy_decorator(self.objective_function.compute_objective_function, **kwargs),
+            x0=self.objective_function.current_point.flatten(),
+            cons=self.domain.get_constraint_list(),
+            rhobeg=self.optimization_parameters.rhobeg,
+            rhoend=self.optimization_parameters.rhoend,
+            maxfun=self.optimization_parameters.maxfun,
+            catol=self.optimization_parameters.catol,
+            disp=0,  # Suppresses output from the routine.
+        )
+        if self._num_points == 1:
+            shaped_point = unshaped_point
+        else:
+            shaped_point = unshaped_point.reshape(self._num_points, self.domain.dim)
+        self.objective_function.current_point = shaped_point

moe/optimal_learning/python/repeated_domain.py

@@ -79,6 +79,13 @@ def get_bounding_box(self):
         """Return a list of ClosedIntervals representing a bounding box for this domain."""
         return self._domain.get_bounding_box()
 
+    def get_constraint_list(self, start_index=0):
+        """Return a list of lambda functions expressing the domain bounds as linear constraints. Used by COBYLA."""
+        constraints = []
+        for i in xrange(self.num_repeats):
+            constraints.extend(self._domain.get_constraint_list(start_index=self._domain.dim * i))
+        return constraints
+
     def generate_random_point_in_domain(self, random_source=None):
         """Generate ``point`` uniformly at random such that ``self.check_point_inside(point)`` is True.
 

moe/tests/optimal_learning/python/python_version/optimization_test.py

@@ -7,7 +7,7 @@
 from moe.optimal_learning.python.geometry_utils import ClosedInterval
 from moe.optimal_learning.python.interfaces.optimization_interface import OptimizableInterface
 from moe.optimal_learning.python.python_version.domain import TensorProductDomain
-from moe.optimal_learning.python.python_version.optimization import multistart_optimize, LBFGSBParameters, GradientDescentParameters, NullOptimizer, GradientDescentOptimizer, MultistartOptimizer, LBFGSBOptimizer
+from moe.optimal_learning.python.python_version.optimization import multistart_optimize, LBFGSBParameters, GradientDescentParameters, NullOptimizer, GradientDescentOptimizer, MultistartOptimizer, LBFGSBOptimizer, ConstrainedDFOOptimizer, ConstrainedDFOParameters
 from moe.tests.optimal_learning.python.optimal_learning_test_case import OptimalLearningTestCase
 
 
@@ -163,6 +163,9 @@ def base_setup(self):
         domain_bounds = [ClosedInterval(-1.0, 1.0)] * self.dim
         self.domain = TensorProductDomain(domain_bounds)
 
+        large_domain_bounds = [ClosedInterval(-1.0, 1.0)] * self.dim
+        self.large_domain = TensorProductDomain(large_domain_bounds)
+
         maxima_point = numpy.full(self.dim, 0.5)
         current_point = numpy.zeros(self.dim)
         self.polynomial = QuadraticFunction(maxima_point, current_point)
@@ -199,36 +202,16 @@ def base_setup(self):
             epsilon,
         )
 
-    def test_gradient_descent_optimizer(self):
-        """Check that gradient descent can find the optimum of the quadratic test objective."""
-        # Check the claimed optima is an optima
-        optimum_point = self.polynomial.optimum_point
-        self.polynomial.current_point = optimum_point
-        gradient = self.polynomial.compute_grad_objective_function()
-        self.assert_vector_within_relative(gradient, numpy.zeros(self.polynomial.dim), 0.0)
-
-        # Verify that gradient descent does not move from the optima if we start it there.
-        gradient_descent_optimizer = GradientDescentOptimizer(self.domain, self.polynomial, self.gd_parameters)
-        gradient_descent_optimizer.optimize()
-        output = gradient_descent_optimizer.objective_function.current_point
-        self.assert_vector_within_relative(output, optimum_point, 0.0)
-
-        # Start at a wrong point and check optimization
-        tolerance = 2.0e-13
-        initial_guess = numpy.full(self.polynomial.dim, 0.2)
-        gradient_descent_optimizer.objective_function.current_point = initial_guess
-        gradient_descent_optimizer.optimize()
-        output = gradient_descent_optimizer.objective_function.current_point
-        # Verify coordinates
-        self.assert_vector_within_relative(output, optimum_point, tolerance)
-
-        # Verify function value
-        value = self.polynomial.compute_objective_function()
-        self.assert_scalar_within_relative(value, self.polynomial.optimum_value, tolerance)
-
-        # Verify derivative
-        gradient = self.polynomial.compute_grad_objective_function()
-        self.assert_vector_within_relative(gradient, numpy.zeros(self.polynomial.dim), tolerance)
+        maxfun = 1000
+        rhobeg = 1.0
+        rhoend = 1.0e-13
+        catol = 2.0e-13
+        self.COBYLA_parameters = ConstrainedDFOParameters(
+            rhobeg,
+            rhoend,
+            maxfun,
+            catol,
+        )
 
     def test_get_averaging_range(self):
         """Test the method used to produce what interval to average over in Polyak-Ruppert averaging."""
@@ -351,49 +334,25 @@ def test_multistarted_gradient_descent_optimizer_crippled_start(self):
         for value in (test_best_point - self.polynomial.optimum_point):
             T.assert_equal(value, 0.0)
 
-    def test_multistarted_gradient_descent_optimizer(self):
-        """Check that multistarted GD can find the optimum in a 'very' large domain."""
-        # Set a large domain: a single GD run is unlikely to reach the optimum
-        domain_bounds = [ClosedInterval(-10.0, 10.0)] * self.dim
-        domain = TensorProductDomain(domain_bounds)
-
-        tolerance = 2.0e-10
-        num_points = 10
-        gradient_descent_optimizer = GradientDescentOptimizer(domain, self.polynomial, self.gd_parameters)
-        multistart_optimizer = MultistartOptimizer(gradient_descent_optimizer, num_points)
-
-        output, _ = multistart_optimizer.optimize()
-        # Verify coordinates
-        self.assert_vector_within_relative(output, self.polynomial.optimum_point, tolerance)
-
-        # Verify function value
-        value = self.polynomial.compute_objective_function()
-        self.assert_scalar_within_relative(value, self.polynomial.optimum_value, tolerance)
-
-        # Verify derivative
-        gradient = self.polynomial.compute_grad_objective_function()
-        self.assert_vector_within_relative(gradient, numpy.zeros(self.polynomial.dim), tolerance)
-
-    def test_bfgs_optimizer(self):
-        """Check that BFGS can find the optimum of the quadratic test objective."""
+    def optimizer_test(self, optimizer):
+        """Check that the optimizer can find the optimum of the quadratic test objective."""
         # Check the claimed optima is an optima
         optimum_point = self.polynomial.optimum_point
         self.polynomial.current_point = optimum_point
         gradient = self.polynomial.compute_grad_objective_function()
         self.assert_vector_within_relative(gradient, numpy.zeros(self.polynomial.dim), 0.0)
 
-        # Verify that gradient descent does not move from the optima if we start it there.
-        bfgs_optimizer = LBFGSBOptimizer(self.domain, self.polynomial, self.BFGS_parameters)
-        bfgs_optimizer.optimize()
-        output = bfgs_optimizer.objective_function.current_point
-        self.assert_vector_within_relative(output, optimum_point, 0.0)
+        # Verify that the optimizer does not move from the optima if we start it there.
+        tolerance = 2.0e-13
+        optimizer.optimize()
+        output = optimizer.objective_function.current_point
+        self.assert_vector_within_relative(output, optimum_point, tolerance)
 
         # Start at a wrong point and check optimization
-        tolerance = 2.0e-13
         initial_guess = numpy.full(self.polynomial.dim, 0.2)
-        bfgs_optimizer.objective_function.current_point = initial_guess
-        bfgs_optimizer.optimize()
-        output = bfgs_optimizer.objective_function.current_point
+        optimizer.objective_function.current_point = initial_guess
+        optimizer.optimize()
+        output = optimizer.objective_function.current_point
         # Verify coordinates
         self.assert_vector_within_relative(output, optimum_point, tolerance)
 
@@ -405,16 +364,11 @@ def test_bfgs_optimizer(self):
         gradient = self.polynomial.compute_grad_objective_function()
         self.assert_vector_within_relative(gradient, numpy.zeros(self.polynomial.dim), tolerance)
 
-    def test_multistarted_bfgs_optimizer(self):
-        """Check that multistarted GD can find the optimum in a 'very' large domain."""
-        # Set a large domain: a single GD run is unlikely to reach the optimum
-        domain_bounds = [ClosedInterval(-10.0, 10.0)] * self.dim
-        domain = TensorProductDomain(domain_bounds)
-
+    def multistarted_optimizer_test(self, optimizer):
+        """Check that the multistarted optimizer can find the optimum in a 'very' large domain."""
         tolerance = 2.0e-10
         num_points = 10
-        bfgs_optimizer = LBFGSBOptimizer(domain, self.polynomial, self.BFGS_parameters)
-        multistart_optimizer = MultistartOptimizer(bfgs_optimizer, num_points)
+        multistart_optimizer = MultistartOptimizer(optimizer, num_points)
 
         output, _ = multistart_optimizer.optimize()
         # Verify coordinates
@@ -427,3 +381,33 @@ def test_multistarted_bfgs_optimizer(self):
         # Verify derivative
         gradient = self.polynomial.compute_grad_objective_function()
         self.assert_vector_within_relative(gradient, numpy.zeros(self.polynomial.dim), tolerance)
+
+    def test_cobyla_optimizer(self):
+        """Test if COBYLA can optimize a simple objective function."""
+        constrained_optimizer = ConstrainedDFOOptimizer(self.domain, self.polynomial, self.COBYLA_parameters)
+        self.optimizer_test(constrained_optimizer)
+
+    def test_bfgs_optimizer(self):
+        """Test if BFGS can optimize a simple objective function."""
+        bfgs_optimizer = LBFGSBOptimizer(self.domain, self.polynomial, self.BFGS_parameters)
+        self.optimizer_test(bfgs_optimizer)
+
+    def test_gradient_descent_optimizer(self):
+        """Test if Gradient Descent can optimize a simple objective function."""
+        gradient_descent_optimizer = GradientDescentOptimizer(self.domain, self.polynomial, self.gd_parameters)
+        self.optimizer_test(gradient_descent_optimizer)
+
+    def test_cobyla_multistarted_optimizer(self):
+        """Test if COBYLA can optimize a "hard" objective function with multistarts."""
+        constrained_optimizer = ConstrainedDFOOptimizer(self.large_domain, self.polynomial, self.COBYLA_parameters)
+        self.multistarted_optimizer_test(constrained_optimizer)
+
+    def test_bfgs_multistarted_optimizer(self):
+        """Test if BFGS can optimize a "hard" objective function with multistarts."""
+        bfgs_optimizer = LBFGSBOptimizer(self.large_domain, self.polynomial, self.BFGS_parameters)
+        self.multistarted_optimizer_test(bfgs_optimizer)
+
+    def test_gradient_descent_multistarted_optimizer(self):
+        """Test if Gradient Descent can optimize a "hard" objective function with multistarts."""
+        gradient_descent_optimizer = GradientDescentOptimizer(self.large_domain, self.polynomial, self.gd_parameters)
+        self.multistarted_optimizer_test(gradient_descent_optimizer)