ENH: seed Monte Carlo per simulation index for worker-invariant inputs

.github/copilot-instructions.md

@@ -160,12 +160,7 @@ docs/
 
 ### Function Definition
 ```python
-def calculate_drag_force(
-    velocity,
-    air_density,
-    drag_coefficient,
-    reference_area
-):
+def calculate_drag_force(velocity, air_density, drag_coefficient, reference_area):
     """Calculate drag force using the standard drag equation.
 
     Parameters
@@ -211,7 +206,9 @@ def test_calculate_drag_force_returns_correct_value():
     expected_force = 30.625  # N
 
     # Act
-    result = calculate_drag_force(velocity, air_density, drag_coefficient, reference_area)
+    result = calculate_drag_force(
+        velocity, air_density, drag_coefficient, reference_area
+    )
 
     # Assert
     assert abs(result - expected_force) < 1e-6

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,30 @@ 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"``.
+        """
+        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 +129,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 +206,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 +277,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 +303,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 +330,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 +528,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 +567,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]
@@ -550,6 +585,11 @@ def visualize_attributes(self):
         Model object. The report includes the variable name, the nominal value,
         the standard deviation, and the distribution function used to generate
         the random attributes.
+
+        Returns
+        -------
+        str
+            The formatted report. It is also printed for interactive use.
         """
 
         def format_attribute(attr, value):
@@ -630,4 +670,10 @@ def format_attribute(attr, value):
                 format_attribute(attr, attributes[attr]) for attr in custom_attributes
             )
 
-        print("\n".join(filter(None, report)))
+        # This is an explicit, user-invoked display method, so it prints
+        # unconditionally (like ``info``/``all_info`` elsewhere) rather than
+        # logging at INFO level, which is silenced by default. The report is
+        # also returned so it can be used programmatically.
+        report_str = "\n".join(filter(None, report))
+        print(report_str)
+        return report_str

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
@@ -22,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,
@@ -161,6 +163,13 @@ def __init__(  # pylint: disable=too-many-arguments
         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,
@@ -180,25 +189,47 @@ def __init__(  # pylint: disable=too-many-arguments
             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.
 
@@ -207,8 +238,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
         -------
@@ -220,7 +252,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),
@@ -472,8 +504,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):
@@ -498,14 +531,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.
 
@@ -666,7 +709,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):
@@ -676,8 +719,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
         ----------