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

CHANGELOG.md

@@ -45,6 +45,7 @@ Attention: The newest changes should be on top -->
 - CI: make the Gemini PR reviewer actually review [#1140](https://github.com/RocketPy-Team/RocketPy/pull/1140)
 - MNT: declare dependency floors the package can actually run on [#1108](https://github.com/RocketPy-Team/RocketPy/pull/1108)
 - CI: build the docs for pull requests into develop as well [#1104](https://github.com/RocketPy-Team/RocketPy/pull/1104)
+- ENH: Make Monte Carlo input sampling reproducible per simulation index via a `random_seed` argument. Fixed-seed samples change, serial log indices are now zero-based to match the parallel path, and checkpoints written by the previous scheme cannot be resumed. [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054)
 - CI: make changelog automation LLM-based (Gemini) and race-safe [#1082](https://github.com/RocketPy-Team/RocketPy/pull/1082)
 - ENH: Resolve pressure_ISA discretization bounds TODO [#1056](https://github.com/RocketPy-Team/RocketPy/pull/1056)
 

rocketpy/stochastic/stochastic_model.py

@@ -3,6 +3,8 @@
 Stochastic classes.
 """
 
+from copy import deepcopy
+
 import numpy as np
 
 from rocketpy.mathutils.function import Function
@@ -11,6 +13,21 @@
 from ..tools import get_distribution
 
 
+def _snapshot_of(value):
+    """A nominal that writing through the wrapped object cannot reach.
+
+    Containers are copied, since ``obj.outline[:] = ...`` would otherwise move
+    the value this is meant to hold still. Anything else is kept by reference:
+    nothing here mutates a ``Function`` or a motor in place, and copying one per
+    model would cost more than it protects.
+    """
+    if isinstance(value, np.ndarray):
+        return value.copy()
+    if isinstance(value, (list, dict, set)):
+        return deepcopy(value)
+    return value
+
+
 def _names_as_spawn_key(input_names):
     """Encode names into spawn-key words that no other set of names produces.
 
@@ -95,6 +112,11 @@ class StochasticModel:
     # are, instead of each deciding for itself.
     array_valued_inputs = ()
 
+    # Collections a child class builds itself and overwrites in its own
+    # ``dict_generator``. Drawing from them here would only move the stream by
+    # however many components happen to be attached, and the draw is discarded.
+    component_collections = ()
+
     def __init__(self, obj, seed=None, **kwargs):
         """
         Initialize the StochasticModel class with validated input arguments.
@@ -122,8 +144,32 @@ 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 built.
+
+        Read once and kept, because ``StochasticEnvironment.create_object``
+        writes the randomised value back onto ``self.obj``. Re-reading it on a
+        reseed would take one simulation's output as the next one's nominal, so
+        a factor would multiply the factor before it. Containers are copied on
+        the way in, so rebinding the attribute and writing through it both leave
+        this where it was; anything else is held by reference and follows the
+        object (see ``_snapshot_of``).
+
+        A custom ``getter`` reads a component's own attribute, which nothing
+        writes back to, so it passes straight through. Caching those would be
+        wrong anyway: every component's position arrives under one name.
+        """
+        if getter is not getattr:
+            return getter(self.obj, input_name)
+        if input_name not in self.__nominal_values:
+            self.__nominal_values[input_name] = _snapshot_of(
+                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.
@@ -170,7 +216,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):
@@ -305,7 +351,7 @@ def _validate_tuple_length_two(self, input_name, input_value, getattr=getattr):
             # object passed.
             dist_func = get_distribution(input_value[1], self.__random_number_generator)
             return (
-                self._nominal_value(input_name, getattr(self.obj, input_name)),
+                self._nominal_value(input_name, self._nominal(input_name, getattr)),
                 input_value[0],
                 dist_func,
             )
@@ -381,7 +427,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
 
@@ -407,7 +453,7 @@ def _validate_scalar(self, input_name, input_value, getattr=getattr):  # pylint:
                 distribution function).
         """
         return (
-            self._nominal_value(input_name, getattr(self.obj, input_name)),
+            self._nominal_value(input_name, self._nominal(input_name, getattr)),
             input_value,
             get_distribution("normal", self.__random_number_generator),
         )
@@ -434,7 +480,7 @@ def _validate_factors(self, input_name, input_value):
             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)
@@ -697,6 +743,8 @@ def dict_generator(self):
         """
         generated_dict = {}
         for arg, value in self.__dict__.items():
+            if arg in self.component_collections:
+                continue
             if isinstance(value, tuple):
                 dist_sampler = value[-1]
                 generated_dict[arg] = dist_sampler(value[0], value[1])

rocketpy/stochastic/stochastic_rocket.py

@@ -2,6 +2,8 @@
 
 import warnings
 
+import numpy as np
+
 from rocketpy.control import _Controller
 from rocketpy.mathutils.vector_matrix import Vector
 from rocketpy.motors.empty_motor import EmptyMotor
@@ -21,6 +23,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,
@@ -90,6 +93,16 @@ class StochasticRocket(StochasticModel):
         can not be a randomized.
     """
 
+    # Overwritten in ``dict_generator`` below, so the base class must not draw
+    # from them.
+    component_collections = (
+        "motors",
+        "aerodynamic_surfaces",
+        "rail_buttons",
+        "air_brakes",
+        "parachutes",
+    )
+
     def __init__(
         self,
         rocket,
@@ -155,6 +168,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,
@@ -173,25 +193,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.
 
@@ -200,8 +242,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
         -------
@@ -213,7 +256,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),
@@ -456,8 +499,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):
@@ -482,14 +526,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.
 
@@ -740,6 +794,11 @@ def create_object(self):
         """Creates and returns a Rocket object from the randomly generated input
         arguments.
 
+        Motors, aerodynamic surfaces, air brakes and their controllers, rail
+        buttons and parachutes are carried onto the new rocket. Sensors are not:
+        one added to the wrapped rocket does not reach a Monte Carlo flight, and
+        so cannot affect it either.
+
         Returns
         -------
         rocket : Rocket

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)