diff --git a/docs/user/stochastic.rst b/docs/user/stochastic.rst index 6e3376236..904737b0d 100644 --- a/docs/user/stochastic.rst +++ b/docs/user/stochastic.rst @@ -278,6 +278,45 @@ reliability of your simulations over time. .. which parameters most significantly impact your simulation results. +Seeding a rocket's components +----------------------------- + +A ``StochasticRocket`` holds nested stochastic objects: the motors, the +aerodynamic surfaces, the rail buttons, the parachutes and the air brakes. Each +time the rocket is reset, every component attached to it at that moment is given +a stream of its own, spawned from the seed the reset was given. A main and a +drogue parachute built from the same ``cd_s`` and ``lag`` are then no longer +made to consume the same draws as each other. Independent streams can still land +on equal values, and a specification with no spread always will. + +The reset is what builds the tree, not ``add_parachute`` or ``add_nose``. A +rocket resets itself once while being constructed, when it has no components +yet, so anything attached afterwards keeps the generator it was built with until +something resets it again. A serial ``MonteCarlo`` run does not reset it at +all. A parallel one tries to, once per worker, but hands the model a +``SeedSequence`` where an integer is wanted, so that path does not get as far +as building the tree either. Resetting per simulation, from an integer seed the +caller chooses, is what the Monte Carlo seeding work adds. + +Each collection is spawned separately, so adding an aerodynamic surface does not +move what the parachutes draw. Within a collection the stream follows insertion +order, so adding a component ahead of another does change what the later one +draws under a fixed seed. The rocket's own inputs, such as ``mass`` and +``radius``, use the seed exactly as given. + +.. note:: + A component's *position* is a property of the rocket rather than of the + component, so it is drawn from the rocket's own stream. + +.. note:: + A stream belongs to one stochastic wrapper. Storing the same wrapper twice, + or sharing one between two rockets, is not supported: the second reset + replaces the first, and the two entries end up drawing from one generator. + + A shared ``CustomSampler.seed_group`` keeps its own rule: a group belongs to + one model. Sharing one between two components leaves each of them seeding it + from their own child, and the last one to be reset decides what both draw. + Conclusion ---------- diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 65cfb5ebe..9bf9d5d36 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -2,6 +2,7 @@ import warnings + from rocketpy.control import _Controller from rocketpy.mathutils.vector_matrix import Vector from rocketpy.motors.empty_motor import EmptyMotor @@ -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_from, _seed_sequence_to_int from .stochastic_aero_surfaces import ( StochasticAirBrakes, @@ -173,25 +175,45 @@ def __init__( coordinate_system_orientation=None, ) + # Nested stochastic objects, in spawn order. A collection's stream is + # addressed by where its name falls in the two tuples read end to end, so + # appending to the first one moves every name in the second. + _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. + Each component takes its own child of its collection's root, so + components stay independent and one seed still reproduces the rocket. + The body keeps the seed as given, and a collection has a root of its + own, so adding a fin does not move every parachute. + 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) + names = self._stochastic_collections() + roots = dict(zip(names, _seed_sequence_from(seed).spawn(len(names)))) + for name in self._POSITIONED_COLLECTIONS: + setattr( + self, name, self.__reset_components(getattr(self, name), roots[name]) + ) + for name in self._PLAIN_COLLECTIONS: + for component in getattr(self, name): + component._set_stochastic( + _seed_sequence_to_int(roots[name].spawn(1)[0]) + ) - def __reset_components(self, components, seed): + def __reset_components(self, components, root): """Creates a new Components whose stochastic structures and their positions are reset. @@ -200,8 +222,8 @@ 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 reseed's root. Each component takes its own spawned child. Returns ------- @@ -213,7 +235,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), diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 0d7f1a74e..c074be9e8 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -1467,6 +1467,26 @@ def find_obj_from_hash(obj, hash_, depth_limit=None): return None +def _seed_sequence_to_int(seed_sequence): + """Returns a ``SeedSequence`` as the 128-bit ``int`` seed ``random.Random`` + accepts, combined by value so it does not depend on byte order. + """ + words = seed_sequence.generate_state(4, dtype=np.uint32) + return sum(int(word) << (32 * position) for position, word in enumerate(words)) + + +def _seed_sequence_from(seed): + """Returns a ``SeedSequence`` of the caller's own to spawn from. + + A parallel run is handed one that ``SeedSequence`` will not take as + entropy, and spawning from it directly would advance the counter of an + object the caller still holds, so it is copied from its full state. + """ + if isinstance(seed, np.random.SeedSequence): + return np.random.SeedSequence(**seed.state) + return np.random.SeedSequence(seed) + + if __name__ == "__main__": # pragma: no cover import doctest diff --git a/tests/unit/stochastic/test_stochastic_rocket_seeding.py b/tests/unit/stochastic/test_stochastic_rocket_seeding.py new file mode 100644 index 000000000..a3430762b --- /dev/null +++ b/tests/unit/stochastic/test_stochastic_rocket_seeding.py @@ -0,0 +1,313 @@ +"""Every nested component of a StochasticRocket is reseeded from its own child +of the run's seed, so components that sample the same distribution do not draw +the same values, and one seed still reproduces the whole rocket. +""" + +import ast +import inspect + +import numpy as np + +from rocketpy.rocket.components import Components +from rocketpy.tools import _seed_sequence_from +from rocketpy.stochastic import ( + StochasticAirBrakes, + StochasticParachute, + StochasticRocket, +) +from rocketpy.stochastic.stochastic_model import StochasticModel + +# Captured before any patching, so wrapping it twice in one test does not stack. +_REAL_SET_STOCHASTIC = StochasticModel._set_stochastic + + +def _seeds_handed_out(monkeypatch, rocket, seed): + """The seeds every nested component received during one reseed.""" + recorded = [] + + def recording(self, seed=None): + recorded.append(seed) + return _REAL_SET_STOCHASTIC(self, seed) + + monkeypatch.setattr(StochasticModel, "_set_stochastic", recording) + rocket._set_stochastic(seed) + return recorded + + +def _drawn(component): + return next(component.dict_generator()) + + +def _members_of(collection): + """Components yields (component, position) pairs; a plain list does not.""" + if isinstance(collection, Components): + return [component for component, _ in collection] + return list(collection) + + +def test_two_components_with_one_spec_do_not_share_one_stream( + stochastic_calisto, calisto_main_chute +): + """The whole rocket shared one seed, so two parachutes built from the same + spec drew the same ``cd_s`` and the same ``lag``, every time. A study of a + main and a drogue was really a study of one chute counted twice. + """ + stochastic_calisto.parachutes = [] + for _ in range(2): + stochastic_calisto.add_parachute( + StochasticParachute(parachute=calisto_main_chute, cd_s=0.1, lag=0.2) + ) + + stochastic_calisto._set_stochastic(99) + first, second = (_drawn(chute) for chute in stochastic_calisto.parachutes) + + assert first["cd_s"] != second["cd_s"] + assert first["lag"] != second["lag"] + + +def test_one_seed_reproduces_every_component(stochastic_calisto, calisto_main_chute): + """Independent is not enough on its own; it still has to follow the seed.""" + stochastic_calisto.parachutes = [] + stochastic_calisto.add_parachute( + StochasticParachute(parachute=calisto_main_chute, cd_s=0.1, lag=0.2) + ) + + def drawn_with(seed): + stochastic_calisto._set_stochastic(seed) + return [_drawn(chute) for chute in stochastic_calisto.parachutes] + [ + _drawn(stochastic_calisto)["mass"] + ] + + assert drawn_with(2718) == drawn_with(2718) + assert drawn_with(2718) != drawn_with(2719) + + +def test_component_seeds_do_not_collide(monkeypatch, stochastic_calisto): + """The same statement across every component type, not only the parachutes: + the body, each aerodynamic surface, the motor and the rail buttons. + """ + seeds = _seeds_handed_out(monkeypatch, stochastic_calisto, 42) + + assert len(seeds) > 3, "expected the rocket body and several components" + assert len(seeds) == len(set(seeds)), ( + "components share a seed, so they draw perfectly correlated samples" + ) + + +def test_the_reseed_covers_every_collection_create_object_uses(stochastic_calisto): + """Whatever ``create_object`` iterates has to be reseeded too. + + Read off the source rather than off a fixture, because the collection no + fixture populates is exactly the one that gets missed: air brakes were built + and sampled and never reseeded, and every seeding test passed anyway. + """ + rocket = stochastic_calisto + tree = ast.parse(inspect.getsource(type(rocket).create_object).lstrip()) + iterated = { + node.iter.attr + for node in ast.walk(tree) + # Comprehensions too. This scan exists to catch a collection added + # later, and one written as a comprehension would slip past a For walk. + if isinstance(node, (ast.For, ast.comprehension)) + and isinstance(node.iter, ast.Attribute) + and isinstance(node.iter.value, ast.Name) + and node.iter.value.id == "self" + and not node.iter.attr.startswith("_") + } + declared = set(type(rocket)._stochastic_collections()) + + assert iterated, "found no collections in create_object, so the scan is broken" + # Both directions. One left in the reseed after create_object stopped using + # it still spawns a child and moves every stream that follows. + assert iterated == declared, { + "sampled but never reseeded": sorted(iterated - declared), + "reseeded but never sampled": sorted(declared - iterated), + } + + +def test_an_air_brake_answers_to_the_seed( + stochastic_calisto, calisto_air_brakes_clamp_on +): + """Air brakes were in ``create_object`` and not in the reseed, so a fixed + seed did not reproduce them: 0.683, then 0.586, then 0.488 for one seed + asked three times. + + Built here rather than taken from the fixture, since wrapping an + ``AirBrakes`` with no arguments gives every parameter a zero standard + deviation, and that draws the same value under any seed whether or not + anything reseeds it. + """ + stochastic_calisto.add_air_brakes( + StochasticAirBrakes( + air_brakes=calisto_air_brakes_clamp_on.air_brakes[0], + deployment_level=(0.5, 0.1), + ), + calisto_air_brakes_clamp_on._controllers[0], + ) + air_brake = stochastic_calisto.air_brakes[0] + + def drawn(seed): + stochastic_calisto._set_stochastic(seed) + return _drawn(air_brake)["deployment_level"] + + first = drawn(1234) + + assert drawn(1234) == first, "the same seed drew a different air brake" + assert drawn(1235) != first, "a different seed drew the same air brake" + + +def test_adding_a_surface_leaves_the_other_collections_alone( + stochastic_calisto, calisto_main_chute, stochastic_nose_cone +): + """Each collection has a root of its own, so an unrelated component in one + of them does not move the streams in the others.""" + stochastic_calisto.parachutes = [] + stochastic_calisto.add_parachute( + StochasticParachute(parachute=calisto_main_chute, cd_s=0.1, lag=0.2) + ) + + def parachute_draw(): + stochastic_calisto._set_stochastic(5) + return _drawn(stochastic_calisto.parachutes[0]) + + before = parachute_draw() + # The deterministic nose, so add_nose builds a wrapper of its own. Adding + # the fixture again would store one wrapper twice, which is #1172. + stochastic_calisto.add_nose(stochastic_nose_cone.obj, position=1.1) + + assert parachute_draw() == before + + +def test_every_entry_is_reseeded_exactly_once( + monkeypatch, stochastic_calisto, calisto_main_chute, calisto_air_brakes_clamp_on +): + """Counted rather than read off the source. + + The scan above reads ``create_object`` for ``for x in self.collection``, so + a helper, a local alias or a ``getattr`` would hide a collection from it. + This counts what actually happens. + """ + stochastic_calisto.add_parachute( + StochasticParachute(parachute=calisto_main_chute, cd_s=0.1, lag=0.2) + ) + stochastic_calisto.add_air_brakes( + StochasticAirBrakes( + air_brakes=calisto_air_brakes_clamp_on.air_brakes[0], + deployment_level=(0.5, 0.1), + ), + calisto_air_brakes_clamp_on._controllers[0], + ) + + counted = {} + + def recording(self, seed=None): + counted[id(self)] = counted.get(id(self), 0) + 1 + return _REAL_SET_STOCHASTIC(self, seed) + + monkeypatch.setattr(StochasticModel, "_set_stochastic", recording) + stochastic_calisto._set_stochastic(3) + + entries = [ + component + for name in type(stochastic_calisto)._stochastic_collections() + for component in _members_of(getattr(stochastic_calisto, name)) + ] + + assert entries, "no components to count" + assert all(counted.get(id(entry)) == 1 for entry in entries), { + type(entry).__name__: counted.get(id(entry)) for entry in entries + } + + +def test_two_air_brakes_with_one_spec_stay_independent( + stochastic_calisto, calisto_air_brakes_clamp_on +): + """The air brakes are a plain list, so they take a different route through + the reseed than the positioned collections do. + + About the streams only. Both are added with one controller because the + rocket keeps a single one, which is a separate problem recorded in #1172. + """ + for _ in range(2): + stochastic_calisto.add_air_brakes( + StochasticAirBrakes( + air_brakes=calisto_air_brakes_clamp_on.air_brakes[0], + deployment_level=(0.5, 0.1), + ), + calisto_air_brakes_clamp_on._controllers[0], + ) + + def drawn(seed): + stochastic_calisto._set_stochastic(seed) + return [ + _drawn(brake)["deployment_level"] for brake in stochastic_calisto.air_brakes + ] + + first = drawn(808) + + assert first[0] != first[1], "two air brakes drew the same value" + assert drawn(808) == first + assert drawn(809) != first + + +def test_the_rocket_body_keeps_the_seed_as_given(monkeypatch, stochastic_calisto): + """Fixing the nested components did not need the body's stream to move. + + Reproducibility and seed uniqueness both hold with the body on a spawned + child, so neither of them would notice it going back there and taking every + fixed-seed mass and radius baseline with it. + """ + seeds = _seeds_handed_out(monkeypatch, stochastic_calisto, 42) + + assert seeds[0] == 42 + + +def test_the_tree_is_built_by_the_reset_not_by_the_add(calisto, calisto_main_chute): + """A rocket resets itself while being constructed, with nothing attached. + + So a component added afterwards keeps the generator it was built with until + the next reset, which is the guarantee the documentation states and the one + a Monte Carlo relies on. + """ + rocket = StochasticRocket(rocket=calisto, radius=0.0127 / 2) + chute = StochasticParachute(parachute=calisto_main_chute, cd_s=0.1, lag=0.2) + rocket.add_parachute(chute) + + assert getattr(chute, "_seed", None) is None + + rocket._set_stochastic(99) + + assert chute._seed is not None + + +def test_a_worker_seed_sequence_is_copied_rather_than_spawned_from(): + # A parallel run hands each worker a SeedSequence. Spawning from it would + # advance a counter the caller still holds, so the next use of the same + # object would build a different tree. + worker = np.random.SeedSequence(7).spawn(2)[0] + + _seed_sequence_from(worker).spawn(3) + + assert worker.n_children_spawned == 0 + + +def test_a_worker_seed_sequence_keeps_its_place_in_the_tree(): + worker = np.random.SeedSequence(7).spawn(2)[0] + + children = _seed_sequence_from(worker).spawn(2) + + assert [child.spawn_key for child in children] == [(0, 0), (0, 1)] + + +def test_two_workers_do_not_get_the_same_collection_roots(): + first, second = np.random.SeedSequence(7).spawn(2) + + one = _seed_sequence_from(first).spawn(1)[0].generate_state(4) + other = _seed_sequence_from(second).spawn(1)[0].generate_state(4) + + assert not np.array_equal(one, other) + + +def test_an_integer_seed_still_roots_the_collections(): + assert _seed_sequence_from(42).spawn(1)[0].spawn_key == (0,) + assert _seed_sequence_from(None).spawn(1)[0].spawn_key == (0,) diff --git a/tests/unit/test_tools.py b/tests/unit/test_tools.py index 3b8df37a3..59f04c6eb 100644 --- a/tests/unit/test_tools.py +++ b/tests/unit/test_tools.py @@ -5,6 +5,7 @@ from rocketpy import Environment from rocketpy.tools import ( + _seed_sequence_to_int, calculate_confidence_ellipse, calculate_cubic_hermite_coefficients, convert_local_extent_to_wgs84, @@ -347,3 +348,25 @@ def test_mercator_extent_to_local_preserves_offset_sign( assert local_extent[0] < local_extent[1] assert local_extent[2] < local_extent[3] assert all(expected_sign * value > 0 for value in local_extent) + + +def test_seed_sequence_to_int_keeps_the_full_width(): + """All four words have to reach the seed. + + Taking only the first one would still hand every component a different + number, so every seeding test would pass over a 32-bit collapse that puts + two streams back together near 2**16 of them. + """ + root = np.random.SeedSequence(12345) + a, b = root.spawn(2) + words = a.generate_state(4, dtype=np.uint32) + # Every word, in one comparison. Asserting only that the high bits are not + # zero leaves a 64 or 96 bit truncation passing. + expected = sum(int(word) << (32 * at) for at, word in enumerate(words)) + + seed = _seed_sequence_to_int(a) + + assert seed == expected + assert seed.bit_length() <= 128 + assert seed != _seed_sequence_to_int(b) + assert seed == _seed_sequence_to_int(a), "reading it twice moved the seed"