diff --git a/CHANGELOG.md b/CHANGELOG.md index e68153f6e..fe6332fdd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,7 @@ Attention: The newest changes should be on top --> ### Added - ENH: update master with develop [#1081](https://github.com/RocketPy-Team/RocketPy/pull/1081) +- ENH: reproducible Monte Carlo runs via a random_seed argument [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054) ### Changed diff --git a/requirements.txt b/requirements.txt index 61a594320..2cc8a61bc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy>=1.13 +numpy>=1.17 scipy>=1.0 matplotlib>=3.9.0 # Released May 15th 2024 netCDF4>=1.6.4 diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index 21c665d01..345b580a6 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -19,7 +19,7 @@ import traceback import warnings from pathlib import Path -from time import time +from time import monotonic, time import numpy as np import simplekml @@ -30,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, @@ -170,6 +171,8 @@ def simulate( append=False, parallel=False, n_workers=None, + *, + random_seed=None, **kwargs, ): """ @@ -189,6 +192,22 @@ def simulate( 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 the seed it was built from. 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: @@ -219,9 +238,35 @@ 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) + # multiprocess is an optional extra. Imported here, an install + # without rocketpy[monte-carlo] raised only after __setup_files had + # already emptied the previous run's results. + _import_multiprocess() + self._export_config = kwargs self.number_of_simulations = number_of_simulations self._initial_sim_idx = self.num_of_loaded_sims if append else 0 + if append: + _check_the_checkpoint_supports_appending( + self.input_file, self.output_file, self._initial_sim_idx + ) + # 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 + + # 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") @@ -232,6 +277,7 @@ def simulate( else: self.__run_in_serial() + self.__check_each_index_was_recorded_once() self.__terminate_simulation() def __setup_files(self, append): @@ -267,10 +313,168 @@ def __setup_files(self, append): except OSError as error: raise OSError(f"Error creating files: {error}") from error + @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 the seed the generator was built from. + ``rng.bit_generator.seed_seq`` also works, but only on NumPy 1.25 and + above, which is later than this package's floor. + """ + 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, a sequence of non-negative " + "integers, or a numpy.random.SeedSequence, not a " + f"{type(random_seed).__name__}. Pass the seed the generator " + "was built from; rng.bit_generator.seed_seq also works on " + "NumPy 1.25 and above." + ) + 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. + + Only over the indices this run claimed. An ``append`` run exists to + carry on from a file some earlier run left behind, and the documented + way to reach one is to interrupt a run, so that file can hold a torn + row or a pair that disagrees. Judging this run on that damage would + make the very files ``append`` is for the ones it refuses, so anything + below ``_initial_sim_idx`` is reported and not raised on. + """ + inputs, damaged = _recorded_indices("inputs", self.input_file) + outputs, damaged_outputs = _recorded_indices("outputs", self.output_file) + damaged += damaged_outputs + + # First, because a torn row is the root cause and the checks below are + # its symptoms: a row that will not parse also makes the two files + # disagree, and "files disagree" points at the wrong thing. + # + # Always this run's doing. An append only gets here past a preflight + # that read the checkpoint and found it whole, so anything unreadable + # now was written during this run. + if damaged: + raise RuntimeError( + f"{len(damaged)} row(s) this run wrote cannot be read: " + f"{damaged[:5]}. The results are wrong, so they are not " + f"reported as a successful run." + ) + + ours = lambda counts: { # noqa: E731 + index: count + for index, count in counts.items() + if index >= self._initial_sim_idx + } + mine, theirs = ours(inputs), ours(outputs) + if mine != theirs: + 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(mine, theirs)[:5]} have inputs and no " + f"outputs, {only_in(theirs, mine)[: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 mine.items() if count > 1) + beyond = sorted(index for index in mine if index >= self.number_of_simulations) + missing = ( + [] + if self._interrupted + # The whole range, not this run's share of it. Appending is only + # allowed onto a checkpoint the preflight found complete, so what + # ends up on disk has to be every simulation that was asked for. + else sorted(set(range(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. + The root seed state is captured by ``simulate`` before this runs, so each + simulation index derives its child seed from ``self.__root_state``. + Returns ------- None @@ -281,38 +485,51 @@ def __run_in_serial(self): start_time=time(), ) try: - while sim_monitor.keep_simulating(): - sim_monitor.increment() - inputs_json, outputs_json = "", "" + 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_simulation(self.__child_seed(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) - - 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) + inputs_json = self.__evaluate_flight_inputs(sim_idx) + outputs_json = self.__evaluate_flight_outputs(flight, sim_idx) + _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._interrupted = True print("Keyboard interrupt received. Files saved.") - with open(self._error_file, "a", encoding="utf-8") as f: - f.write(inputs_json) + self.__keep_the_inputs_that_did_not_finish(inputs_json) except Exception as error: print(f"Error on iteration {sim_monitor.count}: {error}") - with open(self._error_file, "a", encoding="utf-8") as f: - f.write(inputs_json) + self.__keep_the_inputs_that_did_not_finish(inputs_json) raise error + 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, 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 ---------- n_workers: int, optional @@ -338,62 +555,73 @@ def __run_in_parallel(self, n_workers=None): start_time=time(), ) - processes = [] - seeds = np.random.SeedSequence().spawn(n_workers) - - for seed in seeds: - sim_producer = multiprocess.Process( - target=self.__sim_producer, - args=( - seed, - 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) + # The event asks them to stop, it does not stop them. Without + # this window a worker part way through a write is cut off and + # leaves exactly the torn row the check below would report. + # Not _bring_the_fleet_down: that sets the event, which on a run + # that finished cleanly is what the crash check reads next. + _wait_for_workers(started_processes, timeout=_WORKER_SHUTDOWN_GRACE) + _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, seed, 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 ---------- - seed : int - The seed to set the random number generator. sim_monitor : _SimMonitor The simulation monitor object to keep track of the simulations. mutex : multiprocess.Lock @@ -402,21 +630,26 @@ 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 = "", "" - + 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(): # Runs in a worker process spawned via multiprocessing: # logging handlers configured in the main process are @@ -431,27 +664,55 @@ def __sim_producer(self, seed, sim_monitor, mutex, error_event): # pylint: disa 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() + if acquired: + 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) + 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() - # See note above: must use print() to remain visible from a - # multiprocessing worker process. - _SimMonitor.reprint( - f"Error on iteration {sim_idx}:\n{traceback.format_exc()}" - ) - error_event.set() - mutex.release() + # 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" + + 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() + + # 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. @@ -1377,7 +1638,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, @@ -1607,6 +1868,243 @@ def export_errors_to_json(self, filename): self._write_log_to_json(self.errors_log, filename) +def _recorded_indices(label, path): + """``({index: how many rows carry it}, [rows that carry no usable index])``. + + Damage is returned rather than raised on. Whether a torn row matters + depends on which run wrote it, and only the caller knows the range this + run claimed: an ``append`` run is recovering from a file some earlier run + damaged, which is the whole reason it is appending. + + ``type(...) is int`` and not ``isinstance``: ``True`` and ``1.0`` both + compare equal to ``1`` and would otherwise pass for it. + """ + written, damaged = {}, [] + 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: + damaged.append(f"{label} row {number} is not readable JSON") + continue + 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 + damaged.append(f"{label} row {number} carries no simulation index") + continue + written[index] = written.get(index, 0) + 1 + return written, damaged + + +def _check_the_checkpoint_supports_appending(input_file, output_file, resume_at): + """Everything that can be judged from the files, before a worker starts. + + ``num_of_loaded_sims`` counts lines rather than indices, so a blank line or + a torn row moves the resume point past an index that was never run. The run + then skips it, and a check scoped to the new range calls that a success. + Measured: two rows plus one blank line resume at 3, plus two blanks at 4, + while the file holds only 0 and 1 either way. + + Held here rather than after the run so a checkpoint that cannot be resumed + costs no simulations and is left exactly as it was found. + + A file with a hole in it is refused rather than repaired. Filling holes + needs the workers to claim from a plan instead of counting on from the end, + which is #1075; until then, refusing loudly beats resuming in the wrong + place quietly. + """ + for label, path in (("inputs", input_file), ("outputs", output_file)): + written, damaged = _recorded_indices(label, path) + if damaged: + raise ValueError( + f"cannot append to {path}: {len(damaged)} row(s) cannot be " + f"read, so the simulations they held cannot be accounted for: " + f"{damaged[:3]}." + ) + _refuse_a_checkpoint_that_does_not_line_up(label, path, written, resume_at) + + inputs, _ = _recorded_indices("inputs", input_file) + outputs, _ = _recorded_indices("outputs", output_file) + if inputs != outputs: + raise ValueError( + f"cannot append: the input and output files hold different " + f"simulations, {sorted(set(inputs) - set(outputs))[:5]} against " + f"{sorted(set(outputs) - set(inputs))[:5]}. Appending would build " + f"on a checkpoint that is already inconsistent." + ) + + +def _refuse_a_checkpoint_that_does_not_line_up(label, path, written, resume_at): + """One file's indices have to be 0..resume_at-1, with nothing repeated.""" + repeated = sorted(index for index, count in written.items() if count > 1) + if repeated: + raise ValueError( + f"cannot append to {path}: {label} hold {len(repeated)} index(es) " + f"more than once {repeated[:5]}." + ) + + indices = set(written) + if indices == set(range(1, len(indices) + 1)) and indices: + # The serial path used to number from 1. Named rather than reported as + # an off-by-one, because the fix is to re-baseline, not to retry. + raise ValueError( + f"cannot append to {path}: the {label} are numbered from 1, which " + f"is how versions before per-index seeding wrote serial runs. This " + f"release numbers from 0, so the two cannot be continued into each " + f"other. Re-run the study, or renumber the file down by one." + ) + if indices != set(range(resume_at)): + missing = sorted(set(range(resume_at)) - indices) + extra = sorted(indices - set(range(resume_at))) + raise ValueError( + f"cannot append to {path}: the run would start at index " + f"{resume_at}, but the {label} are not the {resume_at} before it. " + f"Missing {missing[:5]}, unexpected {extra[:5]}." + ) + + +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 multiprocess library. diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index ca26f6578..c650097ed 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,37 @@ def __init__(self, obj, seed=None, **kwargs): self.obj = obj self.last_rnd_dict = {} self.__stochastic_dict = kwargs + self.__nominal_values = {} self._set_stochastic(seed) + def _nominal(self, input_name, getter=getattr): + """``self.obj``'s value for ``input_name``, as it was when this model + was built. + + Read once and remembered, because ``StochasticEnvironment`` has + ``create_object`` write the randomised value back onto ``self.obj`` + instead of building a copy. Re-reading it on a reseed would take one + simulation's output as the next one's nominal, and a factor would + multiply the factor before it rather than the original value. + + A custom ``getter`` reads a component's own attribute rather than one + of ``self.obj``'s, and nothing writes back to those, so it is passed + straight through. Caching it here would be wrong as well: every + component's position arrives under the one name ``"position"``. + + This applies to every stochastic model, not only the environment: what + a model samples around is what the wrapped object held when the model + was built. Changing the object afterwards does not move it. Only + ``StochasticEnvironment.create_object`` writes back today, but the rule + is stated for all of them rather than special-cased for one, so a model + means the same thing whichever object it wraps. + """ + if getter is not getattr: + return getter(self.obj, input_name) + if input_name not in self.__nominal_values: + self.__nominal_values[input_name] = getattr(self.obj, input_name) + return self.__nominal_values[input_name] + def _set_stochastic(self, seed=None): """Set the stochastic attributes from the input dictionary. This method is useful to reset or reseed the attributes of the instance. @@ -109,7 +136,7 @@ def _set_stochastic(self, seed=None): "or a custom sampler" ) else: - attr_value = [getattr(self.obj, input_name)] + attr_value = [self._nominal(input_name)] setattr(self, input_name, attr_value) def __repr__(self): @@ -186,7 +213,7 @@ def _validate_tuple_length_two(self, input_name, input_value, getattr=getattr): # function. In this case, the nominal value will be taken from the # object passed. dist_func = get_distribution(input_value[1], self.__random_number_generator) - return (getattr(self.obj, input_name), input_value[0], dist_func) + return (self._nominal(input_name, getattr), input_value[0], dist_func) else: # if second item is an int or float, then it is assumed that the # first item is the nominal value and the second item is the @@ -257,7 +284,7 @@ def _validate_list(self, input_name, input_value, getattr=getattr): # pylint: d If the input is not in a valid format. """ if not input_value: - return [getattr(self.obj, input_name)] + return [self._nominal(input_name, getattr)] else: return input_value @@ -283,7 +310,7 @@ def _validate_scalar(self, input_name, input_value, getattr=getattr): # pylint: distribution function). """ return ( - getattr(self.obj, input_name), + self._nominal(input_name, getattr), input_value, get_distribution("normal", self.__random_number_generator), ) @@ -310,7 +337,7 @@ def _validate_factors(self, input_name, input_value, seed): If the input is not in a valid format. """ attribute_name = input_name.replace("_factor", "") - setattr(self, f"_{attribute_name}", getattr(self.obj, attribute_name)) + setattr(self, f"_{attribute_name}", self._nominal(attribute_name)) if isinstance(input_value, tuple): return self._validate_tuple_factor(input_name, input_value) @@ -508,6 +535,21 @@ def _validate_airfoil(self, airfoil): "the first item" ) + def _random_choice(self, values): + """Choose one value from a list using this model's seeded generator. + + The index is drawn from the seeded generator, not the stdlib global + ``random.choice`` (an unseeded shared instance), so the choice is + governed by ``random_seed``. Indexing rather than ``numpy.random.choice`` + keeps a heterogeneous list -- ``Function`` objects, paths, arrays -- + returned as itself instead of coerced to a common dtype. An empty + ``values`` is returned unchanged. + """ + if not values: + return values + index = int(self.__random_number_generator.integers(len(values))) + return values[index] + def dict_generator(self): """ Generate a dictionary with randomly generated input arguments. @@ -532,7 +574,7 @@ def dict_generator(self): dist_sampler = value[-1] generated_dict[arg] = dist_sampler(value[0], value[1]) elif isinstance(value, list): - generated_dict[arg] = choice(value) if value else value + generated_dict[arg] = self._random_choice(value) elif isinstance(value, CustomSampler): try: generated_dict[arg] = value.sample(n_samples=1)[0] diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 794a66c85..007c52adc 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 @@ -21,6 +22,7 @@ from rocketpy.rocket.rocket import Rocket from rocketpy.stochastic.stochastic_generic_motor import StochasticGenericMotor from rocketpy.stochastic.stochastic_motor_model import StochasticMotorModel +from rocketpy.tools import _seed_sequence_to_int from .stochastic_aero_surfaces import ( StochasticAirBrakes, @@ -154,6 +156,13 @@ def __init__( self.air_brakes = [] self.parachutes = [] self.__components_map = {} + # Raw eccentricity arguments, kept as the caller gave them. + # ``add_cp_eccentricity`` and ``add_thrust_eccentricity`` run after + # ``__init__``, so their values are not in the dict the base class + # re-validates on a reseed. Validating them once would leave the + # distribution bound to the Generator of whichever simulation happened + # to come first, so the raw form is kept and validated again each time. + self.__eccentricity_specs = {} super().__init__( obj=rocket, radius=radius, @@ -172,25 +181,47 @@ def __init__( coordinate_system_orientation=None, ) + # Every collection of nested stochastic objects, in the order their child + # seeds are spawned. Listed here rather than written out inline so that a + # component type cannot end up in ``create_object`` and not in the reseed: + # air brakes were, and their sampling depended on which worker ran the + # index instead of on the index. ``_stochastic_collections`` is asserted + # against the rocket's own attributes in the tests. + _POSITIONED_COLLECTIONS = ("aerodynamic_surfaces", "motors", "rail_buttons") + _PLAIN_COLLECTIONS = ("parachutes", "air_brakes") + + @classmethod + def _stochastic_collections(cls): + """The names of every attribute holding nested stochastic objects.""" + return cls._POSITIONED_COLLECTIONS + cls._PLAIN_COLLECTIONS + def _set_stochastic(self, seed=None): """Set the stochastic attributes for Components, positions and inputs. + Every nested component -- the rocket body, each aerodynamic surface, + motor, rail button, parachute and air brake -- is reseeded from its own + child of a ``SeedSequence`` root, so components that sample the same + distribution do not draw identical values (a main and a drogue parachute + get independent ``cd_s`` and ``lag`` samples, not the same one). Children + are spawned in a fixed order, so the result stays reproducible under + ``random_seed``. + Parameters ---------- seed : int, optional Seed for the random number generator. """ - super()._set_stochastic(seed) - self.aerodynamic_surfaces = self.__reset_components( - self.aerodynamic_surfaces, seed - ) - self.motors = self.__reset_components(self.motors, seed) - self.rail_buttons = self.__reset_components(self.rail_buttons, seed) - for parachute in self.parachutes: - parachute._set_stochastic(seed) - - def __reset_components(self, components, seed): + root = np.random.SeedSequence(seed) + super()._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0])) + self.__apply_eccentricity_specs() + for name in self._POSITIONED_COLLECTIONS: + setattr(self, name, self.__reset_components(getattr(self, name), root)) + for name in self._PLAIN_COLLECTIONS: + for child in getattr(self, name): + child._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0])) + + def __reset_components(self, components, root): """Creates a new Components whose stochastic structures and their positions are reset. @@ -199,8 +230,9 @@ def __reset_components(self, components, seed): components : Components The components which contains the stochastic structure that will be used to create the new components. - seed : int, optional - Seed for the random number generator. + root : numpy.random.SeedSequence + The run's seed root. Each component is reseeded from its own spawned + child, so components sampling the same distribution stay decorrelated. Returns ------- @@ -212,7 +244,7 @@ def __reset_components(self, components, seed): new_components = Components() for stochastic_obj, _ in components: stochastic_obj_position_info = self.__components_map[stochastic_obj] - stochastic_obj._set_stochastic(seed) + stochastic_obj._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0])) new_components.add( stochastic_obj, self._validate_position(stochastic_obj, stochastic_obj_position_info), @@ -434,8 +466,9 @@ def add_cp_eccentricity(self, x=None, y=None): self : StochasticRocket Object of the StochasticRocket class. """ - self.cp_eccentricity_x = self._validate_eccentricity("cp_eccentricity_x", x) - self.cp_eccentricity_y = self._validate_eccentricity("cp_eccentricity_y", y) + self.__eccentricity_specs["cp_eccentricity_x"] = x + self.__eccentricity_specs["cp_eccentricity_y"] = y + self.__apply_eccentricity_specs() return self def add_thrust_eccentricity(self, x=None, y=None): @@ -460,14 +493,24 @@ def add_thrust_eccentricity(self, x=None, y=None): self : StochasticRocket Object of the StochasticRocket class. """ - self.thrust_eccentricity_x = self._validate_eccentricity( - "thrust_eccentricity_x", x - ) - self.thrust_eccentricity_y = self._validate_eccentricity( - "thrust_eccentricity_y", y - ) + self.__eccentricity_specs["thrust_eccentricity_x"] = x + self.__eccentricity_specs["thrust_eccentricity_y"] = y + self.__apply_eccentricity_specs() return self + def __apply_eccentricity_specs(self): + """Re-validate the eccentricities against the current Generator. + + Validation stores a distribution as a method bound to the Generator + that was live at the time, so a tuple validated once keeps sampling + from that one. Re-running it after every reseed is what ties the draw + to the simulation index rather than to whichever index the worker + happened to run first. ``get_distribution`` only binds a method, so + this consumes no randomness and does not shift any other draw. + """ + for name, spec in self.__eccentricity_specs.items(): + setattr(self, name, self._validate_eccentricity(name, spec)) + def _validate_eccentricity(self, eccentricity, position): """Validate the eccentricity argument. @@ -628,7 +671,7 @@ def _randomize_position(self, position): return position[-1](position[0].z, position[1]) return position[-1](position[0], position[1]) elif isinstance(position, list): - return choice(position) if position else position + return self._random_choice(position) # pylint: disable=stop-iteration-return def dict_generator(self): @@ -638,8 +681,8 @@ def dict_generator(self): all attributes of the class and generating a random value for each attribute. The random values are generated according to the format of each attribute. Tuples are generated using the distribution function - specified in the tuple. Lists are generated using the random.choice - function. + specified in the tuple. Lists are sampled through the model's seeded + generator so the choice is governed by ``random_seed``. Parameters ---------- diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 0d7f1a74e..9df900eb5 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -1467,6 +1467,29 @@ def find_obj_from_hash(obj, hash_, depth_limit=None): return None +def _seed_sequence_to_int(seed_sequence): + """Collapse a ``SeedSequence`` into a 128-bit Python ``int``. + + A plain ``int`` is what ``numpy.random.default_rng`` and the stdlib + ``random.Random`` both accept (``random.Random`` rejects a ``SeedSequence`` + with a ``TypeError`` since Python 3.11), so a custom sampler whose + ``reset_seed`` documents an ``int`` and builds a modern generator keeps + working. The legacy ``numpy.random.RandomState`` is the exception: it caps a + single-integer seed at ``2**32 - 1``, so a sampler still built on it would + have to reduce the value (``RandomState`` is a frozen legacy API NumPy steers + new code away from). All four ``uint32`` words are combined to keep the full + 128-bit pool, so sub-streams stay decorrelated instead of collapsing to a + single 32-bit word. + + The words are combined by value (little-endian word order), not via + ``tobytes()``, so the seed is the same on big- and little-endian machines -- + a byte-order-dependent seed would break the cross-platform reproducibility + this exists to provide. + """ + words = seed_sequence.generate_state(4, dtype=np.uint32) + return sum(int(word) << (32 * position) for position, word in enumerate(words)) + + if __name__ == "__main__": # pragma: no cover import doctest diff --git a/tests/fixtures/monte_carlo/custom_sampler_fixtures.py b/tests/fixtures/monte_carlo/custom_sampler_fixtures.py index 8a4ff497d..e3ad85d50 100644 --- a/tests/fixtures/monte_carlo/custom_sampler_fixtures.py +++ b/tests/fixtures/monte_carlo/custom_sampler_fixtures.py @@ -33,7 +33,7 @@ def __init__(self, means_tuple, sd_tuple, prob_tuple, seed=None): 2-Tuple that contains the probability of each normal distribution of the mixture. Its entries should be non-negative and sum up to 1. """ - np.random.default_rng(seed) + self.reset_seed(seed) self.means_tuple = means_tuple self.sd_tuple = sd_tuple self.prob_tuple = prob_tuple @@ -52,16 +52,12 @@ def sample(self, n_samples=1): List containing n_samples samples """ samples_list = [0] * n_samples - mixture_id_list = np.random.binomial(1, self.prob_tuple[0], n_samples) + mixture_id_list = self.rng.binomial(1, self.prob_tuple[0], n_samples) for i, mixture_id in enumerate(mixture_id_list): if mixture_id: - samples_list[i] = np.random.normal( - self.means_tuple[0], self.sd_tuple[0] - ) + samples_list[i] = self.rng.normal(self.means_tuple[0], self.sd_tuple[0]) else: - samples_list[i] = np.random.normal( - self.means_tuple[1], self.sd_tuple[1] - ) + samples_list[i] = self.rng.normal(self.means_tuple[1], self.sd_tuple[1]) return samples_list @@ -73,4 +69,7 @@ def reset_seed(self, seed=None): seed : int, optional Seed for the random number generator. """ - np.random.default_rng(seed) + # Kept on the instance. Building a generator and dropping it made this + # a no-op, and sample() went on drawing from the process-global + # np.random, so nothing here answered to the seed at all. + self.rng = np.random.default_rng(seed) 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..87e0724ba --- /dev/null +++ b/tests/integration/simulation/test_monte_carlo_determinism.py @@ -0,0 +1,904 @@ +"""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) + + +def test_the_real_parallel_path_is_worker_invariant_on_this_platform( + tmp_path, + stochastic_environment_with_wind, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """The same property as the test below, on whatever start method this + platform uses, and without the ``slow`` marker. + + The thorough version covers fork, spawn and forkserver, but it is marked + slow and pull-request CI skips slow tests, so the path this change exists + to support gated nothing. This one is small enough to run every time, and + because it takes the platform default, each CI job ends up gating the start + method it actually uses: spawn on Windows and macOS, forkserver on Python + 3.14's POSIX default, fork below that. + """ + count = 2 + common = {"number_of_simulations": count, "random_seed": 24680} + models = ( + stochastic_environment_with_wind, + stochastic_calisto_numpy_only, + stochastic_flight, + ) + serial = _real_run_inputs(tmp_path, *models, "here-serial", **common)[1] + parallel = _real_run_inputs( + tmp_path, *models, "here-p2", parallel=True, n_workers=2, **common + )[1] + + assert sorted(serial) == list(range(count)) + assert sorted(parallel) == list(range(count)) + for index in range(count): + expected = _sampled_only(json.loads(serial[index])) + actual = _sampled_only(json.loads(parallel[index])) + assert len(expected) > 20, f"only {len(expected)} fields left to compare" + assert actual == expected, ( + f"{multiprocessing.get_start_method()}: serial and parallel(2) " + f"differ at index {index}" + ) + + +@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. + + This run wrote the whole file, so the damage is its own. A run appending + onto a file an earlier run damaged is judged only on what it added, which + ``test_an_append_run_is_not_judged_on_the_damage_it_inherited`` covers. + """ + 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="cannot be read"): + 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_continues_a_checkpoint_and_leaves_the_whole_range( + tmp_path, stochastic_environment, stochastic_calisto, stochastic_flight +): + """A second run carries on from the first, and the pair ends up whole. + + This test used to empty both logs before appending and then assert that a + four-simulation result holding only indices 2 and 3 was a success. That is + the shape of the bug it was meant to guard: the resume point came from a + row count rather than the indices actually on disk, so nothing noticed the + first two were gone. The run is judged on the whole range now. + """ + 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 + + 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" + ) + for label, path in ( + ("inputs", montecarlo.input_file), + ("outputs", montecarlo.output_file), + ): + assert sorted(_read_inputs_by_index(path)) == [0, 1, 2, 3], ( + f"the {label} do not hold every simulation that was asked for" + ) + + +def test_appending_onto_a_checkpoint_with_a_hole_is_refused( + tmp_path, stochastic_environment, stochastic_calisto, stochastic_flight +): + """The other half, and the reason the resume point cannot be a row count. + + Two rows plus a blank line load as three simulations, so the next run would + start at index 2 and leave index 1 missing for good while reporting + success. Refused before it runs, with both files left as they were found. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / "holed"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + montecarlo.simulate(number_of_simulations=2, random_seed=606) + + with open(montecarlo.output_file, "a", encoding="utf-8") as log: + log.write("\n") + montecarlo.set_num_of_loaded_sims() + assert montecarlo.num_of_loaded_sims == 3, "the blank line was not counted" + before = ( + montecarlo.input_file.read_bytes(), + montecarlo.output_file.read_bytes(), + ) + + with pytest.raises(ValueError): + montecarlo.simulate(number_of_simulations=5, append=True, random_seed=606) + + assert ( + montecarlo.input_file.read_bytes(), + montecarlo.output_file.read_bytes(), + ) == before, "a refused checkpoint was modified on the way out" + + +def test_a_missing_parallel_dependency_does_not_cost_the_previous_run( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """``multiprocess`` is an optional extra, so an install without + ``rocketpy[monte-carlo]`` cannot run in parallel at all. + + It used to be imported inside the parallel path, which runs after + ``__setup_files`` has opened both logs "w+" and emptied them, so asking for + a parallel run on such an install destroyed the previous results on the way + to the ImportError. + """ + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + montecarlo = MonteCarlo( + filename=str(tmp_path / "kept"), + 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") + + def no_multiprocess(): + raise ImportError("No module named 'multiprocess'") + + monkeypatch.setattr(mc_module, "_import_multiprocess", no_multiprocess) + + with pytest.raises(ImportError): + montecarlo.simulate( + number_of_simulations=2, parallel=True, n_workers=2, random_seed=7 + ) + + with open(montecarlo.input_file, encoding="utf-8") as kept: + assert kept.read() == "previous results\n" 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..01b20f41b --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_log_integrity.py @@ -0,0 +1,466 @@ +"""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", "cannot be read"), + "a row that is not an object": (COMPLETE + "[]\n", "cannot be read"), + "a row with no index": (COMPLETE + '{"foo": 1}\n', "cannot be read"), + "a boolean index": ('{"index": 0}\n{"index": true}\n', "cannot be read"), + "a float index": ('{"index": 0}\n{"index": 1.0}\n', "cannot be read"), + "a negative index": (COMPLETE + '{"index": -1}\n', "cannot be read"), + "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_checkpoint_that_cannot_be_read_is_refused_before_the_run(tmp_path): + """Where the history is judged: before anything runs, not after. + + An earlier round tolerated inherited damage at the end of the run, on the + grounds that the documented way to reach an append is to interrupt a run. + That was the wrong place for it. A torn row holds an index nobody can + recover, so the resume point cannot be trusted either, and resuming at the + wrong one silently skips a simulation. The preflight refuses instead, with + both files left exactly as they were found. + """ + rows = '{"index": 0}\n{not json\n' + inputs, outputs = tmp_path / "i.txt", tmp_path / "o.txt" + inputs.write_text(rows) + outputs.write_text(rows) + before = inputs.read_bytes(), outputs.read_bytes() + + with pytest.raises(ValueError, match="cannot be read"): + mc._check_the_checkpoint_supports_appending(inputs, outputs, 2) + + assert (inputs.read_bytes(), outputs.read_bytes()) == before, ( + "a refused checkpoint was modified on the way out" + ) + + +def test_this_run_is_still_judged_strictly_while_appending(tmp_path): + """The other half. Tolerating the history must not tolerate the rows this + run added, or appending would become a way to launder a bad run.""" + runner = _runner( + tmp_path, + '{"index": 0}\n{"index": 2}\n{"index": 2}\n', + count=4, + initial=2, + ) + with pytest.raises(RuntimeError, match="more than once"): + _check(runner) + + +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" + ) + + +def test_a_history_that_went_missing_is_still_caught_at_the_end(tmp_path): + """The run is judged on every index asked for, not on its own share. + + Appending normally reaches this past a preflight that found the checkpoint + whole, so the two questions have the same answer there. They do not when + the check is asked directly, and the invariant worth stating is the one + about the whole file: a four-simulation result holds four simulations. + """ + runner = _runner(tmp_path, '{"index": 2}\n{"index": 3}\n', count=4, initial=2) + + with pytest.raises(RuntimeError, match="never written"): + _check(runner) + + +def test_a_checkpoint_numbered_from_one_is_named_as_such(tmp_path): + """Serial runs used to number from 1. Appending onto one would rewrite the + last index rather than continue, so it is refused by name: the fix is to + re-baseline, not to retry, and an off-by-one message would not say that. + """ + rows = "".join('{"index": %d}\n' % index for index in (1, 2, 3)) + inputs, outputs = tmp_path / "i.txt", tmp_path / "o.txt" + inputs.write_text(rows) + outputs.write_text(rows) + + with pytest.raises(ValueError, match="numbered from 1"): + mc._check_the_checkpoint_supports_appending(inputs, outputs, 3) 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..51330dfbe --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_worker_exit.py @@ -0,0 +1,240 @@ +"""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 +from time import monotonic + +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_a_worker_still_writing_gets_a_window_before_it_is_signalled( + parallel_runner, monkeypatch +): + """The event asks the fleet to stop, it does not stop it. + + Cutting a worker off the moment another one reports an error truncates + whatever row it was part way through, which is the corruption the + completeness check then reports. Measured before the fix: the main error + path gave a running worker 0.0 ms, while the interrupt path gave it the + full grace period. + """ + grace = 0.3 + monkeypatch.setattr(mc, "_WORKER_SHUTDOWN_GRACE", grace) + signalled = [] + original_terminate = _Process.terminate + + def note_when(self): + signalled.append(monotonic()) + self._alive = False + original_terminate(self) + + monkeypatch.setattr(_Process, "start", lambda self: setattr(self, "_alive", True)) + monkeypatch.setattr( + _Process, "is_alive", lambda self: getattr(self, "_alive", False) + ) + monkeypatch.setattr(_Process, "terminate", note_when) + + # Another worker has already reported an error while this one is going. + class _AlreadyFailed(_Event): + def __init__(self): + super().__init__() + self.flag = True + + class _FailedManager: # pylint: disable=invalid-name + def Lock(self): # noqa: N802 + return types.SimpleNamespace(acquire=lambda: None, release=lambda: None) + + def Event(self): # noqa: N802 + return _AlreadyFailed() + + def _SimMonitor(self, **kwargs): # noqa: N802 + return _Monitor(**kwargs) + + @contextmanager + def failed_manager(*_a, **_k): + yield _FailedManager() + + monkeypatch.setattr(mc, "_create_multiprocess_manager", failed_manager) + + began = monotonic() + with pytest.raises(RuntimeError): + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=2) + + assert signalled, "the worker was never signalled at all" + assert signalled[0] - began >= grace, ( + f"the worker was cut off after {(signalled[0] - began) * 1000:.1f} ms, " + f"before the {grace * 1000:.0f} ms window it is meant to get" + ) + + +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/stochastic/test_stochastic_model.py b/tests/unit/stochastic/test_stochastic_model.py index 9e35a5330..e41a62751 100644 --- a/tests/unit/stochastic/test_stochastic_model.py +++ b/tests/unit/stochastic/test_stochastic_model.py @@ -1,5 +1,42 @@ +from types import SimpleNamespace + +import numpy as np + 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", @@ -21,3 +58,134 @@ def test_visualize_attributes(request, fixture_name): 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}" + + +def test_a_custom_sampler_answers_to_the_seed_it_is_given(elevation_sampler): + """The 128-bit int this package hands a sampler has to reach its draws. + + The fixture used to build a generator in ``reset_seed`` and drop it, while + ``sample`` drew from the process-global ``np.random``, so nothing in it + answered to a seed and the guarantee went untested. + """ + wide = 271828182845904523536028747135266249775 + + elevation_sampler.reset_seed(wide) + first = elevation_sampler.sample(5) + elevation_sampler.reset_seed(wide) + again = elevation_sampler.sample(5) + + assert first == again, "the same seed gave a different sample" + + elevation_sampler.reset_seed(wide + 1) + other = elevation_sampler.sample(5) + + assert other != first, "a different seed gave the same sample" + + +def test_a_custom_sampler_is_not_moved_by_the_global_generator(elevation_sampler): + """The control for the test above. Drawing from the global stream in + between must not change what the seeded sampler produces, or the sampler is + still reading from somewhere this package does not seed.""" + seed = 12345678901234567890123456789012345678 + + elevation_sampler.reset_seed(seed) + expected = elevation_sampler.sample(5) + + elevation_sampler.reset_seed(seed) + np.random.random(100) + assert elevation_sampler.sample(5) == expected + + +def test_the_nominal_is_the_one_the_model_was_built_with(example_plain_env): + """Snapshot semantics, stated once and pinned here. + + A model samples around what the wrapped object held when it was built. + This exists because ``StochasticEnvironment.create_object`` writes the + sampled value back onto that object on purpose, and reading the nominal + back off it made a factor compound from one simulation to the next. The + rule is the same for every model, so a change to the wrapped object after + construction deliberately does not move what is sampled around. + """ + example_plain_env.elevation = 1000 + # A scalar is a spread around the object's own value, so this is the form + # that reads the nominal. A tuple carries its own centre and would not. + model = StochasticEnvironment(environment=example_plain_env, elevation=5) + + model._set_stochastic(4242) + around_first = model.elevation[0] + + example_plain_env.elevation = 9000 + model._set_stochastic(4242) + + assert model.elevation[0] == around_first == 1000, ( + "the model followed the object instead of the value it was built with" + ) 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"