Update photonics with better mutations and recombinations (#477)

nevergrad/benchmark/experiments.py

@@ -12,6 +12,7 @@
 from nevergrad.functions import FarOptimumFunction
 from nevergrad.functions import MultiobjectiveFunction
 from nevergrad.functions import mlda as _mlda
+from nevergrad.functions.photonics import Photonics
 from nevergrad.functions.arcoating import ARCoating
 from nevergrad.functions.powersystems import PowerSystem
 from nevergrad.functions.stsp import STSP
@@ -173,7 +174,7 @@ def multimodal(seed: Optional[int] = None, para: bool = False) -> Iterator[Exper
               "MultiCMA", "TripleCMA", "MultiScaleCMA"]
     if not para:
         optims += ["RSQP", "RCobyla", "RPowell", "SQPCMA", "SQP", "Cobyla", "Powell"]
-    #+ list(sorted(x for x, y in ng.optimizers.registry.items() if "chain" in x or "BO" in x))
+    # + list(sorted(x for x, y in ng.optimizers.registry.items() if "chain" in x or "BO" in x))
     functions = [
         ArtificialFunction(name, block_dimension=bd, useless_variables=bd * uv_factor)
         for name in names
@@ -206,7 +207,7 @@ def yabbob(seed: Optional[int] = None, parallel: bool = False, big: bool = False
               "TwoPointsDE", "OnePointDE", "AlmostRotationInvariantDE", "RotationInvariantDE", "CMandAS2", "CMandAS"]
     if not parallel:
         optims += ["SQP", "Cobyla", "Powell", "chainCMASQP", "chainCMAPowell"]
-    #optims += [x for x, y in ng.optimizers.registry.items() if "chain" in x]
+    # optims += [x for x, y in ng.optimizers.registry.items() if "chain" in x]
     names = ["hm", "rastrigin", "griewank", "rosenbrock", "ackley", "lunacek", "deceptivemultimodal", "bucherastrigin", "multipeak"]
     names += ["sphere", "doublelinearslope", "stepdoublelinearslope"]
     names += ["cigar", "altcigar", "ellipsoid", "altellipsoid", "stepellipsoid", "discus", "bentcigar"]
@@ -737,3 +738,23 @@ def far_optimum_es(seed: tp.Optional[int] = None) -> Iterator[Experiment]:
         for optim in optimizers:
             for budget in [100, 400, 1000, 4000, 10000]:
                 yield Experiment(func, optim, budget=budget, seed=next(seedg))
+
+
+@registry.register
+def photonics(seed: tp.Optional[int] = None) -> Iterator[Experiment]:
+    seedg = create_seed_generator(seed)
+    popsizes = [10, 40, 100]
+    es = [ng.families.EvolutionStrategy(recombination_ratio=recomb, only_offsprings=False, popsize=pop)
+          for recomb in [0.1, 1] for pop in popsizes]
+    es += [ng.families.EvolutionStrategy(recombination_ratio=recomb, only_offsprings=only, popsize=pop,
+                                         offsprings={10: 20, 40: 60, 100: 150}[pop])
+           for only in [True, False] for recomb in [0.1, 1] for pop in popsizes]
+    algos = ["TwoPointsDE", "DE", "PSO", "OnePlusOne", "ParametrizationDE", "NaiveTBPSA"] + es  # type: ignore
+    for method in ["clipping", "tanh", "arctan"]:
+        # , "chirped"]]:  # , "morpho"]]:
+        for func in [Photonics(x, 60 if x == "morpho" else 80, bounding_method=method) for x in ["bragg"]]:
+            for budget in [1e2, 1e3, 1e4, 1e5, 1e6]:
+                for algo in algos:
+                    xp = Experiment(func, algo, int(budget), num_workers=1, seed=next(seedg))
+                    if not xp.is_incoherent:
+                        yield xp

nevergrad/functions/photonics/core.py

@@ -19,22 +19,22 @@
 #   Mihailovic, M., Centeno, E., Ciracì, C., Smith, D.R. and Moreau, A., 2016.
 #   Moosh: A Numerical Swiss Army Knife for the Optics of Multilayers in Octave/Matlab. Journal of Open Research Software, 4(1), p.e13.
 
-import typing as tp
 import numpy as np
 from nevergrad.parametrization import parameter as p
+from nevergrad.parametrization.utils import Crossover
 from . import photonics
 from ..base import ExperimentFunction
 
 
-def _make_instrumentation(name: str, dimension: int, transform: str = "tanh") -> p.Instrumentation:
+def _make_instrumentation(name: str, dimension: int, bounding_method: str = "clipping") -> p.Array:
     """Creates appropriate instrumentation for a Photonics problem
 
     Parameters
     name: str
         problem name, among bragg, chirped and morpho
     dimension: int
         size of the problem among 16, 40 and 60 (morpho) or 80 (bragg and chirped)
-    transform: str
+    bounding_method: str
         transform type for the bounding ("arctan", "tanh" or "clipping", see `Array.bounded`)
 
     Returns
@@ -43,30 +43,29 @@ def _make_instrumentation(name: str, dimension: int, transform: str = "tanh") ->
         the instrumentation for the problem
     """
     assert not dimension % 4, f"points length should be a multiple of 4, got {dimension}"
-    n = dimension // 4
-    arrays: tp.List[p.Array] = []
-    ones = np.ones((n,), dtype=float)
     if name == "bragg":
-        # n multiple of 2, from 16 to 80
-        # main (n=60): [2,3]^30 x [0,300]^30
-        arrays.extend([p.Array(init=2.5 * ones).set_bounds(2, 3, method=transform) for _ in range(2)])
-        arrays.extend([p.Array(init=150 * ones).set_bounds(0, 300, method=transform) for _ in range(2)])
+        shape = (2, dimension // 2)
+        bounds = [(2, 3), (30, 180)]
     elif name == "chirped":
-        # n multiple of 2, from 10 to 80
-        # domain (n=60): [0,300]^60
-        arrays = [p.Array(init=150 * ones).set_bounds(0, 300, method=transform) for _ in range(4)]
+        shape = (1, dimension)
+        bounds = [(0, 300)]
     elif name == "morpho":
-        # n multiple of 4, from 16 to 60
-        # domain (n=60): [0,300]^15 x [0,600]^15 x [30,600]^15 x [0,300]^15
-        arrays.extend([p.Array(init=150 * ones).set_bounds(0, 300, method=transform),
-                       p.Array(init=300 * ones).set_bounds(0, 600, method=transform),
-                       p.Array(init=315 * ones).set_bounds(30, 600, method=transform),
-                       p.Array(init=150 * ones).set_bounds(0, 300, method=transform)])
+        shape = (4, dimension // 4)
+        bounds = [(0, 300), (0, 600), (30, 600), (0, 300)]
     else:
         raise NotImplementedError(f"Transform for {name} is not implemented")
-    instrumentation = p.Instrumentation(*arrays)
-    assert instrumentation.dimension == dimension
-    return instrumentation
+    assert shape[0] * shape[1] == dimension, f"Cannot work with dimension {dimension} for {name}: not divisible by {shape[0]}."
+    b_array = np.array(bounds)
+    assert b_array.shape[0] == shape[0]  # pylint: disable=unsubscriptable-object
+    init = np.sum(b_array, axis=1, keepdims=True).dot(np.ones((1, shape[1],))) / 2
+    array = p.Array(init=init)
+    if bounding_method not in ("arctan", "tanh"):
+        # sigma must be adapted for clipping and constraint methods
+        array.set_mutation(sigma=p.Array(init=[[10.0]] if name != "bragg" else [[0.03], [10.0]]).set_mutation(exponent=2.0))  # type: ignore
+    array.set_bounds(b_array[:, [0]], b_array[:, [1]], method=bounding_method, full_range_sampling=True)
+    array.set_recombination(Crossover(2, structured_dimensions=(0,))).set_name("")
+    assert array.dimension == dimension, f"Unexpected {array} for dimension {dimension}"
+    return array
 
 
 class Photonics(ExperimentFunction):
@@ -113,17 +112,16 @@ class Photonics(ExperimentFunction):
       Moosh: A Numerical Swiss Army Knife for the Optics of Multilayers in Octave/Matlab. Journal of Open Research Software, 4(1), p.e13.
     """
 
-    def __init__(self, name: str, dimension: int, transform: str = "tanh") -> None:
-        assert dimension in [8, 16, 40, 60 if name == "morpho" else 80]
+    def __init__(self, name: str, dimension: int, bounding_method: str = "clipping") -> None:
         assert name in ["bragg", "morpho", "chirped"]
         self.name = name
         self._base_func = {"morpho": photonics.morpho, "bragg": photonics.bragg, "chirped": photonics.chirped}[name]
-        super().__init__(self._compute, _make_instrumentation(name=name, dimension=dimension, transform=transform))
-        self.register_initialization(name=name, dimension=dimension, transform=transform)
-        self._descriptors.update(name=name)
+        super().__init__(self._compute, _make_instrumentation(name=name, dimension=dimension, bounding_method=bounding_method))
+        self.register_initialization(name=name, dimension=dimension, bounding_method=bounding_method)
+        self._descriptors.update(name=name, bounding_method=bounding_method)
 
-    def _compute(self, *x: np.ndarray) -> float:
-        x_cat = np.concatenate(x)
+    def _compute(self, x: np.ndarray) -> float:
+        x_cat = x.ravel()
         assert x_cat.size == self.dimension
         try:
             output = self._base_func(x_cat)

nevergrad/functions/photonics/test_core.py

@@ -3,83 +3,114 @@
 # This source code is licensed under the MIT license found in the
 # LICENSE file in the root directory of this source tree.
 
+import os
+from unittest import SkipTest
 import warnings
-from typing import List
-from unittest.mock import patch
+import typing as tp
+import pytest
 import numpy as np
-from ...common import testing
+import nevergrad as ng
+from nevergrad.common import testing
 from . import core
 
 
 @testing.parametrized(
-    bragg=("bragg", [2.93, 2.18, 2.35, 2.12, 31.53, 15.98, 226.69, 193.11]),
+    bragg=("bragg", [2.93, 2.18, 2.35, 2.12, 45.77, 37.99, 143.34, 126.55]),
     morpho=("morpho", [280.36, 52.96, 208.16, 72.69, 89.92, 60.37, 226.69, 193.11]),
     chirped=("chirped", [280.36, 52.96, 104.08, 36.34, 31.53, 15.98, 226.69, 193.11]),
 )
-def test_photonics_transforms(pb: str, expected: List[float]) -> None:
-    with patch("shutil.which", return_value="here"):
-        # dim 8 is easier to test... but it is actually not allowed. Nevermind here, HACK IT NEXT LINE
-        with patch("nevergrad.functions.photonics.core._make_instrumentation", return_value=core._make_instrumentation(pb, 8)):
-            func = core.Photonics(pb, 16)
-    assert func.dimension == 8
+def test_photonics_bounding_methods(pb: str, expected: tp.List[float]) -> None:
+    func = core.Photonics(pb, 8, bounding_method="tanh")
     np.random.seed(24)
     x = np.random.normal(0, 1, size=8)
-    all_x = func.parametrization.spawn_child().set_standardized_data(x).args
-    output = np.concatenate(all_x)
+    output = func.parametrization.spawn_child().set_standardized_data(x).value.ravel()
     np.testing.assert_almost_equal(output, expected, decimal=2)
 
 
 @testing.parametrized(
     # bragg domain (n=60): [2,3]^30 x [0,300]^30
-    bragg_tanh=("bragg", "tanh", [2.5, 2.5, 2.5, 2.5, 150., 150., 150., 150.]),
-    bragg_arctan=("bragg", "arctan", [2.5, 2.5, 2.5, 2.5, 150., 150., 150., 150.]),
+    bragg_tanh=("bragg", "tanh", [2.5, 2.5, 2.5, 2.5, 105., 105., 105., 105.]),
+    bragg_arctan=("bragg", "arctan", [2.5, 2.5, 2.5, 2.5, 105., 105., 105., 105.]),
     # chirped domain (n=60): [0,300]^60
     chirped_tanh=("chirped", "tanh", [150., 150., 150., 150., 150., 150., 150., 150.]),
     chirped_arctan=("chirped", "arctan", [150., 150., 150., 150., 150., 150., 150., 150.]),
     # morpho domain (n=60): [0,300]^15 x [0,600]^15 x [30,600]^15 x [0,300]^15
     morpho_tanh=("morpho", "tanh", [150., 150., 300., 300., 315., 315., 150., 150.]),
     morpho_arctan=("morpho", "arctan", [150., 150., 300., 300., 315., 315., 150., 150.]),
 )
-def test_photonics_transforms_mean(pb: str, transform: str, expected: List[float]) -> None:
-    with patch("shutil.which", return_value="here"):
-        # dim 8 is easier to test... but it is actually not allowed. Nevermind here, HACK IT NEXT LINE
-        with patch("nevergrad.functions.photonics.core._make_instrumentation", return_value=core._make_instrumentation(pb, 8, transform)):
-            func = core.Photonics(pb, 16, transform=transform)
-    all_x = func.parametrization.args
-    output = np.concatenate(all_x)
+def test_photonics_bounding_methods_mean(pb: str, bounding_method: str, expected: tp.List[float]) -> None:
+    func = core.Photonics(pb, 8, bounding_method=bounding_method)
+    all_x = func.parametrization.value
+    output = all_x.ravel()
     np.testing.assert_almost_equal(output, expected, decimal=2)
 
 
-def test_morpho_transform_constraints() -> None:
-    with patch("shutil.which", return_value="here"):
-        func = core.Photonics("morpho", 60)
+def test_morpho_bounding_method_constraints() -> None:
+    func = core.Photonics("morpho", 60, bounding_method="arctan")
     x = np.random.normal(0, 5, size=60)  # std 5 to play with boundaries
-    all_x = func.parametrization.spawn_child().set_standardized_data(x).args
-    output = np.concatenate(all_x)
-    assert np.all(output >= 0)
-    q = len(x) // 4
-    assert np.all(output[:q] <= 300)
-    assert np.all(output[q: 3 * q] <= 600)
-    assert np.all(output[2 * q: 3 * q] >= 30)
-    assert np.all(output[3 * q:] <= 300)
+    output1 = func.parametrization.spawn_child().set_standardized_data(x)
+    output2 = output1.sample()
+    for output in (output1, output2):
+        assert np.all(output.value >= 0)
+        assert np.all(output.value[[0, 3], :] <= 300)
+        assert np.all(output.value[[1, 2], :] <= 600)
+        assert np.all(output.value[2, :] >= 30)
+
+
+def test_photonics_recombination() -> None:
+    func = core.Photonics("morpho", 16)
+    func.parametrization.random_state.seed(24)
+    arrays: tp.List[ng.p.Array] = []
+    for num in [50, 71]:
+        arrays.append(func.parametrization.spawn_child())  # type: ignore
+        arrays[-1].value = num * np.ones(arrays[0].value.shape)
+    arrays[0].recombine(arrays[1])
+    expected = [71, 50, 71, 71]
+    np.testing.assert_array_equal(arrays[0].value, np.ones((4, 1)).dot(np.array(expected)[None, :]))
 
 
 def test_photonics_error() -> None:
     # check error
     photo = core.Photonics("bragg", 16)
-    np.testing.assert_raises(AssertionError, photo, np.zeros(12).tolist())
+    np.testing.assert_raises(AssertionError, photo, np.zeros(12))
     with warnings.catch_warnings(record=True) as ws:
         output = photo(np.zeros(16))
         # one warning on Ubuntu, two warnings with Windows
         assert any(isinstance(w.message, RuntimeWarning) for w in ws)
     np.testing.assert_almost_equal(output, float("inf"))
 
 
+@pytest.mark.parametrize("method", ["clipping", "arctan", "tanh", "constraint"])  # type: ignore
+@pytest.mark.parametrize("name", ["bragg", "morpho", "chirped"])  # type: ignore
+def test_no_warning(name: str, method: str) -> None:
+    with warnings.catch_warnings(record=True) as w:
+        core.Photonics(name, 24, bounding_method=method)
+        assert not w, f"Got a warning at initialization: {w[0]}"
+
+
 @testing.parametrized(
-    morpho=("morpho", 100, 1.1647),
+    # morpho=("morpho", 100, 1.1647),  # too slow
     chirped=("chirped", 150, 0.94439),
     bragg=("bragg", 2.5, 0.93216),
 )
 def test_photonics_values(name: str, value: float, expected: float) -> None:
+    if name == "morpho" and os.environ.get("CIRCLECI", False):
+        raise SkipTest("Too slow in CircleCI")
     photo = core.Photonics(name, 16)
     np.testing.assert_almost_equal(photo(value * np.ones(16)), expected, decimal=4)
+
+
+@testing.parametrized(
+    morpho=("morpho", 1.127904),
+    chirped=("chirped", 0.937039),
+    bragg=("bragg", 0.96776)
+)
+def test_photonics_values_random(name: str, expected: float) -> None:
+    if name == "morpho" and os.environ.get("CIRCLECI", False):
+        raise SkipTest("Too slow in CircleCI")
+    size = 16 if name != "morpho" else 4
+    photo = core.Photonics(name, size)
+    np.random.seed(12)
+    x = np.random.normal(0, 1, size=size)
+    candidate = photo.parametrization.spawn_child().set_standardized_data(x)
+    np.testing.assert_almost_equal(photo(candidate.value), expected, decimal=4)