Merge dev (#24)

.github/workflows/tests.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, macos-latest, windows-latest]
+        os: [ubuntu-latest, windows-latest]
         python-version: ["3.10"]
         include:
           - os: ubuntu-latest

docs/algo/segomoe.md

@@ -60,3 +60,40 @@ g = interface.g  # (n, ng)
 pop = interface.pop  # Population containing all design points
 opt = interface.opt  # Population containing optimal point(s)
 ```
+
+### pymoo API
+
+It is also possible to use the pymoo API to run an optimization:
+```python
+from pymoo.optimize import minimize
+from sb_arch_opt.algo.segomoe_interface import SEGOMOEInterface, SEGOMOEAlgorithm
+
+problem = ...  # Subclass of ArchOptProblemBase
+
+# Define folder to store results in
+results_folder = ...
+
+# Use Mixture of Experts: automatically identifies clusters in the design space
+# with different best surrogates ("experts"). Can be more accurate, however
+# also greatly increases the cost of finding new infill points.
+use_moe = True
+
+# Options passed to the Sego class and to model generation, respectively
+sego_options = {}
+model_options = {}
+
+# Get the interface (will be initialized if the results folder has results)
+interface = SEGOMOEInterface(problem, results_folder, n_init=100, n_infill=50,
+                             use_moe=use_moe, sego_options=sego_options,
+                             model_options=model_options)
+
+# Define the pymoo Algorithm
+algo = SEGOMOEAlgorithm(interface)
+
+# Initialize from other results if you want
+algo.initialize_from_previous_results(problem, results_folder='/optional/other/result/folder')
+
+# Run the optimization
+# Note: no need to give a termination, as that is already defined by the SEGOMOEInterface object (n_init + n_infill)
+result = minimize(problem, algo, seed=42)  # Remove seed in production
+```

sb_arch_opt/__init__.py

@@ -1 +1 @@
-__version__ = '1.5.0'
+__version__ = '1.5.1'

sb_arch_opt/algo/segomoe_interface/algo.py

@@ -25,7 +25,7 @@
 import os
 import logging
 import numpy as np
-from typing import Tuple
+from typing import Tuple, Optional
 from sb_arch_opt.sampling import *
 from sb_arch_opt.algo.arch_sbo.models import *
 
@@ -136,13 +136,17 @@ def h(self) -> np.ndarray:
     @property
     def pop(self) -> Population:
         """Population of all evaluated points"""
-        return self.get_population(self.x, self.y)
+        return self.get_population(self.x, self.y, self.x_failed)
 
     @property
     def opt(self) -> Population:
         """Optimal points (Pareto front if multi-objective)"""
         return self._get_pareto_front(self.pop)
 
+    @property
+    def results_folder(self):
+        return self._results_folder
+
     def initialize_from_previous(self, results_folder: str = None):
         capture_log()
         if results_folder is None:
@@ -175,29 +179,65 @@ def initialize_from_previous(self, results_folder: str = None):
 
         log.info('No previous results found')
 
+    def set_pop(self, pop: Population = None):
+        if pop is None:
+            self._x = self._x_failed = self._y = None
+        else:
+            self._x, self._x_failed, self._y = self._get_xy(pop)
+
     def run_optimization(self):
         capture_log()
 
+        n_doe, n_infills = self._optimization_step()
+        if n_doe is not None:
+            log.info(f'Running DOE of {n_doe} points ({self.n_init} total)')
+            self.run_doe(n_doe)
+
+        if n_infills is not None:
+            log.info(f'Running optimization: {n_infills} infill points (ok DOE points: {self.n})')
+            self.run_infills(n_infills)
+
+        # Save final results and return Pareto front
+        self._save_results()
+        return self.opt
+
+    def _optimization_step(self):
         # Automatically initialize from previous results if reusing the same storage folder
         if self._x is None:
             self.initialize_from_previous()
 
         # Run DOE if needed
         n_available = self.n_tried
+        n_doe = None
         if n_available < self.n_init:
-            log.info(f'Running DOE of {self.n_init-n_available} points ({self.n_init} total)')
-            self.run_doe(self.n_init-n_available)
+            n_doe = self.n_init-n_available
 
-        # Run optimization
-        n_available = self.n_tried
+        # Run optimization (infills)
+        n_available = self.n_tried + (n_doe or 0)
+        n_infills = None
         if n_available < self.n_init+self.n_infill:
             n_infills = self.n_infill - (n_available-self.n_init)
-            log.info(f'Running optimization: {n_infills} infill points (ok DOE points: {self.n})')
-            self.run_infills(n_infills)
 
-        # Save final results and return Pareto front
-        self._save_results()
-        return self.opt
+        return n_doe, n_infills
+
+    def optimization_has_ask(self):
+        n_doe, n_infills = self._optimization_step()
+        return n_doe is not None or n_infills is not None
+
+    def optimization_ask(self) -> Optional[np.ndarray]:
+        n_doe, n_infills = self._optimization_step()
+
+        if n_doe is not None:
+            return self._sample_doe(n_doe)
+
+        if n_infills is not None:
+            return np.array([self._ask_infill()])
+
+    def optimization_tell_pop(self, pop: Population):
+        self._tell_infill(*self._get_xy(pop))
+
+    def optimization_tell(self, x, x_failed, y):
+        self._tell_infill(x, x_failed, y)
 
     def run_doe(self, n: int = None):
         if n is None:
@@ -227,10 +267,7 @@ def _grouped_eval(x):
 
              x, x_failed, y = self._get_xy(self._evaluate(np.array([x])))
 
-             self._x = np.row_stack([self._x, x])
-             self._y = np.row_stack([self._y, y])
-             self._x_failed = np.row_stack([self._x_failed, x_failed])
-             self._save_results()
+             self._tell_infill(x, x_failed, y)
 
              if len(x_failed) > 0:
                  return [], True
@@ -256,10 +293,7 @@ def run_infills_ask_tell(self, n_infills: int = None):
             x, x_failed, y = self._get_xy(self._evaluate(np.array([x])))
 
             # Update and save DOE
-            self._x = np.row_stack([self._x, x])
-            self._y = np.row_stack([self._y, y])
-            self._x_failed = np.row_stack([self._x_failed, x_failed])
-            self._save_results()
+            self._tell_infill(x, x_failed, y)
 
     def _ask_infill(self) -> np.ndarray:
         """
@@ -280,6 +314,12 @@ def _dummy_f_grouped(_):
         # Return latest point as suggested infill point
         return sego.get_x(i=-1)
 
+    def _tell_infill(self, x, x_failed, y):
+        self._x = np.row_stack([self._x, x]) if self._x is not None else x
+        self._y = np.row_stack([self._y, y]) if self._y is not None else y
+        self._x_failed = np.row_stack([self._x_failed, x_failed]) if self._x_failed is not None else x_failed
+        self._save_results()
+
     def _get_sego(self, f_grouped):
         design_space_spec = self._get_design_space()
 
@@ -419,9 +459,16 @@ def _flip_g(self, y: np.ndarray):
         g = -g
         return np.column_stack([f, g, h])
 
-    def get_population(self, x: np.ndarray, y: np.ndarray) -> Population:
+    def get_population(self, x: np.ndarray, y: np.ndarray, x_failed: np.ndarray = None) -> Population:
         # Inequality constraint values are flipped to correctly calculate constraint violation values in pymoo
         f, g, h = self._split_y(y)
+
+        if x_failed is not None and len(x_failed) > 0:
+            x = np.row_stack([x, x_failed])
+            f = np.row_stack([f, np.zeros((x_failed.shape[0], f.shape[1]))*np.inf])
+            g = np.row_stack([g, np.zeros((x_failed.shape[0], g.shape[1]))*np.inf])
+            h = np.row_stack([h, np.zeros((x_failed.shape[0], h.shape[1]))*np.inf])
+
         kwargs = {'X': x, 'F': f, 'G': g, 'H': h}
         pop = Population.new(**kwargs)
         return pop

sb_arch_opt/algo/segomoe_interface/api.py

@@ -23,4 +23,6 @@
 SOFTWARE.
 """
 from sb_arch_opt.algo.segomoe_interface.algo import *
-__all__ = ['HAS_SEGOMOE', 'HAS_SMT', 'SEGOMOEInterface']
+from sb_arch_opt.algo.segomoe_interface.pymoo_algo import *
+
+__all__ = ['HAS_SEGOMOE', 'HAS_SMT', 'SEGOMOEInterface', 'SEGOMOEAlgorithm']

sb_arch_opt/algo/segomoe_interface/pymoo_algo.py

@@ -0,0 +1,125 @@
+"""
+MIT License
+
+Copyright: (c) 2024, Deutsches Zentrum fuer Luft- und Raumfahrt e.V.
+Contact: jasper.bussemaker@dlr.de
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+"""
+import logging
+from sb_arch_opt.problem import *
+from sb_arch_opt.algo.pymoo_interface import *
+from sb_arch_opt.algo.pymoo_interface.metrics import EHVMultiObjectiveOutput
+
+from pymoo.core.algorithm import Algorithm
+from pymoo.core.population import Population
+from pymoo.util.optimum import filter_optimum
+from pymoo.termination.max_eval import MaximumFunctionCallTermination
+
+from sb_arch_opt.algo.segomoe_interface.algo import SEGOMOEInterface, check_dependencies
+
+__all__ = ['SEGOMOEAlgorithm']
+
+log = logging.getLogger('sb_arch_opt.segomoe')
+
+
+class SEGOMOEAlgorithm(Algorithm):
+    """
+    Algorithm that wraps the SEGOMOE interface.
+
+    The population state is managed here, and each time infill points are asked for the SEGOMOE population is updated
+    from the algorithm population.
+    """
+
+    def __init__(self, segomoe: SEGOMOEInterface, output=EHVMultiObjectiveOutput(), **kwargs):
+        check_dependencies()
+        super().__init__(output=output, **kwargs)
+        self.segomoe = segomoe
+
+        self.termination = MaximumFunctionCallTermination(self.segomoe.n_init + self.segomoe.n_infill)
+        self._store_intermediate_results()
+
+        self.initialization = None  # Enable DOE override
+
+    def _initialize_infill(self):
+        if self.initialization is not None:
+            return self.initialization.do(self.problem, self.segomoe.n_init, algorithm=self)
+        return self._infill()
+
+    def _initialize_advance(self, infills=None, **kwargs):
+        self._advance(infills=infills, **kwargs)
+
+    def has_next(self):
+        if not super().has_next():
+            return False
+
+        self._infill_set_pop()
+        if not self.segomoe.optimization_has_ask():
+            return False
+        return True
+
+    def _infill(self):
+        self._infill_set_pop()
+        x_infill = self.segomoe.optimization_ask()
+        off = Population.new(X=x_infill) if x_infill is not None else Population.new()
+
+        # Stop if no new offspring is generated
+        if len(off) == 0:
+            self.termination.force_termination = True
+
+        return off
+
+    def _infill_set_pop(self):
+        if self.pop is None or len(self.pop) == 0:
+            self.segomoe.set_pop(pop=None)
+        else:
+            self.segomoe.set_pop(self.pop)
+
+    def _advance(self, infills=None, **kwargs):
+        if infills is not None:
+            self.segomoe.optimization_tell_pop(infills)
+        self.pop = self.segomoe.pop
+
+    def _set_optimum(self):
+        pop = self.pop
+        i_failed = ArchOptProblemBase.get_failed_points(pop)
+        valid_pop = pop[~i_failed]
+        if len(valid_pop) == 0:
+            self.opt = Population.new(X=[None])
+        else:
+            self.opt = filter_optimum(valid_pop, least_infeasible=True)
+
+    def _store_intermediate_results(self):
+        """Enable intermediate results storage to support restarting"""
+        results_folder = self.segomoe.results_folder
+        self.evaluator = ArchOptEvaluator(results_folder=results_folder)
+        self.callback = ResultsStorageCallback(results_folder, callback=self.callback)
+
+    def initialize_from_previous_results(self, problem: ArchOptProblemBase, results_folder: str = None) -> bool:
+        """Initialize the SBO algorithm from previously stored results"""
+        if results_folder is None:
+            results_folder = self.segomoe.results_folder
+
+        population = load_from_previous_results(problem, results_folder)
+        if population is None:
+            return False
+
+        self.pop = population
+        self._set_optimum()
+        return True