From 9d5164807637a7ac3e27bd504a5c64dc5b65ef80 Mon Sep 17 00:00:00 2001 From: Lucas de Oliveira Prates Date: Sun, 13 Apr 2025 21:14:19 +0200 Subject: [PATCH 1/4] ENH: first implementation of CustomSampler --- rocketpy/__init__.py | 1 + rocketpy/stochastic/__init__.py | 1 + rocketpy/stochastic/custom_sampler.py | 38 ++++ rocketpy/stochastic/stochastic_model.py | 56 +++++- test_custom_sampler.ipynb | 232 ++++++++++++++++++++++++ 5 files changed, 327 insertions(+), 1 deletion(-) create mode 100644 rocketpy/stochastic/custom_sampler.py create mode 100644 test_custom_sampler.ipynb diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index 1e7731073..f99a70f28 100644 --- a/rocketpy/__init__.py +++ b/rocketpy/__init__.py @@ -44,6 +44,7 @@ from .sensors import Accelerometer, Barometer, GnssReceiver, Gyroscope from .simulation import Flight, MonteCarlo, MultivariateRejectionSampler from .stochastic import ( + CustomSampler, StochasticAirBrakes, StochasticEllipticalFins, StochasticEnvironment, diff --git a/rocketpy/stochastic/__init__.py b/rocketpy/stochastic/__init__.py index b1e146246..ffadfaaaf 100644 --- a/rocketpy/stochastic/__init__.py +++ b/rocketpy/stochastic/__init__.py @@ -5,6 +5,7 @@ associated with each input parameter. """ +from .custom_sampler import CustomSampler from .stochastic_aero_surfaces import ( StochasticAirBrakes, StochasticEllipticalFins, diff --git a/rocketpy/stochastic/custom_sampler.py b/rocketpy/stochastic/custom_sampler.py new file mode 100644 index 000000000..c13887060 --- /dev/null +++ b/rocketpy/stochastic/custom_sampler.py @@ -0,0 +1,38 @@ +""" +Provides an abstract class so that users can build custom samplers upon +""" + +from abc import ABC, abstractmethod + + +class CustomSampler(ABC): + """Abstract subclass for user defined samplers""" + + @abstractmethod + def sample(self, n_samples=1): + """Generates n samples from the custom distribution + + Parameters + ---------- + n_samples : int, optional + Numbers of samples to be generated + + Returns + ------- + sample_list : list + A list with n_samples elements, each of which is a valid sample + """ + + @abstractmethod + def reset_seed(self, seed=None): + """Resets the seeds of all associated stochastic generators + + Parameters + ---------- + seed : int, optional + Seed for the random number generator. The default is None + + Returns + ------- + None + """ diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index fe47a252a..103ec7175 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -8,6 +8,7 @@ import numpy as np from rocketpy.mathutils.function import Function +from rocketpy.stochastic.custom_sampler import CustomSampler from ..tools import get_distribution @@ -96,6 +97,10 @@ def _set_stochastic(self, seed=None): attr_value = self._validate_list(input_name, input_value) elif isinstance(input_value, (int, float)): attr_value = self._validate_scalar(input_name, input_value) + elif isinstance(input_value, CustomSampler): + attr_value = self._validate_custom_sampler( + input_name, input_value, seed + ) else: raise AssertionError( f"'{input_name}' must be a tuple, list, int, or float" @@ -436,6 +441,33 @@ def _validate_positive_int_list(self, input_name, input_value): isinstance(member, int) and member >= 0 for member in input_value ), f"`{input_name}` must be a list of positive integers" + def _validate_custom_sampler(self, input_name, sampler, seed=None): + """ + Validate a custom sampler. + + Parameters + ---------- + input_name : str + Name of the input argument. + sampler : CustomSampler object + Custom sampler provided by the user + seed : int, optional + Seed for the random number generator. The default is None + + Raises + ------ + AssertionError + If the input is not in a valid format. + """ + try: + sampler.reset_seed(seed) + except RuntimeError as e: + raise RuntimeError( + f"An error occurred in the 'reset_seed' of {input_name} CustomSampler" + ) from e + + return sampler + def _validate_airfoil(self, airfoil): """ Validate airfoil input. @@ -490,9 +522,17 @@ def dict_generator(self): generated_dict = {} for arg, value in self.__dict__.items(): if isinstance(value, tuple): - generated_dict[arg] = value[-1](value[0], value[1]) + 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 + elif isinstance(value, CustomSampler): + try: + generated_dict[arg] = value.sample(n_samples=1)[0] + except RuntimeError as e: + raise RuntimeError( + f"An error occurred in the 'sample' of {arg} CustomSampler" + ) from e self.last_rnd_dict = generated_dict yield generated_dict @@ -527,6 +567,12 @@ def format_attribute(attr, value): f"{nominal_value:.5f} ± " f"{std_dev:.5f} ({dist_func.__name__})" ) + elif isinstance(value, CustomSampler): + sampler_name = type(value).__name__ + return ( + f"\t{attr.ljust(max_str_length)} " + f"\t{sampler_name.ljust(max_str_length)} " + ) return None attributes = {k: v for k, v in self.__dict__.items() if not k.startswith("_")} @@ -550,6 +596,9 @@ def format_attribute(attr, value): list_attributes = [ attr for attr, val in items if isinstance(val, list) and len(val) > 1 ] + custom_attributes = [ + attr for attr, val in items if isinstance(val, CustomSampler) + ] if constant_attributes: report.append("\nConstant Attributes:") @@ -568,5 +617,10 @@ def format_attribute(attr, value): report.extend( format_attribute(attr, attributes[attr]) for attr in list_attributes ) + if custom_attributes: + report.append("\nStochastic Attributes with Custom user samplers:") + report.extend( + format_attribute(attr, attributes[attr]) for attr in custom_attributes + ) print("\n".join(filter(None, report))) diff --git a/test_custom_sampler.ipynb b/test_custom_sampler.ipynb new file mode 100644 index 000000000..c831f0cf9 --- /dev/null +++ b/test_custom_sampler.ipynb @@ -0,0 +1,232 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "# We import these lines for debugging purposes, only works on Jupyter Notebook\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "import numpy as np\n", + "from rocketpy import SolidMotor, StochasticSolidMotor, CustomSampler" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "class GaussianMixture(CustomSampler):\n", + "\n", + " def __init__(self, means_tuple, sd_tuple, prob_tuple, seed = None):\n", + " \"\"\" Creates a sampler for a mixture of two Gaussians\n", + "\n", + " Parameters\n", + " ----------\n", + " means_tuple : 2-tuple\n", + " 2-Tuple that contains the mean of each normal distribution of the\n", + " mixture\n", + " sd_tuple : 2-tuple\n", + " 2-Tuple that contains the sd of each normal distribution of the\n", + " mixture\n", + " prob_tuple : 2-tuple\n", + " 2-Tuple that contains the probability of each normal distribution of the\n", + " mixture\n", + " \"\"\"\n", + " np.random.default_rng(seed)\n", + " self.means_tuple = means_tuple\n", + " self.sd_tuple = sd_tuple\n", + " self.prob_tuple = prob_tuple\n", + " \n", + " def sample(self, n_samples = 1):\n", + " sample_list = [0] * n_samples\n", + " mixture_id_list = np.random.binomial(1, self.prob_tuple[0], n_samples)\n", + " for i, mixture_id in enumerate(mixture_id_list):\n", + " if mixture_id:\n", + " sample_list[i] = np.random.normal(self.means_tuple[0], self.sd_tuple[0])\n", + " else:\n", + " sample_list[i] = np.random.normal(self.means_tuple[1], self.sd_tuple[1])\n", + "\n", + " return sample_list\n", + " \n", + " def reset_seed(self, seed=None):\n", + " np.random.default_rng(seed)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "total_impulse_sampler = GaussianMixture([6000, 7000], [200, 200], [0.5, 0.5])" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "motor = SolidMotor(\n", + " thrust_source=\"data/motors/cesaroni/Cesaroni_M1670.eng\",\n", + " dry_mass=1.815,\n", + " dry_inertia=(0.125, 0.125, 0.002),\n", + " nozzle_radius=33 / 1000,\n", + " grain_number=5,\n", + " grain_density=1815,\n", + " grain_outer_radius=33 / 1000,\n", + " grain_initial_inner_radius=15 / 1000,\n", + " grain_initial_height=120 / 1000,\n", + " grain_separation=5 / 1000,\n", + " grains_center_of_mass_position=0.397,\n", + " center_of_dry_mass_position=0.317,\n", + " nozzle_position=0,\n", + " burn_time=3.9,\n", + " throat_radius=11 / 1000,\n", + " coordinate_system_orientation=\"nozzle_to_combustion_chamber\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "ename": "AssertionError", + "evalue": "'total_impulse' must be a tuple, list, int, or float", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[26], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m stochastic_motor \u001b[38;5;241m=\u001b[39m \u001b[43mStochasticSolidMotor\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 2\u001b[0m \u001b[43m \u001b[49m\u001b[43msolid_motor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmotor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 3\u001b[0m \u001b[43m \u001b[49m\u001b[43mburn_start_time\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m0.1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mbinomial\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 4\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrains_center_of_mass_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.001\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_density\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m50\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 6\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_separation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 7\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_height\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 8\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_inner_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.375\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 9\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_outer_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.375\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 10\u001b[0m \u001b[43m \u001b[49m\u001b[43mtotal_impulse\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtotal_impulse_sampler\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 11\u001b[0m \u001b[43m \u001b[49m\u001b[43mthroat_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.5\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 12\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.5\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 13\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.001\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 14\u001b[0m \u001b[43m)\u001b[49m\n\u001b[1;32m 15\u001b[0m \u001b[38;5;66;03m#stochastic_motor.visualize_attributes()\u001b[39;00m\n", + "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_solid_motor.py:149\u001b[0m, in \u001b[0;36mStochasticSolidMotor.__init__\u001b[0;34m(self, solid_motor, thrust_source, total_impulse, burn_start_time, burn_out_time, dry_mass, dry_inertia_11, dry_inertia_22, dry_inertia_33, dry_inertia_12, dry_inertia_13, dry_inertia_23, nozzle_radius, grain_number, grain_density, grain_outer_radius, grain_initial_inner_radius, grain_initial_height, grain_separation, grains_center_of_mass_position, center_of_dry_mass_position, nozzle_position, throat_radius)\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\n\u001b[1;32m 69\u001b[0m \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m 70\u001b[0m solid_motor,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 92\u001b[0m throat_radius\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 93\u001b[0m ):\n\u001b[1;32m 94\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Initializes the Stochastic Solid Motor class.\u001b[39;00m\n\u001b[1;32m 95\u001b[0m \n\u001b[1;32m 96\u001b[0m \u001b[38;5;124;03m See Also\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 147\u001b[0m \u001b[38;5;124;03m Radius of the throat in the motor in meters.\u001b[39;00m\n\u001b[1;32m 148\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 149\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[1;32m 150\u001b[0m \u001b[43m \u001b[49m\u001b[43msolid_motor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 151\u001b[0m \u001b[43m \u001b[49m\u001b[43mthrust_source\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mthrust_source\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 152\u001b[0m \u001b[43m \u001b[49m\u001b[43mtotal_impulse\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtotal_impulse\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 153\u001b[0m \u001b[43m \u001b[49m\u001b[43mburn_start_time\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mburn_start_time\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 154\u001b[0m \u001b[43m \u001b[49m\u001b[43mburn_out_time\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mburn_out_time\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 155\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_mass\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_mass\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 156\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_11\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_11\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 157\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_22\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_22\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 158\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_33\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_33\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 159\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_12\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_12\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 160\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_13\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_13\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 161\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_23\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_23\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 162\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnozzle_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 163\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_number\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_number\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 164\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_density\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_density\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 165\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_outer_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_outer_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 166\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_inner_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_initial_inner_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 167\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_height\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_initial_height\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 168\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_separation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_separation\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 169\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrains_center_of_mass_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrains_center_of_mass_position\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 170\u001b[0m \u001b[43m \u001b[49m\u001b[43mcenter_of_dry_mass_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcenter_of_dry_mass_position\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 171\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnozzle_position\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 172\u001b[0m \u001b[43m \u001b[49m\u001b[43mthroat_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mthroat_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 173\u001b[0m \u001b[43m \u001b[49m\u001b[43minterpolate\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m 174\u001b[0m \u001b[43m \u001b[49m\u001b[43mcoordinate_system_orientation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m 175\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_motor_model.py:19\u001b[0m, in \u001b[0;36mStochasticMotorModel.__init__\u001b[0;34m(self, obj, **kwargs)\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[38;5;66;03m# TODO: never vary the grain_number\u001b[39;00m\n\u001b[1;32m 18\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_validate_positive_int_list(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrain_number\u001b[39m\u001b[38;5;124m\"\u001b[39m, kwargs\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrain_number\u001b[39m\u001b[38;5;124m\"\u001b[39m))\n\u001b[0;32m---> 19\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mobj\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_model.py:71\u001b[0m, in \u001b[0;36mStochasticModel.__init__\u001b[0;34m(self, obj, seed, **kwargs)\u001b[0m\n\u001b[1;32m 69\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlast_rnd_dict \u001b[38;5;241m=\u001b[39m {}\n\u001b[1;32m 70\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m__stochastic_dict \u001b[38;5;241m=\u001b[39m kwargs\n\u001b[0;32m---> 71\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_set_stochastic\u001b[49m\u001b[43m(\u001b[49m\u001b[43mseed\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_model.py:105\u001b[0m, in \u001b[0;36mStochasticModel._set_stochastic\u001b[0;34m(self, seed)\u001b[0m\n\u001b[1;32m 101\u001b[0m attr_value \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_validate_custom_sampler(\n\u001b[1;32m 102\u001b[0m input_value, seed\n\u001b[1;32m 103\u001b[0m )\n\u001b[1;32m 104\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 105\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAssertionError\u001b[39;00m(\n\u001b[1;32m 106\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00minput_name\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m must be a tuple, list, int, or float\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 107\u001b[0m )\n\u001b[1;32m 108\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 109\u001b[0m attr_value \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mgetattr\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mobj, input_name)]\n", + "\u001b[0;31mAssertionError\u001b[0m: 'total_impulse' must be a tuple, list, int, or float" + ] + } + ], + "source": [ + "stochastic_motor = StochasticSolidMotor(\n", + " solid_motor=motor,\n", + " burn_start_time=(0, 0.1, \"binomial\"),\n", + " grains_center_of_mass_position=0.001,\n", + " grain_density=50,\n", + " grain_separation=1 / 1000,\n", + " grain_initial_height=1 / 1000,\n", + " grain_initial_inner_radius=0.375 / 1000,\n", + " grain_outer_radius=0.375 / 1000,\n", + " total_impulse=total_impulse_sampler,\n", + " throat_radius=0.5 / 1000,\n", + " nozzle_radius=0.5 / 1000,\n", + " nozzle_position=0.001,\n", + ")\n", + "stochastic_motor.visualize_attributes()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "total_impulse_samples = [\n", + " stochastic_motor.create_object().total_impulse for _ in range(300)\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([3.08336103e-05, 3.08336103e-05, 3.39169713e-04, 6.16672205e-04,\n", + " 7.70840257e-04, 9.55841918e-04, 6.16672205e-04, 4.93337764e-04,\n", + " 5.24171375e-04, 1.54168051e-04, 1.85001662e-04, 1.54168051e-04,\n", + " 4.31670544e-04, 7.70840257e-04, 1.04834275e-03, 9.55841918e-04,\n", + " 5.85838595e-04, 2.77502492e-04, 1.85001662e-04, 1.23334441e-04]),\n", + " array([5385.36421608, 5493.47135319, 5601.5784903 , 5709.68562741,\n", + " 5817.79276452, 5925.89990164, 6034.00703875, 6142.11417586,\n", + " 6250.22131297, 6358.32845008, 6466.43558719, 6574.54272431,\n", + " 6682.64986142, 6790.75699853, 6898.86413564, 7006.97127275,\n", + " 7115.07840987, 7223.18554698, 7331.29268409, 7439.3998212 ,\n", + " 7547.50695831]),\n", + " )" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGdCAYAAAD5ZcJyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALbRJREFUeJzt3X90VPWd//FXQswP0ElEZJKwEaJGAUGCBIahCJ4ya9B0l7i0J6GpIE3h6AoLJyACi6Hr0QVTOeuhsGZxW6MVCrJ7DmqkcdMAspVpgAAqiBQ1NlA7AUwzA1HCj/l8//DL1ZHwY+gGSD7Pxzn3hNzP+/78cOe+zs29d2KMMUYAAACdXOyVXgEAAIDLgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALBC3JVegatJOBzWZ599puuuu04xMTFXenUAAMBFMMbo6NGjSk9PV2zsua/nEHq+4bPPPlNGRsaVXg0AAHAJDhw4oL/5m785Zzuh5xuuu+46SV/tNJfLdYXXBgAAXIxQKKSMjAznPH4uhJ5vOPMnLZfLRegBAKCDudCtKdzIDAAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGCFuCu9AgCAzqXP3Dfbbd6fLs5rt3mj8+NKDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACscEmhZ/ny5erTp48SExPl8Xi0devW89avXbtWffv2VWJiogYOHKj169dHtBtjVFpaqrS0NCUlJcnn82n//v0RNU8//bRGjBihrl27KiUlpc3lNDQ0KC8vT127dlXPnj312GOP6dSpU5eyiQAAoJOJOvSsWbNGJSUlWrhwoXbs2KFBgwYpNzdXhw4darN+y5YtmjBhgoqLi7Vz507l5+crPz9fu3fvdmrKysq0dOlSlZeXq7a2Vt26dVNubq6OHz/u1Jw4cUI/+MEP9Mgjj7S5nNOnTysvL08nTpzQli1b9NJLL6miokKlpaXRbiIAAOiEYowxJpoJPB6Phg4dqmXLlkmSwuGwMjIyNH36dM2dO/es+oKCArW0tKiystIZN3z4cGVnZ6u8vFzGGKWnp2vWrFmaPXu2JCkYDMrtdquiokKFhYUR86uoqNDMmTPV3NwcMf43v/mNvve97+mzzz6T2+2WJJWXl+vxxx/X4cOHFR8ff8FtC4VCSk5OVjAYlMvlima3AAD+P75lHZfbxZ6/o7rSc+LECdXV1cnn8309g9hY+Xw++f3+Nqfx+/0R9ZKUm5vr1NfX1ysQCETUJCcny+PxnHOe51rOwIEDncBzZjmhUEh79uy56PkAAIDOKS6a4iNHjuj06dMRwUKS3G63PvzwwzanCQQCbdYHAgGn/cy4c9VcjHMt55vL+LbW1la1trY6v4dCoYteHgAA6Fisfnpr0aJFSk5OdoaMjIwrvUoAAKCdRBV6evTooS5duqixsTFifGNjo1JTU9ucJjU19bz1Z35GM89olvPNZXzbvHnzFAwGneHAgQMXvTwAANCxRBV64uPjNWTIENXU1DjjwuGwampq5PV625zG6/VG1EtSdXW1U5+ZmanU1NSImlAopNra2nPO81zLef/99yOeIquurpbL5VL//v3bnCYhIUEulytiAAAAnVNU9/RIUklJiSZNmqScnBwNGzZMzz33nFpaWjR58mRJ0sSJE9WrVy8tWrRIkjRjxgyNHj1aS5YsUV5enlavXq3t27drxYoVkqSYmBjNnDlTTz31lLKyspSZmaknnnhC6enpys/Pd5bb0NCgpqYmNTQ06PTp09q1a5ck6dZbb9W1116re++9V/3799eDDz6osrIyBQIBLViwQI8++qgSEhL+yt0EAAA6uqhDT0FBgQ4fPqzS0lIFAgFlZ2erqqrKuWm4oaFBsbFfX0AaMWKEVq1apQULFmj+/PnKysrSunXrNGDAAKdmzpw5amlp0dSpU9Xc3KyRI0eqqqpKiYmJTk1paaleeukl5/fBgwdLkjZu3Kh77rlHXbp0UWVlpR555BF5vV5169ZNkyZN0pNPPhn9XgEAAJ1O1O/p6cx4Tw8A/PV4Tw8ut3Z5Tw8AAEBHRegBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFgh6jcyA7bihWvAlddexyHHoB240gMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALDCJYWe5cuXq0+fPkpMTJTH49HWrVvPW7927Vr17dtXiYmJGjhwoNavXx/RboxRaWmp0tLSlJSUJJ/Pp/3790fUNDU1qaioSC6XSykpKSouLtaxY8ciat566y0NHz5c1113nW688UaNHz9en3766aVsIgAA6GSiDj1r1qxRSUmJFi5cqB07dmjQoEHKzc3VoUOH2qzfsmWLJkyYoOLiYu3cuVP5+fnKz8/X7t27nZqysjItXbpU5eXlqq2tVbdu3ZSbm6vjx487NUVFRdqzZ4+qq6tVWVmpzZs3a+rUqU57fX29xo0bp+9+97vatWuX3nrrLR05ckT/8A//EO0mAgCATijGGGOimcDj8Wjo0KFatmyZJCkcDisjI0PTp0/X3Llzz6ovKChQS0uLKisrnXHDhw9Xdna2ysvLZYxRenq6Zs2apdmzZ0uSgsGg3G63KioqVFhYqL1796p///7atm2bcnJyJElVVVW6//77dfDgQaWnp+u//uu/NGHCBLW2tio29qss98Ybb2jcuHFqbW3VNddcc8FtC4VCSk5OVjAYlMvlima3wAJ95r7ZbvP+dHFeu80buNza81hpLxyDHdvFnr+jutJz4sQJ1dXVyefzfT2D2Fj5fD75/f42p/H7/RH1kpSbm+vU19fXKxAIRNQkJyfL4/E4NX6/XykpKU7gkSSfz6fY2FjV1tZKkoYMGaLY2Fi9+OKLOn36tILBoH71q1/J5/NdVOABAACdW1Sh58iRIzp9+rTcbnfEeLfbrUAg0OY0gUDgvPVnfl6opmfPnhHtcXFx6t69u1OTmZmp//mf/9H8+fOVkJCglJQUHTx4UK+++uo5t6e1tVWhUChiAAAAnVOneXorEAhoypQpmjRpkrZt26a3335b8fHx+v73v69z/QVv0aJFSk5OdoaMjIzLvNYAAOByiSr09OjRQ126dFFjY2PE+MbGRqWmprY5TWpq6nnrz/y8UM23b5Q+deqUmpqanJrly5crOTlZZWVlGjx4sEaNGqVXXnlFNTU1zp/Avm3evHkKBoPOcODAgYvZDQAAoAOKKvTEx8dryJAhqqmpccaFw2HV1NTI6/W2OY3X642ol6Tq6mqnPjMzU6mpqRE1oVBItbW1To3X61Vzc7Pq6uqcmg0bNigcDsvj8UiSvvjiC+cG5jO6dOnirGNbEhIS5HK5IgYAANA5Rf3nrZKSEr3wwgt66aWXtHfvXj3yyCNqaWnR5MmTJUkTJ07UvHnznPoZM2aoqqpKS5Ys0Ycffqif/vSn2r59u6ZNmyZJiomJ0cyZM/XUU0/p9ddf1/vvv6+JEycqPT1d+fn5kqR+/fpp7NixmjJlirZu3ap33nlH06ZNU2FhodLT0yVJeXl52rZtm5588knt379fO3bs0OTJk9W7d28NHjz4r91PAACgg4uLdoKCggIdPnxYpaWlCgQCys7OVlVVlXMjckNDQ8QVlxEjRmjVqlVasGCB5s+fr6ysLK1bt04DBgxwaubMmaOWlhZNnTpVzc3NGjlypKqqqpSYmOjUrFy5UtOmTdOYMWMUGxur8ePHa+nSpU77d7/7Xa1atUplZWUqKytT165d5fV6VVVVpaSkpEvaOQAAoPOI+j09nRnv6cH58J4e4OLwnh5cbu3ynh4AAICOitADAACsQOgBAABWiPpGZuBq1xHvJwCuBI4V2IYrPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsELclV4B2KnP3Dev9CpcVdprf3y6OK9d5gsAHRFXegAAgBUuKfQsX75cffr0UWJiojwej7Zu3Xre+rVr16pv375KTEzUwIEDtX79+oh2Y4xKS0uVlpampKQk+Xw+7d+/P6KmqalJRUVFcrlcSklJUXFxsY4dO3bWfJ599lnddtttSkhIUK9evfT0009fyiYCAIBOJurQs2bNGpWUlGjhwoXasWOHBg0apNzcXB06dKjN+i1btmjChAkqLi7Wzp07lZ+fr/z8fO3evdupKSsr09KlS1VeXq7a2lp169ZNubm5On78uFNTVFSkPXv2qLq6WpWVldq8ebOmTp0asawZM2boP//zP/Xss8/qww8/1Ouvv65hw4ZFu4kAAKATijHGmGgm8Hg8Gjp0qJYtWyZJCofDysjI0PTp0zV37tyz6gsKCtTS0qLKykpn3PDhw5Wdna3y8nIZY5Senq5Zs2Zp9uzZkqRgMCi3262KigoVFhZq79696t+/v7Zt26acnBxJUlVVle6//34dPHhQ6enp2rt3r+68807t3r1bt99++yXtjFAopOTkZAWDQblcrkuaBy4O9/RcHtzTg/PhOPwax0rHdrHn76iu9Jw4cUJ1dXXy+XxfzyA2Vj6fT36/v81p/H5/RL0k5ebmOvX19fUKBAIRNcnJyfJ4PE6N3+9XSkqKE3gkyefzKTY2VrW1tZKkN954QzfffLMqKyuVmZmpPn366Cc/+YmamprOuT2tra0KhUIRAwAA6JyiCj1HjhzR6dOn5Xa7I8a73W4FAoE2pwkEAuetP/PzQjU9e/aMaI+Li1P37t2dmk8++UR//OMftXbtWr388suqqKhQXV2dvv/9759zexYtWqTk5GRnyMjIuNAuAAAAHVSneXorHA6rtbVVL7/8su6++27dc889+sUvfqGNGzdq3759bU4zb948BYNBZzhw4MBlXmsAAHC5RBV6evTooS5duqixsTFifGNjo1JTU9ucJjU19bz1Z35eqObbN0qfOnVKTU1NTk1aWpri4uJ02223OTX9+vWTJDU0NLS5bgkJCXK5XBEDAADonKIKPfHx8RoyZIhqamqcceFwWDU1NfJ6vW1O4/V6I+olqbq62qnPzMxUampqRE0oFFJtba1T4/V61dzcrLq6Oqdmw4YNCofD8ng8kqTvfOc7OnXqlD7++GOn5g9/+IMkqXfv3tFsJgAA6ISifiNzSUmJJk2apJycHA0bNkzPPfecWlpaNHnyZEnSxIkT1atXLy1atEjSV4+Rjx49WkuWLFFeXp5Wr16t7du3a8WKFZKkmJgYzZw5U0899ZSysrKUmZmpJ554Qunp6crPz5f01RWbsWPHasqUKSovL9fJkyc1bdo0FRYWKj09XdJXNzbfdddd+vGPf6znnntO4XBYjz76qP72b/824uoPAACwU9Shp6CgQIcPH1ZpaakCgYCys7NVVVXl3Ijc0NCg2NivLyCNGDFCq1at0oIFCzR//nxlZWVp3bp1GjBggFMzZ84ctbS0aOrUqWpubtbIkSNVVVWlxMREp2blypWaNm2axowZo9jYWI0fP15Lly512mNjY/XGG29o+vTpGjVqlLp166b77rtPS5YsuaQdAwCwR3s+vs/j8FePqN/T05nxnp7Lh/eDXB582OJ8OA4vD47D9tcu7+kBAADoqAg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwApxV3oFAHRMfea+2S7z/XRxXrvMFwC40gMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAVuBrKIBOrL2+KgIAOiKu9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKcVd6BQAA59Zn7ptXehWAToMrPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAVrik0LN8+XL16dNHiYmJ8ng82rp163nr165dq759+yoxMVEDBw7U+vXrI9qNMSotLVVaWpqSkpLk8/m0f//+iJqmpiYVFRXJ5XIpJSVFxcXFOnbsWJvL++ijj3TdddcpJSXlUjYPAAB0QlGHnjVr1qikpEQLFy7Ujh07NGjQIOXm5urQoUNt1m/ZskUTJkxQcXGxdu7cqfz8fOXn52v37t1OTVlZmZYuXary8nLV1taqW7duys3N1fHjx52aoqIi7dmzR9XV1aqsrNTmzZs1derUs5Z38uRJTZgwQXfffXe0mwYAADqxGGOMiWYCj8ejoUOHatmyZZKkcDisjIwMTZ8+XXPnzj2rvqCgQC0tLaqsrHTGDR8+XNnZ2SovL5cxRunp6Zo1a5Zmz54tSQoGg3K73aqoqFBhYaH27t2r/v37a9u2bcrJyZEkVVVV6f7779fBgweVnp7uzPvxxx/XZ599pjFjxmjmzJlqbm6+6G0LhUJKTk5WMBiUy+WKZrcgSrxwDefy6eK8K70KVxWOlY6P/9Pt72LP31Fd6Tlx4oTq6urk8/m+nkFsrHw+n/x+f5vT+P3+iHpJys3Nderr6+sVCAQiapKTk+XxeJwav9+vlJQUJ/BIks/nU2xsrGpra51xGzZs0Nq1a7V8+fKL2p7W1laFQqGIAQAAdE5RhZ4jR47o9OnTcrvdEePdbrcCgUCb0wQCgfPWn/l5oZqePXtGtMfFxal79+5Ozeeff66HHnpIFRUVF32VZtGiRUpOTnaGjIyMi5oOAAB0PJ3m6a0pU6bohz/8oUaNGnXR08ybN0/BYNAZDhw40I5rCAAArqSoQk+PHj3UpUsXNTY2RoxvbGxUampqm9Okpqaet/7MzwvVfPtG6VOnTqmpqcmp2bBhg5599lnFxcUpLi5OxcXFCgaDiouL0y9/+cs21y0hIUEulytiAAAAnVNUoSc+Pl5DhgxRTU2NMy4cDqumpkZer7fNabxeb0S9JFVXVzv1mZmZSk1NjagJhUKqra11arxer5qbm1VXV+fUbNiwQeFwWB6PR9JX9/3s2rXLGZ588kldd9112rVrlx544IFoNhMAAHRCcdFOUFJSokmTJiknJ0fDhg3Tc889p5aWFk2ePFmSNHHiRPXq1UuLFi2SJM2YMUOjR4/WkiVLlJeXp9WrV2v79u1asWKFJCkmJkYzZ87UU089paysLGVmZuqJJ55Qenq68vPzJUn9+vXT2LFjNWXKFJWXl+vkyZOaNm2aCgsLnSe3+vXrF7Ge27dvV2xsrAYMGHDJOwcAAHQeUYeegoICHT58WKWlpQoEAsrOzlZVVZVzI3JDQ4NiY7++gDRixAitWrVKCxYs0Pz585WVlaV169ZFhJE5c+aopaVFU6dOVXNzs0aOHKmqqiolJiY6NStXrtS0adM0ZswYxcbGavz48Vq6dOlfs+0AAMAiUb+npzPjPT2XD+8ewbnwTpNIHCsdH/+n21+7vKcHAACgoyL0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGCFuCu9Ari68Q3PAIDOgis9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWCHuSq8AAACdWZ+5b7bLfD9dnNcu8+3MuNIDAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACpcUepYvX64+ffooMTFRHo9HW7duPW/92rVr1bdvXyUmJmrgwIFav359RLsxRqWlpUpLS1NSUpJ8Pp/2798fUdPU1KSioiK5XC6lpKSouLhYx44dc9o3bdqkcePGKS0tTd26dVN2drZWrlx5KZsHAAA6oahDz5o1a1RSUqKFCxdqx44dGjRokHJzc3Xo0KE267ds2aIJEyaouLhYO3fuVH5+vvLz87V7926npqysTEuXLlV5eblqa2vVrVs35ebm6vjx405NUVGR9uzZo+rqalVWVmrz5s2aOnVqxHLuvPNO/fd//7fee+89TZ48WRMnTlRlZWW0mwgAADqhGGOMiWYCj8ejoUOHatmyZZKkcDisjIwMTZ8+XXPnzj2rvqCgQC0tLRHhY/jw4crOzlZ5ebmMMUpPT9esWbM0e/ZsSVIwGJTb7VZFRYUKCwu1d+9e9e/fX9u2bVNOTo4kqaqqSvfff78OHjyo9PT0Ntc1Ly9Pbrdbv/zlLy9q20KhkJKTkxUMBuVyuaLZLZ1Wn7lvXulVgGU+XZx3pVfhqsIxiHPhWPnaxZ6/o7rSc+LECdXV1cnn8309g9hY+Xw++f3+Nqfx+/0R9ZKUm5vr1NfX1ysQCETUJCcny+PxODV+v18pKSlO4JEkn8+n2NhY1dbWnnN9g8Ggunfvfs721tZWhUKhiAEAAHROUYWeI0eO6PTp03K73RHj3W63AoFAm9MEAoHz1p/5eaGanj17RrTHxcWpe/fu51zuq6++qm3btmny5Mnn3J5FixYpOTnZGTIyMs5ZCwAAOrZO+fTWxo0bNXnyZL3wwgu64447zlk3b948BYNBZzhw4MBlXEsAAHA5RRV6evTooS5duqixsTFifGNjo1JTU9ucJjU19bz1Z35eqObbN0qfOnVKTU1NZy337bff1t/93d/p3/7t3zRx4sTzbk9CQoJcLlfEAAAAOqeoQk98fLyGDBmimpoaZ1w4HFZNTY28Xm+b03i93oh6SaqurnbqMzMzlZqaGlETCoVUW1vr1Hi9XjU3N6uurs6p2bBhg8LhsDwejzNu06ZNysvL0zPPPBPxZBcAAEBctBOUlJRo0qRJysnJ0bBhw/Tcc8+ppaXFuXdm4sSJ6tWrlxYtWiRJmjFjhkaPHq0lS5YoLy9Pq1ev1vbt27VixQpJUkxMjGbOnKmnnnpKWVlZyszM1BNPPKH09HTl5+dLkvr166exY8dqypQpKi8v18mTJzVt2jQVFhY6T25t3LhR3/ve9zRjxgyNHz/eudcnPj7+vDczAwAAO0QdegoKCnT48GGVlpYqEAgoOztbVVVVzo3IDQ0Nio39+gLSiBEjtGrVKi1YsEDz589XVlaW1q1bpwEDBjg1c+bMUUtLi6ZOnarm5maNHDlSVVVVSkxMdGpWrlypadOmacyYMYqNjdX48eO1dOlSp/2ll17SF198oUWLFjmBS5JGjx6tTZs2RbuZADqh9nz8m8eHcbnx/zl6Ub+npzPjPT1n4x0huNza88O2I54kOAZxJXS00NMu7+kBAADoqAg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVov4aCgBoTx31DcQddb0Bm3ClBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgTcyAwCACO31hvFPF+e1y3wvFld6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsEHelV8Amfea+eaVXAQAAa3GlBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFS4p9Cxfvlx9+vRRYmKiPB6Ptm7det76tWvXqm/fvkpMTNTAgQO1fv36iHZjjEpLS5WWlqakpCT5fD7t378/oqapqUlFRUVyuVxKSUlRcXGxjh07FlHz3nvv6e6771ZiYqIyMjJUVlZ2KZsHAAA6oahDz5o1a1RSUqKFCxdqx44dGjRokHJzc3Xo0KE267ds2aIJEyaouLhYO3fuVH5+vvLz87V7926npqysTEuXLlV5eblqa2vVrVs35ebm6vjx405NUVGR9uzZo+rqalVWVmrz5s2aOnWq0x4KhXTvvfeqd+/eqqur089+9jP99Kc/1YoVK6LdRAAA0AnFGGNMNBN4PB4NHTpUy5YtkySFw2FlZGRo+vTpmjt37ln1BQUFamlpUWVlpTNu+PDhys7OVnl5uYwxSk9P16xZszR79mxJUjAYlNvtVkVFhQoLC7V37171799f27ZtU05OjiSpqqpK999/vw4ePKj09HQ9//zz+ud//mcFAgHFx8dLkubOnat169bpww8/vKhtC4VCSk5OVjAYlMvlima3XBTe0wMAsNmni/PaZb4Xe/6O6uWEJ06cUF1dnebNm+eMi42Nlc/nk9/vb3Mav9+vkpKSiHG5ublat26dJKm+vl6BQEA+n89pT05Olsfjkd/vV2Fhofx+v1JSUpzAI0k+n0+xsbGqra3VAw88IL/fr1GjRjmB58xynnnmGf3lL3/R9ddff9a6tba2qrW11fk9GAxK+mrntYdw6xftMl8AADqC9jq/npnvha7jRBV6jhw5otOnT8vtdkeMd7vd57yaEggE2qwPBAJO+5lx56vp2bNn5IrHxal79+4RNZmZmWfN40xbW6Fn0aJF+pd/+ZezxmdkZLS5LQAA4NIlP9e+8z969KiSk5PP2W7111DMmzcv4ipUOBxWU1OTbrjhBsXExFzBNYtOKBRSRkaGDhw40C5/lsNfjz66utE/Vzf65+p3pfvIGKOjR48qPT39vHVRhZ4ePXqoS5cuamxsjBjf2Nio1NTUNqdJTU09b/2Zn42NjUpLS4uoyc7Odmq+faP0qVOn1NTUFDGftpbzzWV8W0JCghISEiLGpaSktFnbEbhcLj4QrnL00dWN/rm60T9XvyvZR+e7wnNGVE9vxcfHa8iQIaqpqXHGhcNh1dTUyOv1tjmN1+uNqJek6upqpz4zM1OpqakRNaFQSLW1tU6N1+tVc3Oz6urqnJoNGzYoHA7L4/E4NZs3b9bJkycjlnP77be3+actAABgGROl1atXm4SEBFNRUWE++OADM3XqVJOSkmICgYAxxpgHH3zQzJ0716l/5513TFxcnHn22WfN3r17zcKFC80111xj3n//fadm8eLFJiUlxbz22mvmvffeM+PGjTOZmZnmyy+/dGrGjh1rBg8ebGpra83vfvc7k5WVZSZMmOC0Nzc3G7fbbR588EGze/dus3r1atO1a1fzH//xH9FuYocTDAaNJBMMBq/0quAc6KOrG/1zdaN/rn4dpY+iDj3GGPPzn//c3HTTTSY+Pt4MGzbM/P73v3faRo8ebSZNmhRR/+qrr5rbbrvNxMfHmzvuuMO8+eabEe3hcNg88cQTxu12m4SEBDNmzBizb9++iJrPP//cTJgwwVx77bXG5XKZyZMnm6NHj0bUvPvuu2bkyJEmISHB9OrVyyxevPhSNq/DOX78uFm4cKE5fvz4lV4VnAN9dHWjf65u9M/Vr6P0UdTv6QEAAOiI+O4tAABgBUIPAACwAqEHAABYgdADAACsQOi5Svz0pz9VTExMxNC3b1+n/Z577jmr/eGHH46YR0NDg/Ly8tS1a1f17NlTjz32mE6dOhVRs2nTJt11111KSEjQrbfeqoqKisuxeZ3Cn/70J/3oRz/SDTfcoKSkJA0cOFDbt2932o0xKi0tVVpampKSkuTz+bR///6IeTQ1NamoqEgul0spKSkqLi7WsWPHImree+893X333UpMTFRGRobKysouy/Z1dBfqn4ceeuisY2js2LER86B/2k+fPn3O2v8xMTF69NFHJUnHjx/Xo48+qhtuuEHXXnutxo8ff9YLZ/mMaz8X6p9Ocw66wk+P4f9buHChueOOO8yf//xnZzh8+LDTPnr0aDNlypSI9m++D+HUqVNmwIABxufzmZ07d5r169ebHj16mHnz5jk1n3zyienataspKSkxH3zwgfn5z39uunTpYqqqqi7rtnZETU1Npnfv3uahhx4ytbW15pNPPjFvvfWW+eijj5yaxYsXm+TkZLNu3Trz7rvvmr//+79v831TgwYNMr///e/N//7v/5pbb7014n1TwWDQuN1uU1RUZHbv3m1+/etfm6SkJCveN/XXuJj+mTRpkhk7dmzEMdTU1BQxH/qn/Rw6dChi31dXVxtJZuPGjcYYYx5++GGTkZFhampqzPbt283w4cPNiBEjnOn5jGtfF+qfznIOIvRcJRYuXGgGDRp0zvbRo0ebGTNmnLN9/fr1JjY21nlJpDHGPP/888blcpnW1lZjjDFz5swxd9xxR8R0BQUFJjc3969adxs8/vjjZuTIkedsD4fDJjU11fzsZz9zxjU3N5uEhATz61//2hhjzAcffGAkmW3btjk1v/nNb0xMTIz505/+ZIwx5t///d/N9ddf7/TZmWXffvvt/9eb1KlcqH+M+Sr0jBs37pzt9M/lNWPGDHPLLbeYcDhsmpubzTXXXGPWrl3rtO/du9dIMn6/3xjDZ9zl9s3+MabznIP489ZVZP/+/UpPT9fNN9+soqIiNTQ0RLSvXLlSPXr00IABAzRv3jx98cUXTpvf79fAgQMjvq0+NzdXoVBIe/bscWp8Pl/EPHNzc+X3+9txqzqH119/XTk5OfrBD36gnj17avDgwXrhhRec9vr6egUCgYj9m5ycLI/H4+xfv9+vlJQU5eTkODU+n0+xsbGqra11akaNGqX4+HinJjc3V/v27dNf/vKX9t7MDutC/XPGpk2b1LNnT91+++165JFH9Pnnnztt9M/lc+LECb3yyiv68Y9/rJiYGNXV1enkyZMRx0/fvn110003RRw/fMZdHt/unzM6wzmI0HOV8Hg8qqioUFVVlZ5//nnV19fr7rvv1tGjRyVJP/zhD/XKK69o48aNmjdvnn71q1/pRz/6kTN9IBCI+M8myfk9EAictyYUCunLL79sz83r8D755BM9//zzysrK0ltvvaVHHnlE//RP/6SXXnpJ0tf7uK39+83937Nnz4j2uLg4de/e/YJ99M1l4GwX6h9JGjt2rF5++WXV1NTomWee0dtvv6377rtPp0+flkT/XE7r1q1Tc3OzHnroIUlf7bv4+PizvvD528cPn3GXx7f7R+o856CovmUd7ee+++5z/n3nnXfK4/God+/eevXVV1VcXKypU6c67QMHDlRaWprGjBmjjz/+WLfccsuVWGWrhMNh5eTk6F//9V8lSYMHD9bu3btVXl6uSZMmXeG1w8X0T2FhoVM/cOBA3Xnnnbrlllu0adMmjRkz5oqst61+8Ytf6L777lN6evqVXhW0oa3+6SznIK70XKVSUlJ022236aOPPmqz/cy3y59pT01NPetJhzO/p6amnrfG5XIpKSnp/3T9O5u0tDT1798/Yly/fv2cP0Ge2cdt7d9v7v9Dhw5FtJ86dUpNTU0X7KNvLgNnu1D/tOXmm29Wjx49Io4h+qf9/fGPf9Rvf/tb/eQnP3HGpaam6sSJE2pubo6o/fbxw2dc+2urf9rSUc9BhJ6r1LFjx/Txxx8rLS2tzfZdu3ZJktPu9Xr1/vvvR3xoV1dXy+VyOScDr9ermpqaiPlUV1fL6/W2wxZ0Lt/5zne0b9++iHF/+MMf1Lt3b0lSZmamUlNTI/ZvKBRSbW2ts3+9Xq+am5tVV1fn1GzYsEHhcNj5APF6vdq8ebNOnjzp1FRXV+v222/X9ddf327b19FdqH/acvDgQX3++ecRxxD90/5efPFF9ezZU3l5ec64IUOG6Jprrok4fvbt26eGhoaI44fPuPbXVv+0pcOegy7bLdM4r1mzZplNmzaZ+vp688477xifz2d69OhhDh06ZD766CPz5JNPmu3bt5v6+nrz2muvmZtvvtmMGjXKmf7M44L33nuv2bVrl6mqqjI33nhjm48LPvbYY2bv3r1m+fLlPM55kbZu3Wri4uLM008/bfbv329Wrlxpunbtal555RWnZvHixSYlJcW89tpr5r333jPjxo1r85H1wYMHm9raWvO73/3OZGVlRTwS3dzcbNxut3nwwQfN7t27zerVq03Xrl15JPoCLtQ/R48eNbNnzzZ+v9/U19eb3/72t+auu+4yWVlZEd8KTf+0r9OnT5ubbrrJPP7442e1Pfzww+amm24yGzZsMNu3bzder9d4vV6nnc+49neu/ulM5yBCz1WioKDApKWlmfj4eNOrVy9TUFDgvGOkoaHBjBo1ynTv3t0kJCSYW2+91Tz22GMR70gwxphPP/3U3HfffSYpKcn06NHDzJo1y5w8eTKiZuPGjSY7O9vEx8ebm2++2bz44ouXaxM7vDfeeMMMGDDAJCQkmL59+5oVK1ZEtIfDYfPEE08Yt9ttEhISzJgxY8y+ffsiaj7//HMzYcIEc+211xqXy2UmT55sjh49GlHz7rvvmpEjR5qEhATTq1cvs3jx4nbfts7gfP3zxRdfmHvvvdfceOON5pprrjG9e/c2U6ZMiXi81hj6p7299dZbRtJZx4Uxxnz55ZfmH//xH831119vunbtah544AHz5z//OaKGz7j2da7+6UznoBhjjLl815UAAACuDO7pAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAK/w9o7PHnvsHhBAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(total_impulse_samples, density = True, bins = 20)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "dev-py3-12", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From f3767055df73c1aa830f9700322f24e607bcc191 Mon Sep 17 00:00:00 2001 From: Lucas de Oliveira Prates Date: Thu, 17 Apr 2025 22:22:17 +0200 Subject: [PATCH 2/4] TST: implementing tests for CustomSampler and making small corrections --- .vscode/settings.json | 1 + rocketpy/stochastic/custom_sampler.py | 4 +- rocketpy/stochastic/stochastic_model.py | 17 +- test_custom_sampler.ipynb | 232 ------------------ tests/conftest.py | 1 + .../monte_carlo/custom_sampler_fixtures.py | 76 ++++++ .../monte_carlo/stochastic_fixtures.py | 30 +++ tests/unit/stochastic/test_custom_sampler.py | 21 ++ tests/unit/test_stochastic_model.py | 1 + 9 files changed, 144 insertions(+), 239 deletions(-) delete mode 100644 test_custom_sampler.ipynb create mode 100644 tests/fixtures/monte_carlo/custom_sampler_fixtures.py create mode 100644 tests/unit/stochastic/test_custom_sampler.py diff --git a/.vscode/settings.json b/.vscode/settings.json index 5cd5a878d..7fa0c1ff8 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -48,6 +48,7 @@ "Bigl", "Bigr", "bijective", + "Bivariate", "bmatrix", "boldsymbol", "boxplot", diff --git a/rocketpy/stochastic/custom_sampler.py b/rocketpy/stochastic/custom_sampler.py index c13887060..82a06dd9f 100644 --- a/rocketpy/stochastic/custom_sampler.py +++ b/rocketpy/stochastic/custom_sampler.py @@ -10,7 +10,7 @@ class CustomSampler(ABC): @abstractmethod def sample(self, n_samples=1): - """Generates n samples from the custom distribution + """Generates samples from the custom distribution Parameters ---------- @@ -19,7 +19,7 @@ def sample(self, n_samples=1): Returns ------- - sample_list : list + samples_list : list A list with n_samples elements, each of which is a valid sample """ diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index 103ec7175..879b61e70 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -89,7 +89,9 @@ def _set_stochastic(self, seed=None): attr_value = None if input_value is not None: if "factor" in input_name: - attr_value = self._validate_factors(input_name, input_value) + attr_value = self._validate_factors( + input_name, input_value, seed + ) elif input_name not in self.exception_list: if isinstance(input_value, tuple): attr_value = self._validate_tuple(input_name, input_value) @@ -104,6 +106,7 @@ def _set_stochastic(self, seed=None): else: raise AssertionError( f"'{input_name}' must be a tuple, list, int, or float" + "or a custom sampler" ) else: attr_value = [getattr(self.obj, input_name)] @@ -285,7 +288,7 @@ def _validate_scalar(self, input_name, input_value, getattr=getattr): # pylint: get_distribution("normal", self.__random_number_generator), ) - def _validate_factors(self, input_name, input_value): + def _validate_factors(self, input_name, input_value, seed): """ Validate factor arguments. @@ -313,8 +316,12 @@ def _validate_factors(self, input_name, input_value): return self._validate_tuple_factor(input_name, input_value) elif isinstance(input_value, list): return self._validate_list_factor(input_name, input_value) + elif isinstance(input_value, CustomSampler): + return self._validate_custom_sampler(input_name, input_value, seed) else: - raise AssertionError(f"`{input_name}`: must be either a tuple or list") + raise AssertionError( + f"`{input_name}`: must be either a tuple or listor a custom sampler" + ) def _validate_tuple_factor(self, input_name, factor_tuple): """ @@ -463,7 +470,7 @@ def _validate_custom_sampler(self, input_name, sampler, seed=None): sampler.reset_seed(seed) except RuntimeError as e: raise RuntimeError( - f"An error occurred in the 'reset_seed' of {input_name} CustomSampler" + f"An error occurred in the 'reset_seed' method of {input_name} CustomSampler" ) from e return sampler @@ -531,7 +538,7 @@ def dict_generator(self): generated_dict[arg] = value.sample(n_samples=1)[0] except RuntimeError as e: raise RuntimeError( - f"An error occurred in the 'sample' of {arg} CustomSampler" + f"An error occurred in the 'sample' method of {arg} CustomSampler" ) from e self.last_rnd_dict = generated_dict yield generated_dict diff --git a/test_custom_sampler.ipynb b/test_custom_sampler.ipynb deleted file mode 100644 index c831f0cf9..000000000 --- a/test_custom_sampler.ipynb +++ /dev/null @@ -1,232 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 20, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The autoreload extension is already loaded. To reload it, use:\n", - " %reload_ext autoreload\n" - ] - } - ], - "source": [ - "# We import these lines for debugging purposes, only works on Jupyter Notebook\n", - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [], - "source": [ - "from matplotlib import pyplot as plt\n", - "import numpy as np\n", - "from rocketpy import SolidMotor, StochasticSolidMotor, CustomSampler" - ] - }, - { - "cell_type": "code", - "execution_count": 22, - "metadata": {}, - "outputs": [], - "source": [ - "class GaussianMixture(CustomSampler):\n", - "\n", - " def __init__(self, means_tuple, sd_tuple, prob_tuple, seed = None):\n", - " \"\"\" Creates a sampler for a mixture of two Gaussians\n", - "\n", - " Parameters\n", - " ----------\n", - " means_tuple : 2-tuple\n", - " 2-Tuple that contains the mean of each normal distribution of the\n", - " mixture\n", - " sd_tuple : 2-tuple\n", - " 2-Tuple that contains the sd of each normal distribution of the\n", - " mixture\n", - " prob_tuple : 2-tuple\n", - " 2-Tuple that contains the probability of each normal distribution of the\n", - " mixture\n", - " \"\"\"\n", - " np.random.default_rng(seed)\n", - " self.means_tuple = means_tuple\n", - " self.sd_tuple = sd_tuple\n", - " self.prob_tuple = prob_tuple\n", - " \n", - " def sample(self, n_samples = 1):\n", - " sample_list = [0] * n_samples\n", - " mixture_id_list = np.random.binomial(1, self.prob_tuple[0], n_samples)\n", - " for i, mixture_id in enumerate(mixture_id_list):\n", - " if mixture_id:\n", - " sample_list[i] = np.random.normal(self.means_tuple[0], self.sd_tuple[0])\n", - " else:\n", - " sample_list[i] = np.random.normal(self.means_tuple[1], self.sd_tuple[1])\n", - "\n", - " return sample_list\n", - " \n", - " def reset_seed(self, seed=None):\n", - " np.random.default_rng(seed)" - ] - }, - { - "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [], - "source": [ - "total_impulse_sampler = GaussianMixture([6000, 7000], [200, 200], [0.5, 0.5])" - ] - }, - { - "cell_type": "code", - "execution_count": 24, - "metadata": {}, - "outputs": [], - "source": [ - "motor = SolidMotor(\n", - " thrust_source=\"data/motors/cesaroni/Cesaroni_M1670.eng\",\n", - " dry_mass=1.815,\n", - " dry_inertia=(0.125, 0.125, 0.002),\n", - " nozzle_radius=33 / 1000,\n", - " grain_number=5,\n", - " grain_density=1815,\n", - " grain_outer_radius=33 / 1000,\n", - " grain_initial_inner_radius=15 / 1000,\n", - " grain_initial_height=120 / 1000,\n", - " grain_separation=5 / 1000,\n", - " grains_center_of_mass_position=0.397,\n", - " center_of_dry_mass_position=0.317,\n", - " nozzle_position=0,\n", - " burn_time=3.9,\n", - " throat_radius=11 / 1000,\n", - " coordinate_system_orientation=\"nozzle_to_combustion_chamber\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [ - { - "ename": "AssertionError", - "evalue": "'total_impulse' must be a tuple, list, int, or float", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[26], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m stochastic_motor \u001b[38;5;241m=\u001b[39m \u001b[43mStochasticSolidMotor\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 2\u001b[0m \u001b[43m \u001b[49m\u001b[43msolid_motor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmotor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 3\u001b[0m \u001b[43m \u001b[49m\u001b[43mburn_start_time\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m0.1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mbinomial\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 4\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrains_center_of_mass_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.001\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_density\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m50\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 6\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_separation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 7\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_height\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 8\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_inner_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.375\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 9\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_outer_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.375\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 10\u001b[0m \u001b[43m \u001b[49m\u001b[43mtotal_impulse\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtotal_impulse_sampler\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 11\u001b[0m \u001b[43m \u001b[49m\u001b[43mthroat_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.5\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 12\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.5\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 13\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.001\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 14\u001b[0m \u001b[43m)\u001b[49m\n\u001b[1;32m 15\u001b[0m \u001b[38;5;66;03m#stochastic_motor.visualize_attributes()\u001b[39;00m\n", - "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_solid_motor.py:149\u001b[0m, in \u001b[0;36mStochasticSolidMotor.__init__\u001b[0;34m(self, solid_motor, thrust_source, total_impulse, burn_start_time, burn_out_time, dry_mass, dry_inertia_11, dry_inertia_22, dry_inertia_33, dry_inertia_12, dry_inertia_13, dry_inertia_23, nozzle_radius, grain_number, grain_density, grain_outer_radius, grain_initial_inner_radius, grain_initial_height, grain_separation, grains_center_of_mass_position, center_of_dry_mass_position, nozzle_position, throat_radius)\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\n\u001b[1;32m 69\u001b[0m \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m 70\u001b[0m solid_motor,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 92\u001b[0m throat_radius\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 93\u001b[0m ):\n\u001b[1;32m 94\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Initializes the Stochastic Solid Motor class.\u001b[39;00m\n\u001b[1;32m 95\u001b[0m \n\u001b[1;32m 96\u001b[0m \u001b[38;5;124;03m See Also\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 147\u001b[0m \u001b[38;5;124;03m Radius of the throat in the motor in meters.\u001b[39;00m\n\u001b[1;32m 148\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 149\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[1;32m 150\u001b[0m \u001b[43m \u001b[49m\u001b[43msolid_motor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 151\u001b[0m \u001b[43m \u001b[49m\u001b[43mthrust_source\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mthrust_source\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 152\u001b[0m \u001b[43m \u001b[49m\u001b[43mtotal_impulse\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtotal_impulse\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 153\u001b[0m \u001b[43m \u001b[49m\u001b[43mburn_start_time\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mburn_start_time\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 154\u001b[0m \u001b[43m \u001b[49m\u001b[43mburn_out_time\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mburn_out_time\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 155\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_mass\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_mass\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 156\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_11\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_11\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 157\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_22\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_22\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 158\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_33\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_33\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 159\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_12\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_12\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 160\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_13\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_13\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 161\u001b[0m \u001b[43m \u001b[49m\u001b[43mdry_I_23\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdry_inertia_23\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 162\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnozzle_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 163\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_number\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_number\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 164\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_density\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_density\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 165\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_outer_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_outer_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 166\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_inner_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_initial_inner_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 167\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_initial_height\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_initial_height\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 168\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrain_separation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrain_separation\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 169\u001b[0m \u001b[43m \u001b[49m\u001b[43mgrains_center_of_mass_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mgrains_center_of_mass_position\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 170\u001b[0m \u001b[43m \u001b[49m\u001b[43mcenter_of_dry_mass_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcenter_of_dry_mass_position\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 171\u001b[0m \u001b[43m \u001b[49m\u001b[43mnozzle_position\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnozzle_position\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 172\u001b[0m \u001b[43m \u001b[49m\u001b[43mthroat_radius\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mthroat_radius\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 173\u001b[0m \u001b[43m \u001b[49m\u001b[43minterpolate\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m 174\u001b[0m \u001b[43m \u001b[49m\u001b[43mcoordinate_system_orientation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m 175\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_motor_model.py:19\u001b[0m, in \u001b[0;36mStochasticMotorModel.__init__\u001b[0;34m(self, obj, **kwargs)\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[38;5;66;03m# TODO: never vary the grain_number\u001b[39;00m\n\u001b[1;32m 18\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_validate_positive_int_list(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrain_number\u001b[39m\u001b[38;5;124m\"\u001b[39m, kwargs\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrain_number\u001b[39m\u001b[38;5;124m\"\u001b[39m))\n\u001b[0;32m---> 19\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mobj\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_model.py:71\u001b[0m, in \u001b[0;36mStochasticModel.__init__\u001b[0;34m(self, obj, seed, **kwargs)\u001b[0m\n\u001b[1;32m 69\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlast_rnd_dict \u001b[38;5;241m=\u001b[39m {}\n\u001b[1;32m 70\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m__stochastic_dict \u001b[38;5;241m=\u001b[39m kwargs\n\u001b[0;32m---> 71\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_set_stochastic\u001b[49m\u001b[43m(\u001b[49m\u001b[43mseed\u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/Desktop/Work/RocketPy/RocketPy/rocketpy/stochastic/stochastic_model.py:105\u001b[0m, in \u001b[0;36mStochasticModel._set_stochastic\u001b[0;34m(self, seed)\u001b[0m\n\u001b[1;32m 101\u001b[0m attr_value \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_validate_custom_sampler(\n\u001b[1;32m 102\u001b[0m input_value, seed\n\u001b[1;32m 103\u001b[0m )\n\u001b[1;32m 104\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 105\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAssertionError\u001b[39;00m(\n\u001b[1;32m 106\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00minput_name\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m must be a tuple, list, int, or float\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 107\u001b[0m )\n\u001b[1;32m 108\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 109\u001b[0m attr_value \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mgetattr\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mobj, input_name)]\n", - "\u001b[0;31mAssertionError\u001b[0m: 'total_impulse' must be a tuple, list, int, or float" - ] - } - ], - "source": [ - "stochastic_motor = StochasticSolidMotor(\n", - " solid_motor=motor,\n", - " burn_start_time=(0, 0.1, \"binomial\"),\n", - " grains_center_of_mass_position=0.001,\n", - " grain_density=50,\n", - " grain_separation=1 / 1000,\n", - " grain_initial_height=1 / 1000,\n", - " grain_initial_inner_radius=0.375 / 1000,\n", - " grain_outer_radius=0.375 / 1000,\n", - " total_impulse=total_impulse_sampler,\n", - " throat_radius=0.5 / 1000,\n", - " nozzle_radius=0.5 / 1000,\n", - " nozzle_position=0.001,\n", - ")\n", - "stochastic_motor.visualize_attributes()" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "total_impulse_samples = [\n", - " stochastic_motor.create_object().total_impulse for _ in range(300)\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(array([3.08336103e-05, 3.08336103e-05, 3.39169713e-04, 6.16672205e-04,\n", - " 7.70840257e-04, 9.55841918e-04, 6.16672205e-04, 4.93337764e-04,\n", - " 5.24171375e-04, 1.54168051e-04, 1.85001662e-04, 1.54168051e-04,\n", - " 4.31670544e-04, 7.70840257e-04, 1.04834275e-03, 9.55841918e-04,\n", - " 5.85838595e-04, 2.77502492e-04, 1.85001662e-04, 1.23334441e-04]),\n", - " array([5385.36421608, 5493.47135319, 5601.5784903 , 5709.68562741,\n", - " 5817.79276452, 5925.89990164, 6034.00703875, 6142.11417586,\n", - " 6250.22131297, 6358.32845008, 6466.43558719, 6574.54272431,\n", - " 6682.64986142, 6790.75699853, 6898.86413564, 7006.97127275,\n", - " 7115.07840987, 7223.18554698, 7331.29268409, 7439.3998212 ,\n", - " 7547.50695831]),\n", - " )" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGdCAYAAAD5ZcJyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALbRJREFUeJzt3X90VPWd//FXQswP0ElEZJKwEaJGAUGCBIahCJ4ya9B0l7i0J6GpIE3h6AoLJyACi6Hr0QVTOeuhsGZxW6MVCrJ7DmqkcdMAspVpgAAqiBQ1NlA7AUwzA1HCj/l8//DL1ZHwY+gGSD7Pxzn3hNzP+/78cOe+zs29d2KMMUYAAACdXOyVXgEAAIDLgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALBC3JVegatJOBzWZ599puuuu04xMTFXenUAAMBFMMbo6NGjSk9PV2zsua/nEHq+4bPPPlNGRsaVXg0AAHAJDhw4oL/5m785Zzuh5xuuu+46SV/tNJfLdYXXBgAAXIxQKKSMjAznPH4uhJ5vOPMnLZfLRegBAKCDudCtKdzIDAAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGCFuCu9AgCAzqXP3Dfbbd6fLs5rt3mj8+NKDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACscEmhZ/ny5erTp48SExPl8Xi0devW89avXbtWffv2VWJiogYOHKj169dHtBtjVFpaqrS0NCUlJcnn82n//v0RNU8//bRGjBihrl27KiUlpc3lNDQ0KC8vT127dlXPnj312GOP6dSpU5eyiQAAoJOJOvSsWbNGJSUlWrhwoXbs2KFBgwYpNzdXhw4darN+y5YtmjBhgoqLi7Vz507l5+crPz9fu3fvdmrKysq0dOlSlZeXq7a2Vt26dVNubq6OHz/u1Jw4cUI/+MEP9Mgjj7S5nNOnTysvL08nTpzQli1b9NJLL6miokKlpaXRbiIAAOiEYowxJpoJPB6Phg4dqmXLlkmSwuGwMjIyNH36dM2dO/es+oKCArW0tKiystIZN3z4cGVnZ6u8vFzGGKWnp2vWrFmaPXu2JCkYDMrtdquiokKFhYUR86uoqNDMmTPV3NwcMf43v/mNvve97+mzzz6T2+2WJJWXl+vxxx/X4cOHFR8ff8FtC4VCSk5OVjAYlMvlima3AAD+P75lHZfbxZ6/o7rSc+LECdXV1cnn8309g9hY+Xw++f3+Nqfx+/0R9ZKUm5vr1NfX1ysQCETUJCcny+PxnHOe51rOwIEDncBzZjmhUEh79uy56PkAAIDOKS6a4iNHjuj06dMRwUKS3G63PvzwwzanCQQCbdYHAgGn/cy4c9VcjHMt55vL+LbW1la1trY6v4dCoYteHgAA6Fisfnpr0aJFSk5OdoaMjIwrvUoAAKCdRBV6evTooS5duqixsTFifGNjo1JTU9ucJjU19bz1Z35GM89olvPNZXzbvHnzFAwGneHAgQMXvTwAANCxRBV64uPjNWTIENXU1DjjwuGwampq5PV625zG6/VG1EtSdXW1U5+ZmanU1NSImlAopNra2nPO81zLef/99yOeIquurpbL5VL//v3bnCYhIUEulytiAAAAnVNU9/RIUklJiSZNmqScnBwNGzZMzz33nFpaWjR58mRJ0sSJE9WrVy8tWrRIkjRjxgyNHj1aS5YsUV5enlavXq3t27drxYoVkqSYmBjNnDlTTz31lLKyspSZmaknnnhC6enpys/Pd5bb0NCgpqYmNTQ06PTp09q1a5ck6dZbb9W1116re++9V/3799eDDz6osrIyBQIBLViwQI8++qgSEhL+yt0EAAA6uqhDT0FBgQ4fPqzS0lIFAgFlZ2erqqrKuWm4oaFBsbFfX0AaMWKEVq1apQULFmj+/PnKysrSunXrNGDAAKdmzpw5amlp0dSpU9Xc3KyRI0eqqqpKiYmJTk1paaleeukl5/fBgwdLkjZu3Kh77rlHXbp0UWVlpR555BF5vV5169ZNkyZN0pNPPhn9XgEAAJ1O1O/p6cx4Tw8A/PV4Tw8ut3Z5Tw8AAEBHRegBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFgh6jcyA7bihWvAlddexyHHoB240gMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALDCJYWe5cuXq0+fPkpMTJTH49HWrVvPW7927Vr17dtXiYmJGjhwoNavXx/RboxRaWmp0tLSlJSUJJ/Pp/3790fUNDU1qaioSC6XSykpKSouLtaxY8ciat566y0NHz5c1113nW688UaNHz9en3766aVsIgAA6GSiDj1r1qxRSUmJFi5cqB07dmjQoEHKzc3VoUOH2qzfsmWLJkyYoOLiYu3cuVP5+fnKz8/X7t27nZqysjItXbpU5eXlqq2tVbdu3ZSbm6vjx487NUVFRdqzZ4+qq6tVWVmpzZs3a+rUqU57fX29xo0bp+9+97vatWuX3nrrLR05ckT/8A//EO0mAgCATijGGGOimcDj8Wjo0KFatmyZJCkcDisjI0PTp0/X3Llzz6ovKChQS0uLKisrnXHDhw9Xdna2ysvLZYxRenq6Zs2apdmzZ0uSgsGg3G63KioqVFhYqL1796p///7atm2bcnJyJElVVVW6//77dfDgQaWnp+u//uu/NGHCBLW2tio29qss98Ybb2jcuHFqbW3VNddcc8FtC4VCSk5OVjAYlMvlima3wAJ95r7ZbvP+dHFeu80buNza81hpLxyDHdvFnr+jutJz4sQJ1dXVyefzfT2D2Fj5fD75/f42p/H7/RH1kpSbm+vU19fXKxAIRNQkJyfL4/E4NX6/XykpKU7gkSSfz6fY2FjV1tZKkoYMGaLY2Fi9+OKLOn36tILBoH71q1/J5/NdVOABAACdW1Sh58iRIzp9+rTcbnfEeLfbrUAg0OY0gUDgvPVnfl6opmfPnhHtcXFx6t69u1OTmZmp//mf/9H8+fOVkJCglJQUHTx4UK+++uo5t6e1tVWhUChiAAAAnVOneXorEAhoypQpmjRpkrZt26a3335b8fHx+v73v69z/QVv0aJFSk5OdoaMjIzLvNYAAOByiSr09OjRQ126dFFjY2PE+MbGRqWmprY5TWpq6nnrz/y8UM23b5Q+deqUmpqanJrly5crOTlZZWVlGjx4sEaNGqVXXnlFNTU1zp/Avm3evHkKBoPOcODAgYvZDQAAoAOKKvTEx8dryJAhqqmpccaFw2HV1NTI6/W2OY3X642ol6Tq6mqnPjMzU6mpqRE1oVBItbW1To3X61Vzc7Pq6uqcmg0bNigcDsvj8UiSvvjiC+cG5jO6dOnirGNbEhIS5HK5IgYAANA5Rf3nrZKSEr3wwgt66aWXtHfvXj3yyCNqaWnR5MmTJUkTJ07UvHnznPoZM2aoqqpKS5Ys0Ycffqif/vSn2r59u6ZNmyZJiomJ0cyZM/XUU0/p9ddf1/vvv6+JEycqPT1d+fn5kqR+/fpp7NixmjJlirZu3ap33nlH06ZNU2FhodLT0yVJeXl52rZtm5588knt379fO3bs0OTJk9W7d28NHjz4r91PAACgg4uLdoKCggIdPnxYpaWlCgQCys7OVlVVlXMjckNDQ8QVlxEjRmjVqlVasGCB5s+fr6ysLK1bt04DBgxwaubMmaOWlhZNnTpVzc3NGjlypKqqqpSYmOjUrFy5UtOmTdOYMWMUGxur8ePHa+nSpU77d7/7Xa1atUplZWUqKytT165d5fV6VVVVpaSkpEvaOQAAoPOI+j09nRnv6cH58J4e4OLwnh5cbu3ynh4AAICOitADAACsQOgBAABWiPpGZuBq1xHvJwCuBI4V2IYrPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsELclV4B2KnP3Dev9CpcVdprf3y6OK9d5gsAHRFXegAAgBUuKfQsX75cffr0UWJiojwej7Zu3Xre+rVr16pv375KTEzUwIEDtX79+oh2Y4xKS0uVlpampKQk+Xw+7d+/P6KmqalJRUVFcrlcSklJUXFxsY4dO3bWfJ599lnddtttSkhIUK9evfT0009fyiYCAIBOJurQs2bNGpWUlGjhwoXasWOHBg0apNzcXB06dKjN+i1btmjChAkqLi7Wzp07lZ+fr/z8fO3evdupKSsr09KlS1VeXq7a2lp169ZNubm5On78uFNTVFSkPXv2qLq6WpWVldq8ebOmTp0asawZM2boP//zP/Xss8/qww8/1Ouvv65hw4ZFu4kAAKATijHGmGgm8Hg8Gjp0qJYtWyZJCofDysjI0PTp0zV37tyz6gsKCtTS0qLKykpn3PDhw5Wdna3y8nIZY5Senq5Zs2Zp9uzZkqRgMCi3262KigoVFhZq79696t+/v7Zt26acnBxJUlVVle6//34dPHhQ6enp2rt3r+68807t3r1bt99++yXtjFAopOTkZAWDQblcrkuaBy4O9/RcHtzTg/PhOPwax0rHdrHn76iu9Jw4cUJ1dXXy+XxfzyA2Vj6fT36/v81p/H5/RL0k5ebmOvX19fUKBAIRNcnJyfJ4PE6N3+9XSkqKE3gkyefzKTY2VrW1tZKkN954QzfffLMqKyuVmZmpPn366Cc/+YmamprOuT2tra0KhUIRAwAA6JyiCj1HjhzR6dOn5Xa7I8a73W4FAoE2pwkEAuetP/PzQjU9e/aMaI+Li1P37t2dmk8++UR//OMftXbtWr388suqqKhQXV2dvv/9759zexYtWqTk5GRnyMjIuNAuAAAAHVSneXorHA6rtbVVL7/8su6++27dc889+sUvfqGNGzdq3759bU4zb948BYNBZzhw4MBlXmsAAHC5RBV6evTooS5duqixsTFifGNjo1JTU9ucJjU19bz1Z35eqObbN0qfOnVKTU1NTk1aWpri4uJ02223OTX9+vWTJDU0NLS5bgkJCXK5XBEDAADonKIKPfHx8RoyZIhqamqcceFwWDU1NfJ6vW1O4/V6I+olqbq62qnPzMxUampqRE0oFFJtba1T4/V61dzcrLq6Oqdmw4YNCofD8ng8kqTvfOc7OnXqlD7++GOn5g9/+IMkqXfv3tFsJgAA6ISifiNzSUmJJk2apJycHA0bNkzPPfecWlpaNHnyZEnSxIkT1atXLy1atEjSV4+Rjx49WkuWLFFeXp5Wr16t7du3a8WKFZKkmJgYzZw5U0899ZSysrKUmZmpJ554Qunp6crPz5f01RWbsWPHasqUKSovL9fJkyc1bdo0FRYWKj09XdJXNzbfdddd+vGPf6znnntO4XBYjz76qP72b/824uoPAACwU9Shp6CgQIcPH1ZpaakCgYCys7NVVVXl3Ijc0NCg2NivLyCNGDFCq1at0oIFCzR//nxlZWVp3bp1GjBggFMzZ84ctbS0aOrUqWpubtbIkSNVVVWlxMREp2blypWaNm2axowZo9jYWI0fP15Lly512mNjY/XGG29o+vTpGjVqlLp166b77rtPS5YsuaQdAwCwR3s+vs/j8FePqN/T05nxnp7Lh/eDXB582OJ8OA4vD47D9tcu7+kBAADoqAg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwApxV3oFAHRMfea+2S7z/XRxXrvMFwC40gMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAVuBrKIBOrL2+KgIAOiKu9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKcVd6BQAA59Zn7ptXehWAToMrPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAVrik0LN8+XL16dNHiYmJ8ng82rp163nr165dq759+yoxMVEDBw7U+vXrI9qNMSotLVVaWpqSkpLk8/m0f//+iJqmpiYVFRXJ5XIpJSVFxcXFOnbsWJvL++ijj3TdddcpJSXlUjYPAAB0QlGHnjVr1qikpEQLFy7Ujh07NGjQIOXm5urQoUNt1m/ZskUTJkxQcXGxdu7cqfz8fOXn52v37t1OTVlZmZYuXary8nLV1taqW7duys3N1fHjx52aoqIi7dmzR9XV1aqsrNTmzZs1derUs5Z38uRJTZgwQXfffXe0mwYAADqxGGOMiWYCj8ejoUOHatmyZZKkcDisjIwMTZ8+XXPnzj2rvqCgQC0tLaqsrHTGDR8+XNnZ2SovL5cxRunp6Zo1a5Zmz54tSQoGg3K73aqoqFBhYaH27t2r/v37a9u2bcrJyZEkVVVV6f7779fBgweVnp7uzPvxxx/XZ599pjFjxmjmzJlqbm6+6G0LhUJKTk5WMBiUy+WKZrcgSrxwDefy6eK8K70KVxWOlY6P/9Pt72LP31Fd6Tlx4oTq6urk8/m+nkFsrHw+n/x+f5vT+P3+iHpJys3Nderr6+sVCAQiapKTk+XxeJwav9+vlJQUJ/BIks/nU2xsrGpra51xGzZs0Nq1a7V8+fKL2p7W1laFQqGIAQAAdE5RhZ4jR47o9OnTcrvdEePdbrcCgUCb0wQCgfPWn/l5oZqePXtGtMfFxal79+5Ozeeff66HHnpIFRUVF32VZtGiRUpOTnaGjIyMi5oOAAB0PJ3m6a0pU6bohz/8oUaNGnXR08ybN0/BYNAZDhw40I5rCAAArqSoQk+PHj3UpUsXNTY2RoxvbGxUampqm9Okpqaet/7MzwvVfPtG6VOnTqmpqcmp2bBhg5599lnFxcUpLi5OxcXFCgaDiouL0y9/+cs21y0hIUEulytiAAAAnVNUoSc+Pl5DhgxRTU2NMy4cDqumpkZer7fNabxeb0S9JFVXVzv1mZmZSk1NjagJhUKqra11arxer5qbm1VXV+fUbNiwQeFwWB6PR9JX9/3s2rXLGZ588kldd9112rVrlx544IFoNhMAAHRCcdFOUFJSokmTJiknJ0fDhg3Tc889p5aWFk2ePFmSNHHiRPXq1UuLFi2SJM2YMUOjR4/WkiVLlJeXp9WrV2v79u1asWKFJCkmJkYzZ87UU089paysLGVmZuqJJ55Qenq68vPzJUn9+vXT2LFjNWXKFJWXl+vkyZOaNm2aCgsLnSe3+vXrF7Ge27dvV2xsrAYMGHDJOwcAAHQeUYeegoICHT58WKWlpQoEAsrOzlZVVZVzI3JDQ4NiY7++gDRixAitWrVKCxYs0Pz585WVlaV169ZFhJE5c+aopaVFU6dOVXNzs0aOHKmqqiolJiY6NStXrtS0adM0ZswYxcbGavz48Vq6dOlfs+0AAMAiUb+npzPjPT2XD+8ewbnwTpNIHCsdH/+n21+7vKcHAACgoyL0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGCFuCu9Ari68Q3PAIDOgis9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWCHuSq8AAACdWZ+5b7bLfD9dnNcu8+3MuNIDAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACpcUepYvX64+ffooMTFRHo9HW7duPW/92rVr1bdvXyUmJmrgwIFav359RLsxRqWlpUpLS1NSUpJ8Pp/2798fUdPU1KSioiK5XC6lpKSouLhYx44dc9o3bdqkcePGKS0tTd26dVN2drZWrlx5KZsHAAA6oahDz5o1a1RSUqKFCxdqx44dGjRokHJzc3Xo0KE267ds2aIJEyaouLhYO3fuVH5+vvLz87V7926npqysTEuXLlV5eblqa2vVrVs35ebm6vjx405NUVGR9uzZo+rqalVWVmrz5s2aOnVqxHLuvPNO/fd//7fee+89TZ48WRMnTlRlZWW0mwgAADqhGGOMiWYCj8ejoUOHatmyZZKkcDisjIwMTZ8+XXPnzj2rvqCgQC0tLRHhY/jw4crOzlZ5ebmMMUpPT9esWbM0e/ZsSVIwGJTb7VZFRYUKCwu1d+9e9e/fX9u2bVNOTo4kqaqqSvfff78OHjyo9PT0Ntc1Ly9Pbrdbv/zlLy9q20KhkJKTkxUMBuVyuaLZLZ1Wn7lvXulVgGU+XZx3pVfhqsIxiHPhWPnaxZ6/o7rSc+LECdXV1cnn8309g9hY+Xw++f3+Nqfx+/0R9ZKUm5vr1NfX1ysQCETUJCcny+PxODV+v18pKSlO4JEkn8+n2NhY1dbWnnN9g8Ggunfvfs721tZWhUKhiAEAAHROUYWeI0eO6PTp03K73RHj3W63AoFAm9MEAoHz1p/5eaGanj17RrTHxcWpe/fu51zuq6++qm3btmny5Mnn3J5FixYpOTnZGTIyMs5ZCwAAOrZO+fTWxo0bNXnyZL3wwgu64447zlk3b948BYNBZzhw4MBlXEsAAHA5RRV6evTooS5duqixsTFifGNjo1JTU9ucJjU19bz1Z35eqObbN0qfOnVKTU1NZy337bff1t/93d/p3/7t3zRx4sTzbk9CQoJcLlfEAAAAOqeoQk98fLyGDBmimpoaZ1w4HFZNTY28Xm+b03i93oh6SaqurnbqMzMzlZqaGlETCoVUW1vr1Hi9XjU3N6uurs6p2bBhg8LhsDwejzNu06ZNysvL0zPPPBPxZBcAAEBctBOUlJRo0qRJysnJ0bBhw/Tcc8+ppaXFuXdm4sSJ6tWrlxYtWiRJmjFjhkaPHq0lS5YoLy9Pq1ev1vbt27VixQpJUkxMjGbOnKmnnnpKWVlZyszM1BNPPKH09HTl5+dLkvr166exY8dqypQpKi8v18mTJzVt2jQVFhY6T25t3LhR3/ve9zRjxgyNHz/eudcnPj7+vDczAwAAO0QdegoKCnT48GGVlpYqEAgoOztbVVVVzo3IDQ0Nio39+gLSiBEjtGrVKi1YsEDz589XVlaW1q1bpwEDBjg1c+bMUUtLi6ZOnarm5maNHDlSVVVVSkxMdGpWrlypadOmacyYMYqNjdX48eO1dOlSp/2ll17SF198oUWLFjmBS5JGjx6tTZs2RbuZADqh9nz8m8eHcbnx/zl6Ub+npzPjPT1n4x0huNza88O2I54kOAZxJXS00NMu7+kBAADoqAg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVov4aCgBoTx31DcQddb0Bm3ClBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgTcyAwCACO31hvFPF+e1y3wvFld6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAKhB4AAGAFQg8AALACoQcAAFiB0AMAAKxA6AEAAFYg9AAAACsQegAAgBUIPQAAwAqEHgAAYAVCDwAAsAKhBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsEHelV8Amfea+eaVXAQAAa3GlBwAAWIHQAwAArEDoAQAAViD0AAAAKxB6AACAFS4p9Cxfvlx9+vRRYmKiPB6Ptm7det76tWvXqm/fvkpMTNTAgQO1fv36iHZjjEpLS5WWlqakpCT5fD7t378/oqapqUlFRUVyuVxKSUlRcXGxjh07FlHz3nvv6e6771ZiYqIyMjJUVlZ2KZsHAAA6oahDz5o1a1RSUqKFCxdqx44dGjRokHJzc3Xo0KE267ds2aIJEyaouLhYO3fuVH5+vvLz87V7926npqysTEuXLlV5eblqa2vVrVs35ebm6vjx405NUVGR9uzZo+rqalVWVmrz5s2aOnWq0x4KhXTvvfeqd+/eqqur089+9jP99Kc/1YoVK6LdRAAA0AnFGGNMNBN4PB4NHTpUy5YtkySFw2FlZGRo+vTpmjt37ln1BQUFamlpUWVlpTNu+PDhys7OVnl5uYwxSk9P16xZszR79mxJUjAYlNvtVkVFhQoLC7V37171799f27ZtU05OjiSpqqpK999/vw4ePKj09HQ9//zz+ud//mcFAgHFx8dLkubOnat169bpww8/vKhtC4VCSk5OVjAYlMvlima3XBTe0wMAsNmni/PaZb4Xe/6O6uWEJ06cUF1dnebNm+eMi42Nlc/nk9/vb3Mav9+vkpKSiHG5ublat26dJKm+vl6BQEA+n89pT05Olsfjkd/vV2Fhofx+v1JSUpzAI0k+n0+xsbGqra3VAw88IL/fr1GjRjmB58xynnnmGf3lL3/R9ddff9a6tba2qrW11fk9GAxK+mrntYdw6xftMl8AADqC9jq/npnvha7jRBV6jhw5otOnT8vtdkeMd7vd57yaEggE2qwPBAJO+5lx56vp2bNn5IrHxal79+4RNZmZmWfN40xbW6Fn0aJF+pd/+ZezxmdkZLS5LQAA4NIlP9e+8z969KiSk5PP2W7111DMmzcv4ipUOBxWU1OTbrjhBsXExFzBNYtOKBRSRkaGDhw40C5/lsNfjz66utE/Vzf65+p3pfvIGKOjR48qPT39vHVRhZ4ePXqoS5cuamxsjBjf2Nio1NTUNqdJTU09b/2Zn42NjUpLS4uoyc7Odmq+faP0qVOn1NTUFDGftpbzzWV8W0JCghISEiLGpaSktFnbEbhcLj4QrnL00dWN/rm60T9XvyvZR+e7wnNGVE9vxcfHa8iQIaqpqXHGhcNh1dTUyOv1tjmN1+uNqJek6upqpz4zM1OpqakRNaFQSLW1tU6N1+tVc3Oz6urqnJoNGzYoHA7L4/E4NZs3b9bJkycjlnP77be3+actAABgGROl1atXm4SEBFNRUWE++OADM3XqVJOSkmICgYAxxpgHH3zQzJ0716l/5513TFxcnHn22WfN3r17zcKFC80111xj3n//fadm8eLFJiUlxbz22mvmvffeM+PGjTOZmZnmyy+/dGrGjh1rBg8ebGpra83vfvc7k5WVZSZMmOC0Nzc3G7fbbR588EGze/dus3r1atO1a1fzH//xH9FuYocTDAaNJBMMBq/0quAc6KOrG/1zdaN/rn4dpY+iDj3GGPPzn//c3HTTTSY+Pt4MGzbM/P73v3faRo8ebSZNmhRR/+qrr5rbbrvNxMfHmzvuuMO8+eabEe3hcNg88cQTxu12m4SEBDNmzBizb9++iJrPP//cTJgwwVx77bXG5XKZyZMnm6NHj0bUvPvuu2bkyJEmISHB9OrVyyxevPhSNq/DOX78uFm4cKE5fvz4lV4VnAN9dHWjf65u9M/Vr6P0UdTv6QEAAOiI+O4tAABgBUIPAACwAqEHAABYgdADAACsQOi5Svz0pz9VTExMxNC3b1+n/Z577jmr/eGHH46YR0NDg/Ly8tS1a1f17NlTjz32mE6dOhVRs2nTJt11111KSEjQrbfeqoqKisuxeZ3Cn/70J/3oRz/SDTfcoKSkJA0cOFDbt2932o0xKi0tVVpampKSkuTz+bR///6IeTQ1NamoqEgul0spKSkqLi7WsWPHImree+893X333UpMTFRGRobKysouy/Z1dBfqn4ceeuisY2js2LER86B/2k+fPn3O2v8xMTF69NFHJUnHjx/Xo48+qhtuuEHXXnutxo8ff9YLZ/mMaz8X6p9Ocw66wk+P4f9buHChueOOO8yf//xnZzh8+LDTPnr0aDNlypSI9m++D+HUqVNmwIABxufzmZ07d5r169ebHj16mHnz5jk1n3zyienataspKSkxH3zwgfn5z39uunTpYqqqqi7rtnZETU1Npnfv3uahhx4ytbW15pNPPjFvvfWW+eijj5yaxYsXm+TkZLNu3Trz7rvvmr//+79v831TgwYNMr///e/N//7v/5pbb7014n1TwWDQuN1uU1RUZHbv3m1+/etfm6SkJCveN/XXuJj+mTRpkhk7dmzEMdTU1BQxH/qn/Rw6dChi31dXVxtJZuPGjcYYYx5++GGTkZFhampqzPbt283w4cPNiBEjnOn5jGtfF+qfznIOIvRcJRYuXGgGDRp0zvbRo0ebGTNmnLN9/fr1JjY21nlJpDHGPP/888blcpnW1lZjjDFz5swxd9xxR8R0BQUFJjc3969adxs8/vjjZuTIkedsD4fDJjU11fzsZz9zxjU3N5uEhATz61//2hhjzAcffGAkmW3btjk1v/nNb0xMTIz505/+ZIwx5t///d/N9ddf7/TZmWXffvvt/9eb1KlcqH+M+Sr0jBs37pzt9M/lNWPGDHPLLbeYcDhsmpubzTXXXGPWrl3rtO/du9dIMn6/3xjDZ9zl9s3+MabznIP489ZVZP/+/UpPT9fNN9+soqIiNTQ0RLSvXLlSPXr00IABAzRv3jx98cUXTpvf79fAgQMjvq0+NzdXoVBIe/bscWp8Pl/EPHNzc+X3+9txqzqH119/XTk5OfrBD36gnj17avDgwXrhhRec9vr6egUCgYj9m5ycLI/H4+xfv9+vlJQU5eTkODU+n0+xsbGqra11akaNGqX4+HinJjc3V/v27dNf/vKX9t7MDutC/XPGpk2b1LNnT91+++165JFH9Pnnnztt9M/lc+LECb3yyiv68Y9/rJiYGNXV1enkyZMRx0/fvn110003RRw/fMZdHt/unzM6wzmI0HOV8Hg8qqioUFVVlZ5//nnV19fr7rvv1tGjRyVJP/zhD/XKK69o48aNmjdvnn71q1/pRz/6kTN9IBCI+M8myfk9EAictyYUCunLL79sz83r8D755BM9//zzysrK0ltvvaVHHnlE//RP/6SXXnpJ0tf7uK39+83937Nnz4j2uLg4de/e/YJ99M1l4GwX6h9JGjt2rF5++WXV1NTomWee0dtvv6377rtPp0+flkT/XE7r1q1Tc3OzHnroIUlf7bv4+PizvvD528cPn3GXx7f7R+o856CovmUd7ee+++5z/n3nnXfK4/God+/eevXVV1VcXKypU6c67QMHDlRaWprGjBmjjz/+WLfccsuVWGWrhMNh5eTk6F//9V8lSYMHD9bu3btVXl6uSZMmXeG1w8X0T2FhoVM/cOBA3Xnnnbrlllu0adMmjRkz5oqst61+8Ytf6L777lN6evqVXhW0oa3+6SznIK70XKVSUlJ022236aOPPmqz/cy3y59pT01NPetJhzO/p6amnrfG5XIpKSnp/3T9O5u0tDT1798/Yly/fv2cP0Ge2cdt7d9v7v9Dhw5FtJ86dUpNTU0X7KNvLgNnu1D/tOXmm29Wjx49Io4h+qf9/fGPf9Rvf/tb/eQnP3HGpaam6sSJE2pubo6o/fbxw2dc+2urf9rSUc9BhJ6r1LFjx/Txxx8rLS2tzfZdu3ZJktPu9Xr1/vvvR3xoV1dXy+VyOScDr9ermpqaiPlUV1fL6/W2wxZ0Lt/5zne0b9++iHF/+MMf1Lt3b0lSZmamUlNTI/ZvKBRSbW2ts3+9Xq+am5tVV1fn1GzYsEHhcNj5APF6vdq8ebNOnjzp1FRXV+v222/X9ddf327b19FdqH/acvDgQX3++ecRxxD90/5efPFF9ezZU3l5ec64IUOG6Jprrok4fvbt26eGhoaI44fPuPbXVv+0pcOegy7bLdM4r1mzZplNmzaZ+vp688477xifz2d69OhhDh06ZD766CPz5JNPmu3bt5v6+nrz2muvmZtvvtmMGjXKmf7M44L33nuv2bVrl6mqqjI33nhjm48LPvbYY2bv3r1m+fLlPM55kbZu3Wri4uLM008/bfbv329Wrlxpunbtal555RWnZvHixSYlJcW89tpr5r333jPjxo1r85H1wYMHm9raWvO73/3OZGVlRTwS3dzcbNxut3nwwQfN7t27zerVq03Xrl15JPoCLtQ/R48eNbNnzzZ+v9/U19eb3/72t+auu+4yWVlZEd8KTf+0r9OnT5ubbrrJPP7442e1Pfzww+amm24yGzZsMNu3bzder9d4vV6nnc+49neu/ulM5yBCz1WioKDApKWlmfj4eNOrVy9TUFDgvGOkoaHBjBo1ynTv3t0kJCSYW2+91Tz22GMR70gwxphPP/3U3HfffSYpKcn06NHDzJo1y5w8eTKiZuPGjSY7O9vEx8ebm2++2bz44ouXaxM7vDfeeMMMGDDAJCQkmL59+5oVK1ZEtIfDYfPEE08Yt9ttEhISzJgxY8y+ffsiaj7//HMzYcIEc+211xqXy2UmT55sjh49GlHz7rvvmpEjR5qEhATTq1cvs3jx4nbfts7gfP3zxRdfmHvvvdfceOON5pprrjG9e/c2U6ZMiXi81hj6p7299dZbRtJZx4Uxxnz55ZfmH//xH831119vunbtah544AHz5z//OaKGz7j2da7+6UznoBhjjLl815UAAACuDO7pAQAAViD0AAAAKxB6AACAFQg9AADACoQeAABgBUIPAACwAqEHAABYgdADAACsQOgBAABWIPQAAAArEHoAAIAVCD0AAMAK/w9o7PHnvsHhBAAAAABJRU5ErkJggg==", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.hist(total_impulse_samples, density = True, bins = 20)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "dev-py3-12", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.2" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/tests/conftest.py b/tests/conftest.py index 370721cf8..980e0b6ac 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -18,6 +18,7 @@ "tests.fixtures.units.numerical_fixtures", "tests.fixtures.monte_carlo.monte_carlo_fixtures", "tests.fixtures.monte_carlo.stochastic_fixtures", + "tests.fixtures.monte_carlo.custom_sampler_fixtures", "tests.fixtures.monte_carlo.stochastic_motors_fixtures", "tests.fixtures.sensors.sensors_fixtures", "tests.fixtures.generic_surfaces.generic_surfaces_fixtures", diff --git a/tests/fixtures/monte_carlo/custom_sampler_fixtures.py b/tests/fixtures/monte_carlo/custom_sampler_fixtures.py new file mode 100644 index 000000000..8a4ff497d --- /dev/null +++ b/tests/fixtures/monte_carlo/custom_sampler_fixtures.py @@ -0,0 +1,76 @@ +"""This file contains fixtures of CustomSampler used in stochastic classes.""" + +import numpy as np +import pytest + +from rocketpy import CustomSampler + + +@pytest.fixture +def elevation_sampler(): + """Fixture to create mixture of two gaussian sampler""" + means_tuple = [1400, 1500] + sd_tuple = [40, 50] + prob_tuple = [0.4, 0.6] + return TwoGaussianMixture(means_tuple, sd_tuple, prob_tuple) + + +class TwoGaussianMixture(CustomSampler): + """Class to sample from a mixture of two Gaussian distributions""" + + def __init__(self, means_tuple, sd_tuple, prob_tuple, seed=None): + """Creates a sampler for a mixture of two Gaussian distributions + + Parameters + ---------- + means_tuple : 2-tuple + 2-Tuple that contains the mean of each normal distribution of the + mixture + sd_tuple : 2-tuple + 2-Tuple that contains the sd of each normal distribution of the + mixture + prob_tuple : 2-tuple + 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.means_tuple = means_tuple + self.sd_tuple = sd_tuple + self.prob_tuple = prob_tuple + + def sample(self, n_samples=1): + """Samples from a mixture of two Gaussian + + Parameters + ---------- + n_samples : int, optional + Number of samples, by default 1 + + Returns + ------- + samples_list + List containing n_samples samples + """ + samples_list = [0] * n_samples + mixture_id_list = np.random.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] + ) + else: + samples_list[i] = np.random.normal( + self.means_tuple[1], self.sd_tuple[1] + ) + + return samples_list + + def reset_seed(self, seed=None): + """Resets all associated random number generators + + Parameters + ---------- + seed : int, optional + Seed for the random number generator. + """ + np.random.default_rng(seed) diff --git a/tests/fixtures/monte_carlo/stochastic_fixtures.py b/tests/fixtures/monte_carlo/stochastic_fixtures.py index bf576e5ed..6610666cf 100644 --- a/tests/fixtures/monte_carlo/stochastic_fixtures.py +++ b/tests/fixtures/monte_carlo/stochastic_fixtures.py @@ -43,6 +43,36 @@ def stochastic_environment(example_spaceport_env): ) +@pytest.fixture +def stochastic_environment_custom_sampler(example_spaceport_env, elevation_sampler): + """This fixture is used to create a stochastic environment object for the + Calisto flight using a custom sampler for the elevation. + + Parameters + ---------- + example_spaceport_env : Environment + This is another fixture. + + elevation_sampler: CustomSampler + This is another fixture. + + Returns + ------- + StochasticEnvironment + The stochastic environment object + """ + return StochasticEnvironment( + environment=example_spaceport_env, + elevation=elevation_sampler, + gravity=None, + latitude=None, + longitude=None, + ensemble_member=None, + wind_velocity_x_factor=(1.0, 0.033, "normal"), + wind_velocity_y_factor=(1.0, 0.033, "normal"), + ) + + @pytest.fixture def stochastic_nose_cone(calisto_nose_cone): """This fixture is used to create a StochasticNoseCone object for the diff --git a/tests/unit/stochastic/test_custom_sampler.py b/tests/unit/stochastic/test_custom_sampler.py new file mode 100644 index 000000000..90774ac50 --- /dev/null +++ b/tests/unit/stochastic/test_custom_sampler.py @@ -0,0 +1,21 @@ +from rocketpy.environment.environment import Environment + + +def test_create_object(stochastic_environment_custom_sampler): + """Test create object method of StochasticEnvironment class. + + This test checks if the create_object method of the StochasticEnvironment + class creates a StochasticEnvironment object from the randomly generated + input arguments. + + Parameters + ---------- + stochastic_environment : StochasticEnvironment + StochasticEnvironment object to be tested. + + Returns + ------- + None + """ + obj = stochastic_environment_custom_sampler.create_object() + assert isinstance(obj, Environment) diff --git a/tests/unit/test_stochastic_model.py b/tests/unit/test_stochastic_model.py index 77c94fb40..0d0a13311 100644 --- a/tests/unit/test_stochastic_model.py +++ b/tests/unit/test_stochastic_model.py @@ -7,6 +7,7 @@ "stochastic_rail_buttons", "stochastic_main_parachute", "stochastic_environment", + "stochastic_environment_custom_sampler", "stochastic_tail", "stochastic_calisto", ], From dcdb9618b3431c75b91820ccecf31225b266f6f6 Mon Sep 17 00:00:00 2001 From: Lucas de Oliveira Prates Date: Thu, 17 Apr 2025 22:23:13 +0200 Subject: [PATCH 3/4] DOC: creating documentation for CustomSampler --- .../stochastic_models/custom_sampler.rst | 5 + .../monte_carlo/stochastic_models/index.rst | 1 + docs/user/custom_sampler.rst | 354 ++++++++++++++++++ docs/user/index.rst | 1 + docs/user/stochastic.rst | 5 + 5 files changed, 366 insertions(+) create mode 100644 docs/reference/classes/monte_carlo/stochastic_models/custom_sampler.rst create mode 100644 docs/user/custom_sampler.rst diff --git a/docs/reference/classes/monte_carlo/stochastic_models/custom_sampler.rst b/docs/reference/classes/monte_carlo/stochastic_models/custom_sampler.rst new file mode 100644 index 000000000..54e612c3d --- /dev/null +++ b/docs/reference/classes/monte_carlo/stochastic_models/custom_sampler.rst @@ -0,0 +1,5 @@ +Custom Sampler +--------------------------- + +.. autoclass:: rocketpy.stochastic.CustomSampler + :members: \ No newline at end of file diff --git a/docs/reference/classes/monte_carlo/stochastic_models/index.rst b/docs/reference/classes/monte_carlo/stochastic_models/index.rst index 7bf705347..ca8b2b1e2 100644 --- a/docs/reference/classes/monte_carlo/stochastic_models/index.rst +++ b/docs/reference/classes/monte_carlo/stochastic_models/index.rst @@ -24,3 +24,4 @@ input parameters, enabling robust Monte Carlo simulations. stochastic_rocket stochastic_parachute stochastic_flight + custom_sampler diff --git a/docs/user/custom_sampler.rst b/docs/user/custom_sampler.rst new file mode 100644 index 000000000..8994c9573 --- /dev/null +++ b/docs/user/custom_sampler.rst @@ -0,0 +1,354 @@ +.. _custom_sampler: + +Implementing custom sampler for Stochastic objects +================================================== + +The :ref:`stochastic_usage` documentation teaches how to work with stochastic +objects, discussing the standard initializations, how to create objects and interpret +outputs. Our goal here is to show how to build a custom sampler that gives complete +control of the distributions used. + +Custom Sampler +-------------- + +Rocketpy provides a ``CustomSampler`` abstract class which works as the backbone of +a custom sampler. We begin by first importing it and some other useful modules + +.. jupyter-execute:: + + from rocketpy import CustomSampler + from matplotlib import pyplot as plt + import numpy as np + +In order to use it, we must create a new class that inherits from +it and it **must** implement two methods: *sample* and *reset_seed*. The *sample* +method has one argument, *n_samples*, and must return a list with *n_samples* +entries, each of which is a sample of the desired variable. The *reset_seed* method +has one argument, *seed*, which is used to reset the pseudorandom generators in order +to avoid unwanted dependency across samples. This is especially important when the +``MonteCarlo`` class is used in parallel mode. + +Below, we give an example of how to implement a mixture of two Gaussian +distributions. + +.. jupyter-execute:: + + class TwoGaussianMixture(CustomSampler): + """Class to sample from a mixture of two Gaussian distributions + """ + + def __init__(self, means_tuple, sd_tuple, prob_tuple, seed = None): + """ Creates a sampler for a mixture of two Gaussian distributions + + Parameters + ---------- + means_tuple : 2-tuple + 2-Tuple that contains the mean of each normal distribution of the + mixture + sd_tuple : 2-tuple + 2-Tuple that contains the sd of each normal distribution of the + mixture + prob_tuple : 2-tuple + 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.means_tuple = means_tuple + self.sd_tuple = sd_tuple + self.prob_tuple = prob_tuple + + def sample(self, n_samples = 1): + """Samples from a mixture of two Gaussian + + Parameters + ---------- + n_samples : int, optional + Number of samples, by default 1 + + Returns + ------- + samples_list + List containing n_samples samples + """ + samples_list = [0] * n_samples + mixture_id_list = np.random.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]) + else: + samples_list[i] = np.random.normal(self.means_tuple[1], self.sd_tuple[1]) + + return samples_list + + def reset_seed(self, seed=None): + """Resets all associated random number generators + + Parameters + ---------- + seed : int, optional + Seed for the random number generator. + """ + np.random.default_rng(seed) + +This is an example of a distribution that is not implemented in numpy. Note that it is +a general distribution, so we can use it for many different variables. + +.. note:: + You can check more information about the mixture of Gaussian distributions + `here `. + Intuitively, if you think of a single Gaussian as a bell curve distribution, + the mixture distribution resembles two bell curves superimposed, each with mode at their + respective mean. + +Example: Gaussian Mixture for Total Impulse +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In order to use the new created sampler in a stochastic object, we first need +to build an object. In this example, we will consider a case where the distribution of +the total impulse is a mixture of two gaussian with mean parameters +:math:`(6000, 7000)`, standard deviations :math:`(300, 100)` and mixture probabilities +:math:`(0.7, 0.3)`. + +.. jupyter-execute:: + + means_tuple = (6000, 7000) + sd_tuple = (300, 100) + prob_tuple = (0.7, 0.3) + total_impulse_sampler = TwoGaussianMixture(means_tuple, sd_tuple, prob_tuple) + +Finally, we can create ``StochasticSolidMotor`` object as we did in the example of +:ref:`stochastic_usage`, but we pass the sampler object instead for the *total_impulse* +argument + +.. jupyter-execute:: + + from rocketpy import SolidMotor, StochasticSolidMotor + + motor = SolidMotor( + thrust_source="../data/motors/cesaroni/Cesaroni_M1670.eng", + dry_mass=1.815, + dry_inertia=(0.125, 0.125, 0.002), + nozzle_radius=33 / 1000, + grain_number=5, + grain_density=1815, + grain_outer_radius=33 / 1000, + grain_initial_inner_radius=15 / 1000, + grain_initial_height=120 / 1000, + grain_separation=5 / 1000, + grains_center_of_mass_position=0.397, + center_of_dry_mass_position=0.317, + nozzle_position=0, + burn_time=3.9, + throat_radius=11 / 1000, + coordinate_system_orientation="nozzle_to_combustion_chamber", + ) + + stochastic_motor = StochasticSolidMotor( + solid_motor=motor, + burn_start_time=(0, 0.1, "binomial"), + grains_center_of_mass_position=0.001, + grain_density=10, + grain_separation=1 / 1000, + grain_initial_height=1 / 1000, + grain_initial_inner_radius=0.375 / 1000, + grain_outer_radius=0.375 / 1000, + throat_radius=0.5 / 1000, + nozzle_radius=0.5 / 1000, + nozzle_position=0.001, + total_impulse=total_impulse_sampler, # total impulse using custom sampler! + ) + + stochastic_motor.visualize_attributes() + +Let's generate some random motors and check the distribution of the total impulse + +.. jupyter-execute:: + + total_impulse_samples = [ + stochastic_motor.create_object().total_impulse for _ in range(200) + ] + plt.hist(total_impulse_samples, density = True, bins = 30) + +Introducing dependency between parameters +----------------------------------------- + +Although probabilistic **independency between samples**, i.e. different flight runs, +is desired for Monte Carlo simulations, it is often important to be able to introduce +**dependency between parameters**. A clear example of this is wind speed: if we know +the wind speed in the x-axis, then our forecast model might tells us that the wind +speed y-axis is more likely to be positive than negative, or vice-versa. These +parameters are then correlated, and using samplers that model these correlations make +the Monte Carlo analysis more robust. + +When we use the default numpy samplers, the Monte Carlo analysis samples the +parameters independently from each other. However, using custom samplers, we can +introduce dependency and correlation! It might be a bit tricky, but we will show how +it can be done. First, let us import the modules required + +.. jupyter-execute:: + + from rocketpy import Environment, StochasticEnvironment + from datetime import datetime, timedelta + +Assume we want to model the x and y axis wind speed as a Bivariate Gaussian with +parameters :math:`\mu = (1, 1)` and variance-covariance matrix +:math:`\Sigma = \begin{bmatrix} 0.2 & 0.17 \\ 0.17 & 0.3 \end{bmatrix}`. This will +make the correlation between the speeds be of :math:`0.7`. + +Now, in order to correlate the parameters using different custom samplers, +**the key trick is to create a common generator that will be used by both.** The code +below implements an example of such a generator + +.. jupyter-execute:: + + class BivariateGaussianGenerator: + """Bivariate Gaussian generator used across custom samplers + """ + def __init__(self, mean, cov, seed = None): + """Initializes the generator + + Parameters + ---------- + mean : tuple, list + Tuple or list with mean of bivariate Gaussian + cov : np.array + Variance-Covariance matrix + seed : int, optional + Number to seed random generator, by default None + """ + np.random.default_rng(seed) + self.samples_list = [] + self.samples_generated = 0 + self.used_samples_x = 0 + self.used_samples_y = 0 + self.mean = mean + self.cov = cov + self.generate_samples(1000) + + def generate_samples(self, n_samples = 1): + """Generate samples from bivariate Gaussian and append to sample list + + Parameters + ---------- + n_samples : int, optional + Number of samples to be generated, by default 1 + """ + samples_generated = np.random.multivariate_normal(self.mean, self.cov, n_samples) + self.samples_generated += n_samples + self.samples_list += list(samples_generated) + + def reset_seed(self, seed=None): + np.random.default_rng(seed) + + def get_samples(self, n_samples, axis): + if axis == "x": + if self.samples_generated < self.used_samples_x: + self.generate_samples(n_samples) + samples_list = [ + sample[0] for sample in self.samples_list[self.used_samples_x:(self.used_samples_x + n_samples)] + ] + if axis == "y": + if self.samples_generated < self.used_samples_y: + self.generate_samples(n_samples) + samples_list = [ + sample[1] for sample in self.samples_list[self.used_samples_y:(self.used_samples_y + n_samples)] + ] + self.update_info(n_samples, axis) + return samples_list + + def update_info(self, n_samples, axis): + """Updates the information of the used samples + + Parameters + ---------- + n_samples : int + Number of samples used + axis : str + Which axis was sampled + """ + if axis == "x": + self.used_samples_x += n_samples + if axis == "y": + self.used_samples_y += n_samples + +This generator samples from the bivariate Gaussian and stores then in a *samples_list* +attribute. Then, the idea is to create two samplers for the wind speeds that will share +an object of this class and their sampling methods only get samples from the stored +sample list. + +.. jupyter-execute:: + + class WindXSampler(CustomSampler): + """Samples from X""" + + def __init__(self, bivariate_gaussian_generator): + self.generator = bivariate_gaussian_generator + + def sample(self, n_samples=1): + samples_list = self.generator.get_samples(n_samples, "x") + return samples_list + + def reset_seed(self, seed=None): + self.generator.reset_seed(seed) + + class WindYSampler(CustomSampler): + """Samples from Y""" + + def __init__(self, bivariate_gaussian_generator): + self.generator = bivariate_gaussian_generator + + def sample(self, n_samples=1): + samples_list = self.generator.get_samples(n_samples, "y") + return samples_list + + def reset_seed(self, seed=None): + self.generator.reset_seed(seed) + +Then, we create the objects + +.. jupyter-execute:: + + mean = [1, 2] + cov_mat = [[0.2, 0.171], [0.171, 0.3]] + bivariate_gaussian_generator = BivariateGaussianGenerator(mean, cov_mat) + wind_x_sampler = WindXSampler(bivariate_gaussian_generator) + wind_y_sampler = WindYSampler(bivariate_gaussian_generator) + +With the sample objects created, we only need to create the stochastic objects and +pass them as argument + +.. jupyter-execute:: + + spaceport_env = Environment( + latitude=32.990254, + longitude=-106.974998, + elevation=1400, + datum="WGS84", + ) + spaceport_env.set_atmospheric_model("custom_atmosphere", wind_u = 1, wind_v = 1) + spaceport_env.set_date(datetime.now() + timedelta(days=1)) + + stochastic_environment = StochasticEnvironment( + environment=spaceport_env, + elevation=(1400, 10, "normal"), + gravity=None, + latitude=None, + longitude=None, + ensemble_member=None, + wind_velocity_x_factor=wind_x_sampler, + wind_velocity_y_factor=wind_y_sampler + ) + +Finally, let us check that if there is a correlation between the wind speed in the +two axis + +.. jupyter-execute:: + + wind_velocity_x_samples = [0] * 200 + wind_velocity_y_samples = [0] * 200 + for i in range(200): + stochastic_environment.create_object() + wind_velocity_x_samples[i] = stochastic_environment.obj.wind_velocity_x(0) + wind_velocity_y_samples[i] = stochastic_environment.obj.wind_velocity_y(0) + + np.corrcoef(wind_velocity_x_samples, wind_velocity_y_samples) diff --git a/docs/user/index.rst b/docs/user/index.rst index f23eae25a..a09872ad5 100644 --- a/docs/user/index.rst +++ b/docs/user/index.rst @@ -32,6 +32,7 @@ RocketPy's User Guide :caption: Monte Carlo Simulations Stochastic Classes + Custom Sampler ../notebooks/monte_carlo_analysis/monte_carlo_class_usage.ipynb ../notebooks/monte_carlo_analysis/monte_carlo_analysis.ipynb ../notebooks/monte_carlo_analysis/parachute_drop_from_helicopter.ipynb diff --git a/docs/user/stochastic.rst b/docs/user/stochastic.rst index 0985ca1da..062b034e2 100644 --- a/docs/user/stochastic.rst +++ b/docs/user/stochastic.rst @@ -88,6 +88,11 @@ passed in a few different ways: used as the parameter value during the simulation. You cannot assign standard \ deviations when using lists, nor can you assign different distribution types. +5. **A CustomSampler object**: \ + An object from a class that inherits from ``CustomSampler``. This object \ + gives you the full control of how the samples are generated. See + :ref:`custom_sampler` for more details. + .. note:: In statistics, the terms "Normal" and "Gaussian" refer to the same type of \ distribution. This distribution is commonly used and is the default for the \ From f1f722005eb1b319502fb43367f5abf231444d00 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 28 Apr 2025 15:10:35 +0200 Subject: [PATCH 4/4] DEV: changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 90e4911f1..2a0c9b976 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,7 +31,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/). Attention: The newest changes should be on top --> ### Added - +- ENH: allow users to provide custom samplers [#803](https://github.com/RocketPy-Team/RocketPy/pull/803) - ENH: Implement Multivariate Rejection Sampling (MRS) [#738] (https://github.com/RocketPy-Team/RocketPy/pull/738) - ENH: Create a rocketpy file to store flight simulations [#800](https://github.com/RocketPy-Team/RocketPy/pull/800) - ENH: Support for the RSE file format has been added to the library [#798](https://github.com/RocketPy-Team/RocketPy/pull/798)