from pathlib import Path
from typing import Any, Literal
import numpy as np
import pyvista
import yaml
from matplotlib.axes import Axes
from numpy.random import Generator
from numpy.typing import ArrayLike
from tqdm import tqdm
from cratermaker.components.counting import Counting
from cratermaker.components.crater import Crater
from cratermaker.components.morphology import Morphology
from cratermaker.components.production import Production
from cratermaker.components.projectile import Projectile
from cratermaker.components.scaling import Scaling
from cratermaker.components.surface import Surface
from cratermaker.components.surface.hireslocal import HiResLocalSurface
from cratermaker.components.target import Target
from cratermaker.constants import _COMPONENT_NAMES, _CONFIG_FILE_NAME, _DSTD, FloatLike, PairOfFloats
from cratermaker.core.base import CratermakerBase, _convert_for_yaml
from cratermaker.utils.general_utils import _set_properties, format_large_units, parameter
[docs]
class Simulation(CratermakerBase):
"""
Creates a simulation of a crater population on a target body. It allows for the generation of craters based on a variety of parameters, including the target body, scaling laws, production functions, and morphology models.
Parameters
----------
target: Target or str, optional, default "Moon"
Name target body for the simulation, default is "Moon".
scaling : Scaling or str, optional
The projectile->crater size scaling model to use from the components library. The default is "montecarlo".
production: Production or str, optional
The production function model to use from the components library that defines the production function used to populate the surface with craters. If none provided,
then the default will be based on the target body, with the NeukumProduction crater-based scaling law used if the target
body is the Moon or Mars, the NeukumProduction projectile-based scaling law if the target body is Mercury, Venus, or
Earth, and a simple power law model otherwise.
morphology : str, optional
The model used to generate the morphology of the crater. If none provided, then the default will "basicmoon", which is similar to the one used by CTEM.
projectile : str, optional
The projectile model to use from the components library, which is used to generate the projectile properties for the simulation, such as velocity and density. The default is "asteroids" when target is Mercury, Venus, Earth, Moon, Mars, Ceres, or Vesta, and "comets" otherwise.
surface : str, optional
The name of the surface used for the surface. Default is "icosphere".
counting : Counting or str, optional
The crater counting model to use from the components library. Default is "depthcount".
simdir : str | Path
|simdir|
rng : numpy.random.Generator | None
|rng|
rng_seed : Any type allowed by the rng_seed argument of numpy.random.Generator, optional
|rng_seed|
rng_state : dict, optional
|rng_state|
ask_overwrite : bool, optional
If True, the user will be prompted before overwriting any existing files. Default is True.
reset : bool, optional
Flag to indicate whether to reset the simulation or resume from an old simulation. If False, the simulation will attempt to load the previous state from the config file. Default is False if `ask_overwrite=False` and a config file is detected, otherwise default is True.
do_counting : bool, optional
If True, the counting component will keep track of observable craters during the simulation. If False, emplaced craters that are large enoug to be observable are saved, but the observability of craters on the surface is not evaluated and observed craters are not tracked. Default is True.
save_actions: list[dict[str, dict]], optional
A dictionary of actions to perform when the save method is called. The keys are the names of the actions and the values are dictionaries of keyword arguments to pass to the corresponding component's save method. For example, if you want to automatically generate a hillshade plot every time the simulation is saved, you can pass `save_actions=[{"plot": {"plot_style": "hillshade", "cmap": "pink", "scalebar": True, "label": "Mars region simulation", "show": True, "save": True}}]`. This will call the surface's save method with the specified keyword arguments every time the simulation is saved. Default is to save a hillshade plot of the surface every time the simulation is saved. Set to None to disable or "default" to use the default save actions.
**kwargs : Any
|kwargs|, including those for component function constructors. Refer to the documentation of each component module for details.
"""
def __init__(
self,
*, # Enforce keyword-only arguments
target: Target | str | None = None,
scaling: Scaling | str | None = None,
production: Production | str | None = None,
morphology: Morphology | str | None = None,
projectile: Projectile | str | None = None,
surface: Surface | str | None = None,
counting: Counting | str | None = None,
simdir: str | Path | None = None,
rng: Generator | None = None,
rng_seed: int | None = None,
rng_state: dict | None = None,
reset: bool = None,
ask_overwrite: bool = True,
do_counting: bool = True,
save_actions: list[dict[str, dict]] | None | str = "default",
**kwargs: Any,
):
object.__setattr__(self, "_target", None)
object.__setattr__(self, "_scaling", None)
object.__setattr__(self, "_production", None)
object.__setattr__(self, "_morphology", None)
object.__setattr__(self, "_projectile", None)
object.__setattr__(self, "_surface", None)
object.__setattr__(self, "_counting", None)
object.__setattr__(self, "_interval", None)
object.__setattr__(self, "_elapsed_time", None)
object.__setattr__(self, "_time", None)
object.__setattr__(self, "_elapsed_n1", None)
object.__setattr__(self, "_config_readonly", True)
object.__setattr__(self, "_is_new", False)
object.__setattr__(self, "_do_quasimc", None)
object.__setattr__(self, "_max_crater_diameter_range", None)
super().__init__(simdir=simdir, rng=rng, rng_seed=rng_seed, rng_state=rng_state, ask_overwrite=ask_overwrite, **kwargs)
if reset is None:
if self.ask_overwrite and self.config_file.exists():
response = input(
"Old run detected. Enter y to reset the simulation and n to resume from the previous state. To disable this message and suppress all prompts about overwriting old files, pass `ask_overwrite=False` as an argument to Simulation() or enter 'a' to suppress prompts about overwriting files. (y/[N]/a): "
)
if response.lower() == "a":
print("All prompts about overwriting files will be suppressed for this and future runs.")
self.ask_overwrite = False
reset = True
elif response.lower() == "y":
print("Resetting simulation.")
reset = True
else:
reset = False
else:
reset = False
self.is_new = reset
object.__setattr__(self, "_config_readonly", not reset)
if not self.is_new and self.config_file.exists():
config_file = self.config_file
object.__setattr__(self, "_config_readonly", True)
else:
config_file = None
config_override = {}
for component in _COMPONENT_NAMES:
# Set to true if a local variable from the argument list with the component name is set to something other than None, otherwise false
config_override[component] = getattr(self, f"_{component}") is not None
_, unmatched = _set_properties(
self,
target=target,
scaling=scaling,
production=production,
morphology=morphology,
projectile=projectile,
surface=surface,
counting=counting,
config_file=config_file,
**vars(self.common_args),
)
for component in _COMPONENT_NAMES:
if config_override[component]:
# If the component is set to something other than None, then remove it from the unmatched dictionary
unmatched.pop(f"{component}_config", None)
production_config = unmatched.pop("production_config", {})
scaling_config = unmatched.pop("scaling_config", {})
surface_config = unmatched.pop("surface_config", {})
morphology_config = unmatched.pop("morphology_config", {})
target_config = unmatched.pop("target_config", {})
projectile_config = unmatched.pop("projectile_config", {})
counting_config = unmatched.pop("counting_config", {})
kwargs.update(unmatched)
kwargs = {**kwargs, **vars(self.common_args)}
target_config = {**target_config, **kwargs}
self.target = Target.maker(self.target, **target_config)
production_config = {**production_config, **kwargs}
self.production = Production.maker(self.production, target=self.target, **production_config)
projectile_config = {**projectile_config, **kwargs}
self.projectile = Projectile.maker(self.projectile, target=self.target, **projectile_config)
scaling_config = {**scaling_config, **kwargs}
self.scaling = Scaling.maker(
self.scaling,
target=self.target,
projectile=self.projectile,
**scaling_config,
)
surface_config = {
**surface_config,
**kwargs,
}
if "superdomain_scale_factor" not in surface_config:
surface_config["superdomain_scale_factor"] = (
None # This will trigger setting of the superdomain after the Morphology and Scaling models are set
)
self.surface = Surface.maker(
self.surface,
target=self.target,
reset=self.is_new,
**surface_config,
)
# If the surface had to be regridded, then we will need to reset the simulation. HiResLocal surfaces are deferred, so won't have the flag set
if self.surface.is_new is not None:
self.is_new = self.surface.is_new
counting_config = {**counting_config, **kwargs}
self.counting = Counting.maker(
self.counting,
surface=self.surface,
reset=self.is_new,
**counting_config,
)
morphology_config = {**morphology_config, **kwargs, "do_counting": do_counting}
self.morphology = Morphology.maker(
self.morphology,
surface=self.surface,
production=self.production,
counting=self.counting,
scaling=self.scaling,
**morphology_config,
)
# If this is a variant of the HiResLocalSurface we need to check to see if it has a grid yet.
# This is because when creating a new Surface object of this type, the grid generation is deferred until the Scaling and Morphology objects are initialized in order to set the superdomain properly.
if issubclass(self.surface.__class__, HiResLocalSurface) and self.surface.uxgrid is None:
self.surface.set_superdomain(
morphology=self.morphology,
reset=self.is_new,
**surface_config,
)
# The surface should now hav its is_new attribute set and will be set to True if it had to be regridded
self.is_new = self.surface.is_new
self._max_crater_diameter_range = (min(self.surface.face_size), self.target.diameter)
self.smallest_crater = max(self.smallest_crater, self._max_crater_diameter_range[0])
self.largest_crater = min(self.largest_crater, self._max_crater_diameter_range[1])
if self.production._Crater is None:
self.production._Crater = self.Crater
if self.production.quasimc_craters is not None:
# Trigger a reprocessing to get the craters into the updated Crater class
self.quasimc_craters = self.quasimc_craters
if self.is_new:
object.__setattr__(self, "_config_readonly", False)
# The Surface has already had its reset method called.
skip_components = ["surface"]
self.reset(skip_component=skip_components)
self.save_actions = save_actions
self.to_config()
return
def __str__(self) -> str:
"""
Returns a string representation of the Simulation object.
"""
str_repr = "\n<Simulation>\n"
str_repr += f"simdir : {str(self.simdir)}\n\n"
str_repr += "<Current state>\n"
str_repr += f"Interval : {self.interval}\n"
if self.time is not None:
str_repr += (
f"Current time : {format_large_units(self.time, quantity='time')} before present\n"
f"Elapsed time: {format_large_units(self.elapsed_time, quantity='time')}\n"
f"Elapsed N(1) : {self.elapsed_n1} # per 10⁶ km²\n"
)
str_repr += "\n<Components>\n"
str_repr += (
f"Target: {self.target}\n"
f"Surface: {self.surface}\n"
f"Morphology: {self.morphology}\n"
f"Production: {self.production}\n"
f"Projectile: {self.projectile}\n"
f"Scaling: {self.scaling}\n"
f"Counting: {self.counting}\n"
)
return str_repr
def __repr__(self) -> str:
config = self.to_config(save_to_file=False)
txt = f"{self.__class__.__name__}("
for k, v in config.items():
if isinstance(v, str):
v = f"'{v}'"
txt += f"\n {k}={v},"
txt += "\n)"
return txt
def __setattr__(self, name, value):
super().__setattr__(name, value)
if name not in self._user_defined:
return
# Avoid recursive calls during initialization or early access
if (
not self._config_readonly
and hasattr(self, "to_config")
and callable(getattr(self, "to_config", None))
and _convert_for_yaml(value) is not None
):
self.to_config()
def _validate_run_args(
self,
time_start: FloatLike | None = None,
time_end: FloatLike | None = None,
N_D: PairOfFloats | None = None,
N_D_end: PairOfFloats | None = None,
time_interval: FloatLike | None = None,
ninterval: int | None = None,
**kwargs: Any,
) -> dict:
"""
Validate all the input arguments to the sample method. This function will raise a ValueError if any of the arguments are invalid.
Parameters
----------
kwargs : Any
A dictionary of all the arguments passed to the run method, including age, time_start, time_end, diameter_number, time_interval, ninterval, and any additional arguments passed through **kwargs.
Returns
-------
A dict containing all arguments listed in Parameters above, as well as `is_time_interval`, which is a boolean flag indicating
whether or not the simulation is being run in equal age intervals or equal number intervals.
Raises
------
ValueError
If any of the following conditions are met:
- Neither the age nore the diameter_number argument is provided.
- Both the age and diameter_number arguments are provided.
- Both the age and either time_start or time_end arguments are provided.
- Both the time_start and diameter_number arguments are provided.
- The age, time_start, or time_end arguments are provided but are not a scalar.
- The time_interval is provided but is not a positive scalar.
- The time_interval provided is greater than age or time_start - time_end
- The diameter_number argument is not a pair of values, or any of them are less than 0
- The time_interval and nintervaql arguments are both provided.
- The ninterval is provided but is not an integer or is less than 1.
"""
# Determine whether we are going to do equal time intervals or equal number intervals
diameter_number = kwargs.pop("diameter_number", None)
diameter_number_end = kwargs.pop("diameter_number_end", None)
age = kwargs.pop("age", None)
if age is not None:
if time_start is not None or time_end is not None:
raise ValueError("Cannot specify both age and time_start or time_end")
if diameter_number is not None:
raise ValueError("Cannot specify both age and diameter_number")
if N_D is not None:
raise ValueError("Cannot specify both age and N_D")
if not np.isscalar(age):
raise ValueError("age must be a scalar value")
# as age is just a convenience variable, we replace it with time_start = age and time_end = 0
time_start = age
time_end = 0.0
del age
if time_start is not None:
if time_end is None:
time_end = 0.0
if diameter_number is not None:
raise ValueError("Cannot specify both time_start and diameter_number")
if N_D is not None:
raise ValueError("Cannot specify both time_start and N_D")
if not np.isscalar(time_start):
raise ValueError("time_start must be a scalar value")
if not np.isscalar(time_end):
raise ValueError("time_end must be a scalar value")
elif time_end is not None:
if self._time is None:
raise ValueError("time_end cannot be used without time_start")
time_start = self.time
elif N_D is not None:
if diameter_number is not None or diameter_number_end is not None:
raise ValueError("Cannot specify both N_D and diameter_number")
if not isinstance(N_D, (list, tuple, ArrayLike)) or len(N_D) != 2:
raise ValueError("N_D must be a pair of values (diameter, number)")
if N_D_end is None:
N_D_end = (N_D[0], 0.0)
elif not isinstance(N_D_end, (list, tuple, ArrayLike)) or len(N_D_end) != 2:
raise ValueError("N_D_end must be a pair of values (diameter, number)")
diameter_number = (
N_D[0] * self.production.D_conversion_factor,
self.surface.area * N_D[1] / self.production.N_conversion_factor,
)
diameter_number_end = (
N_D_end[0] * self.production.D_conversion_factor,
self.surface.area * N_D_end[1] / self.production.N_conversion_factor,
)
elif N_D_end is not None:
raise ValueError("N_D_end cannot be used without N_D")
elif diameter_number is None:
raise ValueError("Must provide one of age, time_start, or N_D")
if ninterval is not None:
if not isinstance(ninterval, int):
raise TypeError("ninterval must be an integer")
if ninterval < 1:
raise ValueError("ninterval must be greater than zero")
if time_interval is not None:
raise ValueError("Cannot specify both ninterval and time_interval")
elif time_interval is None:
ninterval = 1
is_time_interval = time_interval is not None
# Validate arguments using the production function validator first, which will convert time-based values to diameter_number-based ones
kwargs["area"] = self.surface.area
kwargs["time_start"] = time_start
kwargs["time_end"] = time_end
kwargs["diameter_number"] = diameter_number
kwargs["diameter_number_end"] = diameter_number_end
kwargs["ninterval"] = ninterval
kwargs["time_interval"] = time_interval
kwargs = self.production._validate_sample_args(**kwargs)
if is_time_interval:
if time_interval is None:
if ninterval is None:
ninterval = 1
time_interval = (time_start - time_end) / ninterval
else:
if time_interval > time_start - time_end:
raise ValueError("time_interval must be less than age or time_start - time_end")
elif time_interval <= 0:
raise ValueError("time_interval must be greater than zero")
ninterval = int(np.ceil((time_start - time_end) / time_interval))
kwargs["time_interval"] = time_interval
kwargs["ninterval"] = ninterval
else:
diameter_number = kwargs.get("diameter_number", None)
diameter_number_end = kwargs.get("diameter_number_end", None)
if diameter_number is None:
raise ValueError("Something went wrong! diameter_number should be set by self.production_validate_sample_args")
if diameter_number_end is None:
raise ValueError("Something went wrong! diameter_number_end should be set by self.production_validate_sample_args")
if ninterval is None:
ninterval = 1
diameter_number_interval = (
diameter_number[0],
(diameter_number[1] - diameter_number_end[1]) / ninterval,
)
kwargs["diameter_number_interval"] = diameter_number_interval
kwargs["ninterval"] = ninterval
kwargs["is_time_interval"] = is_time_interval
# Remove unecessary arguments that came out of the production._validate_sample_args method
kwargs.pop("diameter_range")
kwargs.pop("area")
kwargs.pop("compute_time")
return kwargs
[docs]
def run(
self,
age: FloatLike | None = None,
time_start: FloatLike | None = None,
time_end: FloatLike | None = None,
N_D: PairOfFloats | None = None,
N_D_end: PairOfFloats | None = None,
time_interval: FloatLike | None = None,
ninterval: int | None = None,
**kwargs: Any,
) -> None:
"""
Run the simulation over a specified interval using the current production function.
Parameters
----------
age : FloatLike, optional
Start age in My relative to the present for the simulation, used to compute the starting point of the production function.
Default is None, which requires 'time_start' or `N_D` to be set.
time_start : Floatlike, optional
An alternative to `age` that specifies the starting time in My relative to the present for the simulation, used to compute the starting point of the production function. This is used in conjunction with `time_end` in order to allow for simulations that span a range of time rather than being of a specific age. Default is None, which requires either `age` or `N_D` to be set.
time_end : FloatLike, optional
Ending time in My relative to the present for the simulation, used to compute the ending point of the production function.
Default is 0 (present day) if not provided but `time_start` is provided.
N_D : PairOfFloats, optional
A pair of numbers, (D, N), representing the starting cumulative number density N above a given diameter D using the N(D) convention. By default, D must be in units of km and N is in number of craters per 10⁶ km². The N value conversion factor can be set using :py:attr:`~cratermaker.components.production.Production.N_conversion_factor`. For example (20, 1000) is intepreted as N(20) = 1000, which is number of craters larger than 20 km per 10⁶ km². If provided, the function will convert this value to a corresponding age and use the production function for a given age. By default None, which requires `age` or `time_start` to be set.
N_D_end : PairOfFloats, optional
A pair of numbers, (D, N) representing the ending cumulative number density N above a given diameter D using the N(D) convention. By default, D must be in units of km and N is in number of craters per 10⁶ km². The N value conversion factor can be set using :py:attr:`~cratermaker.components.production.Production.N_conversion_factor`. For example (20, 1000) is intepreted as N(20) = 1000, which is number of craters larger than 20 km per 10⁶ km². If provided, the function will convert this value to a corresponding age and use the production function for a given age. By default None, which will set the ending time to 0 (present day) if N_D is provided.
time_interval : FloatLike, optional
Interval in My for outputting intermediate results. If not provided, calculated as `age` / `ninterval` or (`time_start - time_end`) / `ninterval` if `ninterval` is provided, otherwise set to the total simulation duration (e.g. `ninterval=1`).
ninterval : int, optional
Number of intervals for outputting results. This has a special use case where one can specify age-based inputs but output are computed in equal cumulative number intervals and vice versa.
**kwargs : Any
|kwargs|
Notes
-----
This function allows defining the simulation parameters either in terms of time or crater frequency (cumulative number). The
arguments `age`, `time_start`, `time_end` are mutually exclusive with `diameter_number` and `time_interval` is mutually exclusive with `ninterval`.
The initial state of the simulation (before any craters are emplaced) is always saved automatically.
As a result, the total number of saved states will be `ninterval + 1`, where `ninterval` is the number
of simulation intervals requested.
Examples
--------
.. code-block:: python
# Create a simulation object with default parameters (Moon, NeukumProduction, etc.)
sim = cratermaker.Simulation()
# Run the simulation for 3.8 billion years, saving the results every 100 million years
sim.run(age=3.8e3, time_interval=100.0)
# Run the simulation for 3.8 billion years, saving 100 intervals with equal cumulative number intervals
sim.run(age=3.8e3, ninterval=100)
# Run the simulation to create 80 craters larger than 300 km and output 100 equal cumulative number intervals
sim.run(diameter_number=(300e3, 80), ninterval=100)
# Run the simulation from 3.8 billion years to 3.0 billion years, saving the results every 100 million years
sim.run(time_start=3.8e3, time_end=3.0e3, time_interval=100.0)
"""
arguments = {
"age": age,
"time_start": time_start,
"time_end": time_end,
"time_interval": time_interval,
"N_D": N_D,
"N_D_end": N_D_end,
"ninterval": ninterval,
**kwargs,
}
arguments = self._validate_run_args(**arguments)
time_start = arguments.pop("time_start", None)
time_end = arguments.pop("time_end", None)
time_interval = arguments.pop("time_interval", None)
diameter_number = arguments.pop("diameter_number", None)
diameter_number_interval = arguments.pop("diameter_number_interval", None)
ninterval = arguments.pop("ninterval", None)
is_time_interval = arguments.pop("is_time_interval", None)
validate_inputs = kwargs.pop("validate_inputs", False)
if ninterval == 1:
is_time_interval = True
time_interval = time_start - time_end
if not is_time_interval:
diameter_number_density_start = (
diameter_number[0],
diameter_number[1] / self.surface.area,
)
time_start = self.production.age_from_D_N(*diameter_number_density_start, validate_inputs=validate_inputs)
# If this is a fresh run, we need to set the value of the current time based on the requested starting condition.
if self._time is None:
self.time = time_start
self.interval = 0
# If this is a restarted run, we need to distinguish between a true restart and a continuation with different parameters
if time_start < self.time:
raise RuntimeError(
"Starting time cannot be later than the current time. Choose a starting time value equal to or larger than the current time, or reset this simulation."
)
if is_time_interval:
initial_interval = int((time_start - self.time) / time_interval)
else:
delta_n1_start = self.production.function(
diameter=_DSTD,
time_start=time_start,
time_end=self.time,
validate_inputs=validate_inputs,
).item()
n1_interval = (
self.production.function(
diameter=_DSTD,
time_start=time_start,
time_end=time_end,
validate_inputs=validate_inputs,
).item()
/ ninterval
)
initial_interval = int(delta_n1_start / n1_interval)
if self.is_new:
self.save(**kwargs)
with tqdm(
total=ninterval,
initial=initial_interval,
unit="interval",
position=3,
leave=True,
) as pbar:
for i in range(initial_interval, ninterval):
self.counting._emplaced = []
if is_time_interval:
time = time_start - i * time_interval
current_time_end = time_start - (i + 1) * time_interval
if current_time_end < time_end:
current_time_end = time_end
time_str = format_large_units(time, quantity="time")
elapsed_str = format_large_units(self.elapsed_time, quantity="time")
pbar.set_description(f"Current: {time_str} bp; Elapsed: {elapsed_str}")
self.populate(time_start=time, time_end=current_time_end, **kwargs)
else:
current_diameter_number = (
diameter_number[0],
diameter_number[1] - i * diameter_number_interval[1],
)
current_diameter_number_end = (
diameter_number[0],
diameter_number[1] - (i + 1) * diameter_number_interval[1],
)
self.populate(
diameter_number=current_diameter_number,
diameter_number_end=current_diameter_number_end,
**kwargs,
)
current_diameter_number_density = (
current_diameter_number[0],
current_diameter_number[1] / self.surface.area,
)
time = self.production.age_from_D_N(*current_diameter_number_density, validate_inputs=validate_inputs)
if current_diameter_number_end[1] > 0:
current_diameter_number_density_end = (
current_diameter_number_end[0],
current_diameter_number_end[1] / self.surface.area,
)
current_time_end = self.production.age_from_D_N(
*current_diameter_number_density_end,
validate_inputs=validate_inputs,
)
else:
current_time_end = 0.0
if current_time_end < 0.0:
current_time_end = 0.0
time_interval = time - current_time_end
self.elapsed_time += time_interval
self.elapsed_n1 += (
self.production.function(
diameter=_DSTD,
time_start=time,
time_end=current_time_end,
validate_inputs=validate_inputs,
).item()
* self.production.N_conversion_factor
)
self.time = current_time_end
self.interval += 1
self.save(**kwargs)
pbar.update(1)
return
[docs]
def populate(
self,
age: FloatLike | None = None,
time_start: FloatLike | None = None,
time_end: FloatLike | None = None,
N_D: PairOfFloats | None = None,
N_D_end: PairOfFloats | None = None,
diameter_number: PairOfFloats | None = None,
diameter_number_end: PairOfFloats | None = None,
**kwargs: Any,
) -> list[Crater]:
"""
Populate the surface with craters over a specified interval using the current production function.
Parameters
----------
age : FloatLike, optional
Age in the past in units of My relative to the present, which is used compute the cumulative SFD.
time_start : FloatLike, optional
An alternative to `age` that specifies the starting time in My relative to the present for the simulation, used to compute the starting point of the production function. This is used in conjunction with `time_end` in order to allow for simulations that span a range of time rather than being of a specific age. Default is None, which requires either `age` or `diameter_number`
time_end : FloatLike, optional
The ending time in My relative to the present for the simulation, used to compute the ending point of the production function. Default is 0 (present day).
N_D : PairOfFloats, optional
A pair of numbers, (D, N), representing the starting cumulative number density N above a given diameter D using the N(D) convention. By default, D is in units of km and N is in number of craters per 10⁶ km². The N value conversion factor can be set using :py:attr:`~cratermaker.components.production.Production.N_conversion_factor`. For example (20, 1000) is intepreted as N(20) = 1000, which is number of craters larger than 20 km per 10⁶ km². If provided, the function will convert this value to a corresponding age and use the production function for a given age.
N_D_end : PairOfFloats, optional
A pair of numbers, (D, N), representing the ending cumulative number density N above a given diameter D using the N(D) convention. By default, D is in units of km and N is in number of craters per 10⁶ km². The N value conversion factor can be set using :py:attr:`~cratermaker.components.production.Production.N_conversion_factor`. For example (20, 1000) is intepreted as N(20) = 1000, which is number of craters larger than 20 km per 10⁶ km². If provided, the function will convert this value to a corresponding age and use the production function for a given age.
diameter_number : PairOfFloats, optional
A pair of numbers, (diameter, number), representing the diameter and total number of craters larger than that diameter, where diameter is in units of m and n is in number of craters. If provided, the function will convert this value to a corresponding age and use the production function for a given age.
diameter_number_end : PairOfFloats, optional
A pair of numbers, (diameter, number), representing the diameter and total number of craters in the production function at the end, where diameter is in units of m and n is in number of craters. If provided, the function will convert this value to a corresponding age and use the production function for a given age.
craters : list[Crater] or Crater, optional
A list of Crater objects to include along with the randomly generated craters. The crater list must include either time or time_min, time_max values.
Returns
-------
list[Crater]
A list of the Crater objects that were emplaced. Returns an empty list if no craters were emplaced.
"""
if not hasattr(self, "production"):
raise RuntimeError("No production function defined for this simulation")
elif not hasattr(self.production, "generator_type"):
raise RuntimeError("The production function is not properly defined. Missing 'generator_type' attribute")
elif self.production.generator_type not in ["crater", "projectile"]:
raise RuntimeError(f"Invalid production function type {self.production.generator_type}")
arguments = {
"age": age,
"time_start": time_start,
"time_end": time_end,
"N_D": N_D,
"N_D_end": N_D_end,
"diameter_number": diameter_number,
"diameter_number_end": diameter_number_end,
**kwargs,
}
# This will convert any time-based arguments into diameter_number
arguments = self._validate_run_args(**arguments)
del age
time_start = arguments.pop("time_start", None)
time_end = arguments.pop("time_end", None)
diameter_number = arguments.pop("diameter_number", None)
diameter_number_end = arguments.pop("diameter_number_end", None)
from_projectile = self.production.generator_type == "projectile"
diam_key = "projectile_diameter" if from_projectile else "diameter"
if from_projectile:
diam_max = self.largest_projectile
diam_min = self.smallest_projectile
else:
diam_max = self.largest_crater
diam_min = self.smallest_crater
# Loop over each face in the mesh to build up a population of craters in this interval. This is done because faces may
# not all have the same surface area, the range of crater sizes that can be formed on each face may be different.
impact_diameters = []
impact_times = []
impact_locations = []
# Process each bin
for i, face_indices in enumerate(self.surface.face_bin_indices):
total_bin_area = self.surface.face_bin_area[i]
area_ratio = total_bin_area / self.surface.area
smallest_crater = self.surface.face_bin_min_sizes[i]
if from_projectile:
smallest_projectile = self._get_diameter_limit("smallest", crater_diameter=smallest_crater)
diam_min = max(smallest_projectile, diam_min)
else:
diam_min = max(smallest_crater, diam_min)
diameter_number_local = (diameter_number[0], diameter_number[1] * area_ratio) if diameter_number is not None else None
if diameter_number_end is not None:
diameter_number_end_local = (
diameter_number_end[0],
diameter_number_end[1] * area_ratio,
)
else:
diameter_number_end_local = None
diameters, times = self.production.sample(
diameter_number=diameter_number_local,
diameter_number_end=diameter_number_end_local,
diameter_range=(diam_min, diam_max),
area=total_bin_area,
**kwargs,
)
if diameters.size > 0:
impact_diameters.extend(diameters.tolist())
impact_times.extend(times.tolist())
# Get the relative probability of impact onto any particular face then get the locations of the impacts
p = self.surface.face_area[face_indices] / total_bin_area
face_indices = self.rng.choice(face_indices, size=diameters.shape, p=p)
locations = self.surface.get_random_location_on_face(face_indices)
impact_locations.extend(np.array(locations).T.tolist())
# If quasi-Monte Carlo is enabled, we need to extract the quasimc craters that
craterlist = []
if len(impact_diameters) > 0:
# Sort the times, diameters, and locations so that they are in order of decreasing age
impact_diameters = np.asarray(impact_diameters)
impact_times = np.asarray(impact_times)
impact_locations = np.array(impact_locations)
for diameter, location, time in tqdm(
zip(impact_diameters, impact_locations, impact_times, strict=False),
total=len(impact_diameters),
desc="Generating crater population",
unit="crater",
position=0,
leave=False,
):
diam_arg = {diam_key: diameter}
craterlist.append(
self.Crater.maker(
location=location,
time=time,
**diam_arg,
**vars(self.common_args),
**kwargs,
)
)
craterlist.sort(key=lambda c: c.time, reverse=True)
if self.do_quasimc:
if time_start is None:
N_D = (
diameter_number[0] / self.production.D_conversion_factor,
diameter_number[1] * self.production.N_conversion_factor / self.surface.area,
)
else:
N_D = None
if time_end is None:
N_D_end = (
diameter_number_end[0] / self.production.D_conversion_factor,
diameter_number_end[1] * self.production.N_conversion_factor / self.surface.area,
)
else:
N_D_end = None
craterlist = self.production.quasimc_merge(
craterlist, time_start=time_start, time_end=time_end, N_D=N_D, N_D_end=N_D_end, **kwargs
)
if len(craterlist) > 0:
return self.emplace(craterlist, **kwargs)
else:
return []
[docs]
def emplace(self, craters: list[Crater] | Crater | None = None, **kwargs: Any) -> list[Crater]:
"""
Emplace one or more craters in the simulation.
This method orchestrates the creation and placement of a crater in the
simulation. It can create a crater directly or based on the characteristics
of a projectile.
Parameters
----------
craters : Crater or list of Crater objects, optional
The Crater object(s) to be emplaced. If provided, this will be used directly. Otherwise, a single crat er will be generated based on the keyword arguments.
**kwargs : Any
|kwargs|
Returns
-------
list[Crater]
A list of the Crater objects that were emplaced in the simulation. Returns an empty list if no craters were emplaced.
Notes
-----
The keyword arguments provided are passed down to :py:meth:`Crater.maker`. Refer to its documentation for a detailed description of valid keyword arguments.
Examples
--------
.. code-block:: python
from cratermaker import Simulation, Crater
sim = Simulation()
# Create a crater with specific diameter
sim.emplace(diameter=10.0e3)
# Create a crater based on a projectile with given mass and projectile_velocity
sim.emplace(projectile_mass=1e15, projectile_velocity=20e3)
# Create a crater with a specific transient diameter and location
sim.emplace(transient_diameter=50e3, location=(43.43, -86.92))
# Create multiple craters
craters = [Crater.maker(diameter=20.0e3), Crater.maker(diameter=20.0e3)]
sim.emplace(craters)
"""
self.is_new = False
return self.morphology.emplace(craters=craters, **kwargs)
[docs]
def save(self, **kwargs: Any) -> None:
"""
Save the current simulation state to a file.
Parameters
----------
**kwargs : Any
Additional keyword argumments to pass to the component save methods.
"""
save_args = {"interval": self.interval, "time_variables": self.time_variables, **kwargs}
self.surface.save(**save_args)
self.counting.save(**save_args)
self.to_config(**kwargs)
super().save(**save_args)
return
[docs]
def export(
self,
driver: str = "OpenCraterTool",
interval: int = -1,
ask_overwrite: bool | None = None,
**kwargs: Any,
) -> None:
"""
Export component output to a specified file format.
Parameters
----------
driver : str, optional
The driver to use export the data to. Supported formats are 'OpenCraterTool', 'VTK' or a driver supported by GeoPandas ('GPKG', 'ESRI Shapefile', etc.). This is overridden if either the filename or file_extension parameters are provided. Default is 'OpenCraterTool'.
interval : int, optional
The interval number to export. Default is -1 (the most current interval saved in the simulation).
ask_overwrite : bool, optional
|ask_overwrite_methods|
**kwargs : Any
|kwargs|
Notes
-----
The default driver is 'OpenCraterTool', which is designed to output data into a format that is relatively easy to import into QGIS with the OpenCraterTool plugin. This will create a GeoTIFF file representation of the surface, and a set of SCC files for the crater counting data if counting is enabled.
"""
# Temporarily set the ask_overwrite attribute for the duration of the export, but reset it to its original value afterwards.
ask_overwrite_orig = self.ask_overwrite
if ask_overwrite is not None:
self.ask_overwrite = ask_overwrite
if interval is not None and interval < 0:
interval = self.interval + 1 + interval
self.save(skip_actions=True, ask_overwrite=ask_overwrite, **kwargs)
if driver.lower() == "opencratertool":
surface_driver = "GeoTIFF"
counting_driver = "SCC"
else:
surface_driver = driver
counting_driver = driver
self.surface.export(
driver=surface_driver,
interval=interval,
ask_overwrite=ask_overwrite,
**kwargs,
)
self.counting.export(
driver=counting_driver,
interval=interval,
ask_overwrite=ask_overwrite,
**kwargs,
)
self.ask_overwrite = ask_overwrite_orig
return
[docs]
def labelmaker(
self,
interval: int | None = None,
interval_label: bool = False,
time_label: bool = False,
age_label: bool = True,
N_label: bool = False,
N_diam_val: float | None = None,
compact: bool = False,
**kwargs: Any,
) -> str:
"""
Generates a label for the current state of the simulation based on the time variables and other parameters.
Parameters
----------
interval : int, optional
The interval number to use for generating the label. Default is None, which will use the most current interval saved in the simulation.
interval_label : bool, optional
If True, the interval number will be included in the label. Default is False.
time_label : bool, optional
If True, the time variable will be included in the label if it is available. Default is True.
age_label : bool, optional
If True, the elapsed time variable will be included in the label if it is available. Default is False.
N_label : bool, optional
If True, the elapsed N(D) variable will be included in the label if it is available. Default is False.
N_diam_val : float, optional
The D value in km to use for the N(D) label. If None, a default value is chosen based on the smallest crater size in the simulation. Default is None.
compact : bool, optional
If True, the label will be formatted in a more compact way with newlines separating the different components. If False, the components will be separated by three spaces. Default is False.
**kwargs : Any
|kwargs|
"""
if interval is None:
interval = self.interval
time_variables = self.time_variables
else:
tda = self.surface.get_time_variables(interval=interval).compute()
time_variables = {}
for c in tda.coords:
time_variables[c] = tda[c].data.item()
if compact:
delim = "\n"
else:
delim = " "
# Set a default value of the D for N(D) format based on the smallest pixel in the simulation
if N_diam_val is None:
if self.smallest_crater < 1e2:
N_diam_val = 1.0
elif self.smallest_crater < 10e3:
N_diam_val = 20.0
else:
N_diam_val = 64.0
time = time_variables.get("time")
time_label = time_label and time is not None
elapsed_time = time_variables.get("elapsed_time")
age_label = age_label and elapsed_time is not None
elapsed_n1 = time_variables.get("elapsed_n1")
N_label = N_label and elapsed_n1 is not None
label = []
if interval_label:
label.append(f"Interval: {interval}")
if time_label:
label.append(f"Time: {time:.0f} My bp")
if age_label:
label.append(f"Age : {elapsed_time:.0f} My")
if N_label:
nfac = self.production.csfd(diameter=N_diam_val * self.production.D_conversion_factor) / self.production.csfd(
diameter=_DSTD
)
dtxt = f"{N_diam_val:.0f}"
ndval = elapsed_n1 * nfac
N_unit = self.production.N_unit.replace("⁶", "$^6$")
label.append(f"Elapsed N({dtxt}): {ndval:.4g} {N_unit}")
return delim.join(label)
[docs]
def plot(
self,
interval: int | None = None,
plot_style: str = "hillshade",
label="default",
show=False,
save=True,
ax: Axes | None = None,
close_when_done: bool = True,
scalebar: bool | None = None,
colorbar: bool = True,
include_counting: bool = False,
observed_color: str | None = "white",
observed_original_color: str | None = None,
emplaced_color: str | None = None,
minimum_plot_width: float | None = 800.0,
superdomain: bool = False,
**kwargs: Any,
) -> Axes:
"""
Plot the current state of the surface.
Parameters
----------
include_counting : bool, optional
If True, the counting data will be included in the plot if counting is enabled. Default is False
interval : int, optional
The interval number to plot. Default is None, which will plot the most current interval saved in the simulation.
plot_style : str, optional
The style to use for surface plots. See :py:meth:`Surface.plot` for more details. Default is 'hillshade'.
cmap : str, optional
The colormap to use for the plot. If None, a default colormap will be used ("cividis" by default and "grey" when plot_style=="hillshade" and variable=="face_elevation").
observed_color : str | None, optional
The color to use for observed craters using their measured properties. If None, observed craters will not be plotted. Default is "white".
observed_original_color : str | None, optional
The color to use for observed craters using their original properties. If None, observed craters will not be plotted. Default is None.
emplaced_color : str | None, optional
The color to use for emplaced craters in the interval. If None, emplaced craters will not be plotted. Default is None.
label : str, optional
The label to use for the plot. Default is "default", which will use the labelmaker method to build a label based on the current time and age of the simulation.
scalebar : bool, optional
If True, a scalebar will be added to the plot. Default is False for global surfaces.
colorbar : bool, optional
If True, a colorbar will be added to the plot when using "map" plot_style or "hillshade" with a variable overlay. Default is True.
show : bool, optional
If True, the plot will be displayed. Default for local surfaces is True.
save : bool, optional
If True, the plot will be saved to the default plot directory. Default is False.
ax : matplotlib.axes.Axes, optional
An existing Axes object to plot on. If None, a new figure and axes will be created.
close_when_done : bool, optional
If True, the figure will be closed after plotting. Default is True when save is True and show is False, and False otherwise.
minimum_plot_width : float, optional
Because the width of the plot is determined by the number of faces, small regions will generate small plots with labels that are hard to read. This parameter sets a lower limit to the width of the image that is generated by the plot. By default it is 800. Set to None to turn it off.
superdomain : bool, optional
|HiResLocalSurface| variants only. If True, the full surface including the superdomain will be plotted. If False, only the local region will be plotted. Default is False
**kwargs : Any
|kwargs|
Returns
-------
Axes
The matplotlib Axes object created by the surface plot method.
"""
if label == "default":
compact = kwargs.pop("compact", issubclass(self.surface.__class__, HiResLocalSurface))
label = self.labelmaker(interval=interval, compact=compact, **kwargs)
if interval is None:
interval = self.interval
self.save(**kwargs, skip_actions=True)
plot_args = {
"interval": interval,
"plot_style": plot_style,
"label": label,
"show": show,
"save": save,
"scalebar": scalebar,
"colorbar": colorbar,
"observed_color": observed_color,
"observed_original_color": observed_original_color,
"emplaced_color": emplaced_color,
"ax": ax,
"close_when_done": close_when_done,
"minimum_plot_width": minimum_plot_width,
"superdomain": superdomain,
**kwargs,
}
if include_counting:
ax = self.counting.plot(**plot_args)
else:
ax = self.surface.plot(**plot_args)
return ax
[docs]
def show3d(
self,
**kwargs: Any,
):
"""
Show the component using an interactive 3D plot.
Valid arguments are those that are passed the engine functions (e.g. arguments to pyvista_plotter(), and plotter.show())
Parameters
----------
engine : str, optional
The engine to use for plotting. Currently, only "pyvista" is supported. Default is "pyvista".
interval : int, optional
The interval number to plot. Default is None, which will plot the most current interval saved in the simulation.
label : str, optional
The label to use for the plot. Default is "default", which will use a label based on the current time and elapsed time of the simulation. When "default" is used, additional options for |sim.labelmaker| can be passed as kwargs.
plotter : pv.Plotter, optional
An existing PyVista Plotter object to use for the plot. If None, a new Plotter object will be created. Default is None.
enable_interactive : bool, optional
If True, the key events for the plotter will be updated to include custom events for navigating between intervals. Default is True.
variable_name : str | None, optional
The name of the variable to plot. If the name of the variable is not already stored on the surface mesh, then the `variable` argument must also be passed. Default is None, which will plot a greyscale image of the surface.
variable : (n_face) array, optional
An array face values that will be used to color the surface mesh. This is required if `variable_name` is not a face variable that is already saved in the the uxds dataset. Default is None.
interval : int | None, optional
The interval number of the surface to plot. If None, the currently loaded surface data will be used. Default is None.
theme : str, optional
The PyVista plot theme to use. If None, the default PyVista theme will be used. Default is None.
transparent_background : bool, optional
If True, the background of the plot will be transparent. Default is None, which will use the default background setting for the chosen plot theme.
superdomain : bool, optional
If True, show the full surface including the superdomain (|HiResLocalSurface| variants only). If False, show only the local region. Default is False.
crater_type: Literal["observed","emplaced","both"] | list[Crater].
This is only used if enable_interactive is False, in which case it controls which crater data will be plotted. Default is "both".
crater_color: str | tuple[str] = ("white", "yellow")
The color to use for the plots. When used with enable_interactive or when craters is set to "both", this is a tuple where the first element is the color of observed craters and the second is the color of emplaced craters. Default is ("white,"red").
crater_style : Literal["rings", "points"], optional
Only used when enable_interactive is False. Sets the style of the mesh. Options are "rings", which creates polyline circles over the rim of each crater or "points" which creates a point at the center.
crater_size_scale_factor : FloatLike, optional
A factor to scale the size of the craters in "point", "impacts", or "spheres" styles. Default is 1.0.
**kwargs : Any
|kwargs|
"""
return super().show3d(**kwargs)
[docs]
def pyvista_plotter(
self,
interval: int | None = None,
label: str = "default",
plotter: pyvista.Plotter | None = None,
enable_interactive: bool = True,
**kwargs: Any,
) -> pyvista.Plotter:
"""
Create a PyVista plotter for the current state of the surface.
Parameters
----------
interval : int, optional
The interval number to plot. Default is None, which will plot the most current interval saved in the simulation.
label : str, optional
The label to use for the plot. Default is "default", which will use a label based on the current time and elapsed time of the simulation. When "default" is used, additional options for |sim.labelmaker| can be passed as kwargs.
plotter : pv.Plotter, optional
An existing PyVista Plotter object to use for the plot. If None, a new Plotter object will be created. Default is None.
enable_interactive : bool, optional
If True, the key events for the plotter will be updated to include custom events for navigating between intervals. Default is True.
variable_name : str | None, optional
The name of the variable to plot. If the name of the variable is not already stored on the surface mesh, then the `variable` argument must also be passed. Default is None, which will plot a greyscale image of the surface.
variable : (n_face) array, optional
An array face values that will be used to color the surface mesh. This is required if `variable_name` is not a face variable that is already saved in the the uxds dataset. Default is None.
interval : int | None, optional
The interval number of the surface to plot. If None, the currently loaded surface data will be used. Default is None.
theme : str, optional
The PyVista plot theme to use. If None, the default PyVista theme will be used. Default is None.
transparent_background : bool, optional
If True, the background of the plot will be transparent. Default is None, which will use the default background setting for the chosen plot theme.
superdomain : bool, optional
If True, show the full surface including the superdomain (|HiResLocalSurface| variants only). If False, show only the local region. Default is False.
crater_type: Literal["observed","emplaced","both"] | list[Crater].
This is only used if enable_interactive is False, in which case it controls which crater data will be plotted. Default is "both".
crater_color: str | tuple[str] = ("white", "yellow")
The color to use for the plots. When used with enable_interactive or when craters is set to "both", this is a tuple where the first element is the color of observed craters and the second is the color of emplaced craters. Default is ("white,"red").
crater_style : Literal["rings", "points"], optional
Only used when enable_interactive is False. Sets the style of the mesh. Options are "rings", which creates polyline circles over the rim of each crater or "points" which creates a point at the center.
crater_size_scale_factor : FloatLike, optional
A factor to scale the size of the craters in "point", "impacts", or "spheres" styles. Default is 1.0.
Returns
-------
pyvista.Plotter
A PyVista Plotter object created by the surface pyvista_plotter method.
"""
from cratermaker.utils.general_utils import toggle_pyvista_actor, update_pyvista_help_message
class IntervalState:
def __init__(self, interval, max_interval):
self.interval = interval
self.max_interval = max_interval
def update(self, forward):
if forward:
new_interval = self.interval + 1
if new_interval > self.max_interval:
return
else:
new_interval = self.interval - 1
if new_interval < 0:
return
self.interval = new_interval
return new_interval
def update_interval(plotter, istate, forward):
interval = istate.update(forward)
plotter.clear_actors()
plotter = self.pyvista_plotter(interval=interval, plotter=plotter, enable_interactive=enable_interactive, **kwargs)
plotter.update()
return
text_arg_keys = ["position", "font_size", "color", "font", "shadow", "viewport", "orientation", "font_file", "render"]
new_plotter = plotter is None
if interval is None:
self.save(**kwargs, skip_actions=True)
interval = self.interval
if self.counting is not None:
plotter = self.counting.pyvista_plotter(
interval=interval, plotter=plotter, enable_interactive=enable_interactive, **kwargs
)
else:
plotter = self.surface.pyvista_plotter(
interval=interval, plotter=plotter, enable_interactive=enable_interactive, **kwargs
)
if label is not None:
if label == "default":
compact = kwargs.pop("compact", True)
label = self.labelmaker(interval=interval, compact=compact, **kwargs)
text_args = {k: kwargs[k] for k in text_arg_keys if k in kwargs}
label_actor = plotter.add_text(text=label, name="simulation-label", **text_args)
plotter.remove_actor("simulation-label")
plotter.add_actor(label_actor)
if new_plotter:
istate = IntervalState(interval, max_interval=self.interval)
if enable_interactive and new_plotter:
if label is not None:
plotter = update_pyvista_help_message(plotter, new_message="l: Show/hide label")
plotter.add_key_event(
"l", lambda plotter=plotter, label_actor=label_actor: toggle_pyvista_actor(plotter, label_actor)
)
plotter = update_pyvista_help_message(plotter, new_message="Right: Next interval")
plotter = update_pyvista_help_message(plotter, new_message="Left: Previous interval")
plotter.add_key_event("Right", lambda plotter=plotter: update_interval(plotter, istate, forward=True))
plotter.add_key_event("Left", lambda plotter=plotter: update_interval(plotter, istate, forward=False))
return plotter
[docs]
def to_config(self, save_to_file: bool = True, **kwargs: Any) -> dict:
"""
Converts values to types that can be used in yaml.safe_dump.
This will convert various types into a format that can be saved in a human-readable YAML file. This will consolidate all of the configuration
parameters into a single dictionary that can be saved to a YAML file. This will also remove any common arguments from the individual configurations for each component model to avoid repeating them.
Parameters
----------
save_to_file : bool, optional
If True, the configuration will be saved to a file. Default is True.
**kwargs : Any
|kwargs|
Returns
-------
dict[str, Any]
A dictionary of the object's attributes that can be serialized to YAML.
Notes
-----
- The function will ignore any attributes that are not serializable to human-readable YAML. Therefore, it will ignore anything that cannot be converted into a str, int, float, or bool.
- The function will convert Numpy types to their native Python types.
"""
sim_config = super().to_config(remove_common_args=False)
for component_name in _COMPONENT_NAMES:
component_config = component_name + "_config"
component = getattr(self, component_name, None)
if component is not None:
sim_config[component_name] = component.component_name if hasattr(component, "component_name") else component
if hasattr(component, "to_config") and callable(getattr(component, "to_config", None)):
sim_config[component_config] = component.to_config(remove_common_args=True)
# drop any empty values or {} from either f"{config} or f"{config}_config" if when they are either None or empty
if component_name in sim_config and (sim_config[component_name] is None or sim_config[component_name] == {}):
sim_config.pop(component_name)
if component_config in sim_config and (sim_config[component_config] is None or sim_config[component_config] == {}):
sim_config.pop(component_config)
# Write the combined configuration to a YAML file
if save_to_file and not self._config_readonly:
with Path.open(self.config_file, "w") as f:
yaml.safe_dump(sim_config, f, indent=4)
return sim_config
[docs]
def reset(
self,
skip_component: str | list[str] | None = None,
) -> None:
"""
Reset the simulation by clearing all data and files associated with it.
Parameters
----------
skip_component : str or list of str, optional
List of component names to skip during the reset process. Default is an empty list, which means all components will be reset.
"""
if skip_component is None:
skip_component = []
elif isinstance(skip_component, str):
skip_component = [skip_component]
elif not isinstance(skip_component, list) or not all(isinstance(c, str) for c in skip_component):
raise TypeError("skip_component must be a string or a list of strings")
ask_overwrite_original = (
self.ask_overwrite
) # Saves the value of ask_overwrite in case the user selects 'a' for this operation
if self.ask_overwrite:
files_to_remove = []
for component in _COMPONENT_NAMES:
if component not in skip_component and hasattr(self, component):
files_to_remove += getattr(self, component).saved_output_files()
if len(files_to_remove) > 0:
print("The following files will be deleted:")
for f in files_to_remove:
print(f" {f}")
print("To disable this message, set `ask_overwrite=False` to this instance.")
response = input(f"Are you sure you want to delete {len(files_to_remove)} files? [y/N/a]: ")
if response.lower() == "a":
self.ask_overwrite = False
elif response.lower() != "y":
raise RuntimeError("User aborted the reset operation.")
for component in _COMPONENT_NAMES:
if component not in skip_component and hasattr(self, component):
getattr(self, component).reset()
self.ask_overwrite = ask_overwrite_original
self._interval = 0
self._elapsed_time = None
self._time = None
self._elapsed_n1 = None
self.save(skip_actions=True)
self.is_new = True
return
[docs]
def update_elevation(self, *args: Any, **kwargs: Any) -> None:
"""
Set the elevation on the surface. Delegates to the Surface object.
Parameters
----------
*args: Variable length argument list to pass to self.surface.update_elevation.
|kwargs|
"""
return self.surface.update_elevation(*args, **kwargs)
@property
def max_crater_diameter_range(self) -> tuple[float, float]:
if self._max_crater_diameter_range is None:
return (0.0, np.inf) # Temporarily set until after the Surface is finished initializing
else:
return self._max_crater_diameter_range
def _get_diameter_limit(
self,
limit: Literal["smallest", "largest"],
crater_diameter: FloatLike | None = None,
projectile_diameter: FloatLike | None = None,
) -> float:
"""
Return an estimated smallest or largest crater given a projectile or projectile given a crater.
Parameters
----------
limit : Literal["smallest", "largest"]
A string indicating whether the scaling should be done to calculate the largest crater (vertical impact at 10x the mean velocity of the Projectile component) or the smallest crater (1 degree impact at 0.1x the mean velocity of the Projectile component).
crater_diameter : FloatLike, optional
If passed, this will be used as the crater diameter and the projectile diameter will be computed from it. If None, then the crater will be scaled from the projectile_diameter.
projectile_diameter : FloatLike, optional
If passed, this will be used as the projectile diameter and the crater diameter will be computed from it. If None, then the projectile will be scaled from the crater_diameter.
Returns
-------
float
The smallest crater diameter given a projectile diameter or the smallest projectile diameter given a crater diameter
"""
scaling = Scaling.maker(self.scaling)
if limit == "smallest":
angle = 1.0
projectile_velocity = scaling.projectile.mean_velocity * 0.1
elif limit == "largest":
angle = 90.0
projectile_velocity = scaling.projectile.mean_velocity * 10.0
else:
raise ValueError("limit must either be 'smallest' or 'largest'")
if (crater_diameter is None and projectile_diameter is None) or (
crater_diameter is not None and projectile_diameter is not None
):
raise ValueError("Must provide exactly one of either crater_diameter or projectile_diameter")
if crater_diameter is None:
crater_diameter = self.Crater.maker(
projectile_diameter=projectile_diameter,
angle=angle,
projectile_velocity=projectile_velocity,
**vars(self.common_args),
).diameter
if limit == "smallest":
return float(max(crater_diameter, self.max_crater_diameter_range[0]))
elif limit == "largest":
return float(min(crater_diameter, self.max_crater_diameter_range[1]))
if projectile_diameter is None:
if limit == "smallest":
crater_diameter = max(crater_diameter, self.max_crater_diameter_range[0])
elif limit == "largest":
crater_diameter = min(crater_diameter, self.max_crater_diameter_range[1])
projectile_diameter = self.Crater.maker(
diameter=crater_diameter,
angle=angle,
projectile_velocity=projectile_velocity,
**vars(self.common_args),
).projectile_diameter
return float(projectile_diameter)
@property
def target(self) -> Target | str:
"""
The target body for the impact simulation. Set during initialization.
"""
return self._target
@target.setter
def target(self, value: Target | str):
if not isinstance(value, (Target | str)):
raise TypeError("target must be an instance of Target or str")
self._target = value
@property
def surface(self) -> Surface:
"""
Surface mesh data for the simulation. Set during initialization.
"""
return self._surface
@surface.setter
def surface(self, value):
if not isinstance(value, (Surface | str | type)) or (isinstance(value, type) and not issubclass(value, Surface)):
raise TypeError("surface must be an instance of Surface, a subclass of Surface, or str")
self._surface = value
@property
def production(self):
"""
The Production class instance used for crater production. Set during initialization.
"""
return self._production
@production.setter
def production(self, value):
if not isinstance(value, (Production | str)):
raise TypeError("production must be a subclass of Production or str")
self._production = value
@property
def scaling(self):
"""
The Scaling object that defines the crater scaling relationships model. Set during initialization.
"""
return self._scaling
@scaling.setter
def scaling(self, value):
if not isinstance(value, (Scaling | str)):
raise TypeError("scaling must be of Scaling type or str")
self._scaling = value
@property
def morphology(self):
"""
The crater morphology model. Set during initialization.
"""
return self._morphology
@morphology.setter
def morphology(self, value):
if not isinstance(value, (Morphology | str)):
raise TypeError("morphology must be of Morphology type or str")
self._morphology = value
@property
def projectile(self):
"""
The crater projectile model. Set during initialization.
"""
return self._projectile
@projectile.setter
def projectile(self, value):
if not isinstance(value, (Projectile | str)):
raise TypeError("projectile must be of Projectile type or str")
self._projectile = value
@property
def counting(self):
"""
The crater counting model. Set during initialization.
"""
return self._counting
@counting.setter
def counting(self, value):
if not isinstance(value, (Counting | str)):
raise TypeError("counting must be of Counting type or str")
self._counting = value
@property
def n_node(self):
"""
Number of nodes in the simulation mesh. Dynamically set based on `surface` attribute.
"""
return self.surface.uxgrid.n_node
@property
def n_face(self):
"""
Number of faces in the simulation mesh. Dynamically set based on `surface` attribute.
"""
return self.surface.uxgrid.n_face
@parameter
def interval(self):
"""
The index of the current time step.
"""
if self._interval is None:
return 0
return self._interval
@interval.setter
def interval(self, value):
if not isinstance(value, int):
raise TypeError("interval must be an integer")
if value < 0:
raise ValueError("interval must be greater than or equal to zero")
self._interval = value
@parameter
def elapsed_time(self):
"""
The elapsed time in My since the start of the simulation.
"""
if self._elapsed_time is None:
return 0.0
return self._elapsed_time
@elapsed_time.setter
def elapsed_time(self, value):
self._elapsed_time = float(value)
@parameter
def time(self):
"""
The age of the current time step in My relative to the present from the chronology of the production function.
"""
return self._time
@time.setter
def time(self, value):
self._time = float(value)
@parameter
def elapsed_n1(self):
"""
The elapsed number of craters larger than 1 km in diameter.
"""
if self._elapsed_n1 is None:
return 0.0
return self._elapsed_n1
@elapsed_n1.setter
def elapsed_n1(self, value):
self._elapsed_n1 = float(value)
@parameter
def smallest_crater(self):
"""
The smallest crater diameter in meters. Set during initialization.
"""
if not isinstance(self.production, Production):
return self.max_crater_diameter_range[0]
if self.production.generator_type == "crater":
return self.production.diameter_range[0]
elif self.production.generator_type == "projectile":
projectile_diameter = self.production.diameter_range[0]
if projectile_diameter > 0:
return self._get_diameter_limit("smallest", projectile_diameter=projectile_diameter)
else:
return projectile_diameter
@smallest_crater.setter
def smallest_crater(self, value):
if not isinstance(self.production, Production):
return
if self.production.generator_type == "crater":
minval = max(value, self.max_crater_diameter_range[0])
maxval = min(self.production.diameter_range[1], self.max_crater_diameter_range[1])
elif self.production.generator_type == "projectile":
minval = self._get_diameter_limit("smallest", crater_diameter=value)
maxval = min(
self.production.diameter_range[1],
self._get_diameter_limit("smallest", crater_diameter=self.max_crater_diameter_range[1]),
)
else:
raise ValueError("Unknown Production generator_type. Expecting either 'crater' or 'projectile'")
self.production.diameter_range = (minval, maxval)
@parameter
def largest_crater(self):
"""
The largest crater diameter in meters. Set during initialization.
"""
if not isinstance(self.production, Production):
return self.max_crater_diameter_range[1]
if self.production.generator_type == "crater":
return self.production.diameter_range[1]
elif self.production.generator_type == "projectile":
projectile_diameter = self.production.diameter_range[1]
if not np.isinf(projectile_diameter):
return self._get_diameter_limit("largest", projectile_diameter=projectile_diameter)
else:
return projectile_diameter
@largest_crater.setter
def largest_crater(self, value):
if not isinstance(self.production, Production):
return
if self.production.generator_type == "crater":
minval = max(self.production.diameter_range[0], self.max_crater_diameter_range[0])
maxval = min(value, self.max_crater_diameter_range[1])
elif self.production.generator_type == "projectile":
minval = max(
self.production.diameter_range[0],
self._get_diameter_limit("largest", crater_diameter=self.max_crater_diameter_range[0]),
)
maxval = self._get_diameter_limit("largest", crater_diameter=value)
else:
raise ValueError("Unknown Production generator_type. Expecting either 'crater' or 'projectile'")
self.production.diameter_range = (minval, maxval)
@parameter
def smallest_projectile(self):
"""
The smallest projectile diameter in meters. Set during initialization.
"""
if not isinstance(self.production, Production):
return self._get_diameter_limit("smallest", crater_diameter=self.max_crater_diameter_range[0])
if self.production.generator_type == "projectile":
return self.production.diameter_range[0]
elif self.production.generator_type == "crater":
crater_diameter = self.production.diameter_range[0]
return self._get_diameter_limit("smallest", crater_diameter=crater_diameter)
@smallest_projectile.setter
def smallest_projectile(self, value):
if not isinstance(self.production, Production):
return
if self.production.generator_type == "projectile":
minval = max(value, self._get_diameter_limit("smallest", crater_diameter=self.max_crater_diameter_range[0]))
maxval = min(self.production.diameter_range[1], self._get_diameter_limit("largest", self.max_crater_diameter_range[1]))
elif self.production.generator_type == "crater":
minval = max(value, self.production.diameter_range[0])
maxval = min(self.production.diameter_range[1], self.max_crater_diameter_range[1])
else:
raise ValueError("Unknown Production generator_type. Expecting either 'crater' or 'projectile'")
self.production.diameter_range = (minval, maxval)
@parameter
def largest_projectile(self):
"""
The largest projectile diameter in meters. Set during initialization.
"""
if not isinstance(self.production, Production):
return self._get_diameter_limit("largest", crater_diameter=self.max_crater_diameter_range[1])
if self.production.generator_type == "projectile":
return self.production.diameter_range[1]
elif self.production.generator_type == "crater":
crater_diameter = self.production.diameter_range[1]
return self._get_diameter_limit("largest", crater_diameter=crater_diameter)
@largest_projectile.setter
def largest_projectile(self, value):
if not isinstance(self.production, Production):
return
if self.production.generator_type == "projectile":
minval = max(
self.production.diameter_range[0],
self._get_diameter_limit("smallest", crater_diameter=self.max_crater_diameter_range[0]),
)
maxval = min(value, self._get_diameter_limit("largest", crater_diameter=self.max_crater_diameter_range[1]))
elif self.production.generator_type == "crater":
minval = max(self.production.diameter_range[0], self.max_crater_diameter_range[1])
maxval = self._get_diameter_limit("largest", crater_diameter=value)
else:
raise ValueError("Unknown Production generator_type. Expecting either 'crater' or 'projectile'")
self.production.diameter_range = (minval, maxval)
@property
def name(self):
"""
The name of the simulation.
"""
return "Cratermaker Simulation object"
@property
def config_file(self):
"""
The path to the configuration file for the simulation.
"""
return self.simdir / _CONFIG_FILE_NAME
@property
def config_readonly(self) -> bool:
"""
Flag indicating whether the configuration is read-only.
"""
return self._config_readonly
@property
def time_variables(self) -> dict[str, float]:
"""
A dictionary of time-related variables for the simulation.
Returns
-------
dict[str, float]
A dictionary containing the following keys:
- "time": The current age in My relative to the present.
- "elapsed_time": The elapsed time in My since the start of the simulation.
- "elapsed_n1": The elapsed number of craters larger than 1 km in diameter.
"""
if self.time is None:
return {}
else:
return {
"time": self.time,
"elapsed_time": self.elapsed_time,
"elapsed_n1": self.elapsed_n1,
}
@property
def do_counting(self) -> bool:
"""
A boolean flag indicating whether or not counting is enabled for the simulation. This is determined by whether or not a counting model is present and has counting enabled.
Returns
-------
bool
True if counting is enabled, False otherwise.
"""
return self.counting is not None and self.morphology.do_counting if self.morphology is not None else False
@property
def observed(self) -> dict[Crater] | None:
"""
Pass-through to retrieve the current observed craters from the counting model, if it is enabled.
"""
if self.do_counting:
return self.counting.observed
else:
return None
@property
def emplaced(self) -> list[Crater] | None:
"""
Pass-through to retrieve the current emplaced craters from the morphology model.
"""
return self.counting.emplaced
@property
def n_observed(self) -> int | None:
"""
Pass-through to retrieve the current number of observed craters from the counting model, if it is enabled.
"""
if self.do_counting:
return self.counting.n_observed
else:
return None
@property
def n_emplaced(self) -> int | None:
"""
Pass-through to retrieve the current number of emplaced craters from the counting model.
"""
return self.counting.n_emplaced
[docs]
def output_filename(self, interval=None, **kwargs):
return None
@property
def is_new(self) -> bool:
"""
A boolean flag indicating that this is a new simulation run, which will trigger a reset action when run is called.
"""
return self._is_new
@is_new.setter
def is_new(self, value: bool):
if not isinstance(value, bool):
raise TypeError("is_new must be a boolean value")
self._is_new = value
@property
def Crater(self) -> type[Crater]:
"""
The Crater class used for crater generation in the simulation. Set during initialization.
This is a property that returns |morphology.Crater| if a morphology model is present, or the base |Crater| class if no morphology model is present.
"""
return self.morphology.Crater if self.morphology is not None else Crater
@property
def do_quasimc(self) -> bool:
"""
A boolean flag indicating whether or not quasi-Monte Carlo sampling is being used for crater generation in the simulation. This is determined by whether or not the production function is using quasi-Monte Carlo sampling.
Returns
-------
bool
True if quasi-Monte Carlo sampling is being used, False otherwise.
"""
if self._do_quasimc is None:
self._do_quasimc = (
len(self.production.quasimc_craters) > 0
if self.production is not None and self.production.quasimc_craters is not None
else False
)
return self._do_quasimc
@do_quasimc.setter
def do_quasimc(self, value: bool | None):
if not isinstance(value, (bool, None)):
raise TypeError("do_quasimc must be a boolean value or None")
self._do_quasimc = value
@property
def quasimc_craters(self) -> list[Crater]:
"""
List of craters to be emplaced using quasi-Monte Carlo.
When assigned a list of Crater objects with production metadata (production_time, production_ND, and/or production_sequence), they will be processed to set their |crater.time| value. Any craters without production metadata will be dropped. When accessed, if any of the craters have a None value for time, they will be reprocessed to set their time value. This allows for the quasimc_craters to be dynamically updated.
This is a wrapper for |production.quasimc_craters|.
"""
if self.production is not None:
return self.production.quasimc_craters
@quasimc_craters.setter
def quasimc_craters(self, value: list[Crater] | None):
if self.production is not None:
self.production.quasimc_craters = value
@property
def quasimc_file(self) -> Path:
"""
File containing the quasi-Monte Carlo craters.
This is a wrapper for |production.quasimc_file|.
"""
if self.production is not None:
return self.production.quasimc_file
@quasimc_file.setter
def quasimc_file(self, value: Path | None):
if self.production is not None:
self.production.quasimc_file = value
[docs]
def quasimc_merge(
self,
craters: list[Crater],
time_start: float | None = None,
time_end: float | None = None,
N_D: PairOfFloats | None = None,
N_D_end: PairOfFloats | None = None,
**kwargs,
):
"""
Merge a randomly-generated list of craters with the quasi-Monte Carlo craters over a given time interval.
It will extract any craters from |production.quasimc_craters| that have time values that overlap those of the arguments. Then it will successively check the to see if any in the "craters" list are larger than the largest in
Parameters
----------
craters : list[Crater]
The list of Crater objects to merge with the quasi-Monte Carlo craters.
time_start : float, optional
The starting time in units of My relative to the present in which to extract the quasi-Monte Carlo craters.
time_end : float, optional
The ending time in units of My relative to the present in which to extract the quasi-Monte Carlo craters.
N_D : PairOfFloats, optional
A pair of D, N values that represent the N(D) format for the cumulative number density at the start of the sample, which is used to extract the quasi-Monte Carlo craters. The units given by |production.N_D_units| and can be adjusted by setting |production.D_conversion_factor| and |production.N_conversion_factor|.
N_D_end : PairOfFloats, optional
A pair of D, N values that represent the N(D) format for the cumulative number density at the end of the sample, which is used to extract the quasi-Monte Carlo craters. The units given by |production.N_D_units| and can be adjusted by setting |production.D_conversion_factor| and |production.N_conversion_factor|.
**kwargs: Any
|kwargs|
Returns
-------
list[Crater]
The merged list of Crater objects, sorted in order of decreasing age.
Notes
-----
This is a wrapper for |production.quasimc_merge|
"""
if self.production is not None:
return self.production.quasimc_merge(
craters=craters, time_start=time_start, time_end=time_end, N_D=N_D, N_D_end=N_D_end, **kwargs
)
else:
raise RuntimeError("Production model must be defined to use quasimc_merge")
@property
def face_variables(self) -> list[str]:
"""
Returns a list of the data variables currently being stored on the faces.
"""
return self.surface.face_variables
@CratermakerBase.save_actions.setter
def save_actions(self, value):
if value == "default":
value = [{"plot": {"plot_style": "hillshade", "show": False, "save": True}}]
CratermakerBase.save_actions.fset(self, value)
