diff --git a/CHANGELOG.md b/CHANGELOG.md index 38c9028b4..ce883ebdf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -41,6 +41,7 @@ Attention: The newest changes should be on top --> - MNT: declare dependency floors the package can actually run on [#1108](https://github.com/RocketPy-Team/RocketPy/pull/1108) - CI: build the docs for pull requests into develop as well [#1104](https://github.com/RocketPy-Team/RocketPy/pull/1104) +- ENH: Make Monte Carlo input sampling reproducible per simulation index via a `random_seed` argument. Fixed-seed samples change, serial log indices are now zero-based to match the parallel path, and checkpoints written by the previous scheme cannot be resumed. [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054) - CI: make changelog automation LLM-based (Gemini) and race-safe [#1082](https://github.com/RocketPy-Team/RocketPy/pull/1082) - ENH: Resolve pressure_ISA discretization bounds TODO [#1056](https://github.com/RocketPy-Team/RocketPy/pull/1056) diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index 21c665d01..56625262f 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -18,8 +18,10 @@ import os import traceback import warnings +from contextlib import contextmanager +from copy import deepcopy from pathlib import Path -from time import time +from time import monotonic, sleep, time import numpy as np import simplekml @@ -30,12 +32,16 @@ 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, ) # TODO: Create evolution plots to analyze convergence +# Written by the run itself and used to pair an inputs row with its outputs row. +# A collector that supplies one can relabel a row without tripping any check. +_RESERVED_RECORD_KEYS = frozenset({"index"}) class MonteCarlo: # pylint: disable=too-many-public-methods @@ -170,6 +176,8 @@ def simulate( append=False, parallel=False, n_workers=None, + *, + random_seed=None, **kwargs, ): """ @@ -180,8 +188,11 @@ def simulate( number_of_simulations : int Number of simulations to be run, must be non-negative. append : bool, optional - If True, the results will be appended to the existing files. If - False, the files will be overwritten. Default is False. + If True, resume the existing files. ``number_of_simulations`` is + then the target total rather than a number to add, and a value + below what the files already hold is refused. The root seed is not + stored in them, so pass the same ``random_seed`` to keep the + streams. If False, the files will be overwritten. Default is False. parallel : bool, optional If True, the simulations will be run in parallel. Default is False. n_workers : int, optional @@ -189,6 +200,42 @@ 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 mapping from simulation + index to sampled inputs is 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. The rows themselves are written in completion order, since + a worker takes the log lock once its simulation is done, so the file + order can differ between runs. Compare by the recorded index rather + than byte for byte. + + What is reproduced is the sampled values. With + ``include_function_data=True`` a record also carries a + ``Function``'s signature hash and serialised source, which describe + the object rather than the value drawn for it, so a run under + ``spawn`` or ``forkserver`` writes different ones for the same + inputs. Measured on a real run, six fields differ across that + boundary and all six are these. Pass + ``include_function_data=False`` when the records need to compare + field for field. + + Reproducibility is also scoped to one environment. NumPy promises a + stream only for the same BitGenerator, seed, call sequence, build + and machine, and reserves the right to change what ``default_rng`` + returns. A seed fixes the lineage of a run; it is not an archive + format that survives a version bump. 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: @@ -213,15 +260,60 @@ def simulate( the simulation by running the ``simulate`` method again with the same number of simulations and setting `append=True`. + Not every interruption leaves a checkpoint that can be continued. The + two logs have to hold the same simulations as a complete run of indices + from zero, and a parallel run can stop with one worker's index missing + while a later one is already written. Such a checkpoint is refused + rather than repaired, which is #1075. + Important --------- If you use `append=False` and the files already exist, they will be 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) + # Again rather than only in __init__: the attribute is public and a key + # added after construction would otherwise reach the logs unchecked. + self._check_data_collector(self.data_collector) + 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 + ) + # ``number_of_simulations`` is the target to reach, not a batch to + # add. Below the checkpoint it ran nothing, reported success, and + # left a file with more simulations than the caller had asked for. + if number_of_simulations < self._initial_sim_idx: + raise ValueError( + f"number_of_simulations is the total to reach when " + f"append=True. The checkpoint already holds " + f"{self._initial_sim_idx} simulations, more than the " + f"requested {number_of_simulations}." + ) + # 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 +324,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 +360,172 @@ 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. + + Deep-copied, because ``SeedSequence`` keeps a sequence entropy by + reference. Without it a caller who mutates the list they passed changes + the children this run derives, which is the opposite of a snapshot. + """ + root = self.__root_seed_sequence(random_seed) + self.__root_state = ( + deepcopy(root.entropy), + tuple(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 @@ -280,39 +535,70 @@ def __run_in_serial(self): n_simulations=self.number_of_simulations, start_time=time(), ) + # Bound before the loop: a failure on the very first iteration + # would otherwise reach the error record with no index at all. + sim_idx = self._initial_sim_idx 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 + ) + # The pair is on disk. Cleared before the monitor call so a + # failure there reports itself rather than reporting a row that + # has already been committed as one that never finished. + inputs_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) - raise error + # Captured before reporting, which may fail and must not be what + # gets recorded or raised. + _record_failure( + self._error_file, sim_idx, inputs_json, traceback.format_exc() + ) + # Bare, so the handler's own line does not join the traceback. + raise + + def __keep_the_inputs_that_did_not_finish(self, inputs_json): + """Append the inputs of a simulation that stopped part way through. + + Best effort: an unwritable error file must not turn a clean interrupt + into a crash. + """ + _best_effort( + lambda: _write_unfinished_inputs(self._error_file, inputs_json), + "interrupted simulation inputs", + ) 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 +624,69 @@ 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. + _start_the_fleet( + multiprocess, + self.__sim_producer, + n_workers, + (sim_monitor, mutex, simulation_error_event), + started_processes, + ) + _wait_for_workers(started_processes, simulation_error_event) + _close_the_fleet_down( + 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() + except KeyboardInterrupt: + _bring_the_fleet_down(started_processes, simulation_error_event) + self._interrupted = True - if not isinstance(error, KeyboardInterrupt): - raise error + # Handle error from the main process + except Exception: + _bring_the_fleet_down(started_processes, simulation_error_event) + self._interrupted = False + # Bare, so the handler's own line does not join the traceback. + raise + finally: + # Also best effort: this runs while an exception may be on its + # way out, and a cleanup that raises here would replace it. + _best_effort( + lambda: _stop_any_worker_still_running(started_processes), + "final worker shutdown", + ) 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 not _is_whole_number(n_workers): + 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 +695,23 @@ 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) - try: - mutex.acquire() + with _manager_mutex(mutex): if error_event.is_set(): # Runs in a worker process spawned via multiprocessing: # logging handlers configured in the main process are @@ -431,27 +726,45 @@ 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 + ) + # Same as the serial path: the pair is on disk, so a failure + # in the monitor call below must report itself and not the + # row that has just been committed. + inputs_json, outputs_json = "", "" sim_monitor.print_update_status() - finally: - mutex.release() - - except Exception: # pylint: disable=broad-except - mutex.acquire() - with open(self.error_file, "a", encoding="utf-8") as f: - f.write(inputs_json) - - # 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() + + except Exception: + # Set first, so a parent waiting on the join learns why. Best effort + # like everything below it: this is a manager proxy, the manager may + # already be gone, and reporting must not replace what it reports. + try: + error_event.set() + except Exception: # pylint: disable=broad-exception-caught + pass + details = traceback.format_exc() + + # The failure goes onto the inputs record rather than replacing it, + # the same shape the serial path writes. + record = _build_error_record(sim_idx, inputs_json, details) + + try: + with _manager_mutex(mutex): + 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 + + # 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. @@ -680,18 +993,18 @@ def __evaluate_flight_outputs(self, flight, sim_idx): export_item: getattr(flight, export_item) for export_item in self.export_list } - outputs_dict["index"] = sim_idx - if self.data_collector is not None: - additional_exports = {} for key, callback in self.data_collector.items(): try: - additional_exports[key] = callback(flight) + outputs_dict[key] = callback(flight) except Exception as e: raise ValueError( f"An error was encountered running 'data_collector' callback {key}. " ) from e - outputs_dict = outputs_dict | additional_exports + + # Last, so that the index a row is filed under is the one the run + # assigned even if the collector changed under a validated one. + outputs_dict["index"] = sim_idx return ( json.dumps(outputs_dict, cls=RocketPyEncoder, **self._export_config) + "\n" @@ -835,6 +1148,17 @@ def _check_data_collector(self, data_collector): ) for key, callback in data_collector.items(): + if not isinstance(key, str): + raise ValueError( + "Invalid 'data_collector' key! " + f"Keys must be strings, not {type(key).__name__}." + ) + if key in _RESERVED_RECORD_KEYS: + raise ValueError( + f"Invalid 'data_collector' key '{key}'! " + "That name is reserved for the record metadata that " + "pairs an inputs row with its outputs row." + ) if key in self.export_list: raise ValueError( "Invalid 'data_collector' key! " @@ -1377,7 +1701,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 +1931,425 @@ 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. Renumbering the rows would line the " + f"indices up without lining the seeds up: those rows came from the " + f"old sequential scheme, not from the per-index derivation this " + f"release would use for the same indices." + ) + 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 _is_whole_number(value): + """A Python or NumPy integer, and not a bool. + + ``type(value) in (int, np.integer)`` rejected every NumPy integer, because + ``type(np.int64(2))`` is ``np.int64``. ``True`` still has to go: it is an + ``int`` to ``isinstance`` and would quietly run one simulation. + """ + if isinstance(value, (bool, np.bool_)): + return False + return isinstance(value, (int, np.integer)) + + +def _validate_simulation_count(number_of_simulations): + """A count has to be a whole non-negative number, checked before any file. + + A float ran ``int(count)`` simulations and then failed the completeness + check with a range it could never have satisfied. + """ + if not _is_whole_number(number_of_simulations): + 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 +# One sleep per round, not per worker, so the fleet size does not set how soon +# an error is noticed. +_WORKER_POLL_INTERVAL = 0.05 + + +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 _record_failure(error_file, sim_idx, inputs_json, details): + """Append the failure being handled, the way the workers record theirs. + + Best effort, and it says so rather than going quiet: an unwritable error + file must not become the exception the caller sees in place of the one it + was called to record. + """ + try: + with open(error_file, "a", encoding="utf-8") as handle: + handle.write(_build_error_record(sim_idx, inputs_json, details)) + except Exception as reporting_error: # pylint: disable=broad-exception-caught + _say_so_without_raising( + f"The simulation failed and its error record could not be written: " + f"{reporting_error!r}" + ) + + +def _build_error_record(sim_idx, inputs_json, details): + """One failed simulation as a row: what it drew, and what went wrong. + + The traceback goes onto the inputs rather than replacing them, so the file + the run tells the user to read says both which inputs failed and why. + """ + try: + record = json.loads(inputs_json) if inputs_json else {"index": sim_idx} + except (TypeError, ValueError): + record = {"index": sim_idx} + record["error"] = details + return json.dumps(record) + "\n" + + +def _start_the_fleet(multiprocess, target, n_workers, args, started_processes): + """Start the workers, appending each as it starts. + + Appended one at a time so a ``start()`` that fails part way through leaves + the caller holding exactly those already running. + """ + for _ in range(n_workers): + sim_producer = multiprocess.Process(target=target, args=args) + sim_producer.start() + started_processes.append(sim_producer) + + +def _close_the_fleet_down(started_processes, error_event, error_file): + """Let the fleet finish its writes, stop the rest, then check the logs. + + The event asks workers to stop, it does not stop them, and without this + window one part way through a write is cut off and leaves exactly the torn + row the check reports. Not ``_bring_the_fleet_down``: that sets the event, + which on a clean run is what the crash check reads next. + """ + # Best effort, so that housekeeping cannot report itself in place of the + # worker result below, which is the authoritative verdict on the run. + _best_effort( + lambda: _wait_for_workers(started_processes, timeout=_WORKER_SHUTDOWN_GRACE), + "graceful worker wait", + ) + _best_effort( + lambda: _stop_any_worker_still_running(started_processes), + "forced worker shutdown", + ) + _fail_if_a_worker_did_not_finish(started_processes, error_event, error_file) + + +def _bring_the_fleet_down(started_processes, error_event): + """Stop everything, without raising over the failure being handled. + + Every step is best effort, not only the event: the manager may be the thing + that died, and a shutdown that raises would replace the failure that started + it. Bounded window to notice, then whatever is left gets stopped. + """ + _best_effort(error_event.set, "error notification") + _best_effort( + lambda: _wait_for_workers(started_processes, timeout=_WORKER_SHUTDOWN_GRACE), + "graceful worker wait", + ) + _best_effort( + lambda: _stop_any_worker_still_running(started_processes), + "forced worker shutdown", + ) + + +def _write_unfinished_inputs(error_file, inputs_json): + """Module level for the same reason as ``_record_simulation``: the run paths + are driven by stub objects in the tests, which carry no private methods.""" + with open(error_file, "a", encoding="utf-8") as f: + f.write(inputs_json) + + +def _say_so_without_raising(message): + """Report a secondary failure in a way that cannot become the primary one. + + ``warnings.warn`` raises when the caller has turned ``RuntimeWarning`` into + an error, which is exactly how a diagnostic ends up replacing the failure it + describes. The filter is overridden for this one warning only. + """ + try: + with warnings.catch_warnings(): + warnings.simplefilter("always", RuntimeWarning) + warnings.warn(message, RuntimeWarning, stacklevel=3) + except Exception: # pylint: disable=broad-exception-caught + pass + + +def _best_effort(action, description): + """Run one shutdown step, reporting a failure rather than raising it.""" + try: + action() + except Exception as cleanup_error: # pylint: disable=broad-exception-caught + _say_so_without_raising( + f"Worker cleanup failed during {description}: {cleanup_error!r}" + ) + + +@contextmanager +def _manager_mutex(mutex): + """Hold a manager lock without letting its release replace a failure. + + The proxy can die while the lock is held, and a raw release in ``finally`` + then becomes the exception the caller sees rather than the one already on + its way out. A release that fails with nothing in flight is still raised. + """ + mutex.acquire() + try: + yield + except BaseException: + _best_effort(mutex.release, "manager mutex release") + raise + else: + mutex.release() + + +def _read_error_event(error_event): + """Whether a worker reported an error, and what went wrong asking. + + An unreachable proxy is itself a reason to stop and to fail the run, so it + reads as reported rather than letting the exception past the crash list. + """ + try: + return bool(error_event.is_set()), None + except Exception as event_error: # pylint: disable=broad-exception-caught + return True, event_error + + +def _workers_that_crashed(started_processes): + """Those already known to have exited abnormally. + + ``join(timeout=0)`` first: an unjoined child has no exit code yet, so it + would read as ``None`` and pass for one still running. + """ + crashed = [] + for process in started_processes: + process.join(timeout=0) + if process.exitcode not in (None, 0): + crashed.append(process) + return crashed + + +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. + + The wait itself is one sleep per round rather than a blocking join on each + worker in turn, so how soon the event is noticed does not grow with the + fleet. Blocking 0.1 s per worker meant 100 of them delayed the next check by + 10 s. + """ + 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 _read_error_event(error_event)[0]: + break + # A worker killed outright sets no event. If it died holding the shared + # lock its siblings never return either, and only the unbounded wait + # has nothing else to end it. + if deadline is None and _workers_that_crashed(started_processes): + break + if deadline is None: + sleep(_WORKER_POLL_INTERVAL) + continue + remaining = deadline - monotonic() + if remaining <= 0: + break + sleep(min(_WORKER_POLL_INTERVAL, remaining)) + + # 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 _join_until(processes, deadline): + """Wait on the fleet against one clock rather than one clock each. + + A full grace per worker made the wait scale with the fleet: eight stubborn + ones could hold the parent for eight times what the grace period promised, + and twice over, once for terminate and once for kill. + """ + for process in processes: + process.join(timeout=max(0.0, deadline - monotonic())) + + +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() + _join_until(alive, monotonic() + 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() + _join_until(stubborn, monotonic() + 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 + ] + reported, event_error = _read_error_event(error_event) + if event_error is not None: + crashed.append(f"the worker error event became unavailable: {event_error!r}") + if reported 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). + """ + with _manager_mutex(mutex): + if not sim_monitor.keep_simulating(): + return None + return sim_monitor.increment() - 1 + + 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 5bb0598bf..994f38abc 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 @@ -79,6 +77,11 @@ class StochasticModel: "ensemble_member", ] + # Collections a child class builds itself and overwrites in its own + # ``dict_generator``. Drawing from them here would only move the stream by + # however many components happen to be attached, and the draw is discarded. + component_collections = () + def __init__(self, obj, seed=None, **kwargs): """ Initialize the StochasticModel class with validated input arguments. @@ -106,8 +109,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. @@ -147,7 +179,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): @@ -224,7 +256,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 @@ -295,7 +327,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 @@ -321,7 +353,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), ) @@ -348,7 +380,7 @@ def _validate_factors(self, input_name, input_value): If the input is not in a valid format. """ attribute_name = input_name.replace("_factor", "") - setattr(self, f"_{attribute_name}", getattr(self.obj, attribute_name)) + setattr(self, f"_{attribute_name}", self._nominal(attribute_name)) if isinstance(input_value, tuple): return self._validate_tuple_factor(input_name, input_value) @@ -587,6 +619,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. @@ -607,11 +654,13 @@ def dict_generator(self): """ generated_dict = {} for arg, value in self.__dict__.items(): + if arg in self.component_collections: + continue if isinstance(value, tuple): dist_sampler = value[-1] generated_dict[arg] = dist_sampler(value[0], value[1]) 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 33a364f18..b43e812dc 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, @@ -89,6 +91,16 @@ class StochasticRocket(StochasticModel): can not be a randomized. """ + # Overwritten in ``dict_generator`` below, so the base class must not draw + # from them. + component_collections = ( + "motors", + "aerodynamic_surfaces", + "rail_buttons", + "air_brakes", + "parachutes", + ) + def __init__( self, rocket, @@ -154,6 +166,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 +191,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 +240,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 +254,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 +476,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 +503,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 +681,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 +691,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..2aa2b14b5 --- /dev/null +++ b/tests/integration/simulation/test_monte_carlo_determinism.py @@ -0,0 +1,932 @@ +"""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``. + +Two tests cover the start methods, both unmarked. +``test_seed_derivation_is_start_method_invariant`` checks the derivation with a +top-level picklable target, and +``test_the_real_parallel_path_is_worker_invariant_under_every_start_method`` +drives the shipped parallel path. Both set the method rather than taking the +platform default, which is worth the seconds it costs: ``multiprocess`` +hard-codes ``fork`` on every POSIX platform, macOS and 3.14 included, so a +default-taking test would gate ``spawn`` on Windows and nothing anywhere else. +""" + +import json +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 + +# Parametrizing over start methods runs while tests are collected, and +# `multiprocess` is an optional extra, so skipping there would take the whole +# module down instead of skipping it. The tests themselves still importorskip. +try: + import multiprocess as _start_methods_from +except ImportError: + import multiprocessing as _start_methods_from + + +def _available_start_methods(): + """The start methods this platform supports, as ``multiprocess`` sees them. + + Asked of ``multiprocess`` rather than the standard library, because that is + what creates the workers and the two do not have to agree. + """ + supported = _start_methods_from.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`` + 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) + + multiprocess = pytest.importorskip("multiprocess") + context = multiprocess.get_context(start_method) + chunks = [(root_state, indices[0::2]), (root_state, indices[1::2])] + with context.Pool(2) as pool: + results = pool.starmap(_derive_index_seeds, chunks) + + combined = {} + for result in results: + combined.update(result) + assert combined == expected + assert sorted(combined) == indices + + +def _assert_the_same_environment_was_flown(runs, expected_indices, start_method): + """Every worker count flew index i with the same effective environment. + + This is the half the inputs file cannot show. It records + ``wind_velocity_x_factor``, which is the same for index i however the run + was executed even when the baseline it multiplies has drifted from one + simulation to the next. + """ + effective = { + label: _read_inputs_by_index(montecarlo.output_file) + for label, (montecarlo, _inputs) in runs.items() + } + for label, by_index in effective.items(): + assert sorted(by_index) == expected_indices, f"{label}: outputs are incomplete" + + for index in expected_indices: + reference = json.loads(effective["serial"][index]) + for key in _EFFECTIVE_ENVIRONMENT: + assert key in reference, f"{key} was not recorded" + assert reference["effective_wind_x"] != 0.0, ( + "the wind baseline is zero, so a compounding baseline cannot show" + ) + for label in ("parallel-2", "parallel-4"): + drawn = json.loads(effective[label][index]) + for key in _EFFECTIVE_ENVIRONMENT: + assert drawn[key] == reference[key], ( + f"{start_method}: {label} flew a different {key} at index " + f"{index}: {drawn[key]} against {reference[key]}" + ) + + +def _assert_the_run_is_complete(label, montecarlo, inputs, count): + """Every index written once, to both files, with nothing in the error log. + + The row counts are taken before anything is keyed by index: two workers + claiming the same index write two rows, and the second overwrites the first + in the dict, so a duplicate looks like a complete run. + """ + expected_indices = list(range(count)) + rows = _count_rows(montecarlo.input_file) + + assert sorted(inputs) == expected_indices, ( + f"{label}: indices {sorted(inputs)}, expected {expected_indices}" + ) + assert rows == count, ( + f"{label}: {rows} rows for {count} simulations, so an index was claimed " + f"more than once" + ) + assert _count_rows(montecarlo.output_file) == count, ( + f"{label}: the output rows do not match the simulations run" + ) + assert sorted(_read_inputs_by_index(montecarlo.output_file)) == expected_indices, ( + f"{label}: the outputs do not match the inputs" + ) + assert not os.path.getsize(montecarlo.error_file), ( + f"{label}: the run wrote to its error file" + ) + + +@pytest.fixture +def stochastic_environment_with_wind(example_spaceport_env): + """A stochastic environment whose wind is not zero. + + The shared ``stochastic_environment`` fixture sits on an Environment whose + ``wind_velocity_x`` is 0 at every altitude, and zero times any factor is + zero, so a baseline that compounds from one simulation to the next cannot + show up in it at all. Measured: with the baseline fix reverted, every + assertion in this file still passed. A wind that is actually blowing is + what makes the property testable. + """ + environment = Environment( + latitude=example_spaceport_env.latitude, + longitude=example_spaceport_env.longitude, + elevation=example_spaceport_env.elevation, + ) + environment.set_atmospheric_model( + type="custom_atmosphere", wind_u=12.0, wind_v=-7.0 + ) + return StochasticEnvironment( + environment=environment, + elevation=(1400, 10, "normal"), + wind_velocity_x_factor=(1.0, 0.05, "normal"), + wind_velocity_y_factor=(1.0, 0.05, "normal"), + ) + + +def _wind_x(flight): + """The wind the simulation actually flew with, not the factor drawn for it.""" + return float(flight.env.wind_velocity_x(0)) + + +def _wind_y(flight): + return float(flight.env.wind_velocity_y(0)) + + +def _elevation(flight): + return float(flight.env.elevation) + + +_EFFECTIVE_ENVIRONMENT = { + "effective_wind_x": _wind_x, + "effective_wind_y": _wind_y, + "effective_elevation": _elevation, +} + + +def _sampled_only(record): + """The recorded inputs with object identity stripped out. + + A ``Function``'s ``signature.hash`` and its serialised ``source`` encode the + object, not the value drawn for it, and an object built in another process + has a different one. Under ``fork`` they happen to agree because the child + inherits the parent's objects; under ``spawn`` and ``forkserver`` they + cannot. Measured on a real run: six fields differ across the boundary and + all six are these, while every sampled quantity matches exactly. + """ + flat = {} + + def walk(value, path=""): + if isinstance(value, dict): + for key, item in value.items(): + walk(item, f"{path}.{key}" if path else str(key)) + elif isinstance(value, list): + for position, item in enumerate(value): + walk(item, f"{path}[{position}]") + else: + flat[path] = value + + walk(record) + return { + key: value + for key, value in flat.items() + if "signature" not in key and not key.endswith(".source") + } + + +def _real_run_inputs(tmp_path, environment, rocket, flight, tag, **simulate_kwargs): + """Run a real Monte Carlo, no stub, and return the inputs keyed by index. + + Deliberately without the ``Flight`` stub. Stubbing is what confines the test + above to ``fork``: it replaces a module-level symbol in the parent, and a + ``spawn`` or ``forkserver`` child re-imports the module instead of inheriting + it. A real run has nothing that needs to cross the boundary except the + pickled MonteCarlo, which is the thing worth testing. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / tag), + environment=environment, + rocket=rocket, + flight=flight, + data_collector=_EFFECTIVE_ENVIRONMENT, + ) + montecarlo.simulate(**simulate_kwargs) + return montecarlo, _read_inputs_by_index(montecarlo.input_file) + + +@pytest.fixture +def restore_start_method(): + """Set the start method for one test and put it back afterwards.""" + multiprocess = pytest.importorskip("multiprocess") + original = multiprocess.get_start_method() + yield multiprocess + multiprocess.set_start_method(original, force=True) + + +@pytest.mark.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``. + + Sets each method rather than taking the platform default, because + ``multiprocess`` hard-codes ``fork`` on every POSIX platform including macOS + and 3.14, so a default-taking test gates ``spawn`` on Windows and nothing + else. Unmarked despite the cost: 26 s for the three, against 89 s for the + rest of this directory, and it is the only thing covering the path this + change exists to support. + """ + 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" + + +def test_appending_below_the_checkpoint_is_refused_and_changes_nothing( + tmp_path, stochastic_environment, stochastic_calisto, stochastic_flight +): + """``number_of_simulations`` is the total to reach, not a batch to add. + + Asking for fewer than the checkpoint already holds ran nothing and returned + success: every index it wanted was present, so the completeness check was + satisfied by simulations an earlier run had made. The caller was told three + while the file held five. + """ + montecarlo = MonteCarlo( + filename=str(tmp_path / "shrunk"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + montecarlo.simulate(number_of_simulations=3, random_seed=606) + before = ( + montecarlo.input_file.read_bytes(), + montecarlo.output_file.read_bytes(), + ) + + with pytest.raises(ValueError, match="already holds 3"): + montecarlo.simulate(number_of_simulations=2, append=True, random_seed=606) + + assert ( + montecarlo.input_file.read_bytes(), + montecarlo.output_file.read_bytes(), + ) == before, "the refusal touched the checkpoint it was protecting" + + +def test_appending_to_the_size_it_already_has_is_allowed( + tmp_path, stochastic_environment, stochastic_calisto, stochastic_flight +): + """The control. The guard must refuse only what is below the checkpoint, + not every append that adds no work.""" + montecarlo = MonteCarlo( + filename=str(tmp_path / "unchanged"), + environment=stochastic_environment, + rocket=stochastic_calisto, + flight=stochastic_flight, + ) + montecarlo.simulate(number_of_simulations=2, random_seed=606) + + montecarlo.simulate(number_of_simulations=2, append=True, random_seed=606) + + assert sorted(_read_inputs_by_index(montecarlo.input_file)) == [0, 1] 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..8040449ff --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_determinism.py @@ -0,0 +1,365 @@ +"""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, + _validate_simulation_count, + _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 + + +@pytest.mark.parametrize("count", [3, np.int32(3), np.int64(3), np.uint64(3)], ids=str) +def test_a_numpy_integer_is_a_valid_simulation_count(count): + """``type(count) in (int, np.integer)`` is False for every NumPy integer: + ``type(np.int64(3))`` is ``np.int64``, and ``np.integer`` is only its base. + A count read out of an array or a ``range`` product was refused.""" + _validate_simulation_count(count) + + +@pytest.mark.parametrize( + "count", [True, False, np.bool_(True), 3.0, "3", None], ids=str +) +def test_a_count_that_is_not_a_whole_number_is_still_refused(count): + """The control for the test above. ``True`` is the one that matters: it is + an ``int`` to ``isinstance`` and would quietly run one simulation.""" + with pytest.raises(TypeError): + _validate_simulation_count(count) + + +@pytest.mark.parametrize("wrapped", [False, True], ids=["sequence", "SeedSequence"]) +def test_mutating_the_caller_s_entropy_does_not_move_the_captured_root(wrapped): + """`SeedSequence` keeps a sequence entropy by reference, and so did the + capture, so a caller who reused and edited their list changed the children + of a run that had already read it. An int seed was never exposed to this. + """ + entropy = [1, 2, 3] + seed = np.random.SeedSequence(entropy) if wrapped else entropy + + runner = MonteCarlo.__new__(MonteCarlo) + MonteCarlo._MonteCarlo__capture_root_state(runner, seed) + before = _entropy(_child_seed(runner, 7)) + + entropy[0] = 999999 + + assert _entropy(_child_seed(runner, 7)) == before 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..1e0e6c7de --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_log_integrity.py @@ -0,0 +1,558 @@ +"""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 json +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") as raised: + mc._check_the_checkpoint_supports_appending(inputs, outputs, 3) + + # The message used to offer renumbering as an alternative to re-running, + # which lines the indices up and leaves the seeds behind: those rows came + # from the old sequential scheme, not from this release's per-index one. + assert "Renumbering" in str(raised.value) + assert "without lining the seeds up" in str(raised.value) + + +class _Collector: + """A model with only what the collector checks and the writer touch.""" + + export_list = ("apogee",) + + def __init__(self, data_collector=None): + self.data_collector = data_collector + self._export_config = {} + + +def _outputs(data_collector, sim_idx): + model = _Collector(data_collector) + flight = types.SimpleNamespace(apogee=100.0) + return json.loads( + mc.MonteCarlo._MonteCarlo__evaluate_flight_outputs(model, flight, sim_idx) + ) + + +def test_a_collector_cannot_relabel_the_row_it_is_attached_to(): + """`index` pairs an inputs row with its outputs row. The custom fields used + to be merged over it, so a collector could file a row under a different + simulation than the one that produced it.""" + labels = iter([1, 0]) + collector = {"index": lambda _flight: next(labels)} + + assert _outputs(collector, 0)["index"] == 0 + assert _outputs(collector, 1)["index"] == 1 + + +def test_a_permutation_of_valid_indices_would_pass_every_other_check(): + """Why the one above matters more than a malformed value would. + + A collector returning `-1` writes a row the completeness check rejects. One + returning a permutation writes the same index set with the same counts, so + nothing downstream can tell the outputs are on the wrong simulations. + """ + labels = iter([1, 0]) + rows = [_outputs({"custom": lambda _f: next(labels)}, i) for i in (0, 1)] + + assert sorted(r["index"] for r in rows) == [0, 1] + assert [r["custom"] for r in rows] == [1, 0], "the collector still runs" + + +@pytest.mark.parametrize( + "key, expected", + [("index", "reserved"), (7, "must be strings"), (None, "must be strings")], + ids=["reserved-name", "int-key", "none-key"], +) +def test_a_collector_key_that_cannot_be_written_is_refused(key, expected): + """A non-string key does not survive the round trip these files exist for. + + ``json.dumps`` stringifies it on the way out, so ``7`` comes back as + ``"7"``, ``None`` as ``"null"``. Worse, it can collide: a collector holding + both ``1`` and ``"1"`` writes ``{"1": ..., "1": ...}``, and reading that + back leaves one column where there were two. + """ + with pytest.raises(ValueError, match=expected): + mc.MonteCarlo._check_data_collector(_Collector(), {key: lambda _f: 0}) + + +def test_a_collector_changed_after_construction_is_checked_again(tmp_path): + """`data_collector` is public and mutable, so validating it once in + ``__init__`` is not enough. The check has to run before ``__setup_files`` + opens the logs "w+", or a rejected run destroys the previous one.""" + inputs = tmp_path / "inputs.txt" + outputs = tmp_path / "outputs.txt" + inputs.write_text("previous input\n", encoding="utf-8") + outputs.write_text("previous output\n", encoding="utf-8") + runner = types.SimpleNamespace( + input_file=inputs, + output_file=outputs, + export_list=("apogee",), + data_collector={"index": lambda _flight: 0}, + _check_data_collector=lambda collector: mc.MonteCarlo._check_data_collector( + runner, collector + ), + ) + + with pytest.raises(ValueError, match="reserved"): + mc.MonteCarlo.simulate(runner, number_of_simulations=1, random_seed=42) + + assert inputs.read_text(encoding="utf-8") == "previous input\n" + assert outputs.read_text(encoding="utf-8") == "previous output\n" 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..8904b07f7 --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_worker_exit.py @@ -0,0 +1,505 @@ +"""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 traceback +import types +from contextlib import contextmanager +from time import monotonic, sleep + +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" + ) + + +class _Crashed: + """Killed outright: gone, non-zero exit code, no event set.""" + + def __init__(self): + self.exitcode = 7 + + def join(self, *_a, **_k): + pass + + def is_alive(self): + return False + + +class _BlockedForever: + """Waiting on the lock the crashed worker still holds.""" + + def __init__(self, give_up_after=50): + self.exitcode = None + self.joins = 0 + self._give_up_after = give_up_after + + def join(self, *_a, **_k): + self.joins += 1 + if self.joins > self._give_up_after: + raise AssertionError( + "the parent is still waiting on a worker that a dead sibling " + "has blocked, and nothing else will end this wait" + ) + + def is_alive(self): + return True + + +def test_the_parent_stops_waiting_when_a_worker_dies_holding_the_lock(): + """A worker killed outright sets no event, so the wait had only + ``is_alive`` to end it, and a sibling blocked on the lock it held kept that + true forever. The exit-code check downstream was never reached.""" + crashed, blocked = _Crashed(), _BlockedForever() + + mc._wait_for_workers([crashed, blocked], _Event()) + + assert blocked.is_alive(), "the blocked worker is meant to still be running" + + +def test_the_shutdown_window_is_not_cut_short_by_a_worker_already_known_dead(): + """The grace period exists so the survivors can finish their writes. It is + bounded by its own timeout, so a crash must not end it early. + + On elapsed time rather than on how many times a worker was joined: that + counted the old per-worker blocking join, which is the mechanism and not the + behaviour. The wait can only overshoot the deadline, never undershoot it, so + a floor well under the timeout is safe on a coarse clock. + """ + crashed, blocked = _Crashed(), _BlockedForever(give_up_after=10**6) + + started = monotonic() + mc._wait_for_workers([crashed, blocked], timeout=0.3) + elapsed = monotonic() - started + + assert elapsed >= 0.2, f"the grace period returned after {elapsed:.3f}s" + + +class _Stubborn: + """A worker that sits through terminate and kill, recording its waits.""" + + def __init__(self, waits): + self.exitcode = None + self._waits = waits + + def join(self, timeout=None, **_k): + self._waits.append(timeout) + if timeout: + sleep(timeout) + + def is_alive(self): + return True + + def terminate(self): + pass + + def kill(self): + pass + + +def _granted_shutting_down(fleet_size, grace): + """How long the fleet was granted in total, across both phases.""" + waits = [] + mc._stop_any_worker_still_running( + [_Stubborn(waits) for _ in range(fleet_size)], grace=grace + ) + assert len(waits) == 2 * fleet_size, "every worker is still waited on" + return sum(waits) + + +def test_shutdown_does_not_grow_with_the_fleet(): + """Each worker used to get the full grace to itself, so the wait scaled with + the fleet: six stubborn workers held the parent for six grace periods per + phase rather than one. + + Compares two fleet sizes rather than checking either against the grace. + A single fleet has to be measured against a constant, and the slack that + needs is the machine's timer granularity, which on Windows is 15 ms against + a 200 ms grace. Both measurements carry the same granularity, so comparing + them cancels it. + """ + grace = 0.2 + + alone = _granted_shutting_down(1, grace) + crowd = _granted_shutting_down(6, grace) + + assert crowd < alone * 2, ( + f"six workers were granted {crowd:.2f}s against {alone:.2f}s for one, " + f"so the wait is still scaling with the fleet" + ) + + +def test_the_parent_does_not_repeat_its_own_frame_in_the_traceback( + parallel_runner, monkeypatch +): + """``raise error`` names the exception again, so the handler's line joins the + traceback and the reader walks past it. The serial path already re-raises + bare; this is the parent doing the same. + + On the duplicate rather than the original frame: the failing call survives + either spelling, so asserting it is there passes both. + """ + + def start_then_fail(self): + raise OSError("cannot start") + + monkeypatch.setattr(_Process, "start", start_then_fail) + + with pytest.raises(OSError) as raised: + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=2) + + frames = [frame.name for frame in traceback.extract_tb(raised.value.__traceback__)] + assert "start_then_fail" in frames, frames + assert frames.count("__run_in_parallel") == 1, frames + + +class _NeverFinishes: + """Alive throughout, and costly to join, as a real blocked worker is.""" + + def __init__(self): + self.exitcode = None + + def join(self, timeout=None, **_k): + if timeout: + sleep(timeout) + + def is_alive(self): + return True + + +class _SetMidRound: + """Not set when a round begins, set while the parent is inside it.""" + + def __init__(self): + self.checks = 0 + + def is_set(self): + self.checks += 1 + return self.checks > 1 + + +@pytest.mark.parametrize("fleet", [1, 24], ids=["one", "twenty_four"]) +def test_noticing_an_error_does_not_get_slower_with_a_bigger_fleet(fleet): + """Joining every worker for 0.1s before rechecking the event made the delay + the fleet size times that: 24 workers took 2.4s to notice. One sleep per + round instead, so the cost is the round and not the fleet. + + Both sizes are measured against the same ceiling rather than against each + other, because the point is that neither depends on the count. + """ + started = monotonic() + mc._wait_for_workers([_NeverFinishes() for _ in range(fleet)], _SetMidRound()) + elapsed = monotonic() - started + + assert elapsed < 0.5, f"{fleet} workers delayed the check by {elapsed:.2f}s" + + +class _OriginalFailure(RuntimeError): + """What the caller should end up seeing.""" + + +class _CleanupFailure(RuntimeError): + """What must not take its place.""" + + +def test_the_final_cleanup_does_not_replace_the_failure_on_its_way_out( + parallel_runner, monkeypatch +): + """The handler's own cleanup is best effort, but ``finally`` runs again on + the way out and was calling the same thing raw. A failure there landed on + the caller instead of the one being re-raised.""" + monkeypatch.setattr( + mc, "_start_the_fleet", _raiser(_OriginalFailure, "cannot start") + ) + monkeypatch.setattr( + mc, + "_stop_any_worker_still_running", + _raiser(_CleanupFailure, "cannot clean up"), + ) + + with pytest.warns(RuntimeWarning, match="final worker shutdown"): + with pytest.raises(_OriginalFailure, match="cannot start"): + mc.MonteCarlo._MonteCarlo__run_in_parallel(parallel_runner, n_workers=2) + + +def _raiser(exception, message): + def raise_it(*_args, **_kwargs): + raise exception(message) + + return raise_it + + +class _ExitedWorker: + """A worker that already left, and did not leave cleanly.""" + + name = "worker-0" + exitcode = 7 + + def is_alive(self): + return False + + def join(self, timeout=None): + pass + + +class _SetEvent: + def is_set(self): + return True + + +class _DeadEvent: + def is_set(self): + raise ConnectionResetError("manager is gone") + + +def test_cleanup_failure_does_not_hide_which_worker_crashed(monkeypatch): + """The close path ran the wait and the forced stop raw, so a failure there + reached the caller before anything looked at the exit codes. The worker + result is the verdict on the run; housekeeping is not.""" + monkeypatch.setattr( + mc, "_stop_any_worker_still_running", _raiser(OSError, "cleanup failed") + ) + + with pytest.warns(RuntimeWarning, match="forced worker shutdown"): + with pytest.raises(RuntimeError, match="exited with 7"): + mc._close_the_fleet_down([_ExitedWorker()], _SetEvent(), "errors.txt") + + +def test_an_unreachable_event_is_reported_rather_than_raised(): + """``error_event.is_set() or crashed`` asked the proxy first, so a dead + manager threw before the crash list was ever read.""" + with pytest.raises(RuntimeError) as raised: + mc._fail_if_a_worker_did_not_finish( + [_ExitedWorker()], _DeadEvent(), "errors.txt" + ) + + assert "exited with 7" in str(raised.value) + assert "became unavailable" in str(raised.value) 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..9aa4e0f03 --- /dev/null +++ b/tests/unit/simulation/test_monte_carlo_worker_failures.py @@ -0,0 +1,518 @@ +"""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 traceback +import warnings +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" + + +class _OneThenStop: + """A monitor that allows exactly one simulation, then whatever is asked.""" + + def __init__(self, on_update=None): + self.count = 0 + self._on_update = on_update + + def keep_simulating(self): + return self.count < 1 + + def increment(self): + self.count += 1 + return self.count + + def print_update_status(self): + if self._on_update is not None: + self._on_update() + + def print_final_status(self): + pass + + +def _serial_runner(tmp_path, monkeypatch, monitor, **overrides): + """A stand-in carrying only what ``__run_in_serial`` touches.""" + runner = _worker(tmp_path, **overrides) + runner._error_file = runner.error_file + runner._initial_sim_idx = 0 + runner.number_of_simulations = 1 + runner._interrupted = False + runner._MonteCarlo__keep_the_inputs_that_did_not_finish = lambda payload: ( + runner.error_file.open("a", encoding="utf-8").write(payload) + ) + monkeypatch.setattr(mc, "_SimMonitor", lambda **_kwargs: monitor) + return runner + + +def test_a_serial_failure_records_the_traceback_not_only_the_inputs( + tmp_path, monkeypatch +): + """The worker path writes `{index, ...inputs, error: traceback}`. Serial + wrote the inputs alone, so a failure after sampling named which inputs + failed and never why, in the file the run points the user at.""" + runner = _serial_runner( + tmp_path, + monkeypatch, + _OneThenStop(), + _MonteCarlo__evaluate_flight_outputs=_raise, + ) + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + record = json.loads(runner.error_file.read_text(encoding="utf-8").splitlines()[0]) + assert "_Boom" in record["error"] + assert "injected" in record["error"] + + +def test_a_serial_failure_does_not_repeat_the_handler_frame(tmp_path, monkeypatch): + """`raise error` names the exception again, so the handler's own line joins + the traceback and the reader walks past it to reach the real one. A bare + `raise` leaves the frame it came from. + + On the duplicate rather than on the original frame: the failing call + survives either way, so asserting it is there passes both spellings. + """ + runner = _serial_runner( + tmp_path, + monkeypatch, + _OneThenStop(), + _MonteCarlo__evaluate_flight_outputs=_raise, + ) + + with pytest.raises(_Boom) as raised: + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + frames = [frame.name for frame in traceback.extract_tb(raised.value.__traceback__)] + assert "_raise" in frames, frames + assert frames.count("__run_in_serial") == 1, frames + + +def test_a_committed_row_is_not_reported_as_unfinished(tmp_path, monkeypatch): + """The pair is on disk before the progress call. A failure there used to + append those same inputs to the error file, so one simulation appeared in + both the inputs log and the failures.""" + runner = _serial_runner(tmp_path, monkeypatch, _OneThenStop(on_update=_raise)) + + with pytest.raises(_Boom): + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + assert runner.input_file.read_text(encoding="utf-8").strip() == "{}" + record = json.loads(runner.error_file.read_text(encoding="utf-8").splitlines()[0]) + assert record == {"index": 0, "error": record["error"]} + assert "_Boom" in record["error"] + + +class _ClaimOnceThenFail: + """Lets one simulation through, then fails the progress call.""" + + def __init__(self): + self.count = 0 + + def keep_simulating(self): + return self.count < 1 + + def increment(self): + self.count += 1 + return self.count + + def print_update_status(self): + raise _Boom("progress failed") + + +def test_a_worker_progress_failure_does_not_repeat_committed_inputs(tmp_path): + """The serial path clears its payload once the pair is on disk. The worker + kept it, so a failure in the progress call wrote the same sampled inputs to + the error file and one simulation appeared in the logs and the failures.""" + worker = _worker( + tmp_path, + _MonteCarlo__evaluate_flight_inputs=lambda index: '{"index": 0, "drew": 42}\n', + _MonteCarlo__evaluate_flight_outputs=lambda flight, index: '{"index": 0}\n', + ) + mutex, event = _RecordingMutex(), _Event() + + with pytest.raises(_Boom, match="progress failed"): + mc.MonteCarlo._MonteCarlo__sim_producer( + worker, _ClaimOnceThenFail(), mutex, event + ) + + committed = json.loads( + worker.input_file.read_text(encoding="utf-8").splitlines()[0] + ) + failure = json.loads(worker.error_file.read_text(encoding="utf-8").splitlines()[0]) + + assert committed == {"index": 0, "drew": 42} + assert "drew" not in failure, failure + assert "_Boom" in failure["error"] + + +def test_a_serial_reporting_failure_does_not_replace_the_original( + tmp_path, monkeypatch +): + """The worker guards its own error write. Serial did not, so an unwritable + error file raised OSError in place of the failure it was recording.""" + runner = _serial_runner( + tmp_path, + monkeypatch, + _OneThenStop(), + _MonteCarlo__evaluate_flight_outputs=_raise, + ) + real_open = open + + def refuse_the_error_file(path, *args, **kwargs): + if str(path) == str(runner.error_file): + raise OSError("error disk unavailable") + return real_open(path, *args, **kwargs) + + monkeypatch.setattr("builtins.open", refuse_the_error_file) + + with pytest.warns(RuntimeWarning, match="could not be written"): + with pytest.raises(_Boom, match="injected"): + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + +def test_a_shutdown_failure_does_not_replace_the_error_it_is_handling(monkeypatch): + """``_bring_the_fleet_down`` promised not to raise over the failure being + handled, but only the event was guarded. The two waits below it were not.""" + monkeypatch.setattr(mc, "_wait_for_workers", _raise) + + with pytest.warns(RuntimeWarning, match="graceful worker wait"): + mc._bring_the_fleet_down([], _Event()) + + +def test_reporting_a_failure_cannot_escape_when_warnings_are_errors( + tmp_path, monkeypatch +): + """The guard around the error write is only best effort under the default + filter. Turn RuntimeWarning into an error, as a strict application or test + run does, and the diagnostic becomes the exception it was describing.""" + runner = _serial_runner( + tmp_path, + monkeypatch, + _OneThenStop(), + _MonteCarlo__evaluate_flight_outputs=_raise, + ) + real_open = open + + def refuse_the_error_file(path, *args, **kwargs): + if str(path) == str(runner.error_file): + raise OSError("error disk unavailable") + return real_open(path, *args, **kwargs) + + monkeypatch.setattr("builtins.open", refuse_the_error_file) + + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + + with pytest.raises(_Boom, match="injected"): + mc.MonteCarlo._MonteCarlo__run_in_serial(runner) + + +class _MutexDiesAfter(_RecordingMutex): + """Releases cleanly a few times, then the manager connection goes. + + A mutex that fails on every release is the wrong model: the first clean + claim would fail on its own release, before there is anything to mask. + """ + + def __init__(self, healthy_releases=1): + super().__init__() + self._left = healthy_releases + + def release(self): + super().release() + self._left -= 1 + if self._left < 0: + raise BrokenPipeError("[Errno 32] Broken pipe") + + +def test_a_dying_mutex_does_not_replace_the_simulation_failure(tmp_path): + """The release in the critical section used to run raw in ``finally``, so a + manager that died while the lock was held handed the caller its own + BrokenPipeError and left the real failure as context.""" + worker = _worker(tmp_path) + monitor = _ClaimOnceThenFail() + event = _Event() + + with pytest.warns(RuntimeWarning, match="manager mutex release"): + with pytest.raises(_Boom, match="progress failed"): + mc.MonteCarlo._MonteCarlo__sim_producer( + worker, monitor, _MutexDiesAfter(healthy_releases=1), event + ) + + +def test_a_dying_mutex_does_not_replace_the_claim_failure(): + """Same shape one level down, where the claim itself is what fails.""" + + class _ClaimRaises: + def keep_simulating(self): + raise _Boom("claim failed") + + with pytest.warns(RuntimeWarning, match="manager mutex release"): + with pytest.raises(_Boom, match="claim failed"): + mc._claim_next_index(_ClaimRaises(), _MutexDiesAfter(healthy_releases=0)) + + +def test_a_release_that_fails_on_its_own_is_still_raised(): + """The other half. With nothing in flight a broken release is the failure, + not a warning, so the run does not carry on against a dead manager.""" + + class _Fine: + def keep_simulating(self): + return True + + def increment(self): + return 1 + + with pytest.raises(BrokenPipeError): + mc._claim_next_index(_Fine(), _MutexDiesAfter(healthy_releases=0)) diff --git a/tests/unit/stochastic/test_custom_sampler.py b/tests/unit/stochastic/test_custom_sampler.py index ae3d906ba..286d97cd3 100644 --- a/tests/unit/stochastic/test_custom_sampler.py +++ b/tests/unit/stochastic/test_custom_sampler.py @@ -1,3 +1,5 @@ +import subprocess +import sys from types import SimpleNamespace import numpy as np @@ -361,9 +363,6 @@ def test_a_non_sampler_is_refused_even_under_optimisation(): def test_the_refusal_survives_python_dash_o(): """The mechanism, not just the behaviour: run it in a child with -O and check the exception still arrives.""" - import subprocess - import sys - program = ( "from types import SimpleNamespace;" "from rocketpy.stochastic.stochastic_model import StochasticModel;" 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..e68720622 --- /dev/null +++ b/tests/unit/stochastic/test_stochastic_rocket_seeding.py @@ -0,0 +1,222 @@ +"""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" + + +def _mass_drawn_with(rocket_factory, stochastic_parachutes, seed=42): + rocket = rocket_factory() + for parachute in stochastic_parachutes: + rocket.add_parachute(parachute) + rocket._set_stochastic(seed) + return next(rocket.dict_generator())["mass"] + + +@pytest.mark.parametrize("attached", [0, 1, 2]) +def test_attaching_components_does_not_move_the_rocket_body_stream( + attached, + stochastic_calisto, + stochastic_main_parachute, + stochastic_drogue_parachute, +): + """The base generator walked the whole instance and drew from + ``parachutes`` too, using the model's own generator since this branch + started seeding list choices. The subclass then discards that draw, so the + only thing it did was shift every later draw by however many components + happened to be attached. Two chutes changed the sampled mass. + + One is not enough to catch it: ``integers(1)`` has a single outcome and + NumPy returns it without consuming any state. + """ + available = [stochastic_main_parachute, stochastic_drogue_parachute] + + def bare(): + stochastic_calisto.parachutes = [] + stochastic_calisto.air_brakes = [] + return stochastic_calisto + + alone = _mass_drawn_with(bare, []) + with_components = _mass_drawn_with(bare, available[:attached]) + + assert with_components == alone + + +def test_the_discarded_component_lists_are_still_reported_empty( + stochastic_calisto, stochastic_main_parachute +): + """Skipping them must not change what ``dict_generator`` yields.""" + stochastic_calisto.parachutes = [] + stochastic_calisto.add_parachute(stochastic_main_parachute) + stochastic_calisto._set_stochastic(42) + + generated = next(stochastic_calisto.dict_generator()) + + assert generated["parachutes"] == [] + assert generated["air_brakes"] == []