ENH: reproducible Monte Carlo via per-simulation-index seeding

CHANGELOG.md

@@ -33,6 +33,7 @@ Attention: The newest changes should be on top -->
 ### Added
 
 - ENH: update master with develop [#1081](https://github.com/RocketPy-Team/RocketPy/pull/1081)
+- ENH: reproducible Monte Carlo runs via a random_seed argument [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054)
 
 ### Changed
 

requirements.txt

@@ -1,4 +1,4 @@
-numpy>=1.13
+numpy>=1.17
 scipy>=1.0
 matplotlib>=3.9.0  # Released May 15th 2024
 netCDF4>=1.6.4

rocketpy/stochastic/stochastic_model.py

@@ -3,8 +3,6 @@
 Stochastic classes.
 """
 
-from random import choice
-
 import numpy as np
 
 from rocketpy.mathutils.function import Function
@@ -68,8 +66,37 @@ def __init__(self, obj, seed=None, **kwargs):
         self.obj = obj
         self.last_rnd_dict = {}
         self.__stochastic_dict = kwargs
+        self.__nominal_values = {}
         self._set_stochastic(seed)
 
+    def _nominal(self, input_name, getter=getattr):
+        """``self.obj``'s value for ``input_name``, as it was when this model
+        was built.
+
+        Read once and remembered, because ``StochasticEnvironment`` has
+        ``create_object`` write the randomised value back onto ``self.obj``
+        instead of building a copy. Re-reading it on a reseed would take one
+        simulation's output as the next one's nominal, and a factor would
+        multiply the factor before it rather than the original value.
+
+        A custom ``getter`` reads a component's own attribute rather than one
+        of ``self.obj``'s, and nothing writes back to those, so it is passed
+        straight through. Caching it here would be wrong as well: every
+        component's position arrives under the one name ``"position"``.
+
+        This applies to every stochastic model, not only the environment: what
+        a model samples around is what the wrapped object held when the model
+        was built. Changing the object afterwards does not move it. Only
+        ``StochasticEnvironment.create_object`` writes back today, but the rule
+        is stated for all of them rather than special-cased for one, so a model
+        means the same thing whichever object it wraps.
+        """
+        if getter is not getattr:
+            return getter(self.obj, input_name)
+        if input_name not in self.__nominal_values:
+            self.__nominal_values[input_name] = getattr(self.obj, input_name)
+        return self.__nominal_values[input_name]
+
     def _set_stochastic(self, seed=None):
         """Set the stochastic attributes from the input dictionary.
         This method is useful to reset or reseed the attributes of the instance.
@@ -109,7 +136,7 @@ def _set_stochastic(self, seed=None):
                                 "or a custom sampler"
                             )
                 else:
-                    attr_value = [getattr(self.obj, input_name)]
+                    attr_value = [self._nominal(input_name)]
                 setattr(self, input_name, attr_value)
 
     def __repr__(self):
@@ -186,7 +213,7 @@ def _validate_tuple_length_two(self, input_name, input_value, getattr=getattr):
             # function. In this case, the nominal value will be taken from the
             # object passed.
             dist_func = get_distribution(input_value[1], self.__random_number_generator)
-            return (getattr(self.obj, input_name), input_value[0], dist_func)
+            return (self._nominal(input_name, getattr), input_value[0], dist_func)
         else:
             # if second item is an int or float, then it is assumed that the
             # first item is the nominal value and the second item is the
@@ -257,7 +284,7 @@ def _validate_list(self, input_name, input_value, getattr=getattr):  # pylint: d
             If the input is not in a valid format.
         """
         if not input_value:
-            return [getattr(self.obj, input_name)]
+            return [self._nominal(input_name, getattr)]
         else:
             return input_value
 
@@ -283,7 +310,7 @@ def _validate_scalar(self, input_name, input_value, getattr=getattr):  # pylint:
                 distribution function).
         """
         return (
-            getattr(self.obj, input_name),
+            self._nominal(input_name, getattr),
             input_value,
             get_distribution("normal", self.__random_number_generator),
         )
@@ -310,7 +337,7 @@ def _validate_factors(self, input_name, input_value, seed):
             If the input is not in a valid format.
         """
         attribute_name = input_name.replace("_factor", "")
-        setattr(self, f"_{attribute_name}", getattr(self.obj, attribute_name))
+        setattr(self, f"_{attribute_name}", self._nominal(attribute_name))
 
         if isinstance(input_value, tuple):
             return self._validate_tuple_factor(input_name, input_value)
@@ -508,6 +535,21 @@ def _validate_airfoil(self, airfoil):
                         "the first item"
                     )
 
+    def _random_choice(self, values):
+        """Choose one value from a list using this model's seeded generator.
+
+        The index is drawn from the seeded generator, not the stdlib global
+        ``random.choice`` (an unseeded shared instance), so the choice is
+        governed by ``random_seed``. Indexing rather than ``numpy.random.choice``
+        keeps a heterogeneous list -- ``Function`` objects, paths, arrays --
+        returned as itself instead of coerced to a common dtype. An empty
+        ``values`` is returned unchanged.
+        """
+        if not values:
+            return values
+        index = int(self.__random_number_generator.integers(len(values)))
+        return values[index]
+
     def dict_generator(self):
         """
         Generate a dictionary with randomly generated input arguments.
@@ -532,7 +574,7 @@ def dict_generator(self):
                 dist_sampler = value[-1]
                 generated_dict[arg] = dist_sampler(value[0], value[1])
             elif isinstance(value, list):
-                generated_dict[arg] = choice(value) if value else value
+                generated_dict[arg] = self._random_choice(value)
             elif isinstance(value, CustomSampler):
                 try:
                     generated_dict[arg] = value.sample(n_samples=1)[0]

rocketpy/stochastic/stochastic_rocket.py

@@ -1,7 +1,8 @@
 """Defines the StochasticRocket class."""
 
 import warnings
-from random import choice
+
+import numpy as np
 
 from rocketpy.control import _Controller
 from rocketpy.mathutils.vector_matrix import Vector
@@ -21,6 +22,7 @@
 from rocketpy.rocket.rocket import Rocket
 from rocketpy.stochastic.stochastic_generic_motor import StochasticGenericMotor
 from rocketpy.stochastic.stochastic_motor_model import StochasticMotorModel
+from rocketpy.tools import _seed_sequence_to_int
 
 from .stochastic_aero_surfaces import (
     StochasticAirBrakes,
@@ -154,6 +156,13 @@ def __init__(
         self.air_brakes = []
         self.parachutes = []
         self.__components_map = {}
+        # Raw eccentricity arguments, kept as the caller gave them.
+        # ``add_cp_eccentricity`` and ``add_thrust_eccentricity`` run after
+        # ``__init__``, so their values are not in the dict the base class
+        # re-validates on a reseed. Validating them once would leave the
+        # distribution bound to the Generator of whichever simulation happened
+        # to come first, so the raw form is kept and validated again each time.
+        self.__eccentricity_specs = {}
         super().__init__(
             obj=rocket,
             radius=radius,
@@ -172,25 +181,47 @@ def __init__(
             coordinate_system_orientation=None,
         )
 
+    # Every collection of nested stochastic objects, in the order their child
+    # seeds are spawned. Listed here rather than written out inline so that a
+    # component type cannot end up in ``create_object`` and not in the reseed:
+    # air brakes were, and their sampling depended on which worker ran the
+    # index instead of on the index. ``_stochastic_collections`` is asserted
+    # against the rocket's own attributes in the tests.
+    _POSITIONED_COLLECTIONS = ("aerodynamic_surfaces", "motors", "rail_buttons")
+    _PLAIN_COLLECTIONS = ("parachutes", "air_brakes")
+
+    @classmethod
+    def _stochastic_collections(cls):
+        """The names of every attribute holding nested stochastic objects."""
+        return cls._POSITIONED_COLLECTIONS + cls._PLAIN_COLLECTIONS
+
     def _set_stochastic(self, seed=None):
         """Set the stochastic attributes for Components, positions and
         inputs.
 
+        Every nested component -- the rocket body, each aerodynamic surface,
+        motor, rail button, parachute and air brake -- is reseeded from its own
+        child of a ``SeedSequence`` root, so components that sample the same
+        distribution do not draw identical values (a main and a drogue parachute
+        get independent ``cd_s`` and ``lag`` samples, not the same one). Children
+        are spawned in a fixed order, so the result stays reproducible under
+        ``random_seed``.
+
         Parameters
         ----------
         seed : int, optional
             Seed for the random number generator.
         """
-        super()._set_stochastic(seed)
-        self.aerodynamic_surfaces = self.__reset_components(
-            self.aerodynamic_surfaces, seed
-        )
-        self.motors = self.__reset_components(self.motors, seed)
-        self.rail_buttons = self.__reset_components(self.rail_buttons, seed)
-        for parachute in self.parachutes:
-            parachute._set_stochastic(seed)
-
-    def __reset_components(self, components, seed):
+        root = np.random.SeedSequence(seed)
+        super()._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0]))
+        self.__apply_eccentricity_specs()
+        for name in self._POSITIONED_COLLECTIONS:
+            setattr(self, name, self.__reset_components(getattr(self, name), root))
+        for name in self._PLAIN_COLLECTIONS:
+            for child in getattr(self, name):
+                child._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0]))
+
+    def __reset_components(self, components, root):
         """Creates a new Components whose stochastic structures
         and their positions are reset.
 
@@ -199,8 +230,9 @@ def __reset_components(self, components, seed):
         components : Components
             The components which contains the stochastic structure that
             will be used to create the new components.
-        seed : int, optional
-            Seed for the random number generator.
+        root : numpy.random.SeedSequence
+            The run's seed root. Each component is reseeded from its own spawned
+            child, so components sampling the same distribution stay decorrelated.
 
         Returns
         -------
@@ -212,7 +244,7 @@ def __reset_components(self, components, seed):
         new_components = Components()
         for stochastic_obj, _ in components:
             stochastic_obj_position_info = self.__components_map[stochastic_obj]
-            stochastic_obj._set_stochastic(seed)
+            stochastic_obj._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0]))
             new_components.add(
                 stochastic_obj,
                 self._validate_position(stochastic_obj, stochastic_obj_position_info),
@@ -434,8 +466,9 @@ def add_cp_eccentricity(self, x=None, y=None):
         self : StochasticRocket
             Object of the StochasticRocket class.
         """
-        self.cp_eccentricity_x = self._validate_eccentricity("cp_eccentricity_x", x)
-        self.cp_eccentricity_y = self._validate_eccentricity("cp_eccentricity_y", y)
+        self.__eccentricity_specs["cp_eccentricity_x"] = x
+        self.__eccentricity_specs["cp_eccentricity_y"] = y
+        self.__apply_eccentricity_specs()
         return self
 
     def add_thrust_eccentricity(self, x=None, y=None):
@@ -460,14 +493,24 @@ def add_thrust_eccentricity(self, x=None, y=None):
         self : StochasticRocket
             Object of the StochasticRocket class.
         """
-        self.thrust_eccentricity_x = self._validate_eccentricity(
-            "thrust_eccentricity_x", x
-        )
-        self.thrust_eccentricity_y = self._validate_eccentricity(
-            "thrust_eccentricity_y", y
-        )
+        self.__eccentricity_specs["thrust_eccentricity_x"] = x
+        self.__eccentricity_specs["thrust_eccentricity_y"] = y
+        self.__apply_eccentricity_specs()
         return self
 
+    def __apply_eccentricity_specs(self):
+        """Re-validate the eccentricities against the current Generator.
+
+        Validation stores a distribution as a method bound to the Generator
+        that was live at the time, so a tuple validated once keeps sampling
+        from that one. Re-running it after every reseed is what ties the draw
+        to the simulation index rather than to whichever index the worker
+        happened to run first. ``get_distribution`` only binds a method, so
+        this consumes no randomness and does not shift any other draw.
+        """
+        for name, spec in self.__eccentricity_specs.items():
+            setattr(self, name, self._validate_eccentricity(name, spec))
+
     def _validate_eccentricity(self, eccentricity, position):
         """Validate the eccentricity argument.
 
@@ -628,7 +671,7 @@ def _randomize_position(self, position):
                 return position[-1](position[0].z, position[1])
             return position[-1](position[0], position[1])
         elif isinstance(position, list):
-            return choice(position) if position else position
+            return self._random_choice(position)
 
     # pylint: disable=stop-iteration-return
     def dict_generator(self):
@@ -638,8 +681,8 @@ def dict_generator(self):
         all attributes of the class and generating a random value for each
         attribute. The random values are generated according to the format of
         each attribute. Tuples are generated using the distribution function
-        specified in the tuple. Lists are generated using the random.choice
-        function.
+        specified in the tuple. Lists are sampled through the model's seeded
+        generator so the choice is governed by ``random_seed``.
 
         Parameters
         ----------

rocketpy/tools.py

@@ -1467,6 +1467,29 @@ def find_obj_from_hash(obj, hash_, depth_limit=None):
     return None
 
 
+def _seed_sequence_to_int(seed_sequence):
+    """Collapse a ``SeedSequence`` into a 128-bit Python ``int``.
+
+    A plain ``int`` is what ``numpy.random.default_rng`` and the stdlib
+    ``random.Random`` both accept (``random.Random`` rejects a ``SeedSequence``
+    with a ``TypeError`` since Python 3.11), so a custom sampler whose
+    ``reset_seed`` documents an ``int`` and builds a modern generator keeps
+    working. The legacy ``numpy.random.RandomState`` is the exception: it caps a
+    single-integer seed at ``2**32 - 1``, so a sampler still built on it would
+    have to reduce the value (``RandomState`` is a frozen legacy API NumPy steers
+    new code away from). All four ``uint32`` words are combined to keep the full
+    128-bit pool, so sub-streams stay decorrelated instead of collapsing to a
+    single 32-bit word.
+
+    The words are combined by value (little-endian word order), not via
+    ``tobytes()``, so the seed is the same on big- and little-endian machines --
+    a byte-order-dependent seed would break the cross-platform reproducibility
+    this exists to provide.
+    """
+    words = seed_sequence.generate_state(4, dtype=np.uint32)
+    return sum(int(word) << (32 * position) for position, word in enumerate(words))
+
+
 if __name__ == "__main__":  # pragma: no cover
     import doctest
 

tests/fixtures/monte_carlo/custom_sampler_fixtures.py

@@ -33,7 +33,7 @@ def __init__(self, means_tuple, sd_tuple, prob_tuple, seed=None):
             2-Tuple that contains the probability of each normal distribution of the
             mixture. Its entries should be non-negative and sum up to 1.
         """
-        np.random.default_rng(seed)
+        self.reset_seed(seed)
         self.means_tuple = means_tuple
         self.sd_tuple = sd_tuple
         self.prob_tuple = prob_tuple
@@ -52,16 +52,12 @@ def sample(self, n_samples=1):
             List containing n_samples samples
         """
         samples_list = [0] * n_samples
-        mixture_id_list = np.random.binomial(1, self.prob_tuple[0], n_samples)
+        mixture_id_list = self.rng.binomial(1, self.prob_tuple[0], n_samples)
         for i, mixture_id in enumerate(mixture_id_list):
             if mixture_id:
-                samples_list[i] = np.random.normal(
-                    self.means_tuple[0], self.sd_tuple[0]
-                )
+                samples_list[i] = self.rng.normal(self.means_tuple[0], self.sd_tuple[0])
             else:
-                samples_list[i] = np.random.normal(
-                    self.means_tuple[1], self.sd_tuple[1]
-                )
+                samples_list[i] = self.rng.normal(self.means_tuple[1], self.sd_tuple[1])
 
         return samples_list
 
@@ -73,4 +69,7 @@ def reset_seed(self, seed=None):
         seed : int, optional
             Seed for the random number generator.
         """
-        np.random.default_rng(seed)
+        # Kept on the instance. Building a generator and dropping it made this
+        # a no-op, and sample() went on drawing from the process-global
+        # np.random, so nothing here answered to the seed at all.
+        self.rng = np.random.default_rng(seed)