Stop using triangle wave transform in optimisations

CHANGELOG.md

@@ -5,13 +5,20 @@ All notable changes to this project will be documented in this file.
 ## Unreleased
 
 ### Added
-- [#1420](https://github.com/pints-team/pints/pull/1420) The `Optimisation` objects now distinguish between a best-visited point (`x_best`, with score `f_best`) and a best-guessed point (`x_guessed`, with approximate score `f_guessed`). For most optimisers, the two values are equivalent. As before, the `OptimisationController` now tracks `x_best` and `f_best` by default, but this can be modified using the methods `set_f_guessed_tracking` and `f_guessed_tracking`.
-- 
+- [#1420](https://github.com/pints-team/pints/pull/1420) The `Optimiser` class now distinguishes between a best-visited point (`x_best`, with score `f_best`) and a best-guessed point (`x_guessed`, with approximate score `f_guessed`). For most optimisers, the two values are equivalent. The `OptimisationController` still tracks `x_best` and `f_best` by default, but this can be modified using the methods `set_f_guessed_tracking` and `f_guessed_tracking`.
+
 ### Changed
+- [#1424](https://github.com/pints-team/pints/pull/1424) Fixed a bug in PSO that caused it to use more particles than advertised.
+- [#1424](https://github.com/pints-team/pints/pull/1424) xNES, SNES, PSO, and BareCMAES no longer use a `TriangleWaveTransform` to handle rectangular boundaries (this was found to lead to optimisers diverging in some cases).
+
 ### Deprecated
+
 ### Removed
+- [#1424](https://github.com/pints-team/pints/pull/1424) Removed the `TriangleWaveTransform` class previously used in some optimisers.
+
 ### Fixed
 
+
 ## [0.4.0] - 2021-12-07
 
 ### Added

docs/source/optimisers/index.rst

@@ -16,7 +16,6 @@ or the :class:`OptimisationController` class.
     running
     base_classes
     convenience_methods
-    boundary_transformations
     cmaes_bare
     cmaes
     gradient_descent

pints/__init__.py

@@ -170,7 +170,6 @@ def version(formatted=False):
     optimise,
     Optimiser,
     PopulationBasedOptimiser,
-    TriangleWaveTransform,
 )
 from ._optimisers._cmaes import CMAES
 from ._optimisers._cmaes_bare import BareCMAES

pints/_optimisers/__init__.py

@@ -973,43 +973,6 @@ def optimise(
         function, x0, sigma0, boundaries, transformation, method).run()
 
 
-class TriangleWaveTransform(object):
-    """
-    Transforms from unbounded to (rectangular) bounded parameter space using a
-    periodic triangle-wave transform.
-
-    Note: The transform is applied _inside_ optimisation methods, there is no
-    need to wrap this around your own problem or score function.
-
-    This can be applied as a transformation on ``x`` to implement _rectangular_
-    boundaries in methods with no natural boundary mechanism. It effectively
-    mirrors the search space at every boundary, leading to a continuous (but
-    non-smooth) periodic landscape. While this effectively creates an infinite
-    number of minima/maxima, each one maps to the same point in parameter
-    space.
-
-    It should work well for methods that maintain a single search position or a
-    single search distribution (e.g. :class:`CMAES`, :class:`xNES`,
-    :class:`SNES`), which will end up in one of the many mirror images.
-    However, for methods that use independent search particles (e.g.
-    :class:`PSO`) it could lead to a scattered population, with different
-    particles exploring different mirror images. Other strategies should be
-    used for such problems.
-    """
-
-    def __init__(self, boundaries):
-        self._lower = boundaries.lower()
-        self._upper = boundaries.upper()
-        self._range = self._upper - self._lower
-        self._range2 = 2 * self._range
-
-    def __call__(self, x):
-        y = np.remainder(x - self._lower, self._range2)
-        z = np.remainder(y, self._range)
-        return ((self._lower + z) * (y < self._range)
-                + (self._upper - z) * (y >= self._range))
-
-
 def curve_fit(f, x, y, p0, boundaries=None, threshold=None, max_iter=None,
               max_unchanged=200, verbose=False, parallel=False, method=None):
     """

pints/_optimisers/_cmaes_bare.py

@@ -83,9 +83,6 @@ def __init__(self, x0, sigma0=0.1, boundaries=None):
             / gamma(self._n_parameters / 2)
         )
 
-        # Optional transformation to within-the-boundaries
-        self._boundary_transform = None
-
     def ask(self):
         """ See :meth:`Optimiser.ask()`. """
         # Initialise on first call
@@ -106,19 +103,14 @@ def ask(self):
         # Samples from N(mu, eta**2 * C)
         self._xs = np.array([self._mu + self._eta * y for y in self._ys])
 
-        # Apply boundaries; creating safe points for evaluation
-        # Rectangular boundaries? Then perform boundary transform
-        if self._boundary_transform is not None:
-            self._xs = self._boundary_transform(self._xs)
-
-        # Manual boundaries? Then pass only xs that are within bounds
-        if self._manual_boundaries:
+        # Boundaries? Then only pass user xs that are within bounds
+        if self._boundaries is not None:
             self._user_ids = np.nonzero(
                 [self._boundaries.check(x) for x in self._xs])
             self._user_xs = self._xs[self._user_ids]
             if len(self._user_xs) == 0:     # pragma: no cover
-                warnings.warn('All points requested by CMA-ES are outside the'
-                              ' boundaries.')
+                warnings.warn('All points requested by BareCMAES are outside'
+                              ' the boundaries.')
         else:
             self._user_xs = self._xs
 
@@ -161,14 +153,6 @@ def _initialise(self):
         """
         assert (not self._running)
 
-        # Create boundary transform, or use manual boundary checking
-        self._manual_boundaries = False
-        if isinstance(self._boundaries, pints.RectangularBoundaries):
-            self._boundary_transform = pints.TriangleWaveTransform(
-                self._boundaries)
-        elif self._boundaries is not None:
-            self._manual_boundaries = True
-
         # Parent generation population size
         # The parameter parent_pop_size is the mu in the papers. It represents
         # the size of a parent population used to update our paramters.
@@ -280,8 +264,8 @@ def tell(self, fx):
         npo = self._population_size
         npa = self._parent_pop_size
 
-        # Manual boundaries? Then reconstruct full fx vector
-        if self._manual_boundaries and len(fx) < npo:
+        # Boundaries? Then reconstruct full fx vector
+        if self._boundaries is not None and len(fx) < npo:
             user_fx = fx
             fx = np.ones((npo, )) * float('inf')
             fx[self._user_ids] = user_fx
@@ -377,6 +361,4 @@ def x_best(self):
 
     def x_guessed(self):
         """ See :meth:`Optimiser.x_guessed()`. """
-        if self._boundary_transform is not None:
-            return self._boundary_transform(self._mu)
         return np.array(self._mu, copy=True)

pints/_optimisers/_pso.py

@@ -113,19 +113,22 @@ def _initialise(self):
 
         # Set initial positions
         self._xs.append(np.array(self._x0, copy=True))
+
+        # Attempt to sample n - 1 points from the boundaries
         if self._boundaries is not None:
-            # Attempt to sample n - 1 points from the boundaries
             try:
                 self._xs.extend(
                     self._boundaries.sample(self._population_size - 1))
             except NotImplementedError:
-                # Not all boundaries implement sampling
+                # Not all boundaries implement sampling.
                 pass
+
         # If we couldn't sample from the boundaries, use gaussian sampling
         # around x0.
-        for i in range(1, self._population_size):
-            self._xs.append(np.random.normal(self._x0, self._sigma0))
-        self._xs = np.array(self._xs, copy=True)
+        if len(self._xs) < self._population_size:
+            for i in range(self._population_size - 1):
+                self._xs.append(np.random.normal(self._x0, self._sigma0))
+        self._xs = np.array(self._xs)
 
         # Set initial velocities
         for i in range(self._population_size):
@@ -141,29 +144,15 @@ def _initialise(self):
         self._fg = float('inf')
         self._pg = self._xs[0]
 
-        # Create boundary transform, or use manual boundary checking
-        self._manual_boundaries = False
-        self._boundary_transform = None
-        if isinstance(self._boundaries, pints.RectangularBoundaries):
-            self._boundary_transform = pints.TriangleWaveTransform(
-                self._boundaries)
-        elif self._boundaries is not None:
-            self._manual_boundaries = True
-
-        # Create safe xs to pass to user
-        if self._boundary_transform is not None:
-            # Rectangular boundaries? Then apply transform to xs
-            self._xs = self._boundary_transform(self._xs)
-        if self._manual_boundaries:
-            # Manual boundaries? Then filter out out-of-bounds points from xs
+        # Boundaries? Then filter out out-of-bounds points from xs
+        self._user_xs = self._xs
+        if self._boundaries is not None:
             self._user_ids = np.nonzero(
                 [self._boundaries.check(x) for x in self._xs])
             self._user_xs = self._xs[self._user_ids]
             if len(self._user_xs) == 0:     # pragma: no cover
                 warnings.warn(
                     'All initial PSO particles are outside the boundaries.')
-        else:
-            self._user_xs = self._xs
 
         # Set local/global exploration balance
         self.set_local_global_balance()
@@ -194,8 +183,8 @@ def running(self):
     def set_local_global_balance(self, r=0.5):
         """
         Set the balance between local and global exploration for each particle,
-        using a parameter `r` such that `r = 1` is a fully local search and
-        `r = 0` is a fully global search.
+        using a parameter ``r`` such that ``r = 1`` is a fully local search and
+        ``r = 0`` is a fully global search.
         """
         if self._running:
             raise Exception('Cannot change settings during run.')
@@ -234,8 +223,8 @@ def tell(self, fx):
             raise Exception('ask() not called before tell()')
         self._ready_for_tell = False
 
-        # Manual boundaries? Then reconstruct full fx vector
-        if self._manual_boundaries and len(fx) < self._population_size:
+        # Boundaries? Then reconstruct full fx vector
+        if self._boundaries is not None and len(fx) < self._population_size:
             user_fx = fx
             fx = np.ones((self._population_size, )) * float('inf')
             fx[self._user_ids] = user_fx
@@ -265,19 +254,14 @@ def tell(self, fx):
             # Update position
             self._xs[i] += self._vs[i]
 
-        # Create safe xs to pass to user
-        if self._boundary_transform is not None:
-            # Rectangular boundaries? Then apply transform to xs
-            self._user_xs = self._xs = self._boundary_transform(self._xs)
-        elif self._manual_boundaries:
-            # Manual boundaries? Then filter out out-of-bounds points from xs
+        # Boundaries? Then filter out out-of-bounds points from xs
+        self._user_xs = self._xs
+        if self._boundaries is not None:
             self._user_ids = np.nonzero(
                 [self._boundaries.check(x) for x in self._xs])
             self._user_xs = self._xs[self._user_ids]
             if len(self._user_xs) == 0:     # pragma: no cover
                 warnings.warn('All PSO particles are outside the boundaries.')
-        else:
-            self._user_xs = self._xs
 
         # Update global best score
         i = np.argmin(self._fl)

pints/_optimisers/_snes.py

@@ -53,9 +53,6 @@ def __init__(self, x0, sigma0=None, boundaries=None):
         # We don't have f(mu), so we approximate it by max f(sample)
         self._f_guessed = float('inf')
 
-        # Optional transformation to within-the-boundaries
-        self._boundary_transform = None
-
     def ask(self):
         """ See :meth:`Optimiser.ask()`. """
         # Initialise on first call
@@ -70,12 +67,8 @@ def ask(self):
                             for i in range(self._population_size)])
         self._xs = self._mu + self._sigmas * self._ss
 
-        # Create safe xs to pass to user
-        if self._boundary_transform is not None:
-            # Rectangular boundaries? Then perform boundary transform
-            self._xs = self._boundary_transform(self._xs)
-        if self._manual_boundaries:
-            # Manual boundaries? Then pass only xs that are within bounds
+        # Boundaries? Then only pass user xs that are within bounds
+        if self._boundaries is not None:
             self._user_ids = np.nonzero(
                 [self._boundaries.check(x) for x in self._xs])
             self._user_xs = self._xs[self._user_ids]
@@ -103,15 +96,6 @@ def _initialise(self):
         """
         assert(not self._running)
 
-        # Create boundary transform, or use manual boundary checking
-        self._manual_boundaries = False
-        self._boundary_transform = None
-        if isinstance(self._boundaries, pints.RectangularBoundaries):
-            self._boundary_transform = pints.TriangleWaveTransform(
-                self._boundaries)
-        elif self._boundaries is not None:
-            self._manual_boundaries = True
-
         # Shorthands
         d = self._n_parameters
         n = self._population_size
@@ -154,8 +138,8 @@ def tell(self, fx):
             raise Exception('ask() not called before tell()')
         self._ready_for_tell = False
 
-        # Manual boundaries? Then reconstruct full fx vector
-        if self._manual_boundaries and len(fx) < self._population_size:
+        # Boundaries? Then reconstruct full fx vector
+        if self._boundaries is not None and len(fx) < self._population_size:
             user_fx = fx
             fx = np.ones((self._population_size, )) * float('inf')
             fx[self._user_ids] = user_fx
@@ -186,7 +170,5 @@ def x_best(self):
 
     def x_guessed(self):
         """ See :meth:`Optimiser.x_guessed()`. """
-        if self._boundary_transform is not None:
-            return self._boundary_transform(self._mu)
         return np.array(self._mu, copy=True)
 

pints/_optimisers/_xnes.py

@@ -58,9 +58,6 @@ def __init__(self, x0, sigma0=None, boundaries=None):
         # We don't have f(mu), so we approximate it by min f(sample)
         self._f_guessed = float('inf')
 
-        # Optional transformation to within-the-boundaries
-        self._boundary_transform = None
-
     def ask(self):
         """ See :meth:`Optimiser.ask()`. """
         # Initialise on first call
@@ -76,12 +73,8 @@ def ask(self):
         self._xs = np.array([self._mu + np.dot(self._A, self._zs[i])
                              for i in range(self._population_size)])
 
-        # Create safe xs to pass to user
-        if self._boundary_transform is not None:
-            # Rectangular boundaries? Then perform boundary transform
-            self._xs = self._boundary_transform(self._xs)
-        if self._manual_boundaries:
-            # Manual boundaries? Then pass only xs that are within bounds
+        # Boundaries? Then only pass user xs that are within bounds
+        if self._boundaries is not None:
             self._bounded_ids = np.nonzero(
                 [self._boundaries.check(x) for x in self._xs])
             self._bounded_xs = self._xs[self._bounded_ids]
@@ -109,15 +102,6 @@ def _initialise(self):
         """
         assert(not self._running)
 
-        # Create boundary transform, or use manual boundary checking
-        self._manual_boundaries = False
-        self._boundary_transform = None
-        if isinstance(self._boundaries, pints.RectangularBoundaries):
-            self._boundary_transform = pints.TriangleWaveTransform(
-                self._boundaries)
-        elif self._boundaries is not None:
-            self._manual_boundaries = True
-
         # Shorthands
         d = self._n_parameters
         n = self._population_size
@@ -163,8 +147,8 @@ def tell(self, fx):
             raise Exception('ask() not called before tell()')
         self._ready_for_tell = False
 
-        # Manual boundaries? Then reconstruct full fx vector
-        if self._manual_boundaries and len(fx) < self._population_size:
+        # Boundaries? Then reconstruct full fx vector
+        if self._boundaries is not None and len(fx) < self._population_size:
             bounded_fx = fx
             fx = np.ones((self._population_size, )) * float('inf')
             fx[self._bounded_ids] = bounded_fx
@@ -198,7 +182,5 @@ def x_best(self):
 
     def x_guessed(self):
         """ See :meth:`Optimiser.x_guessed()`. """
-        if self._boundary_transform is not None:
-            return self._boundary_transform(self._mu)
         return np.array(self._mu, copy=True)