From ff11f4591f773cbae1873acf8a8f39509c4e9ad7 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Sun, 28 Jun 2026 22:39:33 +0800 Subject: [PATCH 1/4] ENH: seed Monte Carlo per simulation index for worker-invariant inputs MonteCarlo seeded the stochastic models once per worker (parallel) or once before the loop (serial), so the generated inputs depended on how many workers ran and which worker drew which simulation. Two runs of the same seed with different worker counts produced different inputs. Seed per simulation index instead: spawn one SeedSequence child per index from the run's root seed, and reseed the stochastic models from child_seeds[i] before simulation i. SeedSequence.spawn is invariant to the spawn count, so index i always maps to the same child, including in append mode. Each index seed is split three ways so the environment, rocket and flight do not share a stream. Inputs are now identical across serial, parallel(2) and parallel(N) for a fixed seed, and a split append run reproduces a single run. This changes the numbers a fixed seed produces (the per-model decorrelation, plus a serial index that now counts from 0 like the parallel path already did), so stored baselines regenerate. Adds tests/unit/simulation/test_montecarlo_determinism.py covering serial reproducibility, worker-invariance, append reproducibility, and the None-seed path. Marked slow (each simulation rebuilds a full rocket) per the existing test_monte_carlo_simulate convention. --- rocketpy/simulation/monte_carlo.py | 77 +++- .../simulation/test_montecarlo_determinism.py | 391 ++++++++++++++++++ 2 files changed, 450 insertions(+), 18 deletions(-) create mode 100644 tests/unit/simulation/test_montecarlo_determinism.py diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index b6b71e0c3..ae45f85d9 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -291,16 +291,23 @@ def __run_in_serial(self, random_seed=None): start_time=time(), ) inputs_json = "" + sim_idx = self._initial_sim_idx + # One independent seed per simulation index makes the generated inputs + # invariant to the execution mode and to the worker count (see + # ``__seed_components``). The ``SeedSequence`` is created fresh here so + # that index ``i`` always maps to the same child seed, including in + # ``append`` mode (``spawn`` is invariant to the spawn count). + child_seeds = np.random.SeedSequence(random_seed).spawn( + self.number_of_simulations + ) try: - self.environment._set_stochastic(random_seed) - self.rocket._set_stochastic(random_seed) - self.flight._set_stochastic(random_seed) while sim_monitor.keep_simulating(): - sim_monitor.increment() + sim_idx = sim_monitor.increment() - 1 + self.__seed_components(child_seeds[sim_idx]) flight = self.__run_single_simulation() - inputs_json = self.__evaluate_flight_inputs(sim_monitor.count) - outputs_json = self.__evaluate_flight_outputs(flight, sim_monitor.count) + inputs_json = self.__evaluate_flight_inputs(sim_idx) + outputs_json = self.__evaluate_flight_outputs(flight, sim_idx) self.__append_serial_results(inputs_json, outputs_json) sim_monitor.print_update_status() @@ -312,7 +319,7 @@ def __run_in_serial(self, random_seed=None): except Exception as error: self.__save_serial_error( - inputs_json, f"Error on iteration {sim_monitor.count}: {error}" + inputs_json, f"Error on iteration {sim_idx}: {error}" ) raise error @@ -359,13 +366,20 @@ def __run_in_parallel(self, random_seed=None, n_workers=None): ) processes = [] - seeds = np.random.SeedSequence(random_seed).spawn(n_workers) + # One independent seed per simulation index (not per worker) makes + # the generated inputs invariant to ``n_workers``: the worker that + # runs simulation ``i`` always seeds from ``child_seeds[i]``, + # regardless of how many workers there are. The full list is shared + # with every worker; the shared atomic counter assigns indices. + child_seeds = np.random.SeedSequence(random_seed).spawn( + self.number_of_simulations + ) - for seed in seeds: + for _ in range(n_workers): sim_producer = multiprocess.Process( target=self.__sim_producer, args=( - seed, + child_seeds, sim_monitor, mutex, simulation_error_event, @@ -407,13 +421,16 @@ def __validate_number_of_workers(self, n_workers): raise ValueError("Number of workers must be at least 2 for parallel mode.") return n_workers - def __sim_producer(self, seed, sim_monitor, mutex, error_event): # pylint: disable=too-many-statements + def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event): # pylint: disable=too-many-statements """Simulation producer to be used in parallel by multiprocessing. Parameters ---------- - seed : int - The seed to set the random number generator. + child_seeds : list[numpy.random.SeedSequence] + One seed sequence per simulation index. Before each simulation the + worker seeds the stochastic models from ``child_seeds[sim_idx]``, + where ``sim_idx`` is pulled from the shared atomic counter. This + keeps the generated inputs invariant to the number of workers. sim_monitor : _SimMonitor The simulation monitor object to keep track of the simulations. mutex : multiprocess.Lock @@ -422,15 +439,14 @@ def __sim_producer(self, seed, sim_monitor, mutex, error_event): # pylint: disa Event signaling an error occurred during the simulation. """ try: - # Ensure Processes generate different random numbers - self.environment._set_stochastic(seed) - self.rocket._set_stochastic(seed) - self.flight._set_stochastic(seed) - while sim_monitor.keep_simulating(): sim_idx = sim_monitor.increment() - 1 inputs_json, outputs_json = "", "" + # Seed per simulation index so the inputs are reproducible and + # independent of which worker happens to run this index. + self.__seed_components(child_seeds[sim_idx]) + flight = self.__run_single_simulation() inputs_json = self.__evaluate_flight_inputs(sim_idx) outputs_json = self.__evaluate_flight_outputs(flight, sim_idx) @@ -466,6 +482,31 @@ def __sim_producer(self, seed, sim_monitor, mutex, error_event): # pylint: disa error_event.set() mutex.release() + def __seed_components(self, simulation_seed): + """Seed the stochastic models for a single simulation index. + + The given seed sequence is split into three independent sub-streams so + that the environment, rocket and flight do not share the same random + draws (sharing a single seed would correlate their first sampled + values). Seeding per simulation index -- rather than once per worker -- + is what makes the Monte Carlo inputs invariant to the execution mode + (serial vs parallel) and to the number of workers. + + Parameters + ---------- + simulation_seed : numpy.random.SeedSequence + The seed sequence assigned to the current simulation index. It is + spawned into three child sequences, one per stochastic model. + + Returns + ------- + None + """ + env_seed, rocket_seed, flight_seed = simulation_seed.spawn(3) + self.environment._set_stochastic(env_seed) + self.rocket._set_stochastic(rocket_seed) + self.flight._set_stochastic(flight_seed) + def __run_single_simulation(self): """Runs a single simulation and returns the inputs and outputs. diff --git a/tests/unit/simulation/test_montecarlo_determinism.py b/tests/unit/simulation/test_montecarlo_determinism.py new file mode 100644 index 000000000..2c4629db4 --- /dev/null +++ b/tests/unit/simulation/test_montecarlo_determinism.py @@ -0,0 +1,391 @@ +"""Determinism tests for :class:`rocketpy.simulation.MonteCarlo`. + +These tests verify that the random *inputs* generated by a Monte Carlo run are +invariant to the execution mode (serial vs parallel) and to the number of +parallel workers when a fixed ``random_seed`` is provided. This is the +worker-invariance guarantee introduced by per-simulation-index seeding: the +seed of simulation ``i`` is derived from ``i`` only, never from the worker that +happens to run it. + +The trajectory integration (:class:`rocketpy.Flight`) is replaced by a +lightweight stub so the tests run fast. Worker-invariance is a property of the +*input sampling*, which happens in ``create_object``/``_randomize_*`` before the +``Flight`` object is built and is independent of the (expensive) physics +integration. Stubbing the module-level ``Flight`` symbol propagates to the +parallel workers because the parallel backend uses the ``fork`` start method, so +forked workers inherit the patched module. + +A dedicated ``stochastic_calisto_numpy_only`` rocket is used so that *all* +randomness flows through the seeded numpy generator. List-valued stochastic +attributes are sampled with the standard-library ``random.choice`` (an unseeded +*global* generator) which the per-index seeding does not control; the fixture +removes the only such attribute (a multi-element ``thrust_source``) so the +generated inputs are byte-for-byte reproducible from the seed alone. + +These tests are additive and isolated; they do not modify any inherited test. +""" + +import json + +import multiprocess +import pytest + +import rocketpy.simulation.monte_carlo as mc_module +from rocketpy.simulation import MonteCarlo +from rocketpy.stochastic import StochasticRocket, StochasticSolidMotor + +# Every test here drives a real ``MonteCarlo.simulate`` run (the ``Flight`` +# physics is stubbed, but a full Calisto rocket is still rebuilt per simulation), +# so they are gated behind ``--runslow`` to keep the default suite fast, matching +# the existing ``test_monte_carlo_simulate`` convention. +pytestmark = pytest.mark.slow + +# The stub-based parallel tests rely on workers inheriting the monkeypatched +# ``Flight`` symbol, which only happens with the ``fork`` start method. +requires_fork = pytest.mark.skipif( + multiprocess.get_start_method() != "fork", + reason="stub-based parallel determinism test requires the 'fork' start method", +) + + +class _StubFlight: + """Minimal stand-in for :class:`rocketpy.Flight` that skips integration. + + The Monte Carlo input sampling happens in + ``StochasticRocket.create_object``, ``StochasticEnvironment.create_object`` + and the ``StochasticFlight`` randomize helpers *before* the ``Flight`` + object is built. Replacing ``Flight`` with this stub therefore preserves the + exact random draws while skipping the costly trajectory integration. Any + output attribute requested by the exporter resolves to ``0.0``. + """ + + def __init__(self, **kwargs): + # Accepts the keyword arguments that ``MonteCarlo`` passes to ``Flight`` + # and intentionally ignores them. + pass + + def __getattr__(self, name): + return 0.0 + + +@pytest.fixture +def stochastic_calisto_numpy_only( + cesaroni_m1670, + calisto_robust, + stochastic_nose_cone, + stochastic_trapezoidal_fins, + stochastic_tail, + stochastic_rail_buttons, + stochastic_main_parachute, + stochastic_drogue_parachute, +): + """A ``StochasticRocket`` whose randomness flows entirely through numpy. + + It mirrors the shared ``stochastic_calisto`` fixture but gives the solid + motor a single fixed ``thrust_source`` instead of a multi-element list. + List-valued stochastic attributes are drawn with the standard-library + ``random.choice`` (an unseeded *global* generator) which per-simulation-index + seeding does not govern. Removing the only such attribute makes every + generated input reproducible from the seed alone, so serial and parallel + runs can be compared byte-for-byte. + + Returns + ------- + StochasticRocket + A stochastic Calisto rocket with no global-``random`` dependence. + """ + motor = StochasticSolidMotor( + solid_motor=cesaroni_m1670, + burn_out_time=(4, 0.1), + grains_center_of_mass_position=0.001, + grain_density=50, + grain_separation=1 / 1000, + grain_initial_height=1 / 1000, + grain_initial_inner_radius=0.375 / 1000, + grain_outer_radius=0.375 / 1000, + total_impulse=(6500, 1000), + throat_radius=0.5 / 1000, + nozzle_radius=0.5 / 1000, + nozzle_position=0.001, + ) + rocket = StochasticRocket( + rocket=calisto_robust, + radius=0.0127 / 2000, + mass=(15.426, 0.5, "normal"), + inertia_11=(6.321, 0), + inertia_22=0.01, + inertia_33=0.01, + center_of_mass_without_motor=0, + ) + rocket.add_motor(motor, position=0.001) + rocket.add_nose(stochastic_nose_cone, position=(1.134, 0.001)) + rocket.add_trapezoidal_fins(stochastic_trapezoidal_fins, position=(0.001, "normal")) + rocket.add_tail(stochastic_tail) + rocket.set_rail_buttons( + stochastic_rail_buttons, lower_button_position=(-0.618, 0.001, "normal") + ) + rocket.add_parachute(stochastic_main_parachute) + rocket.add_parachute(stochastic_drogue_parachute) + return rocket + + +def _read_inputs_by_index(input_file): + """Read a Monte Carlo ``.inputs.txt`` file into ``{index: raw_json_line}``. + + Keying by the ``index`` field makes the comparison robust to the order in + which parallel workers append their results. + + Parameters + ---------- + input_file : str or pathlib.Path + Path to the Monte Carlo inputs file. + + Returns + ------- + dict[int, str] + Mapping of simulation index to the raw (stripped) JSON line written for + that simulation. + """ + by_index = {} + with open(input_file, mode="r", encoding="utf-8") as rows: + for line in rows: + line = line.strip() + if not line: + continue + record = json.loads(line) + by_index[record["index"]] = line + return by_index + + +def _simulate_inputs( + monkeypatch, + tmp_path, + environment, + rocket, + flight, + tag, + *, + n_simulations, + random_seed, + parallel, + n_workers=None, +): + """Run a Monte Carlo simulation with a stubbed ``Flight`` and return inputs. + + Parameters + ---------- + monkeypatch : pytest.MonkeyPatch + Fixture used to swap ``Flight`` for :class:`_StubFlight`. + tmp_path : pathlib.Path + Temporary directory for the Monte Carlo output files. + environment, rocket, flight : StochasticModel + Stochastic models shared across runs. + tag : str + Unique filename stem so successive runs do not clobber each other. + n_simulations : int + Number of simulations to run. + random_seed : int or None + Seed forwarded to ``MonteCarlo.simulate``. + parallel : bool + Whether to run in parallel mode. + n_workers : int, optional + Number of workers for parallel mode. + + Returns + ------- + dict[int, str] + Mapping of simulation index to the raw input JSON line. + """ + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + montecarlo = MonteCarlo( + filename=str(tmp_path / tag), + environment=environment, + rocket=rocket, + flight=flight, + ) + montecarlo.simulate( + number_of_simulations=n_simulations, + append=False, + parallel=parallel, + random_seed=random_seed, + n_workers=n_workers, + ) + return _read_inputs_by_index(montecarlo.input_file) + + +def test_montecarlo_inputs_serial_reproducible( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """Two serial runs with the same seed yield identical inputs per index.""" + models = ( + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, + ) + run_a = _simulate_inputs( + monkeypatch, + tmp_path, + *models, + "serial_a", + n_simulations=6, + random_seed=20240617, + parallel=False, + ) + run_b = _simulate_inputs( + monkeypatch, + tmp_path, + *models, + "serial_b", + n_simulations=6, + random_seed=20240617, + parallel=False, + ) + + assert sorted(run_a) == list(range(6)) + assert run_a == run_b + + +@requires_fork +def test_montecarlo_inputs_worker_invariant( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """Serial == parallel(2) == parallel(N): inputs are bit-identical per index. + + The per-iteration input JSON for a given simulation index must be byte-for- + byte identical regardless of whether the run was serial, parallel with two + workers, or parallel with a larger number of workers. + """ + n_simulations = 8 + random_seed = 314159 + models = ( + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, + ) + + serial = _simulate_inputs( + monkeypatch, + tmp_path, + *models, + "serial", + n_simulations=n_simulations, + random_seed=random_seed, + parallel=False, + ) + parallel_2 = _simulate_inputs( + monkeypatch, + tmp_path, + *models, + "parallel_2", + n_simulations=n_simulations, + random_seed=random_seed, + parallel=True, + n_workers=2, + ) + parallel_n = _simulate_inputs( + monkeypatch, + tmp_path, + *models, + "parallel_n", + n_simulations=n_simulations, + random_seed=random_seed, + parallel=True, + n_workers=4, + ) + + # Every index must be present exactly once in each run. + expected_indices = list(range(n_simulations)) + assert sorted(serial) == expected_indices + assert sorted(parallel_2) == expected_indices + assert sorted(parallel_n) == expected_indices + + # Worker-invariance: bit-identical input line per index across all modes. + for index in expected_indices: + assert serial[index] == parallel_2[index], ( + f"serial vs parallel(2) inputs differ at index {index}" + ) + assert serial[index] == parallel_n[index], ( + f"serial vs parallel(4) inputs differ at index {index}" + ) + + +def test_montecarlo_none_seed_runs( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """``random_seed=None`` stays functional (non-reproducible but complete). + + ``numpy.random.SeedSequence(None)`` pulls fresh entropy, so the run cannot + be reproduced, but it must still execute and export one record per index. + """ + inputs = _simulate_inputs( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, + "none_seed", + n_simulations=5, + random_seed=None, + parallel=False, + ) + + assert sorted(inputs) == list(range(5)) + + +def test_montecarlo_inputs_append_reproducible( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """Appending to a partial run reproduces a single full run, index for index. + + ``SeedSequence.spawn`` is invariant to the spawn count, so simulation ``i`` + draws the same seed whether it was produced in the first batch or in a later + ``append`` batch. A 3 + 3 (append) run must therefore match a single run of 6 + with the same seed. ``number_of_simulations`` is the cumulative target, so the + append call passes 6, not the 3 additional simulations. + """ + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + seed = 2718281 + environment = stochastic_environment + rocket = stochastic_calisto_numpy_only + flight = stochastic_flight + + split = MonteCarlo( + filename=str(tmp_path / "split"), + environment=environment, + rocket=rocket, + flight=flight, + ) + split.simulate(number_of_simulations=3, append=False, random_seed=seed) + split.simulate(number_of_simulations=6, append=True, random_seed=seed) + split_inputs = _read_inputs_by_index(split.input_file) + + single_inputs = _simulate_inputs( + monkeypatch, + tmp_path, + environment, + rocket, + flight, + "single", + n_simulations=6, + random_seed=seed, + parallel=False, + ) + + assert sorted(split_inputs) == list(range(6)) + assert split_inputs == single_inputs From 03bc924de173474e78edea10abee73424bc1883f Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Thu, 6 Aug 2026 17:17:30 +0800 Subject: [PATCH 2/4] MNT: format the example in the Copilot instructions `ruff format --check .` is one of the Linters jobs and it covers this file, so the unformatted example in it fails that job on develop today. Splitting it out of the seeding change because it has nothing to do with it. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- .github/copilot-instructions.md | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md index f5366cb3b..d735ada89 100644 --- a/.github/copilot-instructions.md +++ b/.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 From 46c814771f1ee822eb27a8c8b00ec2f3da62bf53 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Thu, 6 Aug 2026 17:17:30 +0800 Subject: [PATCH 3/4] ENH: take the per-index Monte Carlo seeding from RocketPy#1054 Replaces the earlier draft of this change with the implementation that is now on RocketPy#1054, so the two carry the same code and a later sync from upstream is a no-op rather than a conflict to resolve. Simulation index i draws its seed from its own child of the run's root, derived before that simulation runs, so index i maps to the same inputs whether the run is serial or parallel and however many workers it uses. The child is built by extending the root's spawn_key, which is how SeedSequence.spawn derives it, so no per-index list is materialised or sent to a worker. Beyond the seeding itself, the upstream version carries the review rounds this draft predates: each stochastic model holds its nominal values steady instead of reading them back off an object that create_object mutates; air brakes are in the reseed and a source scan keeps a later collection from missing it; CP and thrust eccentricity are reapplied after each reseed; the parallel index claim holds the mutex across the check and the increment; a worker preserves the failure that killed it rather than raising over it; the parent bounds its wait and checks exit codes; and the completeness check rejects a log it cannot read rather than skipping past it. Four things stay ours, since they are what the fork is for: simulate returns the results dict, each Flight gets the stochastic flight's max_time, max_time_step and min_time_step, set_processed_results stays disabled with its TODO, and StochasticRocket keeps add_linear_generic_surface and volume. The rocket file was merged three ways against the shared RocketPy base rather than overwritten. Local run: ruff clean, pylint rocketpy/ tests/ docs/ exits 0, 1640 unit tests and 152 integration tests pass, and the start-method gate passes under fork, spawn and forkserver. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- rocketpy/simulation/monte_carlo.py | 1071 ++++++++++++++--- rocketpy/stochastic/stochastic_model.py | 64 +- rocketpy/stochastic/stochastic_rocket.py | 93 +- rocketpy/tools.py | 50 +- .../test_monte_carlo_determinism.py | 787 ++++++++++++ .../test_monte_carlo_determinism.py | 328 +++++ .../test_monte_carlo_log_integrity.py | 401 ++++++ .../test_monte_carlo_worker_exit.py | 179 +++ .../test_monte_carlo_worker_failures.py | 256 ++++ .../simulation/test_montecarlo_determinism.py | 391 ------ .../unit/stochastic/test_stochastic_model.py | 114 +- .../test_stochastic_rocket_seeding.py | 171 +++ 12 files changed, 3289 insertions(+), 616 deletions(-) create mode 100644 tests/integration/simulation/test_monte_carlo_determinism.py create mode 100644 tests/unit/simulation/test_monte_carlo_determinism.py create mode 100644 tests/unit/simulation/test_monte_carlo_log_integrity.py create mode 100644 tests/unit/simulation/test_monte_carlo_worker_exit.py create mode 100644 tests/unit/simulation/test_monte_carlo_worker_failures.py delete mode 100644 tests/unit/simulation/test_montecarlo_determinism.py create mode 100644 tests/unit/stochastic/test_stochastic_rocket_seeding.py diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index ae45f85d9..cf28eeab2 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -13,12 +13,13 @@ latest documentation. """ +import csv import json import os import traceback import warnings from pathlib import Path -from time import time +from time import monotonic, time import numpy as np import simplekml @@ -29,6 +30,7 @@ from rocketpy.prints.monte_carlo_prints import _MonteCarloPrints from rocketpy.simulation.flight import Flight from rocketpy.tools import ( + _seed_sequence_to_int, generate_monte_carlo_ellipses, generate_monte_carlo_ellipses_coordinates, import_optional_dependency, @@ -37,7 +39,7 @@ # TODO: Create evolution plots to analyze convergence -class MonteCarlo: +class MonteCarlo: # pylint: disable=too-many-public-methods """Class to run a Monte Carlo simulation of a rocket flight. Attributes @@ -95,7 +97,7 @@ def __init__( flight, export_list=None, data_collector=None, - ): # pylint: disable=too-many-statements + ): """ Initialize a MonteCarlo object. @@ -168,10 +170,11 @@ def simulate( number_of_simulations, append=False, parallel=False, - random_seed=None, n_workers=None, + *, + random_seed=None, **kwargs, - ): # pylint: disable=too-many-statements + ): """ Runs the Monte Carlo simulation and saves all data. @@ -184,13 +187,27 @@ def simulate( False, the files will be overwritten. Default is False. parallel : bool, optional If True, the simulations will be run in parallel. Default is False. - random_seed : int, optional - The seed to set the random number generator. Default is None. n_workers : int, optional Number of workers to be used if ``parallel=True``. If None, the number of workers will be equal to the number of CPUs available. A minimum of 2 workers is required for parallel mode. Default is None. + random_seed : int, numpy integer, sequence of ints, or SeedSequence, optional + Root seed for the run. When provided, the sampled inputs are + reproducible and identical across serial and parallel execution and + across any number of workers: each simulation index derives its own + decorrelated child stream from this root, so index ``i`` receives the + same inputs no matter which worker runs it. A supplied ``SeedSequence`` + is copied from its full state rather than consumed, so repeated calls + with the same seed reproduce the same inputs. Each model is reseeded + with a 128-bit integer -- the seed type a custom sampler's + ``reset_seed`` accepts. A stateful ``numpy.random.Generator`` or + ``BitGenerator`` is rejected (it is an RNG to draw from, not a fixed + seed); pass ``rng.bit_generator.seed_seq`` to seed from one. Default is + None, which draws fresh entropy on each run -- the previous, + non-reproducible default. This seeding is informed by Scientific Python + SPEC 7 but keeps immutable seed-snapshot semantics rather than sharing a + ``Generator``. kwargs : dict Custom arguments for simulation export of the ``inputs`` file. Options are: @@ -222,19 +239,38 @@ def simulate( overwritten. Make sure to save the files with the results before running the simulation again with `append=False`. """ + # Everything that can be judged from the arguments alone happens before + # __setup_files, which opens both logs "w+" and empties them. Raising + # after that point destroys the previous run on the way out. + _validate_simulation_count(number_of_simulations) + if parallel: + n_workers = self.__validate_number_of_workers(n_workers) + self._export_config = kwargs self.number_of_simulations = number_of_simulations self._initial_sim_idx = self.num_of_loaded_sims if append else 0 + # Both run paths catch Ctrl-C, save what they have and return, so a + # stopped run is incomplete on purpose and the completeness check below + # has to know the difference between that and a worker going missing. + self._interrupted = False - _SimMonitor.reprint("Starting Monte Carlo analysis") + # Capture the small, picklable root seed state once per run (every + # simulation index derives its child seed from it, see __child_seed). + # This validates random_seed *before* __setup_files truncates any + # existing output, so an invalid seed cannot destroy prior results on + # the way to raising. + self.__capture_root_state(random_seed) + + print("Starting Monte Carlo analysis") self.__setup_files(append) if parallel: - self.__run_in_parallel(random_seed, n_workers) + self.__run_in_parallel(n_workers) else: - self.__run_in_serial(random_seed) + self.__run_in_serial() + self.__check_each_index_was_recorded_once() self.__terminate_simulation() return self.results @@ -272,14 +308,141 @@ def __setup_files(self, append): except OSError as error: raise OSError(f"Error creating files: {error}") from error - def __run_in_serial(self, random_seed=None): + @staticmethod + def __root_seed_sequence(random_seed): + """Build a fresh ``SeedSequence`` root from ``random_seed``. + + ``random_seed`` may be an int (or any entropy ``numpy.random.SeedSequence`` + accepts), an existing ``SeedSequence``, or ``None`` for fresh entropy. A + supplied ``SeedSequence`` is copied from its full ``state``, so the + spawning below neither mutates the caller's object nor advances a shared + child counter between calls; repeated ``simulate`` calls with the same + seed then stay reproducible. A stateful ``Generator``/``BitGenerator`` is + not accepted, since using it as an immutable seed would contradict its + consume-on-use semantics; pass ``rng.bit_generator.seed_seq`` to seed + from an existing generator's stream. + """ + if isinstance(random_seed, np.random.SeedSequence): + return np.random.SeedSequence(**random_seed.state) + if isinstance(random_seed, (np.random.Generator, np.random.BitGenerator)): + raise TypeError( + "random_seed must be an int or a numpy.random.SeedSequence, not " + f"a {type(random_seed).__name__}; to seed from an existing " + "generator pass rng.bit_generator.seed_seq." + ) + return np.random.SeedSequence(random_seed) + + def __capture_root_state(self, random_seed): + """Capture the small, picklable root seed state for this run. + + Stored once so serial mode and every parallel worker derive the same + per-index child seeds from it (see ``__child_seed``), instead of + materializing and pickling the full ``spawn(number_of_simulations)`` + list to each process. + """ + root = self.__root_seed_sequence(random_seed) + self.__root_state = ( + root.entropy, + root.spawn_key, + root.pool_size, + root.n_children_spawned, + ) + + def __child_seed(self, sim_idx): + """Return the seed sequence for a single simulation index. + + This equals ``root.spawn(number_of_simulations)[sim_idx]`` but is O(1) + in time and memory: ``SeedSequence.spawn`` derives child ``i`` by + appending ``n_children_spawned + i`` to the parent ``spawn_key``, so + rebuilding that one child directly reproduces it bit-for-bit while + letting a worker reconstruct any index from the small root state alone. + """ + entropy, spawn_key, pool_size, base = self.__root_state + return np.random.SeedSequence( + entropy=entropy, + spawn_key=(*spawn_key, base + sim_idx), + pool_size=pool_size, + ) + + def __seed_simulation(self, child_seed): + """Reseed the stochastic models for a single simulation index. + + The per-index child seed is split three ways so the environment, + rocket and flight draw from independent streams instead of sharing + one. Seeding per simulation index (not per worker) is what makes the + sampled inputs invariant to the execution mode and to the number of + workers. Each sub-stream is handed over as a 128-bit ``int`` (see + ``_seed_sequence_to_int``) so custom samplers keep working. + """ + env_seed, rocket_seed, flight_seed = child_seed.spawn(3) + self.environment._set_stochastic(_seed_sequence_to_int(env_seed)) + self.rocket._set_stochastic(_seed_sequence_to_int(rocket_seed)) + self.flight._set_stochastic(_seed_sequence_to_int(flight_seed)) + + def __check_each_index_was_recorded_once(self): + """Every index this run claimed left exactly one input and one output row. + + The counter hands each index out once, so a missing one means a worker + stopped between claiming and writing, and a repeated one means two + claimed the same index. Neither is visible in the files themselves: the + rows look well formed, and reading them back keyed by index hides the + duplicate behind the row that overwrote it. Both make the results wrong + while the run reports success, which is the thing per-index seeding is + supposed to rule out. + + Only over the range this run produced. ``append=True`` leaves earlier + runs in the same files, and ``number_of_simulations`` is the total to + reach rather than a count to add, so the new indices are + ``_initial_sim_idx`` up to it. + + A run stopped with Ctrl-C is short on purpose, so the indices it never + reached are not an error. What it did write is still held to the rest: + readable rows, one row per index, and nothing outside the range. + + Indices below ``_initial_sim_idx`` are an earlier run's and are left + alone. Reconciling a history with holes in it is a separate job, and + this only answers for the simulations this run claimed. + """ + inputs = _recorded_indices("inputs", self.input_file) + outputs = _recorded_indices("outputs", self.output_file) + if inputs != outputs: + only_in = lambda a, b: sorted(set(a) - set(b)) # noqa: E731 + raise RuntimeError( + f"the input and output files disagree about which simulations " + f"ran: {only_in(inputs, outputs)[:5]} have inputs and no " + f"outputs, {only_in(outputs, inputs)[:5]} the other way round. " + f"A worker stopped between the two writes, so the results are " + f"not reported as a successful run." + ) + + repeated = sorted(index for index, count in inputs.items() if count > 1) + beyond = sorted( + index for index in inputs if index >= self.number_of_simulations + ) + missing = ( + [] + if self._interrupted + else sorted( + set(range(self._initial_sim_idx, self.number_of_simulations)) + - set(inputs) + ) + ) + if missing or repeated or beyond: + raise RuntimeError( + f"the files do not match the simulations that ran: " + f"{len(missing)} never written {missing[:5]}, " + f"{len(repeated)} written more than once {repeated[:5]}, " + f"{len(beyond)} outside the range this run claimed {beyond[:5]}. " + f"The results are wrong, so they are not reported as a " + f"successful run." + ) + + def __run_in_serial(self): """ Runs the monte carlo simulation in serial mode. - Parameters - ---------- - random_seed : int, optional - The seed to set the random number generator in serial mode. Default is None. + The root seed state is captured by ``simulate`` before this runs, so each + simulation index derives its child seed from ``self.__root_state``. Returns ------- @@ -290,58 +453,54 @@ def __run_in_serial(self, random_seed=None): n_simulations=self.number_of_simulations, start_time=time(), ) - inputs_json = "" - sim_idx = self._initial_sim_idx - # One independent seed per simulation index makes the generated inputs - # invariant to the execution mode and to the worker count (see - # ``__seed_components``). The ``SeedSequence`` is created fresh here so - # that index ``i`` always maps to the same child seed, including in - # ``append`` mode (``spawn`` is invariant to the spawn count). - child_seeds = np.random.SeedSequence(random_seed).spawn( - self.number_of_simulations - ) try: - while sim_monitor.keep_simulating(): + while True: + # First statement in the loop, so it is bound before the two + # monitor calls rather than after them. Ctrl-C in either one + # used to leave it unbound, or holding the last completed row. + inputs_json = "" + + if not sim_monitor.keep_simulating(): + break sim_idx = sim_monitor.increment() - 1 - self.__seed_components(child_seeds[sim_idx]) + self.__seed_simulation(self.__child_seed(sim_idx)) flight = self.__run_single_simulation() inputs_json = self.__evaluate_flight_inputs(sim_idx) outputs_json = self.__evaluate_flight_outputs(flight, sim_idx) - self.__append_serial_results(inputs_json, outputs_json) + _record_simulation( + self.input_file, self.output_file, inputs_json, outputs_json + ) sim_monitor.print_update_status() sim_monitor.print_final_status() except KeyboardInterrupt: - self.__save_serial_error(inputs_json, "Keyboard Interrupt, files saved.") + self._interrupted = True + print("Keyboard interrupt received. Files saved.") + self.__keep_the_inputs_that_did_not_finish(inputs_json) except Exception as error: - self.__save_serial_error( - inputs_json, f"Error on iteration {sim_idx}: {error}" - ) + print(f"Error on iteration {sim_monitor.count}: {error}") + self.__keep_the_inputs_that_did_not_finish(inputs_json) raise error - def __append_serial_results(self, inputs_json, outputs_json): - with open(self.input_file, "a", encoding="utf-8") as f: - f.write(inputs_json) - with open(self.output_file, "a", encoding="utf-8") as f: - f.write(outputs_json) - - def __save_serial_error(self, inputs_json, message): - _SimMonitor.reprint(message) + def __keep_the_inputs_that_did_not_finish(self, inputs_json): + """Append the inputs of a simulation that stopped part way through.""" with open(self._error_file, "a", encoding="utf-8") as f: f.write(inputs_json) - def __run_in_parallel(self, random_seed=None, n_workers=None): + def __run_in_parallel(self, n_workers=None): """ Runs the monte carlo simulation in parallel. + The root seed state is captured by ``simulate`` before this runs and + travels with the pickled instance, so every worker derives the same + per-index child seed from ``self.__root_state``. + Parameters ---------- - random_seed : int, optional - The seed to set the random sequence generator in parallel mode. Default is None. n_workers: int, optional Number of workers to be used. If None, the number of workers will be equal to the number of CPUs available. Default is None. @@ -352,7 +511,7 @@ def __run_in_parallel(self, random_seed=None, n_workers=None): """ n_workers = self.__validate_number_of_workers(n_workers) - _SimMonitor.reprint(f"Running Monte Carlo simulation with {n_workers} workers.") + print(f"Running Monte Carlo simulation with {n_workers} workers.") multiprocess, managers = _import_multiprocess() @@ -365,72 +524,67 @@ def __run_in_parallel(self, random_seed=None, n_workers=None): start_time=time(), ) - processes = [] - # One independent seed per simulation index (not per worker) makes - # the generated inputs invariant to ``n_workers``: the worker that - # runs simulation ``i`` always seeds from ``child_seeds[i]``, - # regardless of how many workers there are. The full list is shared - # with every worker; the shared atomic counter assigns indices. - child_seeds = np.random.SeedSequence(random_seed).spawn( - self.number_of_simulations - ) - - for _ in range(n_workers): - sim_producer = multiprocess.Process( - target=self.__sim_producer, - args=( - child_seeds, - sim_monitor, - mutex, - simulation_error_event, - ), - ) - processes.append(sim_producer) - sim_producer.start() - + # Started workers only, and inside the try, so a ``start()`` that + # fails part way through the fleet does not leave the ones already + # running with nobody to clean them up. + started_processes = [] try: - for sim_producer in processes: - sim_producer.join() - - # Handle error from the child processes - if simulation_error_event.is_set(): - raise RuntimeError( - "An error occurred during the simulation. \n" - f"Check the logs and error file {self.error_file} " - "for more information." + # Each worker derives one independent child seed per simulation + # index (not per worker) from the shared root state: the counter + # assigns indices and index i always seeds from __child_seed(i), + # so the sampled inputs do not depend on the number of workers. + # The root state is small and travels with the pickled instance, + # so no per-index seed list is materialized or sent. + for _ in range(n_workers): + sim_producer = multiprocess.Process( + target=self.__sim_producer, + args=( + sim_monitor, + mutex, + simulation_error_event, + ), ) + sim_producer.start() + started_processes.append(sim_producer) + + _wait_for_workers(started_processes, simulation_error_event) + _stop_any_worker_still_running(started_processes) + _fail_if_a_worker_did_not_finish( + started_processes, simulation_error_event, self.error_file + ) sim_monitor.print_final_status() # Handle error from the main process # pylint: disable=broad-except except (Exception, KeyboardInterrupt) as error: - simulation_error_event.set() - - for sim_producer in processes: - sim_producer.join() - - if not isinstance(error, KeyboardInterrupt): + _bring_the_fleet_down(started_processes, simulation_error_event) + self._interrupted = isinstance(error, KeyboardInterrupt) + if not self._interrupted: raise error + finally: + _stop_any_worker_still_running(started_processes) def __validate_number_of_workers(self, n_workers): - if n_workers is None or n_workers > os.cpu_count(): - n_workers = os.cpu_count() + # os.cpu_count() is documented as possibly None, and comparing against + # it then raises rather than falling back to a usable default. + available = os.cpu_count() or 2 + if n_workers is not None and type(n_workers) not in (int, np.integer): # noqa: E721 + raise TypeError( + f"Number of workers must be an integer, not {type(n_workers).__name__}." + ) + if n_workers is None or n_workers > available: + n_workers = available if n_workers < 2: raise ValueError("Number of workers must be at least 2 for parallel mode.") return n_workers - def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event): # pylint: disable=too-many-statements + def __sim_producer(self, sim_monitor, mutex, error_event): # pylint: disable=too-many-statements """Simulation producer to be used in parallel by multiprocessing. Parameters ---------- - child_seeds : list[numpy.random.SeedSequence] - One seed sequence per simulation index. Before each simulation the - worker seeds the stochastic models from ``child_seeds[sim_idx]``, - where ``sim_idx`` is pulled from the shared atomic counter. This - keeps the generated inputs invariant to the number of workers. sim_monitor : _SimMonitor The simulation monitor object to keep track of the simulations. mutex : multiprocess.Lock @@ -439,23 +593,33 @@ def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event): # pylin Event signaling an error occurred during the simulation. """ try: - while sim_monitor.keep_simulating(): - sim_idx = sim_monitor.increment() - 1 - inputs_json, outputs_json = "", "" - - # Seed per simulation index so the inputs are reproducible and - # independent of which worker happens to run this index. - self.__seed_components(child_seeds[sim_idx]) - + while True: + # First statement in the loop, so it is bound before the claim + # rather than after it. A claim that failed left these unassigned + # and the handler raised UnboundLocalError over the real error; + # a claim that failed on a later lap reported the previous row. + sim_idx, inputs_json, outputs_json = None, "", "" + + sim_idx = _claim_next_index(sim_monitor, mutex) + if sim_idx is None: + break + + self.__seed_simulation(self.__child_seed(sim_idx)) flight = self.__run_single_simulation() inputs_json = self.__evaluate_flight_inputs(sim_idx) outputs_json = self.__evaluate_flight_outputs(flight, sim_idx) + acquired = False try: mutex.acquire() + acquired = True if error_event.is_set(): - sim_monitor.reprint( - "Simulation Interrupt, files from simulation " + # Runs in a worker process spawned via multiprocessing: + # logging handlers configured in the main process are + # not guaranteed to be inherited (e.g. Windows "spawn"), + # so this must use print() to remain visible. + _SimMonitor.reprint( + f"Simulation interrupt. Files from simulation " f"{sim_idx} saved." ) with open(self.error_file, "a", encoding="utf-8") as f: @@ -463,49 +627,55 @@ def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event): # pylin break - with open(self.input_file, "a", encoding="utf-8") as f: - f.write(inputs_json) - with open(self.output_file, "a", encoding="utf-8") as f: - f.write(outputs_json) - + _record_simulation( + self.input_file, self.output_file, inputs_json, outputs_json + ) sim_monitor.print_update_status() finally: - mutex.release() - - except Exception: # pylint: disable=broad-except - mutex.acquire() - with open(self.error_file, "a", encoding="utf-8") as f: - f.write(inputs_json) - - sim_monitor.reprint(f"Error on iteration {sim_idx}:") - sim_monitor.reprint(traceback.format_exc()) - error_event.set() - mutex.release() + if acquired: + mutex.release() - def __seed_components(self, simulation_seed): - """Seed the stochastic models for a single simulation index. + except Exception: + # Set first, so a parent waiting on the join learns why. Best effort + # like everything below it: this is a manager proxy, the manager may + # already be gone, and reporting must not replace what it reports. + try: + error_event.set() + except Exception: # pylint: disable=broad-exception-caught + pass + details = traceback.format_exc() - The given seed sequence is split into three independent sub-streams so - that the environment, rocket and flight do not share the same random - draws (sharing a single seed would correlate their first sampled - values). Seeding per simulation index -- rather than once per worker -- - is what makes the Monte Carlo inputs invariant to the execution mode - (serial vs parallel) and to the number of workers. + # The failure goes onto the inputs record rather than replacing it. + # Writing one or the other dropped the traceback for every failure + # after sampling, from the file the run tells the user to read. + try: + record = json.loads(inputs_json) if inputs_json else {"index": sim_idx} + except ValueError: + record = {"index": sim_idx} + record["error"] = details + record = json.dumps(record) + "\n" - Parameters - ---------- - simulation_seed : numpy.random.SeedSequence - The seed sequence assigned to the current simulation index. It is - spawned into three child sequences, one per stochastic model. + acquired = False + try: + mutex.acquire() + acquired = True + with open(self.error_file, "a", encoding="utf-8") as f: + f.write(record) + + # See note above: must use print() to remain visible from a + # multiprocessing worker process. + _SimMonitor.reprint(f"Error on iteration {sim_idx}:\n{details}") + except Exception: # pylint: disable=broad-exception-caught + # The mutex or the error file is unreachable too. Reporting is + # not worth losing the failure that started this. + pass + finally: + if acquired: + mutex.release() - Returns - ------- - None - """ - env_seed, rocket_seed, flight_seed = simulation_seed.spawn(3) - self.environment._set_stochastic(env_seed) - self.rocket._set_stochastic(rocket_seed) - self.flight._set_stochastic(flight_seed) + # The worker exits non-zero, so the parent can tell a crash from a + # clean finish rather than only from the error event. + raise def __run_single_simulation(self): """Runs a single simulation and returns the inputs and outputs. @@ -590,6 +760,107 @@ def estimate_confidence_interval( return res.confidence_interval + def simulate_convergence( + self, + target_attribute="apogee_time", + target_confidence=0.95, + tolerance=0.5, + max_simulations=1000, + batch_size=50, + parallel=False, + n_workers=None, + ): + """Run Monte Carlo simulations in batches until the confidence interval + width converges within the specified tolerance or the maximum number of + simulations is reached. + + Parameters + ---------- + target_attribute : str + The target attribute to track its convergence (e.g., "apogee", "apogee_time", etc.). + target_confidence : float, optional + The confidence level for the interval (between 0 and 1). Default is 0.95. + tolerance : float, optional + The desired width of the confidence interval in seconds, meters, or other units. Default is 0.5. + max_simulations : int, optional + The maximum number of simulations to run to avoid infinite loops. Default is 1000. + batch_size : int, optional + The number of simulations to run in each batch. Default is 50. + parallel : bool, optional + Whether to run simulations in parallel. Default is False. + n_workers : int, optional + The number of worker processes to use if running in parallel. Default is None. + + Returns + ------- + confidence_interval_history : list of float + History of confidence interval widths, one value per batch of simulations. + The last element corresponds to the width when the simulation stopped for + either meeting the tolerance or reaching the maximum number of simulations. + """ + + # Validate inputs up-front. Without this, a non-positive batch_size makes + # the loop run zero new simulations every iteration and spin forever. + if not isinstance(batch_size, (int, np.integer)) or batch_size <= 0: + raise ValueError( + f"'batch_size' must be a positive integer, got {batch_size!r}." + ) + if not isinstance(max_simulations, (int, np.integer)) or max_simulations <= 0: + raise ValueError( + f"'max_simulations' must be a positive integer, got " + f"{max_simulations!r}." + ) + if not isinstance(tolerance, (int, float)) or tolerance <= 0: + raise ValueError( + f"'tolerance' must be a positive number, got {tolerance!r}." + ) + if not 0 < target_confidence < 1: + raise ValueError( + "'target_confidence' must be between 0 and 1 (exclusive), got " + f"{target_confidence!r}." + ) + + self.import_outputs(self.filename.with_suffix(".outputs.txt")) + confidence_interval_history = [] + + while self.num_of_loaded_sims < max_simulations: + total_sims = min(self.num_of_loaded_sims + batch_size, max_simulations) + + self.simulate( + number_of_simulations=total_sims, + append=True, + include_function_data=False, + parallel=parallel, + n_workers=n_workers, + ) + + self.import_outputs(self.filename.with_suffix(".outputs.txt")) + + ci = self.estimate_confidence_interval( + attribute=target_attribute, + confidence_level=target_confidence, + ) + + width = float(ci.high - ci.low) + confidence_interval_history.append(width) + + # A NaN width means the target attribute contains NaN values; the + # tolerance check would never pass, so the loop would run to + # max_simulations and silently return a NaN history. Stop and warn. + if np.isnan(width): + warnings.warn( + f"The confidence interval width for '{target_attribute}' is " + "NaN, likely because the attribute contains NaN values. " + "Stopping convergence early; check the simulation outputs.", + stacklevel=2, + ) + break + + if width <= tolerance: + break + + return confidence_interval_history + def __evaluate_flight_inputs(self, sim_idx): """Evaluates the inputs of a single flight simulation. @@ -666,7 +937,7 @@ def __terminate_simulation(self): self.output_file = self._output_file self.error_file = self._error_file - _SimMonitor.reprint(f"Results saved to {self._output_file}") + print(f"Results saved to {self._output_file}") def __check_export_list(self, export_list): """ @@ -872,58 +1143,232 @@ def error_file(self, value): self._error_file = value self.set_errors_log() + # File format helpers + + @staticmethod + def _detect_file_format(filepath): + """Detect file format from the file extension. + + Parameters + ---------- + filepath : str or Path + Path to the file. + + Returns + ------- + str + One of ``"jsonl"``, ``"csv"``, or ``"json"``. + + Raises + ------ + ValueError + If the file extension is not supported. + """ + suffix = Path(filepath).suffix.lower() + format_map = {".txt": "jsonl", ".csv": "csv", ".json": "json"} + if suffix not in format_map: + raise ValueError( + f"Unsupported file extension '{suffix}'. " + "Expected '.txt', '.csv', or '.json'." + ) + return format_map[suffix] + + @staticmethod + def _parse_csv_value(value): + """Parse a string value from a CSV cell into its appropriate type. + + Parameters + ---------- + value : str + The raw string value from the CSV cell. + + Returns + ------- + int, float, dict, list, or str + The parsed value in its appropriate Python type. + """ + if value == "": + return value + # Try parsing JSON objects/arrays + if value.startswith(("{", "[")): + try: + return json.loads(value) + except (json.JSONDecodeError, ValueError): + pass + # Try numeric types + try: + int_val = int(value) + # Ensure the string was truly an integer (not "1.0") + if str(int_val) == value: + return int_val + except ValueError: + pass + try: + return float(value) + except ValueError: + pass + return value + + def _read_log_file(self, filepath): + """Read a log file in any supported format and return a list of dicts. + + Parameters + ---------- + filepath : str or Path + Path to the log file. Format is detected from the extension. + + Returns + ------- + list of dict + A list of dictionaries, one per simulation record. + """ + fmt = self._detect_file_format(filepath) + result = [] + with open(filepath, mode="r", encoding="utf-8") as f: + if fmt == "jsonl": + for line in f: + line = line.strip() + if line: + result.append(json.loads(line)) + elif fmt == "json": + content = f.read().strip() + if content: + result = json.loads(content) + elif fmt == "csv": + reader = csv.DictReader(f) + for row in reader: + result.append({k: self._parse_csv_value(v) for k, v in row.items()}) + return result + + @staticmethod + def _write_log_to_csv(log_data, filepath, flatten=False): + """Write a list of dicts to a CSV file. + + Parameters + ---------- + log_data : list of dict + The data to write. Each dict is one row. + filepath : str or Path + Output file path. + flatten : bool, optional + If True, non-scalar columns (dicts, lists) are omitted. + If False (default), non-scalar values are serialized as JSON + strings in the CSV cells. + + Raises + ------ + ValueError + If ``log_data`` is empty. + """ + if not log_data: + raise ValueError( + "No data to export. Run a simulation first or import existing data." + ) + # Collect all keys preserving insertion order + all_keys = list(dict.fromkeys(k for row in log_data for k in row)) + + if flatten: + # Identify scalar-only keys + scalar_keys = [] + for key in all_keys: + if all(not isinstance(row.get(key), (dict, list)) for row in log_data): + scalar_keys.append(key) + fieldnames = scalar_keys + else: + fieldnames = all_keys + + with open(filepath, mode="w", newline="", encoding="utf-8") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames, extrasaction="ignore") + writer.writeheader() + for row in log_data: + csv_row = {} + for key in fieldnames: + value = row.get(key, "") + if isinstance(value, (dict, list)): + csv_row[key] = json.dumps(value) + else: + csv_row[key] = value + writer.writerow(csv_row) + + def _write_log_to_json(self, log_data, filepath): + """Write a list of dicts to a JSON file as a proper JSON array. + + Parameters + ---------- + log_data : list of dict + The data to write. Each dict becomes one element of the array. + filepath : str or Path + Output file path. + + Raises + ------ + ValueError + If ``log_data`` is empty. + """ + if not log_data: + raise ValueError( + "No data to export. Run a simulation first or import existing data." + ) + with open(filepath, mode="w", encoding="utf-8") as f: + json.dump(log_data, f, cls=RocketPyEncoder, indent=2) + # Setters for post simulation attributes def set_inputs_log(self): """ Sets inputs_log from a file into an attribute for easy access. + Supports .txt (JSONL), .csv, and .json file formats. Returns ------- None """ - self.inputs_log = [] - with open(self.input_file, mode="r", encoding="utf-8") as rows: - for line in rows: - self.inputs_log.append(json.loads(line)) + self.inputs_log = self._read_log_file(self.input_file) def set_outputs_log(self): """ Sets outputs_log from a file into an attribute for easy access. + Supports .txt (JSONL), .csv, and .json file formats. Returns ------- None """ - self.outputs_log = [] - with open(self.output_file, mode="r", encoding="utf-8") as rows: - for line in rows: - self.outputs_log.append(json.loads(line)) + self.outputs_log = self._read_log_file(self.output_file) def set_errors_log(self): """ Sets errors_log from a file into an attribute for easy access. + Supports .txt (JSONL), .csv, and .json file formats. Returns ------- None """ - self.errors_log = [] - with open(self.error_file, mode="r", encoding="utf-8") as errors: - for line in errors: - self.errors_log.append(json.loads(line)) + self.errors_log = self._read_log_file(self.error_file) def set_num_of_loaded_sims(self): """ Determines the number of simulations loaded from output_file being - currently used. + currently used. Supports .txt (JSONL), .csv, and .json formats. Returns ------- None """ + fmt = self._detect_file_format(self.output_file) with open(self.output_file, mode="r", encoding="utf-8") as outputs: - self.num_of_loaded_sims = sum(1 for _ in outputs) + if fmt == "jsonl": + self.num_of_loaded_sims = sum(1 for _ in outputs) + elif fmt == "csv": + # Subtract 1 for the header row + self.num_of_loaded_sims = max(0, sum(1 for _ in outputs) - 1) + elif fmt == "json": + content = outputs.read().strip() + if content: + self.num_of_loaded_sims = len(json.loads(content)) + else: + self.num_of_loaded_sims = 0 def set_results(self): """ @@ -980,13 +1425,15 @@ def set_processed_results(self): def import_outputs(self, filename=None): """ - Import Monte Carlo results from .txt file and save it into a dictionary. + Import Monte Carlo results from a file and save it into a dictionary. + Supports .txt (JSONL), .csv, and .json file formats. Parameters ---------- filename : str, optional Name or directory path to the file to be imported. If none, - self.filename will be used. + self.filename will be used with the default .outputs.txt suffix. + Files with .csv or .json extensions are also accepted. Returns ------- @@ -994,7 +1441,7 @@ def import_outputs(self, filename=None): Notes ----- - Notice that you can import the outputs, inputs, and errors from the a + Notice that you can import the outputs, inputs, and errors from a file without the need to run simulations. You can use previously saved files to process analyze the results or to continue a simulation. """ @@ -1007,20 +1454,22 @@ def import_outputs(self, filename=None): with open(filepath, "w+", encoding="utf-8"): self.output_file = filepath - _SimMonitor.reprint( - f"A total of {self.num_of_loaded_sims} simulations results were " - f"loaded from the following output file: {self.output_file}\n" + print( + f"A total of {self.num_of_loaded_sims} simulation results were " + f"loaded from: {self.output_file}" ) def import_inputs(self, filename=None): """ - Import Monte Carlo inputs from .txt file and save it into a dictionary. + Import Monte Carlo inputs from a file and save it into a dictionary. + Supports .txt (JSONL), .csv, and .json file formats. Parameters ---------- filename : str, optional Name or directory path to the file to be imported. If none, - self.filename will be used. + self.filename will be used with the default .inputs.txt suffix. + Files with .csv or .json extensions are also accepted. Returns ------- @@ -1035,17 +1484,19 @@ def import_inputs(self, filename=None): with open(filepath, "w+", encoding="utf-8"): self.input_file = filepath - _SimMonitor.reprint(f"The following input file was imported: {self.input_file}") + print(f"The following input file was imported: {self.input_file}") def import_errors(self, filename=None): """ - Import Monte Carlo errors from .txt file and save it into a dictionary. + Import Monte Carlo errors from a file and save it into a dictionary. + Supports .txt (JSONL), .csv, and .json file formats. Parameters ---------- filename : str, optional Name or directory path to the file to be imported. If none, - self.filename will be used. + self.filename will be used with the default .errors.txt suffix. + Files with .csv or .json extensions are also accepted. Returns ------- @@ -1060,11 +1511,11 @@ def import_errors(self, filename=None): with open(filepath, "w+", encoding="utf-8"): self.error_file = filepath - _SimMonitor.reprint(f"The following error file was imported: {self.error_file}") + print(f"The following error file was imported: {self.error_file}") def import_results(self, filename=None): """ - Import Monte Carlo results from .txt file and save it into a dictionary. + Import Monte Carlo results from a file and save it into a dictionary. Parameters ---------- @@ -1153,7 +1604,7 @@ def export_ellipses_to_kml( # pylint: disable=too-many-statements except KeyError as e: raise KeyError("No impact data found. Skipping impact ellipses.") from e - (apogee_ellipses, impact_ellipses) = generate_monte_carlo_ellipses( + apogee_ellipses, impact_ellipses = generate_monte_carlo_ellipses( impact_x, impact_y, apogee_x, @@ -1279,6 +1730,283 @@ def compare_ellipses(self, other_monte_carlo, **kwargs): """ self.plots.ellipses_comparison(other_monte_carlo, **kwargs) + # CSV and JSON export methods + + def export_outputs_to_csv(self, filename): + """Export simulation outputs to a CSV file. + + Each row represents one simulation. All output values are scalar, + so the CSV is directly usable in spreadsheet applications. + + Parameters + ---------- + filename : str + Path to the output CSV file. + + Raises + ------ + ValueError + If no output data is available to export. + """ + self._write_log_to_csv(self.outputs_log, filename) + + def export_outputs_to_json(self, filename): + """Export simulation outputs to a JSON file as an array of objects. + + Parameters + ---------- + filename : str + Path to the output JSON file. + + Raises + ------ + ValueError + If no output data is available to export. + """ + self._write_log_to_json(self.outputs_log, filename) + + def export_inputs_to_csv(self, filename, flatten=False): + """Export simulation inputs to a CSV file. + + Parameters + ---------- + filename : str + Path to the output CSV file. + flatten : bool, optional + If True, columns with non-scalar values (dicts, lists) are + omitted from the CSV. If False (default), non-scalar values + are serialized as JSON strings within the CSV cells. + + Raises + ------ + ValueError + If no input data is available to export. + """ + self._write_log_to_csv(self.inputs_log, filename, flatten=flatten) + + def export_inputs_to_json(self, filename): + """Export simulation inputs to a JSON file as an array of objects. + + Parameters + ---------- + filename : str + Path to the output JSON file. + + Raises + ------ + ValueError + If no input data is available to export. + """ + self._write_log_to_json(self.inputs_log, filename) + + def export_errors_to_csv(self, filename, flatten=False): + """Export simulation errors to a CSV file. + + Parameters + ---------- + filename : str + Path to the output CSV file. + flatten : bool, optional + If True, columns with non-scalar values (dicts, lists) are + omitted from the CSV. If False (default), non-scalar values + are serialized as JSON strings within the CSV cells. + + Raises + ------ + ValueError + If no error data is available to export. + """ + self._write_log_to_csv(self.errors_log, filename, flatten=flatten) + + def export_errors_to_json(self, filename): + """Export simulation errors to a JSON file as an array of objects. + + Parameters + ---------- + filename : str + Path to the output JSON file. + + Raises + ------ + ValueError + If no error data is available to export. + """ + self._write_log_to_json(self.errors_log, filename) + + +def _recorded_indices(label, path): + """``{index: how many rows carry it}`` for one log file. + + Strict about what a row is. A row that will not parse, is not an + object, or carries anything but a non-negative plain ``int`` index is + the corruption this check exists to find, so it is named and raised on + rather than skipped. ``type(...) is int`` and not ``isinstance``: + ``True`` and ``1.0`` both compare equal to ``1`` and would otherwise + pass for it. + """ + written = {} + with open(path, mode="r", encoding="utf-8") as rows: + for number, line in enumerate(rows, start=1): + line = line.strip() + if not line: + continue + try: + record = json.loads(line) + except ValueError as error: + raise RuntimeError( + f"{label} row {number} is not readable JSON, so a " + f"worker was cut off part way through writing it: " + f"{line[:60]!r}" + ) from error + index = record.get("index") if isinstance(record, dict) else None + # isinstance is the wrong tool here, see the docstring: bool is a + # subclass of int, so True would pass for the index 1. + # pylint: disable-next=unidiomatic-typecheck + if type(index) is not int or index < 0: # noqa: E721 + raise RuntimeError( + f"{label} row {number} does not carry a simulation " + f"index: {line[:60]!r}" + ) + written[index] = written.get(index, 0) + 1 + return written + + +def _validate_simulation_count(number_of_simulations): + """A count has to be a whole non-negative number, checked before any file. + + ``type(...) is not int``: ``True`` is an ``int`` to ``isinstance`` and would + quietly run one simulation. A float ran ``int(count)`` of them and then + failed the completeness check with a range it could never have satisfied. + """ + if type(number_of_simulations) not in (int, np.integer): # noqa: E721 + raise TypeError( + f"number_of_simulations must be an integer, not " + f"{type(number_of_simulations).__name__}." + ) + if number_of_simulations < 0: + raise ValueError( + f"number_of_simulations must not be negative, got {number_of_simulations}." + ) + + +_WORKER_SHUTDOWN_GRACE = 5.0 + + +def _record_simulation(input_file, output_file, inputs_json, outputs_json): + """Append one simulation's inputs and outputs to their logs. + + Module level rather than a method: the run paths are driven directly by + stub objects in the tests, and a private method is not reachable on those. + """ + with open(input_file, "a", encoding="utf-8") as f: + f.write(inputs_json) + with open(output_file, "a", encoding="utf-8") as f: + f.write(outputs_json) + + +def _bring_the_fleet_down(started_processes, error_event): + """Stop everything, without raising over the failure being handled. + + Setting the event is best effort like the workers' own reporting: the + manager may be the thing that died. Then a bounded window to notice it and + leave, and whatever is left gets stopped. + """ + try: + error_event.set() + except Exception: # pylint: disable=broad-exception-caught + pass + _wait_for_workers(started_processes, timeout=_WORKER_SHUTDOWN_GRACE) + _stop_any_worker_still_running(started_processes) + + +def _wait_for_workers(started_processes, error_event=None, timeout=None): + """Wait for the fleet, giving up early once one of them reports an error. + + Joining each worker in turn waits on them in the order they were started. A + worker stuck in a native call held the parent on the first join while + another had already set the event, so neither the error nor the cleanup + after it was ever reached. + + No overall deadline on the normal path: a run with no error and one worker + still going is a long simulation, and that is not for this to cut short. + """ + deadline = None if timeout is None else monotonic() + timeout + while any(process.is_alive() for process in started_processes): + if error_event is not None and error_event.is_set(): + break + if deadline is not None and monotonic() >= deadline: + break + for process in started_processes: + process.join(timeout=0.1) + + # Reap whatever has already finished. A worker that was gone before the + # loop started was never joined by it, and an unjoined child has no exit + # code yet, so the crash check downstream would read None and call it one. + for process in started_processes: + process.join(timeout=0) + + +def _stop_any_worker_still_running(started_processes, grace=_WORKER_SHUTDOWN_GRACE): + """Whatever is still going here is not going to stop on its own. + + Signal every worker before waiting on any of them. Terminating one and + joining it before reaching the next let a worker that ignores the signal + keep the rest of the fleet, the manager and the open files alive behind it. + """ + alive = [process for process in started_processes if process.is_alive()] + for process in alive: + process.terminate() + for process in alive: + process.join(timeout=grace) + + # terminate is a request. SIGKILL is not, and a worker that sat through the + # first one would otherwise keep the manager and the files open for good. + stubborn = [process for process in alive if process.is_alive()] + for process in stubborn: + process.kill() + for process in stubborn: + process.join(timeout=grace) + + +def _fail_if_a_worker_did_not_finish(started_processes, error_event, error_file): + """Raise unless every worker finished and none of them reported an error. + + A worker can die without ever setting the event: SystemExit, ``os._exit``, a + segfault in a native extension, a target that will not unpickle under spawn, + or the error handler itself failing. ``join()`` returns None whatever + happened, so the exit status is the only thing that separates a crash from a + clean finish. + """ + crashed = [ + f"{sim_producer.name} exited with {sim_producer.exitcode}" + for sim_producer in started_processes + if sim_producer.exitcode != 0 + ] + if error_event.is_set() or crashed: + raise RuntimeError( + "An error occurred during the simulation. \n" + + (f"Workers that did not exit cleanly: {crashed}. \n" if crashed else "") + + f"Check the logs and error file {error_file} for more information." + ) + + +def _claim_next_index(sim_monitor, mutex): + """Atomically claim the next 0-based simulation index, or ``None`` if done. + + ``keep_simulating()`` and ``increment()`` are two separate manager calls, so + the shared ``mutex`` has to be held across both. Without it, two workers can + each pass the ``count < number_of_simulations`` check at the tail before + either increments, and both then claim an index, running more simulations + than were requested (and duplicating a simulation index). + """ + mutex.acquire() + try: + if not sim_monitor.keep_simulating(): + return None + return sim_monitor.increment() - 1 + finally: + mutex.release() + def _import_multiprocess(): """Import the necessary modules and submodules for the @@ -1373,15 +2101,12 @@ def print_final_status(self): msg = f"Completed {self.count - self.initial_count} iterations." msg += f" In total, {self.count} simulations are exported.\n" msg += f"Total wall time: {time() - self.start_time:.1f} s" - _SimMonitor.reprint(msg, end="\n", flush=True) @staticmethod def reprint(msg, end="\n", flush=True): - """ - Prints a message on the same line as the previous one and replaces the - previous message with the new one, deleting the extra characters from - the previous message. + """Prints a message replacing the previous line to avoid cluttering + the terminal output during concurrent simulation progress updates. Parameters ---------- @@ -1396,12 +2121,8 @@ def reprint(msg, end="\n", flush=True): ------- None """ - padding = "" - if len(msg) < _SimMonitor._last_print_len: padding = " " * (_SimMonitor._last_print_len - len(msg)) - print(msg + padding, end=end, flush=flush) - _SimMonitor._last_print_len = len(msg) diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index 879b61e70..b8af7c516 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/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 diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index c9e0c9ca6..79a397984 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/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 ---------- diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 68ab3404a..9df900eb5 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -11,6 +11,7 @@ import importlib import importlib.metadata import json +import logging import math import re import time @@ -25,6 +26,8 @@ from matplotlib.patches import Ellipse from packaging import version as packaging_version +logger = logging.getLogger(__name__) + # Mapping of module name and the name of the package that should be installed INSTALL_MAPPING = {"IPython": "ipython"} @@ -422,18 +425,18 @@ def inverted_haversine(lat0, lon0, distance, bearing, earth_radius=6.3781e6): lon0_rad = np.deg2rad(lon0) # Apply inverted Haversine formula - lat1_rad = math.asin( - math.sin(lat0_rad) * math.cos(distance / earth_radius) - + math.cos(lat0_rad) - * math.sin(distance / earth_radius) - * math.cos(math.radians(bearing)) + lat1_rad = np.arcsin( + np.sin(lat0_rad) * np.cos(distance / earth_radius) + + np.cos(lat0_rad) + * np.sin(distance / earth_radius) + * np.cos(np.radians(bearing)) ) - lon1_rad = lon0_rad + math.atan2( - math.sin(math.radians(bearing)) - * math.sin(distance / earth_radius) - * math.cos(lat0_rad), - math.cos(distance / earth_radius) - math.sin(lat0_rad) * math.sin(lat1_rad), + lon1_rad = lon0_rad + np.arctan2( + np.sin(np.radians(bearing)) + * np.sin(distance / earth_radius) + * np.cos(lat0_rad), + np.cos(distance / earth_radius) - np.sin(lat0_rad) * np.sin(lat1_rad), ) # Convert back to degrees and then return @@ -1464,11 +1467,34 @@ 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 res = doctest.testmod() if res.failed < 1: - print(f"All the {res.attempted} tests passed!") + logger.info("All the %d tests passed!", res.attempted) else: - print(f"{res.failed} out of {res.attempted} tests failed.") + logger.error("%d out of %d tests failed.", res.failed, res.attempted) diff --git a/tests/integration/simulation/test_monte_carlo_determinism.py b/tests/integration/simulation/test_monte_carlo_determinism.py new file mode 100644 index 000000000..1d65ec78e --- /dev/null +++ b/tests/integration/simulation/test_monte_carlo_determinism.py @@ -0,0 +1,787 @@ +"""End-to-end determinism tests for ``MonteCarlo.simulate(random_seed=...)``. + +With a fixed ``random_seed`` the generated random *inputs* are reproducible and +identical across serial and parallel execution and across any number of workers. +Each simulation index draws from its own child stream spawned from the run's root +seed, and ``SeedSequence.spawn`` is prefix-stable, so index ``i`` maps to the same +seed regardless of the worker that runs it. (The seed-handling helpers themselves +are unit tested in ``tests/unit/simulation/test_monte_carlo_determinism``.) + +The trajectory integration (``Flight``) is stubbed: worker invariance is a +property of the *input sampling*, which happens before ``Flight`` is built, so a +stub keeps the runs fast while still driving the real serial and parallel loops. +Stubbing the module-level ``Flight`` symbol reaches the parallel workers only +under the ``fork`` start method, so the worker-invariance test skips otherwise and +is marked ``slow`` to match the other Monte Carlo multiprocessing tests. + +A dedicated numpy-only rocket keeps the fork-based end-to-end test simple: it +gives the motor a single ``thrust_source`` so the run has no list-valued attribute +at all. List sampling is itself seeded now (it draws through the model generator, +not the stdlib ``random.choice``) and is covered directly in +``tests/unit/stochastic/test_stochastic_model``. + +Seed derivation being independent of the multiprocessing start method (fork, +spawn or forkserver) is verified separately by +``test_seed_derivation_is_start_method_invariant``, which uses a top-level +picklable target so it is safe under ``spawn``/``forkserver`` -- unlike the +``Flight``-stub test above, which reaches workers only under ``fork``. +""" + +import json +import multiprocessing +import os +from types import SimpleNamespace + +import numpy as np +import pytest + +import rocketpy.simulation.monte_carlo as mc_module +from rocketpy import Environment +from rocketpy.simulation import MonteCarlo +from rocketpy.simulation.monte_carlo import _seed_sequence_to_int +from rocketpy.stochastic import ( + StochasticAirBrakes, + StochasticEnvironment, + StochasticRocket, + StochasticSolidMotor, +) + +_child_seed = MonteCarlo._MonteCarlo__child_seed + + +def _available_start_methods(): + """The multiprocessing start methods this platform actually supports.""" + supported = multiprocessing.get_all_start_methods() + return [method for method in ("fork", "spawn", "forkserver") if method in supported] + + +def _derive_index_seeds(root_state, indices): + """Derive the per-index seed fingerprints from ``root_state``. + + Top-level and picklable (only a small tuple and a list of ints cross the + process boundary), so it runs unchanged under every start method -- including + ``spawn``/``forkserver``, which re-import this module rather than inheriting + the parent's memory. It calls the real production helpers (``__child_seed`` + and ``_seed_sequence_to_int``) so the test tracks the shipped derivation. + """ + plan = SimpleNamespace(_MonteCarlo__root_state=root_state) + return {index: _seed_sequence_to_int(_child_seed(plan, index)) for index in indices} + + +class _StubFlight: + """Minimal stand-in for ``Flight`` that skips trajectory integration.""" + + def __init__(self, **kwargs): # accepts and ignores MonteCarlo's Flight kwargs + pass + + def __getattr__(self, name): + return 0.0 + + +@pytest.fixture +def stochastic_calisto_numpy_only( + cesaroni_m1670, + calisto_robust, + stochastic_nose_cone, + stochastic_trapezoidal_fins, + stochastic_tail, + stochastic_rail_buttons, + stochastic_main_parachute, + stochastic_drogue_parachute, +): + """A ``StochasticRocket`` whose randomness flows entirely through numpy. + + Mirrors the shared ``stochastic_calisto`` fixture but gives the solid motor a + single ``thrust_source`` instead of a multi-element list, so no attribute is + sampled through the unseeded standard-library ``random.choice``. + """ + motor = StochasticSolidMotor( + solid_motor=cesaroni_m1670, + burn_out_time=(4, 0.1), + grains_center_of_mass_position=0.001, + grain_density=50, + grain_separation=1 / 1000, + grain_initial_height=1 / 1000, + grain_initial_inner_radius=0.375 / 1000, + grain_outer_radius=0.375 / 1000, + total_impulse=(6500, 1000), + throat_radius=0.5 / 1000, + nozzle_radius=0.5 / 1000, + nozzle_position=0.001, + ) + rocket = StochasticRocket( + rocket=calisto_robust, + radius=0.0127 / 2000, + mass=(15.426, 0.5, "normal"), + inertia_11=(6.321, 0), + inertia_22=0.01, + inertia_33=0.01, + center_of_mass_without_motor=0, + ) + rocket.add_motor(motor, position=0.001) + rocket.add_nose(stochastic_nose_cone, position=(1.134, 0.001)) + rocket.add_trapezoidal_fins(stochastic_trapezoidal_fins, position=(0.001, "normal")) + rocket.add_tail(stochastic_tail) + rocket.set_rail_buttons( + stochastic_rail_buttons, lower_button_position=(-0.618, 0.001, "normal") + ) + rocket.add_parachute(parachute=stochastic_main_parachute) + rocket.add_parachute(parachute=stochastic_drogue_parachute) + return rocket + + +def _read_inputs_by_index(input_file): + """Read a ``.inputs.txt`` file into ``{index: raw_json_line}``.""" + by_index = {} + with open(input_file, mode="r", encoding="utf-8") as rows: + for line in rows: + line = line.strip() + if not line: + continue + by_index[json.loads(line)["index"]] = line + return by_index + + +def _count_rows(log_file): + """How many records were written, before anything is keyed by index. + + Keying by index hides a duplicate: two workers claiming the same index + write two rows and the second overwrites the first in the dict, so the + result looks complete. The claim is meant to be atomic, and the count is + what says so. + """ + with open(log_file, mode="r", encoding="utf-8") as rows: + return sum(1 for line in rows if line.strip()) + + +def _simulate_inputs( + monkeypatch, tmp_path, environment, rocket, flight, tag, **simulate_kwargs +): + """Run a Monte Carlo with a stubbed ``Flight`` and return inputs by index.""" + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + montecarlo = MonteCarlo( + filename=str(tmp_path / tag), + environment=environment, + rocket=rocket, + flight=flight, + ) + montecarlo.simulate(**simulate_kwargs) + return _read_inputs_by_index(montecarlo.input_file) + + +def test_invalid_seed_does_not_truncate_existing_output( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """A rejected seed must fail before any output file is truncated, so passing + an invalid seed cannot destroy the results of a previous run.""" + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + montecarlo = MonteCarlo( + filename=str(tmp_path / "keep"), + environment=stochastic_environment, + rocket=stochastic_calisto_numpy_only, + flight=stochastic_flight, + ) + with open(montecarlo.input_file, "w", encoding="utf-8") as existing: + existing.write("previous results\n") + + # A Generator is not a seed and is rejected; the run must raise before the + # ``w+`` file setup truncates anything. + with pytest.raises(TypeError): + montecarlo.simulate( + number_of_simulations=3, random_seed=np.random.default_rng(0) + ) + + with open(montecarlo.input_file, encoding="utf-8") as kept: + assert kept.read() == "previous results\n" + + +@pytest.mark.parametrize( + ("kwargs", "error"), + [ + ({"number_of_simulations": 2.5}, TypeError), + ({"number_of_simulations": True}, TypeError), + ({"number_of_simulations": -1}, ValueError), + ({"number_of_simulations": 3, "parallel": True, "n_workers": 1}, ValueError), + ], + ids=["float count", "boolean count", "negative count", "one worker"], +) +def test_a_rejected_argument_does_not_truncate_existing_output( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, + kwargs, + error, +): + """Every check that needs only the arguments belongs before the logs open. + + ``__setup_files`` opens both of them "w+", which empties them, and + ``n_workers`` was validated after that. So asking for a worker count the run + cannot use destroyed the previous run's results on the way to raising. + + ``True`` is the one that does not raise on its own: it is an ``int`` to + ``isinstance``, so it would quietly have run one simulation. + """ + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + montecarlo = MonteCarlo( + filename=str(tmp_path / f"keep-{sorted(kwargs.items())}"), + environment=stochastic_environment, + rocket=stochastic_calisto_numpy_only, + flight=stochastic_flight, + ) + with open(montecarlo.input_file, "w", encoding="utf-8") as existing: + existing.write("previous results\n") + + with pytest.raises(error): + montecarlo.simulate(random_seed=11, **kwargs) + + with open(montecarlo.input_file, encoding="utf-8") as kept: + assert kept.read() == "previous results\n" + + +def test_serial_inputs_are_reproducible( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """Two serial runs with the same seed yield byte-identical inputs per index. + + This drives the serial ``simulate`` path end to end; the flexible seed types + are covered by the unit test of ``__root_seed_sequence``. + """ + models = (stochastic_environment, stochastic_calisto_numpy_only, stochastic_flight) + run_a = _simulate_inputs( + monkeypatch, tmp_path, *models, "a", number_of_simulations=3, random_seed=7 + ) + run_b = _simulate_inputs( + monkeypatch, tmp_path, *models, "b", number_of_simulations=3, random_seed=7 + ) + assert sorted(run_a) == list(range(3)) + assert run_a == run_b + + +@pytest.mark.slow +def test_inputs_are_worker_invariant( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """serial == parallel(2) == parallel(4): inputs are bit-identical per index.""" + multiprocess = pytest.importorskip("multiprocess") + if multiprocess.get_start_method() != "fork": + pytest.skip( + "stub-based parallel determinism test requires the 'fork' start method" + ) + + models = (stochastic_environment, stochastic_calisto_numpy_only, stochastic_flight) + common = {"number_of_simulations": 8, "random_seed": 314159} + + serial = _simulate_inputs(monkeypatch, tmp_path, *models, "serial", **common) + par2 = _simulate_inputs( + monkeypatch, tmp_path, *models, "par2", parallel=True, n_workers=2, **common + ) + par4 = _simulate_inputs( + monkeypatch, tmp_path, *models, "par4", parallel=True, n_workers=4, **common + ) + + expected = list(range(8)) + assert sorted(serial) == expected + assert sorted(par2) == expected + assert sorted(par4) == expected + for index in expected: + assert serial[index] == par2[index], f"serial vs parallel(2) differ at {index}" + assert serial[index] == par4[index], f"serial vs parallel(4) differ at {index}" + + +@pytest.mark.parametrize("start_method", _available_start_methods()) +def test_seed_derivation_is_start_method_invariant(start_method): + """Per-index seeds derived in a worker match the main process under every + available start method (fork, spawn, forkserver). + + The full worker-invariance test above stubs the module-level ``Flight`` and so + only reaches workers under ``fork``. This one instead checks the property that + actually has to hold cross-platform -- that a simulation index maps to the same + seed no matter which process derives it -- using a top-level picklable target + and small picklable arguments, so it is valid under ``spawn``/``forkserver`` + (Python 3.14's POSIX default) without relying on any inherited parent state. + Two workers split the indices; their combined result must equal the + single-process derivation. + """ + root = np.random.SeedSequence(2718281828) + root_state = ( + root.entropy, + root.spawn_key, + root.pool_size, + root.n_children_spawned, + ) + indices = list(range(6)) + expected = _derive_index_seeds(root_state, indices) + + context = multiprocessing.get_context(start_method) + chunks = [(root_state, indices[0::2]), (root_state, indices[1::2])] + with context.Pool(2) as pool: + results = pool.starmap(_derive_index_seeds, chunks) + + combined = {} + for result in results: + combined.update(result) + assert combined == expected + assert sorted(combined) == indices + + +def _assert_the_same_environment_was_flown(runs, expected_indices, start_method): + """Every worker count flew index i with the same effective environment. + + This is the half the inputs file cannot show. It records + ``wind_velocity_x_factor``, which is the same for index i however the run + was executed even when the baseline it multiplies has drifted from one + simulation to the next. + """ + effective = { + label: _read_inputs_by_index(montecarlo.output_file) + for label, (montecarlo, _inputs) in runs.items() + } + for label, by_index in effective.items(): + assert sorted(by_index) == expected_indices, f"{label}: outputs are incomplete" + + for index in expected_indices: + reference = json.loads(effective["serial"][index]) + for key in _EFFECTIVE_ENVIRONMENT: + assert key in reference, f"{key} was not recorded" + assert reference["effective_wind_x"] != 0.0, ( + "the wind baseline is zero, so a compounding baseline cannot show" + ) + for label in ("parallel-2", "parallel-4"): + drawn = json.loads(effective[label][index]) + for key in _EFFECTIVE_ENVIRONMENT: + assert drawn[key] == reference[key], ( + f"{start_method}: {label} flew a different {key} at index " + f"{index}: {drawn[key]} against {reference[key]}" + ) + + +def _assert_the_run_is_complete(label, montecarlo, inputs, count): + """Every index written once, to both files, with nothing in the error log. + + The row counts are taken before anything is keyed by index: two workers + claiming the same index write two rows, and the second overwrites the first + in the dict, so a duplicate looks like a complete run. + """ + expected_indices = list(range(count)) + rows = _count_rows(montecarlo.input_file) + + assert sorted(inputs) == expected_indices, ( + f"{label}: indices {sorted(inputs)}, expected {expected_indices}" + ) + assert rows == count, ( + f"{label}: {rows} rows for {count} simulations, so an index was claimed " + f"more than once" + ) + assert _count_rows(montecarlo.output_file) == count, ( + f"{label}: the output rows do not match the simulations run" + ) + assert sorted(_read_inputs_by_index(montecarlo.output_file)) == expected_indices, ( + f"{label}: the outputs do not match the inputs" + ) + assert not os.path.getsize(montecarlo.error_file), ( + f"{label}: the run wrote to its error file" + ) + + +@pytest.fixture +def stochastic_environment_with_wind(example_spaceport_env): + """A stochastic environment whose wind is not zero. + + The shared ``stochastic_environment`` fixture sits on an Environment whose + ``wind_velocity_x`` is 0 at every altitude, and zero times any factor is + zero, so a baseline that compounds from one simulation to the next cannot + show up in it at all. Measured: with the baseline fix reverted, every + assertion in this file still passed. A wind that is actually blowing is + what makes the property testable. + """ + environment = Environment( + latitude=example_spaceport_env.latitude, + longitude=example_spaceport_env.longitude, + elevation=example_spaceport_env.elevation, + ) + environment.set_atmospheric_model( + type="custom_atmosphere", wind_u=12.0, wind_v=-7.0 + ) + return StochasticEnvironment( + environment=environment, + elevation=(1400, 10, "normal"), + wind_velocity_x_factor=(1.0, 0.05, "normal"), + wind_velocity_y_factor=(1.0, 0.05, "normal"), + ) + + +def _wind_x(flight): + """The wind the simulation actually flew with, not the factor drawn for it.""" + return float(flight.env.wind_velocity_x(0)) + + +def _wind_y(flight): + return float(flight.env.wind_velocity_y(0)) + + +def _elevation(flight): + return float(flight.env.elevation) + + +_EFFECTIVE_ENVIRONMENT = { + "effective_wind_x": _wind_x, + "effective_wind_y": _wind_y, + "effective_elevation": _elevation, +} + + +def _sampled_only(record): + """The recorded inputs with object identity stripped out. + + A ``Function``'s ``signature.hash`` and its serialised ``source`` encode the + object, not the value drawn for it, and an object built in another process + has a different one. Under ``fork`` they happen to agree because the child + inherits the parent's objects; under ``spawn`` and ``forkserver`` they + cannot. Measured on a real run: six fields differ across the boundary and + all six are these, while every sampled quantity matches exactly. + """ + flat = {} + + def walk(value, path=""): + if isinstance(value, dict): + for key, item in value.items(): + walk(item, f"{path}.{key}" if path else str(key)) + elif isinstance(value, list): + for position, item in enumerate(value): + walk(item, f"{path}[{position}]") + else: + flat[path] = value + + walk(record) + return { + key: value + for key, value in flat.items() + if "signature" not in key and not key.endswith(".source") + } + + +def _real_run_inputs(tmp_path, environment, rocket, flight, tag, **simulate_kwargs): + """Run a real Monte Carlo, no stub, and return the inputs keyed by index. + + Deliberately without the ``Flight`` stub. Stubbing is what confines the test + above to ``fork``: it replaces a module-level symbol in the parent, and a + ``spawn`` or ``forkserver`` child re-imports the module instead of inheriting + it. A real run has nothing that needs to cross the boundary except the + pickled MonteCarlo, which is the thing worth testing. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / tag), + environment=environment, + rocket=rocket, + flight=flight, + data_collector=_EFFECTIVE_ENVIRONMENT, + ) + montecarlo.simulate(**simulate_kwargs) + return montecarlo, _read_inputs_by_index(montecarlo.input_file) + + +@pytest.fixture +def restore_start_method(): + """Set the start method for one test and put it back afterwards.""" + multiprocess = pytest.importorskip("multiprocess") + original = multiprocess.get_start_method() + yield multiprocess + multiprocess.set_start_method(original, force=True) + + +@pytest.mark.slow +@pytest.mark.parametrize("start_method", _available_start_methods()) +def test_the_real_parallel_path_is_worker_invariant_under_every_start_method( + restore_start_method, + tmp_path, + stochastic_environment_with_wind, + stochastic_calisto_numpy_only, + stochastic_flight, + calisto_air_brakes_clamp_on, + start_method, +): + """The whole parallel path, not just the seed arithmetic. + + ``test_seed_derivation_is_start_method_invariant`` covers the derivation on + every start method, and the stubbed test above covers the real loop on + ``fork``. Neither covers ``multiprocess.Process``, ``__sim_producer``, the + manager proxies or pickling the stochastic object graph anywhere but + ``fork``, and that is what Windows, macOS and Python 3.14's POSIX default + actually run. + """ + multiprocess = restore_start_method + if start_method not in multiprocess.get_all_start_methods(): + pytest.skip(f"{start_method} is not available here") + multiprocess.set_start_method(start_method, force=True) + + # Air brakes and eccentricity are sampled by their own code paths, and each + # one was reseeded from somewhere other than the simulation index. + stochastic_calisto_numpy_only.add_air_brakes( + StochasticAirBrakes( + air_brakes=calisto_air_brakes_clamp_on.air_brakes[0], + drag_coefficient_curve_factor=(1.0, 0.1), + ), + calisto_air_brakes_clamp_on._controllers[0], + ) + stochastic_calisto_numpy_only.add_cp_eccentricity(x=(0.0, 0.001, "normal"), y=0.001) + stochastic_calisto_numpy_only.add_thrust_eccentricity( + x=(0.0, 0.001, "normal"), y=0.001 + ) + + count = 4 + common = {"number_of_simulations": count, "random_seed": 987654321} + models = ( + stochastic_environment_with_wind, + stochastic_calisto_numpy_only, + stochastic_flight, + ) + runs = { + "serial": _real_run_inputs( + tmp_path, *models, f"{start_method}-serial", **common + ), + "parallel-2": _real_run_inputs( + tmp_path, + *models, + f"{start_method}-p2", + parallel=True, + n_workers=2, + **common, + ), + "parallel-4": _real_run_inputs( + tmp_path, + *models, + f"{start_method}-p4", + parallel=True, + n_workers=4, + **common, + ), + } + + expected_indices = list(range(count)) + for label, (montecarlo, inputs) in runs.items(): + _assert_the_run_is_complete(label, montecarlo, inputs, count) + + _assert_the_same_environment_was_flown(runs, expected_indices, start_method) + + serial = runs["serial"][1] + for label in ("parallel-2", "parallel-4"): + for index in expected_indices: + expected = _sampled_only(json.loads(serial[index])) + actual = _sampled_only(json.loads(runs[label][1][index])) + + # Or stripping identity could quietly empty the comparison. + assert len(expected) > 20, f"only {len(expected)} fields left to compare" + assert sum("eccentricity" in key for key in expected) == 4, ( + "the four eccentricities are not among the compared fields" + ) + assert sum("brake" in key for key in expected) >= 1, ( + "the air brake is not among the compared fields" + ) + assert actual == expected, ( + f"{start_method}: serial and {label} differ at index {index} in " + f"{sorted(k for k in set(expected) | set(actual) if expected.get(k) != actual.get(k))}" + ) + + +@pytest.mark.parametrize("parallel", [False, True], ids=["serial", "parallel"]) +def test_a_missing_simulation_is_not_reported_as_a_successful_run( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto, + stochastic_flight, + parallel, +): + """A run that wrote fewer records than it claimed has to fail. + + Neither file shows this on its own: every row is well formed, and reading + them back keyed by index cannot tell four rows from three plus a duplicate. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / f"short-{parallel}"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + kwargs = {"parallel": True, "n_workers": 2} if parallel else {} + + # Lose the simulation from both files, the way a worker that dies between + # the claim and the writes does. Driven through ``simulate`` rather than by + # calling the check afterwards, so this also proves the check is reached. + real_inputs = montecarlo._MonteCarlo__evaluate_flight_inputs + real_outputs = montecarlo._MonteCarlo__evaluate_flight_outputs + + def drop_the_second_inputs(sim_idx): + return "" if sim_idx == 1 else real_inputs(sim_idx) + + def drop_the_second_outputs(flight, sim_idx): + return "" if sim_idx == 1 else real_outputs(flight, sim_idx) + + monkeypatch.setattr( + montecarlo, "_MonteCarlo__evaluate_flight_inputs", drop_the_second_inputs + ) + monkeypatch.setattr( + montecarlo, "_MonteCarlo__evaluate_flight_outputs", drop_the_second_outputs + ) + + with pytest.raises(RuntimeError, match="never written"): + montecarlo.simulate(number_of_simulations=2, random_seed=5150, **kwargs) + + +@pytest.mark.parametrize("parallel", [False, True], ids=["serial", "parallel"]) +def test_a_simulation_whose_outputs_went_missing_also_fails( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto, + stochastic_flight, + parallel, +): + """A worker that wrote its inputs and stopped before its outputs leaves the + two logs disagreeing. Checking each file against the expected range on its + own cannot see that: the inputs file is complete, and it is only complete + because the row it is missing is in the other file. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / f"no-output-{parallel}"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + kwargs = {"parallel": True, "n_workers": 2} if parallel else {} + real = montecarlo._MonteCarlo__evaluate_flight_outputs + + def drop_the_second(flight, sim_idx): + return "" if sim_idx == 1 else real(flight, sim_idx) + + monkeypatch.setattr( + montecarlo, "_MonteCarlo__evaluate_flight_outputs", drop_the_second + ) + + with pytest.raises(RuntimeError, match="disagree about which simulations"): + montecarlo.simulate(number_of_simulations=2, random_seed=5150, **kwargs) + + +@pytest.mark.parametrize("parallel", [False, True], ids=["serial", "parallel"]) +def test_a_row_cut_off_mid_write_is_named_as_unreadable( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto, + stochastic_flight, + parallel, +): + """A worker killed part way through a write leaves a truncated row. + + That row is the corruption this check exists to find, so it is named and + raised on. Skipping it and reporting the index as missing was a worse + answer: with every expected index present, a corrupt file passed. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / f"truncated-{parallel}"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + kwargs = {"parallel": True, "n_workers": 2} if parallel else {} + real = montecarlo._MonteCarlo__evaluate_flight_inputs + + def cut_the_second_short(sim_idx): + row = real(sim_idx) + if sim_idx != 1: + return row + half = row[: len(row) // 2] + "\n" + with pytest.raises(ValueError): + json.loads(half) # the row has to be unparseable for this to test it + return half + + monkeypatch.setattr( + montecarlo, "_MonteCarlo__evaluate_flight_inputs", cut_the_second_short + ) + + with pytest.raises(RuntimeError, match="not readable JSON"): + montecarlo.simulate(number_of_simulations=2, random_seed=5150, **kwargs) + + +def test_a_run_stopped_with_ctrl_c_keeps_what_it_saved( + monkeypatch, tmp_path, stochastic_environment, stochastic_calisto, stochastic_flight +): + """Ctrl-C is a stop, not a fault. + + The run path catches it, prints that the files are saved and returns. The + completeness check then counted the simulations that never ran and called + the run a failure, contradicting the message printed a moment earlier. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / "interrupted"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + real = montecarlo._MonteCarlo__run_single_simulation + finished = [] + + def stop_after_the_first(): + if finished: + raise KeyboardInterrupt("user pressed ctrl-c") + finished.append(1) + return real() + + monkeypatch.setattr( + montecarlo, "_MonteCarlo__run_single_simulation", stop_after_the_first + ) + + montecarlo.simulate(number_of_simulations=3, random_seed=42) + + # Short of the three asked for, so the check really was in a position to + # reject this run, and the one simulation that did finish is still there. + assert _count_rows(montecarlo.input_file) == 1 + + +def test_appending_checks_only_the_simulations_the_run_added( + tmp_path, stochastic_environment, stochastic_calisto, stochastic_flight +): + """``append=True`` leaves the earlier run's records in the same files, and + ``number_of_simulations`` is the total to reach rather than a count to add. + The check has to look at indices ``_initial_sim_idx`` upwards, or a second + run would be judged against records it never wrote. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / "appended"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + montecarlo.simulate(number_of_simulations=2, random_seed=606) + + assert _count_rows(montecarlo.input_file) == 2 + + # Take the first run's records away before appending. The second run is + # judged on what it wrote, so a check counting the whole file would call + # this incomplete even though nothing went wrong. + montecarlo.input_file.write_text("", encoding="utf-8") + montecarlo.output_file.write_text("", encoding="utf-8") + + montecarlo.simulate(number_of_simulations=4, append=True, random_seed=606) + + assert montecarlo._initial_sim_idx == 2, ( + "the second run should have started where the first stopped" + ) + written = _read_inputs_by_index(montecarlo.input_file) + assert sorted(written) == [2, 3], ( + f"the appended run wrote the wrong indices: {sorted(written)}" + ) diff --git a/tests/unit/simulation/test_monte_carlo_determinism.py b/tests/unit/simulation/test_monte_carlo_determinism.py new file mode 100644 index 000000000..546b2c1d1 --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_determinism.py @@ -0,0 +1,328 @@ +"""Unit tests for the Monte Carlo seeding helpers. + +``MonteCarlo.simulate(random_seed=...)`` makes the sampled inputs reproducible by +turning the run's root seed into one independent child stream per simulation +index. Four small helpers do the work: + +* ``__root_seed_sequence`` normalizes the ``random_seed`` argument (an int, a + sequence of ints, a ``SeedSequence`` or None) into a fresh ``SeedSequence``; +* ``__child_seed`` derives the child seed for one simulation index in O(1) by + extending the captured root ``spawn_key`` -- bit-identical to + ``root.spawn(n)[index]`` but without materializing the whole spawned list, so + a worker can rebuild any index from the small root state alone; +* ``_seed_sequence_to_int`` collapses a child into a 128-bit ``int`` (the seed + type a documented ``CustomSampler.reset_seed`` accepts); +* ``__seed_simulation`` splits one per-index child seed three ways so the + environment, rocket and flight draw from independent streams. + +These tests exercise the helpers directly, with no fixtures and no simulation, so +they stay fast. The end-to-end reproducibility of ``simulate`` (serial and across +workers) is covered by ``tests/integration/simulation/test_monte_carlo_determinism``. + +Reaching a name-mangled member is an established pattern in this suite (see +``tests/unit/test_sensitivity.py`` and ``tests/unit/environment/test_environment.py``); +it lets the seeding invariants be asserted without running a Monte Carlo. +""" + +import random as stdlib_random +import sys +import threading +import time +from types import SimpleNamespace + +import numpy as np +import pytest + +from rocketpy.simulation import MonteCarlo +from rocketpy.simulation.monte_carlo import ( + _SimMonitor, + _claim_next_index, + _seed_sequence_to_int, +) + +_root_seed_sequence = MonteCarlo._MonteCarlo__root_seed_sequence +_child_seed = MonteCarlo._MonteCarlo__child_seed +_seed_simulation = MonteCarlo._MonteCarlo__seed_simulation + + +def _entropy(seed_sequence, n=4): + """A stable, comparable fingerprint of a ``SeedSequence``'s stream.""" + return tuple(int(x) for x in seed_sequence.generate_state(n)) + + +def _plan(root): + """A stand-in ``self`` carrying only the root state ``__child_seed`` reads.""" + return SimpleNamespace( + _MonteCarlo__root_state=( + root.entropy, + root.spawn_key, + root.pool_size, + root.n_children_spawned, + ) + ) + + +def _advanced_root(seed, already_spawned): + """A root whose own child counter has advanced (n_children_spawned != 0), + the state a user's already-spawned SeedSequence would arrive in.""" + root = np.random.SeedSequence(seed) + root.spawn(already_spawned) + return root + + +# --------------------------------------------------------------------------- # +# __root_seed_sequence: normalizing the flexible seed argument # +# --------------------------------------------------------------------------- # + + +@pytest.mark.parametrize( + "make_seed", + [ + pytest.param(lambda: 12345, id="int"), + pytest.param(lambda: np.int64(12345), id="numpy-int"), + pytest.param(lambda: [1, 2, 3], id="sequence"), + pytest.param(lambda: np.random.SeedSequence(12345), id="seedsequence"), + ], +) +def test_root_seed_sequence_accepts_seed_like_values(make_seed): + """An int, a numpy integer, a sequence of ints and a SeedSequence are all + accepted, normalize to a SeedSequence, and are reproducible.""" + root = _root_seed_sequence(make_seed()) + assert isinstance(root, np.random.SeedSequence) + assert _entropy(root) == _entropy(_root_seed_sequence(make_seed())) + + +def test_root_seed_sequence_none_draws_fresh_entropy(): + """None yields a SeedSequence seeded from fresh OS entropy (not reproducible).""" + root = _root_seed_sequence(None) + assert isinstance(root, np.random.SeedSequence) + assert root.entropy is not None + + +@pytest.mark.parametrize( + "make_generator", + [ + pytest.param(lambda: np.random.default_rng(999), id="generator"), + pytest.param(lambda: np.random.PCG64(999), id="bitgenerator"), + ], +) +def test_root_seed_sequence_rejects_stateful_generators(make_generator): + """A Generator/BitGenerator is a stateful RNG, not a seed value, so it is + rejected instead of being reduced to its underlying SeedSequence.""" + with pytest.raises(TypeError, match="SeedSequence"): + _root_seed_sequence(make_generator()) + + +def test_root_seed_sequence_copies_full_state_without_mutating_caller(): + """A supplied SeedSequence is copied from its FULL state -- entropy, spawn_key, + pool_size and n_children_spawned -- not just its entropy, and the caller object + is not mutated. Asserting on ``.state`` is what gives this teeth: an + entropy-only copy would silently drop spawn_key/n_children_spawned (making a + spawned-child seed collide with its parent) and fail the state comparison.""" + source = np.random.SeedSequence(2024).spawn(3)[2] # non-empty spawn_key + source.spawn(5) # advance its own child counter, so it is not 0 + assert source.spawn_key == (2,) + assert source.n_children_spawned == 5 + + state_before = dict(source.state) + clone = _root_seed_sequence(source) + + assert clone is not source, "must return a copy, not the caller" + assert clone.state == state_before, "copy must preserve the full seed state" + assert source.state == state_before, "caller must not be mutated" + # The copy reproduces exactly what an independent full-state rebuild produces. + rebuilt = np.random.SeedSequence(**state_before) + assert [_entropy(c) for c in clone.spawn(3)] == [ + _entropy(c) for c in rebuilt.spawn(3) + ] + + +# --------------------------------------------------------------------------- # +# __child_seed: O(1) per-index derivation, bit-identical to spawn(n)[index] # +# --------------------------------------------------------------------------- # + + +@pytest.mark.parametrize( + "make_root", + [ + pytest.param(lambda: np.random.SeedSequence(2024), id="int-root"), + pytest.param(lambda: np.random.SeedSequence([7, 8, 9]), id="sequence-root"), + pytest.param( + lambda: np.random.SeedSequence(2024).spawn(3)[2], id="spawned-root" + ), + pytest.param(lambda: _advanced_root(99, 4), id="advanced-counter-root"), + ], +) +def test_child_seed_matches_spawn_bit_for_bit(make_root): + """Deriving index i by extending the root spawn_key equals ``root.spawn(n)[i]`` + exactly, so the O(1) derivation changes no sampled inputs versus a full spawn. + A fresh, identical root is built on each side so neither run mutates the other. + The ``advanced-counter`` root (n_children_spawned != 0) covers a user passing a + SeedSequence they have already spawned from: the base offset must equal + n_children_spawned or the derived index would collide with those children. + """ + n = 6 + derived = [_entropy(_child_seed(_plan(make_root()), i)) for i in range(n)] + spawned = [_entropy(child) for child in make_root().spawn(n)] + assert derived == spawned + + +def test_child_seed_is_worker_order_independent(): + """Any index maps to the same child regardless of the order indices are asked + for -- the property that makes a run invariant to worker scheduling.""" + plan = _plan(np.random.SeedSequence(2024)) + forward = {i: _entropy(_child_seed(plan, i)) for i in range(5)} + backward = {i: _entropy(_child_seed(plan, i)) for i in reversed(range(5))} + assert forward == backward + + +def test_child_seed_supports_indices_beyond_32_bits(): + """A simulation index past 2**32 is not truncated: it derives a distinct child + from its neighbour and matches the direct spawn_key construction for it.""" + root = np.random.SeedSequence(11) + plan = _plan(root) + big = 2**32 + 5 + assert _entropy(_child_seed(plan, big)) != _entropy(_child_seed(plan, big + 1)) + expected = np.random.SeedSequence( + entropy=root.entropy, spawn_key=(big,), pool_size=root.pool_size + ) + assert _entropy(_child_seed(plan, big)) == _entropy(expected) + + +# --------------------------------------------------------------------------- # +# _seed_sequence_to_int: 128-bit int seed for the samplers # +# --------------------------------------------------------------------------- # + + +def test_seed_sequence_to_int_is_deterministic_128_bit_int(): + def child(): + return np.random.SeedSequence(42).spawn(1)[0] + + seed = _seed_sequence_to_int(child()) + assert isinstance(seed, int) + assert 0 <= seed < 2**128 + assert seed == _seed_sequence_to_int(child()), "must be deterministic" + + +def test_seed_sequence_to_int_uses_all_128_bits(): + """The int combines all four uint32 words, not a single 32-bit word, so it + keeps the full entropy pool rather than collapsing collision risk to n**2 / + 2**32. A single-word reduction would compare unequal here.""" + ss = np.random.SeedSequence(42).spawn(1)[0] + one_word = int(np.random.SeedSequence(42).spawn(1)[0].generate_state(1)[0]) + assert _seed_sequence_to_int(ss) != one_word + assert _seed_sequence_to_int(ss).bit_length() > 32 + + +def test_seed_int_is_accepted_by_the_modern_rng_apis(): + """The 128-bit int a sampler receives works with random.Random and + numpy.random.default_rng -- the paths a CustomSampler uses. Passing a + SeedSequence there instead is unsafe: from Python 3.11 random.Random rejects + it with a TypeError, and before 3.11 it is silently hashed rather than used as + entropy. Either way an int is the right thing to hand a sampler.""" + seed = _seed_sequence_to_int(np.random.SeedSequence(1).spawn(1)[0]) + assert isinstance(stdlib_random.Random(seed).random(), float) + assert np.random.default_rng(seed).random() is not None + if sys.version_info >= (3, 11): + with pytest.raises(TypeError): + stdlib_random.Random(np.random.SeedSequence(1)) + + +# --------------------------------------------------------------------------- # +# __seed_simulation: splitting one child seed across the three models # +# --------------------------------------------------------------------------- # + + +class _RecordingModel: + """Stand-in stochastic model that records the seeds it is handed.""" + + def __init__(self): + self.seeds = [] + + def _set_stochastic(self, seed=None): + self.seeds.append(seed) + + +def _split_seeds(child_seed): + """Run ``__seed_simulation`` against recording models; return the three seeds.""" + models = SimpleNamespace( + environment=_RecordingModel(), + rocket=_RecordingModel(), + flight=_RecordingModel(), + ) + _seed_simulation(models, child_seed) + return models.environment.seeds, models.rocket.seeds, models.flight.seeds + + +def test_seed_simulation_hands_each_model_a_distinct_128_bit_int(): + """The per-index child seed is split three ways, and each model receives a + plain 128-bit int (not a SeedSequence) from an independent stream.""" + env_seeds, rocket_seeds, flight_seeds = _split_seeds(np.random.SeedSequence(2024)) + assert [len(env_seeds), len(rocket_seeds), len(flight_seeds)] == [1, 1, 1] + seeds = [env_seeds[0], rocket_seeds[0], flight_seeds[0]] + assert all(isinstance(s, int) and 0 <= s < 2**128 for s in seeds) + assert len(set(seeds)) == 3, "env/rocket/flight must be decorrelated" + + +def test_seed_simulation_is_deterministic_per_child(): + """A given child seed reseeds the three models identically every time.""" + + def split(child): + env, rocket, flight = _split_seeds(child) + return [env[0], rocket[0], flight[0]] + + assert split(np.random.SeedSequence(2024)) == split(np.random.SeedSequence(2024)) + + +# --------------------------------------------------------------------------- # +# _claim_next_index: atomic hand-out of the next simulation index # +# --------------------------------------------------------------------------- # + + +def test_claim_next_index_hands_out_each_index_once_under_contention(): + """Holding the mutex across keep_simulating() and increment() must hand out + each index exactly once, even when every worker reaches the claim together. + + A barrier releases all workers at once and a widened check-to-increment + window would let an unlocked claim run several workers past the count < n + check before any increments; the lock is what keeps the result to exactly + n_simulations indices (0..n-1, none repeated) and the counter from + overshooting. + """ + n_simulations = 5 + n_workers = 8 + monitor = _SimMonitor(initial_count=0, n_simulations=n_simulations, start_time=0.0) + + # Widen the window between the check and the increment so that, without the + # lock, several workers could pass count < n before any of them increments. + real_keep_simulating = monitor.keep_simulating + + def slow_keep_simulating(): + result = real_keep_simulating() + time.sleep(0.02) + return result + + monitor.keep_simulating = slow_keep_simulating + + mutex = threading.Lock() + barrier = threading.Barrier(n_workers) + claimed = [] + claimed_lock = threading.Lock() + + def worker(): + barrier.wait() + while True: + index = _claim_next_index(monitor, mutex) + if index is None: + break + with claimed_lock: + claimed.append(index) + + workers = [threading.Thread(target=worker) for _ in range(n_workers)] + for thread in workers: + thread.start() + for thread in workers: + thread.join() + + assert sorted(claimed) == list(range(n_simulations)) + assert monitor.count == n_simulations diff --git a/tests/unit/simulation/test_monte_carlo_log_integrity.py b/tests/unit/simulation/test_monte_carlo_log_integrity.py new file mode 100644 index 000000000..80ce38c23 --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_log_integrity.py @@ -0,0 +1,401 @@ +"""What the run is allowed to call a success, and how it comes down when it is not. + +The completeness check reads the two logs back and decides whether the run can +be reported as complete. Everything it accepts is a claim about the results, so +a row it cannot read, an index it cannot trust, or a file that disagrees with +its pair has to stop the run rather than be skipped past. + +The shutdown tests cover the other half: a fleet where one worker is not coming +back has to be brought down in bounded time, and the failure that started it has +to survive that. +""" + +import threading +import types + +import pytest + +import rocketpy.simulation.monte_carlo as mc + + +def _runner(tmp_path, rows, outputs=None, count=2, initial=0, interrupted=False): + """A stand-in carrying only what the completeness check reads.""" + inputs_file = tmp_path / "inputs.txt" + outputs_file = tmp_path / "outputs.txt" + inputs_file.write_text(rows, encoding="utf-8") + outputs_file.write_text(rows if outputs is None else outputs, encoding="utf-8") + return types.SimpleNamespace( + input_file=inputs_file, + output_file=outputs_file, + number_of_simulations=count, + _initial_sim_idx=initial, + _interrupted=interrupted, + ) + + +def _check(runner): + mc.MonteCarlo._MonteCarlo__check_each_index_was_recorded_once(runner) + + +COMPLETE = '{"index": 0}\n{"index": 1}\n' + + +CORRUPT = { + "a row cut off mid-write": (COMPLETE + "{not json\n", "not readable JSON"), + "a row that is not an object": (COMPLETE + "[]\n", "does not carry"), + "a row with no index": (COMPLETE + '{"foo": 1}\n', "does not carry"), + "a boolean index": ('{"index": 0}\n{"index": true}\n', "does not carry"), + "a float index": ('{"index": 0}\n{"index": 1.0}\n', "does not carry"), + "a negative index": (COMPLETE + '{"index": -1}\n', "does not carry"), + "an index past the run": (COMPLETE + '{"index": 99}\n', "outside the range"), + "the same index twice": (COMPLETE + '{"index": 1}\n', "more than once"), + "an index never written": ('{"index": 0}\n', "never written"), +} + + +@pytest.mark.parametrize( + ("rows", "expected"), list(CORRUPT.values()), ids=list(CORRUPT) +) +def test_a_corrupt_log_is_not_a_successful_run(tmp_path, rows, expected): + """Every one of these was accepted as a complete run. + + ``True`` and ``1.0`` are the two that look harmless: both compare equal to + ``1``, so an ``isinstance`` check or a bare dict lookup counts them as the + index they are not. Hence ``type(index) is int``. + """ + with pytest.raises(RuntimeError, match=expected): + _check(_runner(tmp_path, rows)) + + +def test_a_complete_run_is_still_accepted(tmp_path): + """The control. Without this the table above passes on a check that + rejects everything.""" + _check(_runner(tmp_path, COMPLETE)) + + +def test_an_earlier_run_left_in_the_files_is_not_an_error(tmp_path): + """``append=True`` keeps the earlier run's rows, and they are below + ``_initial_sim_idx``. Rejecting anything outside the new range would make + every appended run fail.""" + rows = '{"index": 0}\n{"index": 1}\n{"index": 2}\n{"index": 3}\n' + _check(_runner(tmp_path, rows, count=4, initial=2)) + + +def test_an_interrupted_run_still_has_to_be_readable(tmp_path): + """Being short is allowed after Ctrl-C. Being corrupt is not: skipping the + check entirely meant a duplicate or an unreadable row went unreported.""" + _check(_runner(tmp_path, '{"index": 0}\n', interrupted=True)) + + with pytest.raises(RuntimeError, match="more than once"): + _check(_runner(tmp_path, '{"index": 0}\n{"index": 0}\n', interrupted=True)) + + +def test_the_two_files_have_to_agree(tmp_path): + """A worker that wrote its inputs and stopped before its outputs. Each file + on its own can look complete, because the row one is missing is in the + other.""" + with pytest.raises(RuntimeError, match="disagree about which simulations"): + _check(_runner(tmp_path, COMPLETE, outputs='{"index": 0}\n')) + + +class _Worker: + """A process stub that can be told to ignore termination. + + Every call is appended to a shared ``trace`` so the order across the whole + fleet can be asserted, not just the per-worker counts. A stub join costs no + time, so a test that only counts calls cannot tell "signal everyone, then + wait" from "signal one and wait for it before reaching the next". + """ + + def __init__(self, name="worker", alive=False, deaf=False, exitcode=0, trace=None): + self.name = name + self._alive = alive + self._deaf = deaf + self.exitcode = exitcode + self.terminated = 0 + self.killed = 0 + self.joins = [] + self.trace = [] if trace is None else trace + + def is_alive(self): + return self._alive + + def terminate(self): + self.terminated += 1 + self.trace.append(("terminate", self.name)) + if not self._deaf: + self._alive = False + + def kill(self): + self.killed += 1 + self.trace.append(("kill", self.name)) + self._alive = False + + def join(self, timeout=None): + self.joins.append(timeout) + self.trace.append(("join", self.name)) + + +class _Event: + def __init__(self, flag=False): + self.flag = flag + + def set(self): + self.flag = True + + def is_set(self): + return self.flag + + +def test_the_wait_gives_up_as_soon_as_a_worker_reports_an_error(): + """One worker is not coming back and another has already failed. + + Joining the fleet in order held the parent on the first worker forever, so + the error the second had already reported was never seen and the cleanup + after it never ran. + """ + stuck, failed = _Worker(alive=True), _Worker(exitcode=1) + + mc._wait_for_workers([stuck, failed], _Event(flag=True)) + + assert stuck.is_alive(), "the wait should return, not stop the workers itself" + + +def test_the_wait_is_bounded_when_it_is_given_a_deadline(): + """The interrupt path waits a short while for workers to notice the event. + Without a deadline that wait was the second place Ctrl-C could hang.""" + stuck = _Worker(alive=True) + + mc._wait_for_workers([stuck], timeout=0.2) + + assert stuck.is_alive() + + +def test_the_wait_reaps_workers_that_had_already_finished(): + """A worker gone before the loop starts is never joined by it, and an + unjoined child has no exit code, which the crash check downstream reads as + a crash.""" + done = _Worker(alive=False) + + mc._wait_for_workers([done], _Event()) + + assert done.joins, "a finished worker was never joined, so it was not reaped" + + +def test_every_worker_is_signalled_before_any_of_them_is_waited_on(): + """Signal the whole fleet, then wait on it. + + Terminating one and joining it before reaching the next made every worker + wait out the grace period of the ones ahead of it in the list, so a single + worker that ignores the signal delays the rest by that much each. + """ + trace = [] + deaf = _Worker(name="deaf", alive=True, deaf=True, trace=trace) + ordinary = _Worker(name="ordinary", alive=True, trace=trace) + + mc._stop_any_worker_still_running([deaf, ordinary], grace=0.01) + + first_join = next(i for i, (call, _) in enumerate(trace) if call == "join") + assert not [c for c in trace[first_join:] if c[0] == "terminate"], ( + f"a worker was signalled only after another had been waited on: {trace}" + ) + assert ordinary.terminated == 1, "the second worker was never signalled" + assert deaf.killed == 1, "the worker that sat through terminate was not killed" + assert not deaf.is_alive() + + +def test_a_worker_that_already_exited_is_left_alone(): + """The control: cleanup runs on every path, including the ones where + nothing went wrong.""" + done = _Worker(alive=False) + + mc._stop_any_worker_still_running([done], grace=0.01) + + assert done.terminated == 0 and done.killed == 0 + + +class _RecordingMutex: + def __init__(self, fail_on_acquire=False): + self.acquired = 0 + self.released = 0 + self._fail = fail_on_acquire + self._lock = threading.Lock() + + def acquire(self, *args, **kwargs): + self.acquired += 1 + if self._fail: + raise ConnectionResetError("the manager is gone") + return self._lock.acquire(*args, **kwargs) + + def release(self): + self.released += 1 + return self._lock.release() + + +class _Boom(RuntimeError): + """The failure under test, so it cannot be confused with an incidental one.""" + + +class _Monitor: + """Enough of a monitor for a worker that completes an iteration.""" + + count = 0 + + def print_update_status(self): + pass + + +def _sim_worker(tmp_path, **overrides): + attributes = { + "error_file": tmp_path / "errors.txt", + "input_file": tmp_path / "inputs.txt", + "output_file": tmp_path / "outputs.txt", + "_MonteCarlo__child_seed": lambda index: index, + "_MonteCarlo__seed_simulation": lambda seed: None, + "_MonteCarlo__run_single_simulation": object, + "_MonteCarlo__evaluate_flight_inputs": lambda index: '{"index": 0}\n', + "_MonteCarlo__evaluate_flight_outputs": lambda flight, index: '{"index": 0}\n', + } + attributes.update(overrides) + return types.SimpleNamespace(**attributes) + + +def test_a_claim_that_fails_after_a_completed_run_does_not_report_that_run( + tmp_path, monkeypatch +): + """The state was cleared after the claim rather than before it. + + So a claim that failed on the second lap reached the handler still holding + the row that had just been written successfully, and the error file got a + second copy of a simulation that never failed. + """ + claims = iter([0]) + + def claim_once_then_fail(*_args, **_kwargs): + try: + return next(claims) + except StopIteration: + raise _Boom("the claim failed on the second lap") from None + + monkeypatch.setattr(mc, "_claim_next_index", claim_once_then_fail) + reported = [] + monkeypatch.setattr(mc._SimMonitor, "reprint", staticmethod(reported.append)) + worker = _sim_worker(tmp_path) + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__sim_producer( + worker, _Monitor(), _RecordingMutex(), _Event() + ) + + written = (tmp_path / "errors.txt").read_text() + assert '"index": 0' not in written, ( + f"the completed simulation was written to the error file again: {written!r}" + ) + assert "the claim failed on the second lap" in written + + +def test_a_mutex_that_cannot_be_taken_is_not_then_released(tmp_path, monkeypatch): + """The normal write path released in ``finally`` whether or not it had the + lock, so a manager that died during acquire raised a second error on the way + out and buried the first.""" + monkeypatch.setattr(mc, "_claim_next_index", lambda *a, **k: 0) + mutex = _RecordingMutex(fail_on_acquire=True) + + with pytest.raises(ConnectionResetError): + mc.MonteCarlo._MonteCarlo__sim_producer( + _sim_worker(tmp_path), _Monitor(), mutex, _Event() + ) + + assert mutex.released == 0, "released a lock it never held" + + +def _serial_runner(tmp_path, row=""): + """A stand-in carrying only what ``__run_in_serial`` touches.""" + runner = types.SimpleNamespace( + _initial_sim_idx=0, + number_of_simulations=2, + _interrupted=False, + _error_file=tmp_path / "errors.txt", + input_file=tmp_path / "inputs.txt", + output_file=tmp_path / "outputs.txt", + _MonteCarlo__child_seed=lambda index: index, + _MonteCarlo__seed_simulation=lambda seed: None, + _MonteCarlo__run_single_simulation=object, + _MonteCarlo__evaluate_flight_inputs=lambda index: row, + _MonteCarlo__evaluate_flight_outputs=lambda flight, index: row, + ) + runner._MonteCarlo__keep_the_inputs_that_did_not_finish = lambda payload: ( + mc.MonteCarlo._MonteCarlo__keep_the_inputs_that_did_not_finish(runner, payload) + ) + return runner + + +def test_ctrl_c_before_the_first_row_keeps_the_interrupt(tmp_path, monkeypatch): + """Half of the fix: the payload is bound before the try. + + It was assigned inside the loop body, after the two monitor calls, so Ctrl-C + in either of those reached the handler with it still unbound and the + interrupt came out as an UnboundLocalError instead. + """ + + class _Monitor: + count = 0 + + def __init__(self, **_kwargs): + pass + + def keep_simulating(self): + raise KeyboardInterrupt("ctrl-c before the first simulation") + + monkeypatch.setattr(mc, "_SimMonitor", _Monitor) + runner = _serial_runner(tmp_path) + + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + assert runner._interrupted, "the interrupt was not recorded" + + +def test_ctrl_c_between_rows_does_not_report_the_row_that_succeeded( + tmp_path, monkeypatch +): + """The other half: the payload is cleared before each lap, not after. + + Binding it once before the try stops the UnboundLocalError but leaves it + holding the last completed row, so an interrupt between two iterations + wrote a simulation that had succeeded into the error file. Both halves are + needed, and each one passes the other's test on its own. + """ + + class _Monitor: + count = 0 + + def __init__(self, **_kwargs): + self.laps = 0 + + def keep_simulating(self): + self.laps += 1 + if self.laps > 1: + raise KeyboardInterrupt("ctrl-c after the first simulation") + return True + + def increment(self): + return 1 + + def print_update_status(self): + pass + + def print_final_status(self): + pass + + monkeypatch.setattr(mc, "_SimMonitor", _Monitor) + runner = _serial_runner(tmp_path, row='{"index": 0}\n') + + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + assert runner._interrupted + assert (tmp_path / "inputs.txt").read_text() == '{"index": 0}\n', ( + "the simulation that completed should have been written normally" + ) + assert (tmp_path / "errors.txt").read_text() == "", ( + "a completed simulation was written to the error file as if it failed" + ) diff --git a/tests/unit/simulation/test_monte_carlo_worker_exit.py b/tests/unit/simulation/test_monte_carlo_worker_exit.py new file mode 100644 index 000000000..52cd0bce9 --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_worker_exit.py @@ -0,0 +1,179 @@ +"""The parent has to notice a worker that died without saying so. + +``join()`` returns None however the child ended, so the shared error event was +the only signal the run had. A worker can leave without setting it: ``SystemExit``, +``os._exit``, a segfault in a native extension, a target that will not unpickle +under spawn, or the error handler of the worker itself failing. The exit status +is what separates those from a clean finish. +""" + +import types +from contextlib import contextmanager + +import pytest + +import rocketpy.simulation.monte_carlo as mc + + +class _Process: + """A worker that does not run, and reports the exit code it was given.""" + + instances = [] + + def __init__(self, target=None, args=(), **_kwargs): # pylint: disable=unused-argument + self.name = f"worker-{len(self.instances)}" + self.exitcode = None + self.started = False + self.terminated = False + self._planned_exitcode = 0 + self.instances.append(self) + + def start(self): + self.started = True + + def join(self, *_a, **_k): + self.exitcode = self._planned_exitcode + + def is_alive(self): + return False + + def terminate(self): + self.terminated = True + + +class _Event: + def __init__(self): + self.flag = False + + def set(self): + self.flag = True + + def is_set(self): + return self.flag + + +class _Monitor: + def __init__(self, **_kwargs): + pass + + def print_final_status(self): + pass + + +@pytest.fixture +def parallel_runner(monkeypatch, tmp_path): + """Run ``__run_in_parallel`` over stub workers and hand back the stubs.""" + _Process.instances = [] + fake_multiprocess = types.SimpleNamespace(Process=_Process) + + class _Manager: # pylint: disable=invalid-name + """Method names mirror the multiprocess manager API.""" + + def Lock(self): # noqa: N802 + return types.SimpleNamespace(acquire=lambda: None, release=lambda: None) + + def Event(self): # noqa: N802 + return _Event() + + def _SimMonitor(self, **kwargs): # noqa: N802 + return _Monitor(**kwargs) + + @contextmanager + def fake_manager(*_a, **_k): + yield _Manager() + + monkeypatch.setattr(mc, "_import_multiprocess", lambda: (fake_multiprocess, None)) + monkeypatch.setattr(mc, "_create_multiprocess_manager", fake_manager) + + runner = types.SimpleNamespace( + error_file=tmp_path / "errors.txt", + input_file=tmp_path / "inputs.txt", + output_file=tmp_path / "outputs.txt", + _initial_sim_idx=0, + number_of_simulations=4, + _interrupted=False, + _MonteCarlo__validate_number_of_workers=lambda n: 2, + _MonteCarlo__sim_producer=lambda *a: None, + ) + runner.input_file.write_text("") + runner.output_file.write_text("") + return runner + + +def test_a_worker_that_crashes_without_setting_the_event_fails_the_run( + parallel_runner, +): + """The case the event alone cannot see.""" + original_join = _Process.join + + def crash(self, *a, **k): + original_join(self, *a, **k) + self.exitcode = -11 # SIGSEGV + + _Process.join = crash + try: + with pytest.raises(RuntimeError, match="did not exit cleanly"): + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=2) + finally: + _Process.join = original_join + + +def test_a_clean_run_is_not_reported_as_a_crash(parallel_runner): + """The other half: every worker exits 0, so nothing is raised.""" + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=2) + + assert all(p.exitcode == 0 for p in _Process.instances) + + +def test_a_failed_start_still_cleans_up_the_workers_already_running( + parallel_runner, monkeypatch +): + """The start loop is inside the try for this. A ``start()`` that fails part + way through used to leave the ones already running with nobody to reap + them.""" + started = [] + original_start = _Process.start + + def start_then_fail(self): + if len(started) >= 1: + raise OSError("cannot allocate a process") + original_start(self) + started.append(self) + + monkeypatch.setattr(_Process, "start", start_then_fail) + monkeypatch.setattr(_Process, "is_alive", lambda self: not self.terminated) + + with pytest.raises(OSError): + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=3) + + assert started, "the fixture never started anything" + assert all(p.terminated for p in started), "a started worker was left running" + + +def test_an_interrupted_run_is_not_then_reported_as_incomplete( + parallel_runner, monkeypatch +): + """Ctrl-C in the parent is caught and deliberately not re-raised, so + ``simulate`` carries on to the completeness check with the run unfinished. + The two are composed here in that order, since the check has to be able to + tell a run the user stopped from a worker that went missing. + """ + interrupted = [] + original_join = _Process.join + + def ctrl_c(self, *a, **k): + # Once: the handler joins again on its way out, and that has to work. + if not interrupted: + interrupted.append(True) + raise KeyboardInterrupt("user pressed ctrl-c") + original_join(self, *a, **k) + + monkeypatch.setattr(_Process, "join", ctrl_c) + + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=2) + mc.MonteCarlo._MonteCarlo__check_each_index_was_recorded_once(parallel_runner) + + assert interrupted, "the run was never interrupted" + assert parallel_runner.input_file.read_text() == "", ( + "nothing was written, so the check really was in a position to reject this" + ) diff --git a/tests/unit/simulation/test_monte_carlo_worker_failures.py b/tests/unit/simulation/test_monte_carlo_worker_failures.py new file mode 100644 index 000000000..4646356a9 --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_worker_failures.py @@ -0,0 +1,256 @@ +"""What a worker does when something fails part way through. + +A worker that dies has to leave three things true: the failure that started it +is the one that escapes, the shared mutex is not still held, and the error event +is set. Before this, a failure in the claim itself broke all three at once -- +``sim_idx`` and ``inputs_json`` were only bound inside the loop, so the handler +raised ``UnboundLocalError`` over the real error while holding the mutex, and +``error_event.set()`` sat after the write that never ran. +""" + +import json +import threading +import types + +import pytest + +import rocketpy.simulation.monte_carlo as mc + + +class _RecordingMutex: + """A real lock that counts acquires and releases.""" + + def __init__(self): + self.acquired = 0 + self.released = 0 + self._lock = threading.Lock() + + def acquire(self, *args, **kwargs): + self.acquired += 1 + return self._lock.acquire(*args, **kwargs) + + def release(self): + self.released += 1 + return self._lock.release() + + +class _Event: + def __init__(self): + self.flag = False + + def set(self): + self.flag = True + + def is_set(self): + return self.flag + + +class _Boom(RuntimeError): + """The failure under test, so it cannot be confused with an incidental one.""" + + +def _worker(tmp_path, **overrides): + """A stand-in carrying only the attributes ``__sim_producer`` touches.""" + attributes = { + "error_file": tmp_path / "errors.txt", + "input_file": tmp_path / "inputs.txt", + "output_file": tmp_path / "outputs.txt", + "_MonteCarlo__child_seed": lambda index: index, + "_MonteCarlo__seed_simulation": lambda seed: None, + "_MonteCarlo__run_single_simulation": object, + "_MonteCarlo__evaluate_flight_inputs": lambda index: "{}\n", + "_MonteCarlo__evaluate_flight_outputs": lambda flight, index: "{}\n", + } + attributes.update(overrides) + return types.SimpleNamespace(**attributes) + + +def _raise(*_args, **_kwargs): + raise _Boom("injected") + + +@pytest.mark.parametrize( + "stage", + ["claim", "reseed", "flight", "inputs", "outputs"], + ids=["claim", "reseed", "flight", "input_eval", "output_eval"], +) +def test_a_failure_anywhere_keeps_the_cause_and_frees_the_mutex( + tmp_path, monkeypatch, stage +): + """Whichever stage fails, the same three things have to hold.""" + indices = iter([0, None]) + monkeypatch.setattr(mc, "_claim_next_index", lambda *a, **k: next(indices)) + overrides = {} + if stage == "claim": + monkeypatch.setattr(mc, "_claim_next_index", _raise) + elif stage == "reseed": + overrides["_MonteCarlo__seed_simulation"] = _raise + elif stage == "flight": + overrides["_MonteCarlo__run_single_simulation"] = _raise + elif stage == "inputs": + overrides["_MonteCarlo__evaluate_flight_inputs"] = _raise + elif stage == "outputs": + overrides["_MonteCarlo__evaluate_flight_outputs"] = _raise + + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__sim_producer( + _worker(tmp_path, **overrides), object(), mutex, event + ) + + assert mutex.acquired == mutex.released, "the mutex was left held" + assert event.flag, "the parent was never told an error happened" + + +def test_the_cause_survives_even_when_the_error_report_also_fails( + tmp_path, monkeypatch +): + """Reporting is best effort. If the error file is unwritable too, the + failure that started it is still what comes out, and the mutex is still + released.""" + monkeypatch.setattr(mc, "_claim_next_index", lambda *a, **k: 0) + monkeypatch.setattr( + "builtins.open", lambda *a, **k: (_ for _ in ()).throw(OSError("no disk")) + ) + worker = _worker(tmp_path, _MonteCarlo__run_single_simulation=_raise) + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__sim_producer(worker, object(), mutex, event) + + assert mutex.acquired == mutex.released, "the mutex was left held" + assert event.flag + + +def test_the_failure_is_still_reported_when_the_claim_itself_failed( + tmp_path, monkeypatch +): + """The report has to survive a failure before the loop body ran. + + ``sim_idx`` and ``inputs_json`` are bound before the try for this reason. + Left to the loop, the handler raised ``UnboundLocalError`` at its first + write, so nothing was written and nothing was printed: the run ended with + no record of what went wrong. + """ + monkeypatch.setattr(mc, "_claim_next_index", _raise) + reported = [] + monkeypatch.setattr(mc._SimMonitor, "reprint", staticmethod(reported.append)) + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__sim_producer( + _worker(tmp_path), object(), mutex, event + ) + + assert reported, "the worker died without reporting anything" + assert "injected" in reported[0], f"the report does not name the cause: {reported}" + + +def test_an_interrupt_while_reporting_does_not_leave_the_mutex_held( + tmp_path, monkeypatch +): + """``except Exception`` does not catch ``KeyboardInterrupt``, so the release + has to be in a ``finally``. Ctrl-C between the acquire and the release would + otherwise leave every other worker blocked on it for good.""" + monkeypatch.setattr(mc, "_claim_next_index", lambda *a, **k: 0) + + def interrupt(*_a, **_k): + raise KeyboardInterrupt + + monkeypatch.setattr(mc._SimMonitor, "reprint", staticmethod(interrupt)) + worker = _worker(tmp_path, _MonteCarlo__run_single_simulation=_raise) + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(KeyboardInterrupt): + mc.MonteCarlo._MonteCarlo__sim_producer(worker, object(), mutex, event) + + assert mutex.acquired == mutex.released, "the mutex was left held on interrupt" + + +def test_a_failure_before_the_inputs_exist_still_leaves_a_readable_record( + tmp_path, monkeypatch +): + """The run tells the user to check the error file, so it has to say + something. A failure in the claim has no inputs to write, and the file was + left empty while the traceback went only to a worker's stdout, which under + ``spawn`` on Windows the user may never see. + + It has to stay a JSON line: ``_read_log_file`` parses this file with + ``json.loads`` per line, so free text would make the whole log unreadable. + """ + monkeypatch.setattr(mc, "_claim_next_index", _raise) + worker = _worker(tmp_path) + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__sim_producer(worker, object(), mutex, event) + + lines = [ + line + for line in worker.error_file.read_text(encoding="utf-8").splitlines() + if line.strip() + ] + + assert lines, "the error file was left empty" + record = json.loads(lines[0]) + assert record["index"] is None, "an early failure has no simulation index" + assert "injected" in record["error"], "the record does not carry the cause" + + +@pytest.mark.parametrize("failing_file", ["input_file", "output_file"]) +def test_a_failed_write_is_reported_like_any_other_failure( + tmp_path, monkeypatch, failing_file +): + """A disk that fills up part way through is a failure like any other: the + cause has to escape, the mutex has to come back, and the event has to be + set. These two writes sit inside the loop's own mutex block rather than the + handler, so they are worth exercising separately from the stages above. + """ + monkeypatch.setattr(mc, "_claim_next_index", lambda *a, **k: 0) + worker = _worker(tmp_path) + blocked = str(getattr(worker, failing_file)) + real_open = open + + def selective_open(path, *args, **kwargs): + if str(path) == blocked: + raise OSError("no space left on device") + return real_open(path, *args, **kwargs) + + monkeypatch.setattr("builtins.open", selective_open) + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(OSError, match="no space left"): + mc.MonteCarlo._MonteCarlo__sim_producer(worker, object(), mutex, event) + + assert mutex.acquired == mutex.released, "the mutex was left held" + assert event.flag + + +def test_an_unreachable_error_event_does_not_replace_the_failure_it_reports( + tmp_path, monkeypatch +): + """The event is a manager proxy, so notifying can fail on its own. + + It is set first and every other report is guarded, which left this one + statement able to do the thing the guards exist to prevent: raise over the + failure being reported, so the parent sees a connection error instead. + """ + + class _UnreachableEvent: + def set(self): + raise ConnectionResetError("the manager is gone") + + def is_set(self): + raise ConnectionResetError("the manager is gone") + + monkeypatch.setattr(mc, "_claim_next_index", lambda *a, **k: 0) + worker = _worker(tmp_path, _MonteCarlo__run_single_simulation=_raise) + mutex = _RecordingMutex() + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__sim_producer( + worker, object(), mutex, _UnreachableEvent() + ) + + assert mutex.acquired == mutex.released, "the mutex was left held" diff --git a/tests/unit/simulation/test_montecarlo_determinism.py b/tests/unit/simulation/test_montecarlo_determinism.py deleted file mode 100644 index 2c4629db4..000000000 --- a/tests/unit/simulation/test_montecarlo_determinism.py +++ /dev/null @@ -1,391 +0,0 @@ -"""Determinism tests for :class:`rocketpy.simulation.MonteCarlo`. - -These tests verify that the random *inputs* generated by a Monte Carlo run are -invariant to the execution mode (serial vs parallel) and to the number of -parallel workers when a fixed ``random_seed`` is provided. This is the -worker-invariance guarantee introduced by per-simulation-index seeding: the -seed of simulation ``i`` is derived from ``i`` only, never from the worker that -happens to run it. - -The trajectory integration (:class:`rocketpy.Flight`) is replaced by a -lightweight stub so the tests run fast. Worker-invariance is a property of the -*input sampling*, which happens in ``create_object``/``_randomize_*`` before the -``Flight`` object is built and is independent of the (expensive) physics -integration. Stubbing the module-level ``Flight`` symbol propagates to the -parallel workers because the parallel backend uses the ``fork`` start method, so -forked workers inherit the patched module. - -A dedicated ``stochastic_calisto_numpy_only`` rocket is used so that *all* -randomness flows through the seeded numpy generator. List-valued stochastic -attributes are sampled with the standard-library ``random.choice`` (an unseeded -*global* generator) which the per-index seeding does not control; the fixture -removes the only such attribute (a multi-element ``thrust_source``) so the -generated inputs are byte-for-byte reproducible from the seed alone. - -These tests are additive and isolated; they do not modify any inherited test. -""" - -import json - -import multiprocess -import pytest - -import rocketpy.simulation.monte_carlo as mc_module -from rocketpy.simulation import MonteCarlo -from rocketpy.stochastic import StochasticRocket, StochasticSolidMotor - -# Every test here drives a real ``MonteCarlo.simulate`` run (the ``Flight`` -# physics is stubbed, but a full Calisto rocket is still rebuilt per simulation), -# so they are gated behind ``--runslow`` to keep the default suite fast, matching -# the existing ``test_monte_carlo_simulate`` convention. -pytestmark = pytest.mark.slow - -# The stub-based parallel tests rely on workers inheriting the monkeypatched -# ``Flight`` symbol, which only happens with the ``fork`` start method. -requires_fork = pytest.mark.skipif( - multiprocess.get_start_method() != "fork", - reason="stub-based parallel determinism test requires the 'fork' start method", -) - - -class _StubFlight: - """Minimal stand-in for :class:`rocketpy.Flight` that skips integration. - - The Monte Carlo input sampling happens in - ``StochasticRocket.create_object``, ``StochasticEnvironment.create_object`` - and the ``StochasticFlight`` randomize helpers *before* the ``Flight`` - object is built. Replacing ``Flight`` with this stub therefore preserves the - exact random draws while skipping the costly trajectory integration. Any - output attribute requested by the exporter resolves to ``0.0``. - """ - - def __init__(self, **kwargs): - # Accepts the keyword arguments that ``MonteCarlo`` passes to ``Flight`` - # and intentionally ignores them. - pass - - def __getattr__(self, name): - return 0.0 - - -@pytest.fixture -def stochastic_calisto_numpy_only( - cesaroni_m1670, - calisto_robust, - stochastic_nose_cone, - stochastic_trapezoidal_fins, - stochastic_tail, - stochastic_rail_buttons, - stochastic_main_parachute, - stochastic_drogue_parachute, -): - """A ``StochasticRocket`` whose randomness flows entirely through numpy. - - It mirrors the shared ``stochastic_calisto`` fixture but gives the solid - motor a single fixed ``thrust_source`` instead of a multi-element list. - List-valued stochastic attributes are drawn with the standard-library - ``random.choice`` (an unseeded *global* generator) which per-simulation-index - seeding does not govern. Removing the only such attribute makes every - generated input reproducible from the seed alone, so serial and parallel - runs can be compared byte-for-byte. - - Returns - ------- - StochasticRocket - A stochastic Calisto rocket with no global-``random`` dependence. - """ - motor = StochasticSolidMotor( - solid_motor=cesaroni_m1670, - burn_out_time=(4, 0.1), - grains_center_of_mass_position=0.001, - grain_density=50, - grain_separation=1 / 1000, - grain_initial_height=1 / 1000, - grain_initial_inner_radius=0.375 / 1000, - grain_outer_radius=0.375 / 1000, - total_impulse=(6500, 1000), - throat_radius=0.5 / 1000, - nozzle_radius=0.5 / 1000, - nozzle_position=0.001, - ) - rocket = StochasticRocket( - rocket=calisto_robust, - radius=0.0127 / 2000, - mass=(15.426, 0.5, "normal"), - inertia_11=(6.321, 0), - inertia_22=0.01, - inertia_33=0.01, - center_of_mass_without_motor=0, - ) - rocket.add_motor(motor, position=0.001) - rocket.add_nose(stochastic_nose_cone, position=(1.134, 0.001)) - rocket.add_trapezoidal_fins(stochastic_trapezoidal_fins, position=(0.001, "normal")) - rocket.add_tail(stochastic_tail) - rocket.set_rail_buttons( - stochastic_rail_buttons, lower_button_position=(-0.618, 0.001, "normal") - ) - rocket.add_parachute(stochastic_main_parachute) - rocket.add_parachute(stochastic_drogue_parachute) - return rocket - - -def _read_inputs_by_index(input_file): - """Read a Monte Carlo ``.inputs.txt`` file into ``{index: raw_json_line}``. - - Keying by the ``index`` field makes the comparison robust to the order in - which parallel workers append their results. - - Parameters - ---------- - input_file : str or pathlib.Path - Path to the Monte Carlo inputs file. - - Returns - ------- - dict[int, str] - Mapping of simulation index to the raw (stripped) JSON line written for - that simulation. - """ - by_index = {} - with open(input_file, mode="r", encoding="utf-8") as rows: - for line in rows: - line = line.strip() - if not line: - continue - record = json.loads(line) - by_index[record["index"]] = line - return by_index - - -def _simulate_inputs( - monkeypatch, - tmp_path, - environment, - rocket, - flight, - tag, - *, - n_simulations, - random_seed, - parallel, - n_workers=None, -): - """Run a Monte Carlo simulation with a stubbed ``Flight`` and return inputs. - - Parameters - ---------- - monkeypatch : pytest.MonkeyPatch - Fixture used to swap ``Flight`` for :class:`_StubFlight`. - tmp_path : pathlib.Path - Temporary directory for the Monte Carlo output files. - environment, rocket, flight : StochasticModel - Stochastic models shared across runs. - tag : str - Unique filename stem so successive runs do not clobber each other. - n_simulations : int - Number of simulations to run. - random_seed : int or None - Seed forwarded to ``MonteCarlo.simulate``. - parallel : bool - Whether to run in parallel mode. - n_workers : int, optional - Number of workers for parallel mode. - - Returns - ------- - dict[int, str] - Mapping of simulation index to the raw input JSON line. - """ - monkeypatch.setattr(mc_module, "Flight", _StubFlight) - montecarlo = MonteCarlo( - filename=str(tmp_path / tag), - environment=environment, - rocket=rocket, - flight=flight, - ) - montecarlo.simulate( - number_of_simulations=n_simulations, - append=False, - parallel=parallel, - random_seed=random_seed, - n_workers=n_workers, - ) - return _read_inputs_by_index(montecarlo.input_file) - - -def test_montecarlo_inputs_serial_reproducible( - monkeypatch, - tmp_path, - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, -): - """Two serial runs with the same seed yield identical inputs per index.""" - models = ( - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, - ) - run_a = _simulate_inputs( - monkeypatch, - tmp_path, - *models, - "serial_a", - n_simulations=6, - random_seed=20240617, - parallel=False, - ) - run_b = _simulate_inputs( - monkeypatch, - tmp_path, - *models, - "serial_b", - n_simulations=6, - random_seed=20240617, - parallel=False, - ) - - assert sorted(run_a) == list(range(6)) - assert run_a == run_b - - -@requires_fork -def test_montecarlo_inputs_worker_invariant( - monkeypatch, - tmp_path, - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, -): - """Serial == parallel(2) == parallel(N): inputs are bit-identical per index. - - The per-iteration input JSON for a given simulation index must be byte-for- - byte identical regardless of whether the run was serial, parallel with two - workers, or parallel with a larger number of workers. - """ - n_simulations = 8 - random_seed = 314159 - models = ( - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, - ) - - serial = _simulate_inputs( - monkeypatch, - tmp_path, - *models, - "serial", - n_simulations=n_simulations, - random_seed=random_seed, - parallel=False, - ) - parallel_2 = _simulate_inputs( - monkeypatch, - tmp_path, - *models, - "parallel_2", - n_simulations=n_simulations, - random_seed=random_seed, - parallel=True, - n_workers=2, - ) - parallel_n = _simulate_inputs( - monkeypatch, - tmp_path, - *models, - "parallel_n", - n_simulations=n_simulations, - random_seed=random_seed, - parallel=True, - n_workers=4, - ) - - # Every index must be present exactly once in each run. - expected_indices = list(range(n_simulations)) - assert sorted(serial) == expected_indices - assert sorted(parallel_2) == expected_indices - assert sorted(parallel_n) == expected_indices - - # Worker-invariance: bit-identical input line per index across all modes. - for index in expected_indices: - assert serial[index] == parallel_2[index], ( - f"serial vs parallel(2) inputs differ at index {index}" - ) - assert serial[index] == parallel_n[index], ( - f"serial vs parallel(4) inputs differ at index {index}" - ) - - -def test_montecarlo_none_seed_runs( - monkeypatch, - tmp_path, - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, -): - """``random_seed=None`` stays functional (non-reproducible but complete). - - ``numpy.random.SeedSequence(None)`` pulls fresh entropy, so the run cannot - be reproduced, but it must still execute and export one record per index. - """ - inputs = _simulate_inputs( - monkeypatch, - tmp_path, - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, - "none_seed", - n_simulations=5, - random_seed=None, - parallel=False, - ) - - assert sorted(inputs) == list(range(5)) - - -def test_montecarlo_inputs_append_reproducible( - monkeypatch, - tmp_path, - stochastic_environment, - stochastic_calisto_numpy_only, - stochastic_flight, -): - """Appending to a partial run reproduces a single full run, index for index. - - ``SeedSequence.spawn`` is invariant to the spawn count, so simulation ``i`` - draws the same seed whether it was produced in the first batch or in a later - ``append`` batch. A 3 + 3 (append) run must therefore match a single run of 6 - with the same seed. ``number_of_simulations`` is the cumulative target, so the - append call passes 6, not the 3 additional simulations. - """ - monkeypatch.setattr(mc_module, "Flight", _StubFlight) - seed = 2718281 - environment = stochastic_environment - rocket = stochastic_calisto_numpy_only - flight = stochastic_flight - - split = MonteCarlo( - filename=str(tmp_path / "split"), - environment=environment, - rocket=rocket, - flight=flight, - ) - split.simulate(number_of_simulations=3, append=False, random_seed=seed) - split.simulate(number_of_simulations=6, append=True, random_seed=seed) - split_inputs = _read_inputs_by_index(split.input_file) - - single_inputs = _simulate_inputs( - monkeypatch, - tmp_path, - environment, - rocket, - flight, - "single", - n_simulations=6, - random_seed=seed, - parallel=False, - ) - - assert sorted(split_inputs) == list(range(6)) - assert split_inputs == single_inputs diff --git a/tests/unit/stochastic/test_stochastic_model.py b/tests/unit/stochastic/test_stochastic_model.py index 0d0a13311..35bc86a96 100644 --- a/tests/unit/stochastic/test_stochastic_model.py +++ b/tests/unit/stochastic/test_stochastic_model.py @@ -1,5 +1,40 @@ +from types import SimpleNamespace + import pytest +from rocketpy import Environment +from rocketpy.stochastic import StochasticEnvironment +from rocketpy.stochastic.stochastic_model import StochasticModel + + +def _sampled_option(model): + """Return the value ``dict_generator`` picks for the ``options`` attribute.""" + return next(model.dict_generator())["options"] + + +def test_list_attribute_sampling_is_reproducible_under_seed(): + """A list-valued stochastic attribute is drawn through the model's own seeded + numpy generator, so a fixed seed reproduces the choice. It used to be drawn + with the stdlib ``random.choice`` (an unseeded global instance), which + ``random_seed`` could not govern. Heterogeneous entries (paths, callables, + lists) are returned unchanged rather than coerced to a numpy dtype the way + ``numpy.random.choice`` would. + """ + options = ["/motor/a.eng", "/motor/b.eng", (lambda t: t), [1, 2, 3]] + model = StochasticModel(obj=SimpleNamespace(), options=options) + + model._set_stochastic(42) + first = _sampled_option(model) + model._set_stochastic(42) + assert _sampled_option(model) == first, "same seed must reproduce the choice" + assert any(first is option for option in options), "object returned unchanged" + + chosen_ids = set() + for seed in range(16): + model._set_stochastic(seed) + chosen_ids.add(id(_sampled_option(model))) + assert len(chosen_ids) > 1, "different seeds must be able to pick differently" + @pytest.mark.parametrize( "fixture_name", @@ -13,9 +48,80 @@ ], ) def test_visualize_attributes(request, fixture_name): - """Tests the visualize_attributes method of the StochasticModel class. This - test verifies if the method returns None, which means that the method is - running without breaking. + """Tests the visualize_attributes method of the StochasticModel class. It + must run without breaking and return the formatted report string (which is + also printed), so the report is never silently lost. """ fixture = request.getfixturevalue(fixture_name) - assert fixture.visualize_attributes() is None + report = fixture.visualize_attributes() + assert isinstance(report, str) + assert report + + +def _effective_wind_x(environment): + """The wind the Environment would actually fly with.""" + wind = environment.wind_velocity_x + return float(wind(0)) if callable(wind) else float(wind) + + +def test_reseeding_does_not_take_the_last_run_as_the_next_nominal(): + """Reseeding with the same seed has to give the same inputs. + + ``StochasticEnvironment.create_object`` writes the randomised value back + onto the Environment rather than building a copy, so re-reading the nominal + from it on the next reseed compounded: 10 -> 8.576 -> 7.355 -> 6.308, each + one the last multiplied by the same factor again. + """ + environment = Environment() + environment.set_atmospheric_model(type="standard_atmosphere") + environment.wind_velocity_x = 10.0 + stochastic = StochasticEnvironment( + environment=environment, wind_velocity_x_factor=(1.0, 0.1) + ) + + winds = [] + for _ in range(4): + stochastic._set_stochastic(12345) + winds.append(_effective_wind_x(stochastic.create_object())) + + assert len(set(winds)) == 1, f"the same seed drifted across reseeds: {winds}" + + +def test_a_simulation_index_does_not_depend_on_the_indices_before_it(): + """What the per-index seeding claims: index i gets the same inputs however + it is reached. Running 0, 1, 2 in order has to match running 2 on its own, + which is what a worker that happens to pick up index 2 first would do. + """ + + def wind_for(seeds): + environment = Environment() + environment.set_atmospheric_model(type="standard_atmosphere") + environment.wind_velocity_x = 10.0 + stochastic = StochasticEnvironment( + environment=environment, wind_velocity_x_factor=(1.0, 0.1) + ) + wind = None + for seed in seeds: + stochastic._set_stochastic(seed) + wind = _effective_wind_x(stochastic.create_object()) + return wind + + assert wind_for([101, 102, 103]) == wind_for([103]) + + +def test_a_scalar_nominal_does_not_drift_across_reseeds(): + """Not only the factors. ``_validate_scalar`` and the ``(std, "distribution")`` + tuple both take their nominal from the object, and ``create_object`` writes + the drawn value back onto that same object, so a plain scalar spec drifts + the same way a factor compounds. + """ + environment = Environment() + environment.set_atmospheric_model(type="standard_atmosphere") + stochastic = StochasticEnvironment(environment=environment, elevation=100.0) + + elevations = [] + for _ in range(4): + stochastic._set_stochastic(2024) + elevations.append(float(stochastic.create_object().elevation)) + + assert len(set(elevations)) == 1, f"the nominal elevation drifted: {elevations}" 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..0e3efef8a --- /dev/null +++ b/tests/unit/stochastic/test_stochastic_rocket_seeding.py @@ -0,0 +1,171 @@ +"""Nested StochasticRocket components are reseeded from distinct SeedSequence +children, so components that sample the same distribution (a main and a drogue +parachute, for example) do not draw identical values. Reproducible under a fixed +seed. See the seeding design in ``StochasticRocket._set_stochastic``. +""" + +import ast +import inspect + +import pytest + +from rocketpy.stochastic import StochasticAirBrakes +from rocketpy.stochastic.stochastic_model import StochasticModel + +# Captured once, before any patching, so wrapping it repeatedly in one test does +# not stack (each recorder wraps the real method, not a previous recorder). +_REAL_SET_STOCHASTIC = StochasticModel._set_stochastic + + +def _record_component_seeds(monkeypatch, rocket, seed): + """Return the seeds handed to every nested component for 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 test_rocket_components_receive_distinct_seeds(monkeypatch, stochastic_calisto): + """Every nested component (body, aerodynamic surfaces, motor, rail buttons and + the two parachutes) is reseeded from its own child, so none collide.""" + seeds = _record_component_seeds(monkeypatch, stochastic_calisto, 42) + + assert len(seeds) > 3, "expected the rocket body plus several components" + assert len(seeds) == len(set(seeds)), ( + "components share a seed -- they would draw perfectly correlated samples" + ) + + +def test_rocket_component_seeds_are_reproducible(monkeypatch, stochastic_calisto): + """The same root seed reseeds every component identically; a different root + seed changes them.""" + first = _record_component_seeds(monkeypatch, stochastic_calisto, 42) + again = _record_component_seeds(monkeypatch, stochastic_calisto, 42) + different = _record_component_seeds(monkeypatch, stochastic_calisto, 43) + + assert again == first, "same seed must reproduce every component seed" + assert different != first, "a different seed must change the component seeds" + + +def test_the_reseed_covers_every_collection_create_object_uses(stochastic_calisto): + """Whatever ``create_object`` iterates has to be reseeded too. + + Checked against the source rather than against a fixture, because a + collection that no fixture populates is exactly the one that gets missed: + air brakes were built and sampled and never reseeded, and every seeding + test passed because no fixture had one. + """ + 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 a loop rewritten as one would slip past a For-only 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; the scan is broken" + assert iterated <= declared, ( + f"create_object samples these but the reseed never reaches them: " + f"{sorted(iterated - declared)}" + ) + + +def test_air_brakes_are_reseeded_like_every_other_component( + monkeypatch, stochastic_calisto, calisto_air_brakes_clamp_on +): + """Air brakes were in ``create_object`` and not in the reseed, so their + samples came from wherever the Generator had been left rather than from the + simulation index. Measured before the fix: 3 surfaces, 1 motor, 1 rail + button and 2 parachutes reseeded, air brakes 0 of 1. + """ + stochastic_calisto.add_air_brakes( + calisto_air_brakes_clamp_on.air_brakes[0], + calisto_air_brakes_clamp_on._controllers[0], + ) + air_brake = stochastic_calisto.air_brakes[0] + seen = [] + original = air_brake._set_stochastic + monkeypatch.setattr( + air_brake, + "_set_stochastic", + lambda seed=None: (seen.append(seed), original(seed))[1], + ) + + stochastic_calisto._set_stochastic(42) + + assert seen, "air brakes were not reseeded" + assert seen[0] is not None + + +@pytest.mark.parametrize( + "spec", + [0.001, (0.001, "normal"), (0.0, 0.001, "normal"), [0.0005, 0.001, 0.002]], + ids=["scalar", "tuple2", "tuple3", "list"], +) +def test_eccentricity_is_resampled_from_the_new_generator(stochastic_calisto, spec): + """``add_cp_eccentricity`` and ``add_thrust_eccentricity`` run after + ``__init__``, so their values never reached the dict the base class + re-validates. Validation binds a distribution to the Generator that is live + at the time, so the tuple kept sampling from the one the rocket was built + with: same seed, different eccentricity, while every constructor field + reproduced exactly. + """ + rocket = stochastic_calisto + rocket.add_cp_eccentricity(x=spec, y=spec) + rocket.add_thrust_eccentricity(x=spec, y=spec) + + def sample(): + rocket._set_stochastic(777) + drawn = next(rocket.dict_generator()) + return {k: v for k, v in drawn.items() if "eccentricity" in k} + + first = sample() + + assert len(first) == 4, f"expected four eccentricities, got {sorted(first)}" + assert sample() == first, "the same seed drew a different eccentricity" + + +def test_the_air_brake_sample_follows_the_seed_not_the_call_order( + stochastic_calisto, calisto_air_brakes_clamp_on +): + """That the reseed reaches the air brake is only half of it. + + What matters is the value it draws: the same seed has to give the same + sample, and a different seed a different one. Asserting only that + ``_set_stochastic`` was called would pass over an air brake reseeded with a + constant. + """ + # Built here rather than taken from the fixture: wrapping an AirBrakes with + # no arguments gives every parameter a zero standard deviation, so it draws + # the same values under any seed and the assertions below would hold over an + # air brake that was never reseeded at all. + stochastic_calisto.add_air_brakes( + StochasticAirBrakes( + air_brakes=calisto_air_brakes_clamp_on.air_brakes[0], + drag_coefficient_curve_factor=(1.0, 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 next(air_brake.dict_generator()) + + first = drawn(31337) + + assert first, "the air brake sampled nothing, so this proves nothing" + assert drawn(31337) == first, "the same seed drew a different air brake" + assert drawn(31338) != first, "a different seed drew the same air brake" From 8ed2f5a6b8a2a5222ad9e1e87ca320968acc5e31 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Thu, 6 Aug 2026 17:24:20 +0800 Subject: [PATCH 4/4] Keep the fork's reprint on the serial error paths Missed one of ours when taking the upstream file. The serial run redraws a progress line with a carriage return, so this fork routes its interrupt and error messages through _SimMonitor.reprint rather than print, and a plain print lands on top of that line. The two monitor stubs in the log-integrity tests gain a reprint for the same reason, which is the only place those tests differ from the upstream copy. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- rocketpy/simulation/monte_carlo.py | 6 ++++-- .../unit/simulation/test_monte_carlo_log_integrity.py | 10 ++++++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index cf28eeab2..1f5299db0 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -478,11 +478,13 @@ def __run_in_serial(self): except KeyboardInterrupt: self._interrupted = True - print("Keyboard interrupt received. Files saved.") + # reprint rather than print: the serial run redraws a progress line + # with a carriage return, and a plain print lands on top of it. + _SimMonitor.reprint("Keyboard interrupt received. Files saved.") self.__keep_the_inputs_that_did_not_finish(inputs_json) except Exception as error: - print(f"Error on iteration {sim_monitor.count}: {error}") + _SimMonitor.reprint(f"Error on iteration {sim_monitor.count}: {error}") self.__keep_the_inputs_that_did_not_finish(inputs_json) raise error diff --git a/tests/unit/simulation/test_monte_carlo_log_integrity.py b/tests/unit/simulation/test_monte_carlo_log_integrity.py index 80ce38c23..d9e182414 100644 --- a/tests/unit/simulation/test_monte_carlo_log_integrity.py +++ b/tests/unit/simulation/test_monte_carlo_log_integrity.py @@ -347,6 +347,11 @@ def __init__(self, **_kwargs): def keep_simulating(self): raise KeyboardInterrupt("ctrl-c before the first simulation") + @staticmethod + def reprint(*args, **kwargs): + """This fork routes run messages through the monitor so they do not + land on top of the progress line.""" + monkeypatch.setattr(mc, "_SimMonitor", _Monitor) runner = _serial_runner(tmp_path) @@ -372,6 +377,11 @@ class _Monitor: def __init__(self, **_kwargs): self.laps = 0 + @staticmethod + def reprint(*args, **kwargs): + """This fork routes run messages through the monitor so they do not + land on top of the progress line.""" + def keep_simulating(self): self.laps += 1 if self.laps > 1: