further improved dynamic control surrogate modelling output

moptipyapps/dynamic_control/experiment_surrogate.py

@@ -55,10 +55,9 @@ def make_instances() -> Iterable[Callable[[], SystemModel]]:
     res: list[Callable[[], SystemModel]] = []
     for system in [THREE_COUPLED_OSCILLATORS]:
         controllers = [
-            make_ann(system.state_dims, system.control_dims, [2, 2]),
             make_ann(system.state_dims, system.control_dims, [3, 3])]
         for controller in controllers:
-            for ann_model in [[2, 2], [3, 3]]:
+            for ann_model in [[7, 6], [6, 6, 6]]:
                 res.append(cast(
                     Callable[[], SystemModel],
                     lambda _s=system, _c=controller, _m=make_ann(
@@ -103,20 +102,22 @@ def cmaes_raw(instance: Instance, max_fes: int = MAX_FES) -> Execution:
 def cmaes_surrogate(instance: SystemModel,
                     max_fes: int = MAX_FES,
                     fes_for_training: int = 128,
-                    fes_per_model_run: int = 128) -> Execution:
+                    fes_per_model_run: int = 128,
+                    fancy_logs: bool = True) -> Execution:
     """
     Create the Bi-Pop-CMA-ES setup.
 
     :param instance: the problem instance
     :param max_fes: the maximum FEs
     :param fes_for_training: the FEs for training
     :param fes_per_model_run: the FEs per model run
+    :param fancy_logs: should we do fancy logging?
     :return: the setup
     """
     execution, objective, space = base_setup(instance, max_fes)
     return execution.set_solution_space(space).set_algorithm(SurrogateCmaEs(
         instance, space, objective, max_fes // 4,
-        fes_for_training, fes_per_model_run))
+        fes_for_training, fes_per_model_run, fancy_logs))
 
 
 def on_completion(instance: Any, log_file: Path, process: Process) -> None:
@@ -170,12 +171,15 @@ def run(base_dir: str, n_runs: int = 5) -> None:
         perform_pre_warmup=False,
         on_completion=on_completion)
 
-    for fes in [128, 256, 512, 1024]:
+    for training_fes, run_fes in ((2 ** 13, 2 ** 8),
+                                  (2 ** 16, 2 ** 10)):
         run_experiment(
-            base_dir=use_dir.resolve_inside(f"model_for_{fes}x{fes}_fes"),
+            base_dir=use_dir.resolve_inside(
+                f"model_for_{training_fes}x{run_fes}_fes"),
             instances=instances,
-            setups=[cast(Callable[[Any], Execution], lambda i, __f=fes:
-                         cmaes_surrogate(i, MAX_FES, __f, __f))],
+            setups=[cast(Callable[[Any], Execution],
+                         lambda i, __t=training_fes, __r=run_fes:
+                         cmaes_surrogate(i, MAX_FES, __t, __r))],
             n_runs=n_runs,
             n_threads=Parallelism.ACCURATE_TIME_MEASUREMENTS,
             perform_warmup=False,

moptipyapps/dynamic_control/instance.py

@@ -104,7 +104,7 @@ def log_parameters_to(self, logger: KeyValueLogSection) -> None:
     def describe_parameterization(
             self, title: str | None,
             parameters: np.ndarray, base_name: str,
-            dest_dir: str) -> tuple[Path, Path]:
+            dest_dir: str) -> tuple[Path, ...]:
         """
         Describe the performance of a given system of system.
 

moptipyapps/dynamic_control/model_objective.py

@@ -28,22 +28,22 @@
 The model objective function is used by the
 :mod:`~moptipyapps.dynamic_control.surrogate_cma` algorithm.
 """
-
 from typing import Callable, Final
 
 import numba  # type: ignore
 import numpy as np
 from moptipy.api.objective import Objective
+from moptipy.utils.logger import KeyValueLogSection
 
 from moptipyapps.dynamic_control.controller import Controller
 from moptipyapps.dynamic_control.objective import FigureOfMerit
 
 
-@numba.njit(cache=True, inline="always", fastmath=True, boundscheck=False)
+@numba.njit(cache=False, inline="always", fastmath=True, boundscheck=False)
 def _evaluate(x: np.ndarray, dim: int, pin: np.ndarray,
               pout: np.ndarray, res: np.ndarray,
               eq: Callable[[np.ndarray, float, np.ndarray,
-                            np.ndarray], None]) -> np.ndarray:
+                            np.ndarray], None]) -> float:
     """
     Compute the squared differences between expected and actual model output.
 
@@ -52,14 +52,15 @@ def _evaluate(x: np.ndarray, dim: int, pin: np.ndarray,
     :param pin: the input vectors
     :param pout: the expected output vectors, flattened
     :param eq: the equations
-    :return: the vector of squared differences
+    :return: the mean of the `log(x+1)` of the squared differences `x`
     """
     idx: int = 0
     for row in pin:  # iterate over all row=(s, c) tuples
         nidx: int = idx + dim
         eq(row, 0.0, x, res[idx:nidx])  # store the equation results
         idx = nidx
-    return np.square(np.subtract(res, pout, res), res)  # compute squared diff
+    return np.log1p(np.square(  # compute the log(x+1) of the squared diff x
+        np.subtract(res, pout, res), res), res).mean()
 
 
 class ModelObjective(Objective):
@@ -94,14 +95,6 @@ def __init__(self, real: FigureOfMerit,
         :param real: the objective used for the real optimization problem
         :param model: the model
         """
-        n1 = str(real)
-        n2 = "figureOfMerit"
-        if n1.startswith(n2):
-            n1 = n1[len(n2):]
-        #: the name of this objective
-        self.name: Final[str] = f"model{n1}"
-        #: the result summation method
-        self.__sum: Callable[[np.ndarray], float] = real.sum_up_results
         #: the equations of the model
         self.__equations: Callable[[np.ndarray, float, np.ndarray,
                                     np.ndarray], None] = model.controller
@@ -116,6 +109,8 @@ def __init__(self, real: FigureOfMerit,
             np.ndarray, np.ndarray]]] = real.get_differentials
         #: the dimension of the output
         self.__dim: Final[int] = real.instance.system.state_dims
+        #: the real figure of merit name
+        self.__real_name: Final[str] = str(real)
 
     def begin(self) -> None:
         """
@@ -138,8 +133,8 @@ def evaluate(self, x: np.ndarray) -> float:
         :param x: the model parameterization
         :return: the objective value
         """
-        return self.__sum(_evaluate(x, self.__dim, self.__in, self.__out,
-                                    self.__res, self.__equations))
+        return _evaluate(x, self.__dim, self.__in, self.__out,
+                         self.__res, self.__equations)
 
     def end(self) -> None:
         """End a model optimization run and free the associated memory."""
@@ -153,7 +148,7 @@ def __str__(self) -> str:
 
         :return: the name of this objective
         """
-        return self.name
+        return "model"
 
     def lower_bound(self) -> float:
         """
@@ -162,3 +157,14 @@ def lower_bound(self) -> float:
         :returns: 0.0
         """
         return 0.0
+
+    def log_parameters_to(self, logger: KeyValueLogSection) -> None:
+        """
+        Log all parameters of this component as key-value pairs.
+
+        :param logger: the logger for the parameters
+        """
+        super().log_parameters_to(logger)
+        logger.key_value("figureOfMeritName", self.__real_name)
+        logger.key_value("nSamples",
+                         0 if self.__in is None else len(self.__in))

moptipyapps/dynamic_control/objective.py

@@ -23,6 +23,7 @@
 import numpy as np
 from moptipy.api.objective import Objective
 from moptipy.utils.logger import KeyValueLogSection
+from moptipy.utils.nputils import array_to_str
 from moptipy.utils.types import type_error
 
 from moptipyapps.dynamic_control.instance import Instance
@@ -235,6 +236,13 @@ def log_parameters_to(self, logger: KeyValueLogSection) -> None:
         """
         super().log_parameters_to(logger)
         logger.key_value("modelModeEnabled", self.__collection is not None)
+        logger.key_value("dataCollecting", self.__collect)
+        eq: Final = self.__equations
+        logger.key_value("usingOriginalEquations",
+                         eq is self.instance.system.equations)
+        mp: Final[str] = "modelParameters"
+        if hasattr(self.__equations, mp):
+            logger.key_value(mp, array_to_str(getattr(eq, mp)))
         with logger.scope(SCOPE_INSTANCE) as scope:
             self.instance.log_parameters_to(scope)
 

moptipyapps/dynamic_control/surrogate_cma.py

@@ -66,23 +66,29 @@
 """
 
 
+from copy import copy
+from gc import collect
+from os.path import basename
 from typing import Callable, Final
 
 import numba  # type: ignore
 import numpy as np
 from moptipy.algorithms.so.vector.cmaes_lib import BiPopCMAES
 from moptipy.api.algorithm import Algorithm
 from moptipy.api.execution import Execution
+from moptipy.api.logging import FILE_SUFFIX
 from moptipy.api.process import Process
 from moptipy.api.subprocesses import for_fes
 from moptipy.spaces.vectorspace import VectorSpace
 from moptipy.utils.logger import KeyValueLogSection
-from moptipy.utils.nputils import DEFAULT_FLOAT, rand_seed_generate
+from moptipy.utils.nputils import rand_seed_generate
+from moptipy.utils.path import Path
 from moptipy.utils.types import check_int_range, type_error
 from numpy.random import Generator
 
 from moptipyapps.dynamic_control.model_objective import ModelObjective
 from moptipyapps.dynamic_control.objective import FigureOfMerit
+from moptipyapps.dynamic_control.system import System
 from moptipyapps.dynamic_control.system_model import SystemModel
 
 
@@ -99,7 +105,8 @@ def __init__(self,
                  objective: FigureOfMerit,
                  fes_for_warmup: int,
                  fes_for_training: int,
-                 fes_per_model_run: int) -> None:
+                 fes_per_model_run: int,
+                 fancy_logs: bool = False) -> None:
         """
         Initialize the algorithm.
 
@@ -113,6 +120,7 @@ def __init__(self,
             model
         :param fes_per_model_run: the number of FEs to be applied to each
             optimization run on the model
+        :param fancy_logs: should we perform fancy logging?
         """
         super().__init__()
 
@@ -123,7 +131,11 @@ def __init__(self,
                 controller_space, "controller_space", VectorSpace)
         if not isinstance(objective, FigureOfMerit):
             raise type_error(objective, "objective", FigureOfMerit)
+        if not isinstance(fancy_logs, bool):
+            raise type_error(fancy_logs, "fancy_logs", bool)
 
+        #: should we do fancy logging?
+        self.fancy_logs: Final[bool] = fancy_logs
         #: the number of objective function evaluations to be used for warmup
         self.fes_for_warmup: Final[int] = check_int_range(
             fes_for_warmup, "fes_for_warmup", 1, 1_000_000)
@@ -206,9 +218,6 @@ def solve(self, process: Process) -> None:
         model_equations: Final[Callable[[
             np.ndarray, float, np.ndarray, np.ndarray], None]] =\
             self.system_model.model.controller
-        merged: Final[np.ndarray] = np.empty(
-            self.system_model.system.state_dims
-            + self.system_model.system.control_dims, DEFAULT_FLOAT)
         raw: Final[FigureOfMerit] = self.__control_objective
         random: Final[Generator] = process.get_random()
         training_execute: Final[Execution] = (
@@ -224,15 +233,57 @@ def solve(self, process: Process) -> None:
         result: Final[np.ndarray] = self.__control_cma.space.create()
         orig_init: Callable = raw.initialize
 
+# Get a log dir if logging is enabled and set up all the logging information.
+        log_dir_name: str | None = process.get_log_basename() \
+            if self.fancy_logs else None
+        model_training_dir: Path | None = None
+        model_training_log_name: str | None = None
+        models_dir: Path | None = None
+        models_name: str | None = None
+        tempsys: System | None = None
+        ctrl_dir: Path | None = None
+        ctrl_log_name: str | None = None
+        controllers_dir: Path | None = None
+        controllers_name: str | None = None
+        if log_dir_name is not None:
+            log_dir: Final[Path] = Path.path(log_dir_name)
+            log_dir.ensure_dir_exists()
+            prefix: str = "modelTraining"
+            model_training_dir = log_dir.resolve_inside(prefix)
+            model_training_dir.ensure_dir_exists()
+            base_name: Final[str] = basename(log_dir_name)
+            model_training_log_name = f"{base_name}_{prefix}_"
+            training_execute.set_log_improvements(True)
+            prefix = "controllerOnModel"
+            ctrl_dir = log_dir.resolve_inside(prefix)
+            ctrl_dir.ensure_dir_exists()
+            ctrl_log_name = f"{base_name}_{prefix}_"
+            on_model_execute.set_log_improvements(True)
+            prefix = "model"
+            models_dir = log_dir.resolve_inside(prefix)
+            models_dir.ensure_dir_exists()
+            tempsys = copy(self.system_model.system)
+            models_name = f"{tempsys.name}_{prefix}_"
+            prefix = "controllersOnReal"
+            controllers_dir = log_dir.resolve_inside(prefix)
+            controllers_dir.ensure_dir_exists()
+            controllers_name = f"{base_name}_{prefix}_"
+
 # Now we do the setup run that creates some basic results and
 # gathers the initial information for modelling the system.
         with for_fes(process, self.fes_for_warmup) as prc:
             self.__control_cma.solve(prc)
 
         while not should_terminate():  # until budget exhausted
+            consumed_fes: int = process.get_consumed_fes()
 
             # We now train a model on the data that was gathered.
             training_execute.set_rand_seed(rand_seed_generate(random))
+            if model_training_dir is not None:
+                training_execute.set_log_file(
+                    model_training_dir.resolve_inside(
+                        f"{model_training_log_name}{consumed_fes}"
+                        f"{FILE_SUFFIX}"))
             model_objective.begin()  # get the collected data
             with training_execute.execute() as sub:  # train model
                 sub.get_copy_of_best_y(model)  # get best model
@@ -246,24 +297,50 @@ def solve(self, process: Process) -> None:
                         boundscheck=False)
             def __new_model(state: np.ndarray, time: float,
                             control: np.ndarray, out: np.ndarray,
-                            _merged=merged, _params=model,
-                            _eq=model_equations) -> None:
-                sd: Final[int] = len(state)
-                _merged[0:sd] = state
-                _merged[sd:] = control
-                _eq(_merged, time, _params, out)
+                            _params=model, _eq=model_equations) -> None:
+                _eq(np.hstack((state, control)), time, _params, out)
+
+            setattr(__new_model, "modelParameters", model)  # see objective
+
+            if tempsys is not None:  # plot the model behavior
+                tempsys.equations = __new_model  # type: ignore
+                setattr(tempsys, "name", f"{models_name}{consumed_fes}")
+                tempsys.describe_system_without_control(models_dir)
+
+            collect()  # now we collect all garbage ... there should be much
 
 # OK, now that we got the model, we can perform the model optimization run.
             raw.set_model(__new_model)  # switch to use the model
             on_model_execute.set_rand_seed(rand_seed_generate(random))
+            if ctrl_dir is not None:
+                on_model_execute.set_log_file(ctrl_dir.resolve_inside(
+                    f"{ctrl_log_name}{consumed_fes}{FILE_SUFFIX}"))
             with on_model_execute.execute() as ome:
                 ome.get_copy_of_best_y(result)  # get best controller
             raw.set_raw()  # switch to the actual problem and data collection
             setattr(raw, "initialize", orig_init)  # allow resetting to "raw"
 
+            if tempsys is not None:  # plot the controller on that model
+                setattr(tempsys, "name", f"{models_name}{consumed_fes}")
+                tempsys.describe_system(
+                    f"{models_name}{consumed_fes}",
+                    self.system_model.controller.controller, result,
+                    f"{models_name}{consumed_fes}_synthesized_controller",
+                    models_dir)
+
 # Finally, we re-evaluate the result that we got from the model run on the
 # actual objective function.
-            process.evaluate(result)
+            process.evaluate(result)  # get the real objective value
+
+# plot the actual behavior
+            if controllers_dir is not None:
+                self.system_model.system.describe_system(
+                    f"{self.system_model.system}_{consumed_fes}",
+                    self.system_model.controller.controller,
+                    result, f"{controllers_name}{consumed_fes}",
+                    controllers_dir)
+
+            collect()  # now we collect all garbage ... there should be much
 
     def __str__(self):
         """
@@ -285,6 +362,7 @@ def log_parameters_to(self, logger: KeyValueLogSection) -> None:
         logger.key_value("fesForWarmup", self.fes_for_warmup)
         logger.key_value("fesForTraining", self.fes_for_training)
         logger.key_value("fesPerModelRun", self.fes_per_model_run)
+        logger.key_value("fancyLogs", self.fancy_logs)
         with logger.scope("ctrlCMA") as ccma:
             self.__control_cma.log_parameters_to(ccma)
         with logger.scope("mdlSpace") as mspce: