From e9592b549ea3be62fdf5e991e8b778f119e36442 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 25 Oct 2026 01:03:33 +0000 Subject: [PATCH 1/3] Blend slip noise families for depth-dependent correlation lengths DrawnSlip gains an optional tuple of NoiseFamily, each a stationary von Karman covariance (its own relation and Hurst exponent, or the source's) and a weight over position. Slip's latent then comes from a Blend: every family filters one white noise on a shared padded grid, crossfaded at equal power and normalised by the exact cross-correlations, so the latent stays standard normal at every cell. Rise time and the onset perturbation draw their independent part from the same Blend; the perturbation can keep the source's stationary covariance instead, as genslip v5.6.3 does. Rake is unchanged. CorrelationRelation gains saturation magnitudes, and SUZUKI_2022 is Suzuki et al. (2022)'s relation, the shallow family of Graves & Pitarka's hybrid. Configs take the families as [[slip.noise]]. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01KrujMjfczwcCciJwtErEbs --- rupture_generator/__init__.py | 5 +- rupture_generator/config.py | 28 +++- rupture_generator/rupture/models.py | 120 ++++++++++++-- rupture_generator/rupture/source.py | 33 +++- rupture_generator/sampling/__init__.py | 4 + rupture_generator/sampling/field.py | 166 +++++++++++++++++-- tests/test_blend.py | 217 +++++++++++++++++++++++++ 7 files changed, 536 insertions(+), 37 deletions(-) create mode 100644 tests/test_blend.py diff --git a/rupture_generator/__init__.py b/rupture_generator/__init__.py index 028e3d0..1817bdc 100644 --- a/rupture_generator/__init__.py +++ b/rupture_generator/__init__.py @@ -15,8 +15,9 @@ :mod:`rupture_generator.config` builds the same inputs from a TOML file, and the command line runs that file. :func:`rupture_generator.sampling.sampler` keeps the -last 64 covariance embeddings, two per segment, so drawing a fault system of up to -32 segments again reuses every one. ``sampler.cache_clear()`` releases them. +last 128 covariance embeddings, two per segment and one more per extra slip noise +family, so drawing a fault system of up to 64 segments again reuses every one. +``sampler.cache_clear()`` releases them. """ from importlib.metadata import version diff --git a/rupture_generator/config.py b/rupture_generator/config.py index 3dd6208..0436462 100644 --- a/rupture_generator/config.py +++ b/rupture_generator/config.py @@ -64,6 +64,7 @@ CorrelatedRise, DrawnSlip, GaussianRake, + NoiseFamily, PulseSettings, RuptureModel, ShearFractionSpeed, @@ -241,10 +242,35 @@ def field(self) -> SpatialField: # ----------------------------------------------------------- the default models +@dataclasses.dataclass(frozen=True) +class Noise(NoiseFamily): + """``[[slip.noise]]``: :class:`~rupture_generator.rupture.models.NoiseFamily`. + + Parameters + ---------- + weight : DepthProfile + ``[slip.noise.weight]``, how much of slip's latent the family is. + relation : Correlation, optional + ``[slip.noise.relation]``, the family's correlation lengths. + """ + + weight: DepthProfile + relation: Correlation | None = None + Config = _Strict + + @dataclasses.dataclass(frozen=True) class Slip(DrawnSlip): - """``[slip]``: :class:`~rupture_generator.rupture.models.DrawnSlip`.""" + """``[slip]``: :class:`~rupture_generator.rupture.models.DrawnSlip`. + + Parameters + ---------- + noise : tuple of Noise + ``[[slip.noise]]``, the families slip's latent blends. A segment's list + replaces the global one whole. + """ + noise: tuple[Noise, ...] = () Config = _Strict diff --git a/rupture_generator/rupture/models.py b/rupture_generator/rupture/models.py index 54e322e..74a947c 100644 --- a/rupture_generator/rupture/models.py +++ b/rupture_generator/rupture/models.py @@ -23,7 +23,7 @@ import dataclasses from collections.abc import Callable, Mapping -from typing import Protocol +from typing import Literal, Protocol import numpy as np import scipy as sp @@ -31,15 +31,20 @@ from rupture_generator.errors import RuptureGeneratorError from rupture_generator.geometry import CellArray, CellMask, Geometry from rupture_generator.rupture.medium import SpatialField, constant_field, rigidity_pa -from rupture_generator.rupture.source import SegmentSource +from rupture_generator.rupture.source import ( + CorrelationRelation, + SegmentSource, + magnitude_from_moment, +) from rupture_generator.sampling import ( + Blend, Covariance, Grid, Marginal, MarginalFamily, PreCorrected, + VonKarman, mix, - sampler, standardise, ) @@ -131,7 +136,7 @@ class SlipField: Slip, in metres. The generator zeroes it off the fault. latent : CellArray A standard normal field that slip is monotone in. - latent_covariance : Covariance + latent_covariance : Covariance or Blend The covariance the model drew ``latent`` from, or the nearest one for a slip it didn't draw. A model that mixes its own draw into ``latent`` draws from this. @@ -141,7 +146,7 @@ class SlipField: slip_m: CellArray latent: CellArray - latent_covariance: Covariance + latent_covariance: Covariance | Blend statistics: Mapping[str, float] = dataclasses.field(default_factory=dict) @@ -512,7 +517,7 @@ def _independent( context: SegmentContext, slip: SlipField, rng: np.random.Generator ) -> CellArray: """A standard normal field independent of slip's latent, with its covariance.""" - independent, _ = sampler(context.grid, slip.latent_covariance).draw(rng) + independent, _ = slip.latent_covariance.sampler(context.grid).draw(rng) return independent @@ -532,6 +537,51 @@ def _reach(mask: CellMask, axis: int, *, reverse: bool) -> np.ndarray: return np.flip(counted, axis) if reverse else counted +@dataclasses.dataclass(frozen=True) +class NoiseFamily: + """One stationary noise slip draws from, and where on the fault it applies. + + Several families make Graves & Pitarka's hybrid slip: long correlation lengths near + the surface, short ones at depth, and a crossfade between them. + + Attributes + ---------- + weight : SpatialField + How much of slip's latent this family is at each position, in ``[0, 1]``. + relation : CorrelationRelation, optional + The family's correlation lengths at the segment's magnitude. The source's + lengths when absent. + hurst : float, optional + The family's von Karman Hurst exponent. The source's correlation function when + absent. + """ + + weight: SpatialField + relation: CorrelationRelation | None = None + hurst: float | None = None + + def covariance(self, source: SegmentSource) -> Covariance: + """The family's covariance on a segment. + + Parameters + ---------- + source : SegmentSource + The segment's source, whose moment sets the magnitude. + + Returns + ------- + Covariance + The family's correlation function and lengths. + """ + lengths = source.covariance.lengths_km + if self.relation is not None: + lengths = self.relation.lengths_km(magnitude_from_moment(source.moment_nm)) + correlation = source.covariance.correlation + if self.hurst is not None: + correlation = VonKarman(self.hurst) + return Covariance(correlation, lengths) + + @dataclasses.dataclass(frozen=True) class DrawnSlip: """Slip drawn from the source's covariance, with an exact marginal and tapered edges. @@ -557,6 +607,9 @@ class DrawnSlip: bottom_taper : float The fraction of the fault's width from the bottom edge over which slip ramps to zero. + noise : tuple of NoiseFamily + The families slip's latent blends, each pre-corrected for the marginal. The + source's covariance alone when empty. """ family: MarginalFamily = "truncated_exponential" @@ -565,6 +618,7 @@ class DrawnSlip: side_taper: float = 0.02 top_taper: float = 0.0 bottom_taper: float = 0.0 + noise: tuple[NoiseFamily, ...] = () def __post_init__(self) -> None: """Refuse a taper outside ``[0, 0.5]``, or a spread outside the family's range.""" @@ -605,6 +659,33 @@ def _taper(self, field: CellArray, occupied: CellMask) -> CellArray: ramp *= np.minimum(_reach(occupied, axis, reverse=reverse) / cells, 1.0) return field * ramp * occupied + def _latent_covariance(self, context: SegmentContext) -> Covariance | Blend: + """The covariance to draw the latent from, pre-corrected for the marginal.""" + + def corrected(covariance: Covariance) -> Covariance: + return Covariance( + PreCorrected(covariance.correlation, self.marginal), + covariance.lengths_km, + ) + + if not self.noise: + return corrected(context.source.covariance) + return Blend( + tuple( + corrected(family.covariance(context.source)) for family in self.noise + ), + tuple( + read_at_centres( + context.geometry, + "a slip noise family's weight", + family.weight, + lambda weight: (weight >= 0.0) & (weight <= 1.0), + "weights in [0, 1]", + ) + for family in self.noise + ), + ) + def __call__(self, context: SegmentContext, rng: np.random.Generator) -> SlipField: """Draw slip on a segment. @@ -623,14 +704,12 @@ def __call__(self, context: SegmentContext, rng: np.random.Generator) -> SlipFie Raises ------ RuptureGeneratorError - If the covariance doesn't embed on this chart, or the taper leaves no - slip anywhere. + If the covariance doesn't embed on this chart, a noise family's weight + reads out of range, or the taper leaves no slip anywhere. """ - geometry, covariance = context.geometry, context.source.covariance - latent_covariance = Covariance( - PreCorrected(covariance.correlation, self.marginal), covariance.lengths_km - ) - latent, _ = sampler(context.grid, latent_covariance).draw(rng) + geometry = context.geometry + latent_covariance = self._latent_covariance(context) + latent, _ = latent_covariance.sampler(context.grid).draw(rng) pattern = self._taper(self.marginal.apply(latent), geometry.occupied) total = context.moment_nm(pattern) if not total > 0.0: @@ -806,7 +885,7 @@ def __call__(self, context: SegmentContext, rng: np.random.Generator) -> Drawn: Drawn Rake in degrees. """ - latent, _ = sampler(context.grid, context.source.covariance).draw(rng) + latent, _ = context.source.covariance.sampler(context.grid).draw(rng) return Drawn(context.source.rake_deg + self.sigma_deg * standardise(latent)) @@ -908,12 +987,16 @@ class CorrelatedPerturbation: blend_sigma : float The width of the zone over which the displacement grows in from the seed, in units of the spread. + independent : {"slip", "source"} + The covariance of the displacement's own draw: slip's latent's, which follows + any slip noise families, or the source's, which genslip v5.6.3 keeps. """ offset_s: float = 0.1 coefficient: float = 0.5 correlation: float = 0.8 blend_sigma: float = 4.0 + independent: Literal["slip", "source"] = "slip" def __post_init__(self) -> None: """Refuse a negative spread, or a blend of no width.""" @@ -956,11 +1039,12 @@ def __call__( # genslip's tsfac is negative and this package reads it as a magnitude, so the # displacement subtracts: with a positive correlation, high slip ruptures early. rho = self.correlation + if self.independent == "slip": + independent = _independent(context, slip, rng) + else: + independent, _ = context.source.covariance.sampler(context.grid).draw(rng) displacement_s = -scale_s * standardise( - mix( - (rho, slip.latent), - (np.sqrt(1.0 - rho * rho), _independent(context, slip, rng)), - ) + mix((rho, slip.latent), (np.sqrt(1.0 - rho * rho), independent)) ) return Perturbation( displacement_s, self.blend_sigma * scale_s, {"scale_s": scale_s} diff --git a/rupture_generator/rupture/source.py b/rupture_generator/rupture/source.py index 95f7411..6ed7905 100644 --- a/rupture_generator/rupture/source.py +++ b/rupture_generator/rupture/source.py @@ -12,6 +12,10 @@ Mai, P. M., & Beroza, G. C. (2002). A spatial random field model to characterize complexity in earthquake slip. *Journal of Geophysical Research*, 107(B11), 2308. + +Suzuki, F., Kato, K., & Watanabe, T. (2022). Spatial complexity characteristics of fault +slip from Japanese earthquakes based on self-similar models. *Bulletin of the +Seismological Society of America*, 112(4), 1853-1870. """ import dataclasses @@ -65,8 +69,9 @@ def magnitude_from_moment(moment_nm: float) -> float: class CorrelationRelation: """Slip correlation lengths that grow with magnitude. - ``length = 10 ** (exponent * Mw - offset)`` kilometres, along strike and down dip. - The defaults are Mai & Beroza (2002)'s fit for the von Karman correlation. + ``length = 10 ** (exponent * min(Mw, saturation) - offset)`` kilometres, along strike + and down dip. The defaults are Mai & Beroza (2002)'s fit for the von Karman + correlation, which never saturates. Parameters ---------- @@ -78,12 +83,19 @@ class CorrelationRelation: How fast the length down dip grows with magnitude. dip_offset : float The length down dip's offset, in decades of kilometres. + strike_saturation : float + The magnitude past which the length along strike stops growing. + dip_saturation : float + The magnitude past which the length down dip stops growing, such as the + seismogenic width bounding it. """ strike_exponent: float = 0.5 strike_offset: float = 2.5 dip_exponent: float = 1.0 / 3.0 dip_offset: float = 1.5 + strike_saturation: float = np.inf + dip_saturation: float = np.inf def lengths_km(self, magnitude: float) -> tuple[float, float]: """Compute the correlation lengths for a magnitude. @@ -99,8 +111,16 @@ def lengths_km(self, magnitude: float) -> tuple[float, float]: ``(down dip, along strike)`` in kilometres, the chart's axis order. """ return ( - 10.0 ** (self.dip_exponent * magnitude - self.dip_offset), - 10.0 ** (self.strike_exponent * magnitude - self.strike_offset), + 10.0 + ** ( + self.dip_exponent * min(magnitude, self.dip_saturation) + - self.dip_offset + ), + 10.0 + ** ( + self.strike_exponent * min(magnitude, self.strike_saturation) + - self.strike_offset + ), ) @@ -110,6 +130,10 @@ def lengths_km(self, magnitude: float) -> tuple[float, float]: MAI_BEROZA = CorrelationRelation() """Mai & Beroza (2002): ``10^(Mw/2 - 2.5)`` along strike, ``10^(Mw/3 - 1.5)`` down dip.""" +SUZUKI_2022 = CorrelationRelation(0.5, 1.67, 0.5, 1.69, dip_saturation=6.3) +"""Suzuki et al. (2022): ``10^(Mw/2 - 1.67)`` along strike, ``10^(min(Mw, 6.3)/2 - 1.69)`` +down dip. Smoother than Mai & Beroza, the shallow slip of Graves & Pitarka's hybrid.""" + @dataclasses.dataclass(frozen=True) class SegmentSource: @@ -208,6 +232,7 @@ def split_moment( __all__ = [ "MAI_BEROZA", + "SUZUKI_2022", "VON_KARMAN", "CorrelationRelation", "SegmentSource", diff --git a/rupture_generator/sampling/__init__.py b/rupture_generator/sampling/__init__.py index 39f0d81..7bd418f 100644 --- a/rupture_generator/sampling/__init__.py +++ b/rupture_generator/sampling/__init__.py @@ -1,6 +1,8 @@ """Random fields on the fault: Gaussian machinery, correlation functions, NORTA.""" from rupture_generator.sampling.field import ( + Blend, + BlendedSampler, Correlation, Covariance, FieldArray, @@ -22,6 +24,8 @@ __all__ = [ "HURST", "NORMAL", + "Blend", + "BlendedSampler", "Correlation", "Covariance", "FieldArray", diff --git a/rupture_generator/sampling/field.py b/rupture_generator/sampling/field.py index 284a94b..051e0b1 100644 --- a/rupture_generator/sampling/field.py +++ b/rupture_generator/sampling/field.py @@ -3,7 +3,9 @@ This module embeds any correlation function and doesn't assume anything about what the caller does with the field. A :class:`Covariance` is a :class:`Correlation` and the lengths that scale its argument. :func:`sampler` embeds one on a :class:`Grid` once -and draws from it as often as asked. Every tuple is ``(axis 0, axis 1)``. +and draws from it as often as asked. A :class:`Blend` crossfades several covariances +cell by cell, for a field whose correlation lengths change across the grid. Every tuple +is ``(axis 0, axis 1)``. """ import dataclasses @@ -42,9 +44,10 @@ MAXIMUM_EMBEDDING_CELLS = 1 << 27 """The largest padded grid to transform. """ -SAMPLER_CACHE_SIZE = 64 -"""How many embeddings :func:`sampler` keeps. Each segment takes two, so a fault system -of up to 32 segments draws again without embedding anything.""" +SAMPLER_CACHE_SIZE = 128 +"""How many embeddings :func:`sampler` keeps. Each segment takes two, and one more for +each family a :class:`Blend` adds, so a fault system of up to 64 segments, or 42 with +two families each, draws again without embedding anything.""" @dataclasses.dataclass(frozen=True) @@ -82,6 +85,21 @@ class Covariance: correlation: Correlation lengths_km: tuple[float, float] + def sampler(self, grid: Grid) -> "Sampler": + """Embed the covariance on ``grid``, through :func:`sampler`'s cache. + + Parameters + ---------- + grid : Grid + The grid to draw on. + + Returns + ------- + Sampler + The embedding, ready to draw from. + """ + return sampler(grid, self) + def _decay_length(correlation: Correlation, tolerance: float) -> float: """The lag, in correlation lengths, past which the correlation is below ``tolerance``.""" @@ -126,19 +144,26 @@ def _mirror_weights(extent: int) -> np.ndarray: return weights -def _embed(grid: Grid, covariance: Covariance) -> tuple[FieldArray, tuple[int, int]]: +def _embed( + grid: Grid, covariance: Covariance, floor: tuple[int, int] = (0, 0) +) -> tuple[FieldArray, tuple[int, int]]: """The square-rooted eigenvalues of a circulant embedding, and the padded shape. The margin starts at the correlation's decay length and doubles until the variance the clipped negative eigenvalues drop is within :data:`MAXIMUM_VARIANCE_DEFICIT`. + The padded shape is at least ``floor`` on each axis. """ decay = _decay_length(covariance.correlation, WRAP_TOLERANCE) for doubling in range(MAXIMUM_EMBEDDING_DOUBLINGS): margin = decay * 2**doubling rows, columns = ( - _padded_extent(extent, resolution, length, margin) - for extent, resolution, length in zip( - grid.shape, grid.resolution_km, covariance.lengths_km, strict=True + max(_padded_extent(extent, resolution, length, margin), least) + for extent, resolution, length, least in zip( + grid.shape, + grid.resolution_km, + covariance.lengths_km, + floor, + strict=True, ) ) padded = (rows, columns) @@ -212,12 +237,18 @@ def draw(self, rng: np.random.Generator) -> tuple[FieldArray, FieldArray]: tuple of FieldArray The real and imaginary fields, each shaped like the grid. """ - seed = int(rng.integers(0, np.iinfo(np.uint64).max, dtype=np.uint64)) - return circulant_draw(self.amplitudes, self.padded, self.shape, seed) + return circulant_draw(self.amplitudes, self.padded, self.shape, _seed(rng)) + + +def _seed(rng: np.random.Generator) -> int: + """Take one draw's seed from ``rng``.""" + return int(rng.integers(0, np.iinfo(np.uint64).max, dtype=np.uint64)) @functools.lru_cache(maxsize=SAMPLER_CACHE_SIZE) -def sampler(grid: Grid, covariance: Covariance) -> Sampler: +def sampler( + grid: Grid, covariance: Covariance, floor: tuple[int, int] = (0, 0) +) -> Sampler: """Embed ``covariance`` on ``grid``, or return the embedding already made. Parameters @@ -226,6 +257,8 @@ def sampler(grid: Grid, covariance: Covariance) -> Sampler: The grid to draw on. covariance : Covariance The covariance to embed. + floor : tuple of int, optional + The least padded shape, so covariances that share noise share a padded grid. Returns ------- @@ -238,7 +271,116 @@ def sampler(grid: Grid, covariance: Covariance) -> Sampler: If the covariance doesn't embed within the cell cap, or isn't positive definite at any margin tried. """ - return Sampler(grid.shape, *_embed(grid, covariance)) + return Sampler(grid.shape, *_embed(grid, covariance, floor)) + + +@dataclasses.dataclass(frozen=True, eq=False) +class Blend: + """Stationary covariances that filter one white noise, crossfaded cell by cell. + + It is the covariance of a field whose correlation lengths change across the grid, + such as slip with long lengths near the surface and short ones at depth. Every + component filters the same noise, so a feature changes shape through a crossfade + rather than fading out. A cell's loadings are ``sin(pi w / 2)``, scaled to unit + variance from the components' exact cross-correlations, so the field is standard + normal at every cell however the weights overlap. + + Attributes + ---------- + components : tuple of Covariance + The stationary covariances. + weights : tuple of FieldArray + Each component's weight at each cell, in ``[0, 1]``. Every cell needs a + positive weight. + """ + + components: tuple[Covariance, ...] + weights: tuple[FieldArray, ...] + + def sampler(self, grid: Grid) -> "BlendedSampler": + """Embed every component on ``grid`` at one padded shape, and load them. + + Parameters + ---------- + grid : Grid + The grid to draw on, shaped like the weights. + + Returns + ------- + BlendedSampler + The embeddings and each one's loading, ready to draw from. + + Raises + ------ + RuptureGeneratorError + If a cell has no positive weight, or a component doesn't embed. + """ + gains = np.sin(0.5 * np.pi * np.asarray(self.weights)) + power = np.sum(np.square(gains), axis=0) + if not np.all(power > 0.0): + raise RuptureGeneratorError( + f"{int(np.sum(~(power > 0.0)))} cells have no positive weight in the blend" + ) + floor = (0, 0) + while True: + samplers = [ + sampler(grid, component, floor) for component in self.components + ] + widest = tuple(max(axis) for axis in zip(*(s.padded for s in samplers))) + if all(s.padded == widest for s in samplers): + break + floor = widest + weights = np.outer(*(_mirror_weights(extent) for extent in widest)) + amplitudes = np.array([s.amplitudes for s in samplers]) + cross = np.einsum("ikl,jkl,kl->ij", amplitudes, amplitudes, weights) + cross /= widest[0] * widest[1] + gains /= np.sqrt(power) + variance = np.einsum("ikl,ij,jkl->kl", gains, cross, gains) + return BlendedSampler(tuple(samplers), gains / np.sqrt(variance)) + + +@dataclasses.dataclass(frozen=True, eq=False) +class BlendedSampler: + """A :class:`Blend`'s components, embedded on one padded grid, and their loadings. + + Attributes + ---------- + samplers : tuple of Sampler + The components' embeddings, all at one padded shape. + loadings : numpy.ndarray + Each component's loading at each cell, shaped ``(components, *grid)``. + """ + + samplers: tuple[Sampler, ...] + loadings: np.ndarray + + def draw(self, rng: np.random.Generator) -> tuple[FieldArray, FieldArray]: + """Draw a pair of independent blended fields, standard normal at every cell. + + Every component draws from the one seed, so they filter the same noise. The + imaginary parts are independent of every real part, because the amplitudes + are even. + + Parameters + ---------- + rng : np.random.Generator + The generator the draw's seed comes from. + + Returns + ------- + tuple of FieldArray + The real and imaginary fields, each shaped like the grid. + """ + seed = _seed(rng) + real = np.zeros(self.loadings.shape[1:]) + imaginary = np.zeros(self.loadings.shape[1:]) + for loading, embedding in zip(self.loadings, self.samplers, strict=True): + parts = circulant_draw( + embedding.amplitudes, embedding.padded, embedding.shape, seed + ) + real += loading * parts[0] + imaginary += loading * parts[1] + return real, imaginary def mix(*terms: tuple[float | FieldArray, FieldArray]) -> FieldArray: diff --git a/tests/test_blend.py b/tests/test_blend.py new file mode 100644 index 0000000..a867f7f --- /dev/null +++ b/tests/test_blend.py @@ -0,0 +1,217 @@ +"""Slip noise families: the blended sampler, the models that use it, the config.""" + +import dataclasses + +import numpy as np +import pytest + +from rupture_generator import RuptureGeneratorError, generate +from rupture_generator import config as rupture_config +from rupture_generator.rupture.models import ( + CorrelatedPerturbation, + DrawnSlip, + NoiseFamily, + RuptureModel, + SegmentContext, +) +from rupture_generator.rupture.source import MAI_BEROZA, SUZUKI_2022 +from rupture_generator.sampling import Blend, Covariance, Grid, VonKarman, sampler +from rupture_generator.sampling.field import _mirror_weights + +from .conftest import EXAMPLE + +GRID = Grid((24, 60), (0.5, 0.5)) +SHORT = Covariance(VonKarman(), (2.0, 3.0)) +LONG = Covariance(VonKarman(), (4.0, 12.0)) + + +def _ramp(rows: int = GRID.shape[0]) -> np.ndarray: + """Weight of the deep family: 0 above row 8, 1 below row 16, linear between.""" + column = np.clip((np.arange(rows) - 8.0) / 8.0, 0.0, 1.0) + return np.repeat(column[:, None], GRID.shape[1], axis=1) + + +def _blend() -> Blend: + deep = _ramp() + return Blend((LONG, SHORT), (1.0 - deep, deep)) + + +def _cross(samplers) -> np.ndarray: + padded = samplers[0].padded + weights = np.outer(*(_mirror_weights(extent) for extent in padded)) + return np.array( + [ + [np.sum(weights * a.amplitudes * b.amplitudes) for b in samplers] + for a in samplers + ] + ) / (padded[0] * padded[1]) + + +def test_saturation_stops_one_axis_growing(): + below, above = SUZUKI_2022.lengths_km(6.3), SUZUKI_2022.lengths_km(7.0) + assert above[0] == pytest.approx(below[0]) + assert above[1] > below[1] + assert MAI_BEROZA.lengths_km(7.0) == pytest.approx( + (10.0 ** (7.0 / 3.0 - 1.5), 10.0 ** (3.5 - 2.5)) + ) + + +def test_a_floor_only_ever_widens_the_embedding(): + alone = sampler(GRID, SHORT) + widened = sampler(GRID, SHORT, (alone.padded[0] * 2, alone.padded[1])) + assert widened.padded == (alone.padded[0] * 2, alone.padded[1]) + assert sampler(GRID, SHORT, (2, 2)).padded == alone.padded + + +def test_the_components_share_the_widest_padded_shape(): + across = Covariance(VonKarman(), (12.0, 2.0)) + blended = Blend((across, SHORT, LONG), (np.full(GRID.shape, 1.0),) * 3).sampler( + GRID + ) + own = [sampler(GRID, c).padded for c in (across, SHORT, LONG)] + widest = tuple(max(axis) for axis in zip(*own)) + assert all(s.padded == widest for s in blended.samplers) + + +def test_every_cell_has_unit_variance_exactly(): + blended = _blend().sampler(GRID) + variance = np.einsum( + "ikl,ij,jkl->kl", blended.loadings, _cross(blended.samplers), blended.loadings + ) + np.testing.assert_allclose(variance, 1.0, rtol=1e-12) + + +def test_a_cell_of_one_family_is_that_familys_draw(): + blended = _blend().sampler(GRID) + real, imaginary = blended.draw(np.random.default_rng(3)) + # The same seed draws the same noise for every component. + seed = int( + np.random.default_rng(3).integers(0, np.iinfo(np.uint64).max, dtype=np.uint64) + ) + from rupture_generator._kernels import circulant_draw + + long_draw = circulant_draw( + blended.samplers[0].amplitudes, blended.samplers[0].padded, GRID.shape, seed + ) + scale = np.sqrt(_cross(blended.samplers)[0, 0]) + np.testing.assert_allclose(real[:8], long_draw[0][:8] / scale, rtol=1e-12) + np.testing.assert_allclose(imaginary[:8], long_draw[1][:8] / scale, rtol=1e-12) + + +def test_draws_are_standard_normal_and_the_pair_independent(): + blended = _blend().sampler(GRID) + rng = np.random.default_rng(11) + pairs = np.array([blended.draw(rng) for _ in range(600)]) + real, imaginary = pairs[:, 0], pairs[:, 1] + # Through the crossfade, rows 8 to 16, as well as either side of it. + assert np.abs(real.var(axis=0) - 1.0).max() < 0.25 + assert abs(float(real.var() - 1.0)) < 0.05 + assert abs(float(np.mean(real * imaginary))) < 0.02 + + +def test_a_cell_with_no_weight_is_refused(): + deep = _ramp() + with pytest.raises(RuptureGeneratorError, match="no positive weight"): + Blend((LONG, SHORT), (np.zeros(GRID.shape), deep)).sampler(GRID) + + +# ----------------------------------------------------------------- the models + + +def _hybrid_slip() -> DrawnSlip: + def deep(positions): + return np.clip((positions[..., 2] - 4.0) / 4.0, 0.0, 1.0) + + return DrawnSlip( + noise=( + NoiseFamily(lambda p: 1.0 - deep(p), SUZUKI_2022, hurst=0.6), + NoiseFamily(deep), + ) + ) + + +def _context(scenario, name="bent"): + chart = scenario.realisation[name] + return SegmentContext( + chart, + scenario.sources[name], + scenario.medium.shear_speed_km_s(chart.centres), + scenario.medium.density_g_cm3(chart.centres), + ) + + +def test_slip_draws_a_blend_its_partners_follow(scenario): + context = _context(scenario) + slip = _hybrid_slip()(context, np.random.default_rng(5)) + blend = slip.latent_covariance + assert isinstance(blend, Blend) + long, short = blend.components + assert long.correlation.target == VonKarman(0.6) + assert short.lengths_km == scenario.sources["bent"].covariance.lengths_km + assert context.moment_nm(slip.slip_m) == pytest.approx(context.source.moment_nm) + assert np.all(np.isfinite(slip.latent)) + + +def test_a_weight_out_of_range_is_refused(scenario): + slip = DrawnSlip(noise=(NoiseFamily(lambda p: np.full(p.shape[:-1], 2.0)),)) + with pytest.raises(RuptureGeneratorError, match=r"weights in \[0, 1\]"): + slip(_context(scenario), np.random.default_rng(0)) + + +def test_the_perturbation_can_keep_the_sources_noise(scenario): + def onset(independent): + model = RuptureModel( + slip=_hybrid_slip(), + perturbation=CorrelatedPerturbation(independent=independent), + ) + drawn = generate( + scenario.realisation, scenario.medium, scenario.sources, model, seed=9 + ) + return drawn["bent"].onset_s + + assert not np.allclose(onset("slip"), onset("source")) + + +# ----------------------------------------------------------------- the config + +NOISE = """ +[[slip.noise]] +relation = { strike_offset = 1.67, dip_exponent = 0.5, dip_offset = 1.69, dip_saturation = 6.3 } +weight = { type = "ramp", centre_km = 6.5, half_width_km = 1.5, shallow = 1.0, deep = 0.0 } + +[[slip.noise]] +weight = { type = "ramp", centre_km = 6.5, half_width_km = 1.5, shallow = 0.0, deep = 1.0 } +""" + + +def _hybrid_file(tmp_path, extra=""): + path = tmp_path / "hybrid.toml" + path.write_text( + EXAMPLE.read_text().replace( + "two_faults.geojson", str(EXAMPLE.with_name("two_faults.geojson")) + ) + + NOISE + + extra + ) + return path + + +def test_noise_families_load_and_round_trip(tmp_path): + loaded = rupture_config.load(_hybrid_file(tmp_path)) + shallow, deep = loaded.slip.noise + assert shallow.relation == rupture_config.Correlation( + strike_offset=1.67, dip_exponent=0.5, dip_offset=1.69, dip_saturation=6.3 + ) + assert deep.relation is None and deep.hurst is None + again = tmp_path / "again.toml" + again.write_text(rupture_config.dump(loaded)) + assert rupture_config.load(again) == loaded + + +def test_a_segments_noise_replaces_the_global_list(tmp_path): + extra = '\n[segments.next.slip]\nnoise = [{ weight = { type = "constant", value = 1.0 } }]\n' + loaded = rupture_config.load(_hybrid_file(tmp_path, extra)) + assert len(loaded.slip.noise) == 2 + segment = loaded.segments["next"].slip + assert segment is not None and len(segment.noise) == 1 + assert dataclasses.replace(segment, noise=loaded.slip.noise) == loaded.slip From d84a2da6dbfd2ad3b62934976843ffab2ce27361 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 25 Oct 2026 01:13:32 +0000 Subject: [PATCH 2/3] Give each model its own noise; one sampler for any number of covariances A Noise is a tuple of NoiseFamily (correlation relation, Hurst exponent, weight over position). DrawnSlip, GaussianRake and CorrelatedPerturbation each take one, so making rake hybrid is a config change. embed() builds a noise on a chart, pre-corrected for a marginal, and caches it; SlipField carries the built sampler, so rise time and the onset draw their independent part from it without embedding again. sampling.sampler() now embeds one or more covariances at a shared padded shape and returns one Sampler type; the stationary case is one covariance. Blend, BlendedSampler, Covariance.sampler and the sampler lru_cache go. SegmentSource is a moment and a rake; [source.correlation] and its per-segment override are replaced by [[.noise]] lists. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01KrujMjfczwcCciJwtErEbs --- README.md | 12 +- benchmarks/test_kernels.py | 6 +- examples/alpine_hope.toml | 6 - examples/beavan.toml | 6 - examples/colombia.toml | 6 - examples/crustal.toml | 6 - examples/landers.toml | 6 - examples/northridge.toml | 6 - examples/two_faults_segments.toml | 7 +- examples/venezuela.toml | 6 - rupture_generator/__init__.py | 15 +- rupture_generator/config.py | 85 ++++---- rupture_generator/rupture/models.py | 256 +++++++++++++------------ rupture_generator/rupture/source.py | 62 ++---- rupture_generator/sampling/__init__.py | 4 - rupture_generator/sampling/field.py | 231 +++++++--------------- tests/test_blend.py | 197 ++++++++----------- tests/test_config.py | 16 +- tests/test_models.py | 5 +- tests/test_numerics.py | 4 +- tests/test_pipeline.py | 8 +- 21 files changed, 396 insertions(+), 554 deletions(-) diff --git a/README.md b/README.md index bff9945..39a0b05 100644 --- a/README.md +++ b/README.md @@ -178,7 +178,13 @@ import dataclasses import numpy as np -from rupture_generator import SegmentContext, SlipField, gaussian_anamorphosis +from rupture_generator import ( + STATIONARY, + SegmentContext, + SlipField, + embed, + gaussian_anamorphosis, +) @dataclasses.dataclass(frozen=True, eq=False) @@ -191,9 +197,9 @@ class GivenSlip: return SlipField( slip_m=self.slip_m, latent=gaussian_anamorphosis(self.slip_m, context.geometry.occupied), - # The scores have no covariance of their own; the source's is the + # The scores have no noise of their own; the stationary one is the # nearest, and it is what the other models' draws use. - latent_covariance=context.source.covariance, + sampler=embed(STATIONARY, context.geometry, context.source.magnitude), ) ``` diff --git a/benchmarks/test_kernels.py b/benchmarks/test_kernels.py index 4ecfe8a..8ce9090 100644 --- a/benchmarks/test_kernels.py +++ b/benchmarks/test_kernels.py @@ -13,8 +13,8 @@ eikonal_solve, synthesise_pulses, ) -from rupture_generator.rupture.source import segment_source -from rupture_generator.sampling import Grid +from rupture_generator.rupture.source import MAI_BEROZA +from rupture_generator.sampling import Covariance, Grid, VonKarman from rupture_generator.sampling.field import _embed SPACING_KM = 0.1 @@ -34,7 +34,7 @@ def grid(): @pytest.fixture(scope="module", params=MAGNITUDES, ids=lambda m: f"Mw{m}") def covariance(request): - return segment_source(request.param, 90.0).covariance + return Covariance(VonKarman(), MAI_BEROZA.lengths_km(request.param)) def test_embed(benchmark, grid, covariance): diff --git a/examples/alpine_hope.toml b/examples/alpine_hope.toml index b01c5cb..7bfac88 100644 --- a/examples/alpine_hope.toml +++ b/examples/alpine_hope.toml @@ -143,12 +143,6 @@ Kelly = 180.0 "Te Rapa 1" = 160.0 "Te Rapa 2" = 135.0 -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333333333333333 -dip_offset = 1.5 - [propagation] type = "predetermined" diff --git a/examples/beavan.toml b/examples/beavan.toml index a028207..e3bba26 100644 --- a/examples/beavan.toml +++ b/examples/beavan.toml @@ -156,12 +156,6 @@ Charing_Cross = 66.5 Fault_NNW = 6.23 Fault_NE = 15.46 -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333333333333333 -dip_offset = 1.5 - [propagation] type = "predetermined" diff --git a/examples/colombia.toml b/examples/colombia.toml index 369799c..8860d39 100644 --- a/examples/colombia.toml +++ b/examples/colombia.toml @@ -102,12 +102,6 @@ type = "finite" magnitude = 7.4 rake_deg = 50.54 -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333 -dip_offset = 1.5 - [slip] coefficient_of_variation = 0.75 side_taper = 0.02 diff --git a/examples/crustal.toml b/examples/crustal.toml index e047996..63004b0 100644 --- a/examples/crustal.toml +++ b/examples/crustal.toml @@ -40,12 +40,6 @@ type = "finite" magnitude = 6.2 rake_deg = 175.0 -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333333333333333 -dip_offset = 1.5 - [slip] coefficient_of_variation = 0.9 side_taper = 0.02 diff --git a/examples/landers.toml b/examples/landers.toml index 1fe4580..a936c3d 100644 --- a/examples/landers.toml +++ b/examples/landers.toml @@ -68,12 +68,6 @@ type = "per_fault" "Homestead Valley" = 180.0 "Camp Rock - Emerson" = 180.0 -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333333333333333 -dip_offset = 1.5 - [propagation] type = "predetermined" diff --git a/examples/northridge.toml b/examples/northridge.toml index 029f03c..05de8bc 100644 --- a/examples/northridge.toml +++ b/examples/northridge.toml @@ -63,12 +63,6 @@ type = "finite" magnitude = 6.8 rake_deg = 102.2 -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333333333333333 -dip_offset = 1.5 - [slip] coefficient_of_variation = 0.9 side_taper = 0.02 diff --git a/examples/two_faults_segments.toml b/examples/two_faults_segments.toml index efb8412..876ffa1 100644 --- a/examples/two_faults_segments.toml +++ b/examples/two_faults_segments.toml @@ -1,7 +1,6 @@ # The two faults of two_faults.toml, with settings that differ on one segment. # A [segments.] table starts from the global table of the same name and -# replaces the keys it gives; [segments..correlation] does the same over the -# source's correlation table. +# replaces the keys it gives. A list such as a noise replaces the global one whole. # # The parameters are not tuned to match any observations: this is not a past earthquake. @@ -41,5 +40,5 @@ coefficient_of_variation = 0.6 [segments.next.speed] velocity_fraction = 0.7 -[segments.next.correlation] -strike_offset = 2.2 +[[segments.next.slip.noise]] +relation = { strike_offset = 2.2 } diff --git a/examples/venezuela.toml b/examples/venezuela.toml index 314fdbd..316e6a0 100644 --- a/examples/venezuela.toml +++ b/examples/venezuela.toml @@ -82,12 +82,6 @@ type = "predetermined" [propagation.parents] "Segment 2" = "Segment 1" -[source.correlation] -strike_exponent = 0.5 -strike_offset = 2.5 -dip_exponent = 0.3333333333333333 -dip_offset = 1.5 - [slip] coefficient_of_variation = 0.9 side_taper = 0.02 diff --git a/rupture_generator/__init__.py b/rupture_generator/__init__.py index 1817bdc..6e2fe4d 100644 --- a/rupture_generator/__init__.py +++ b/rupture_generator/__init__.py @@ -14,10 +14,9 @@ later. :mod:`rupture_generator.config` builds the same inputs from a TOML file, and the -command line runs that file. :func:`rupture_generator.sampling.sampler` keeps the -last 128 covariance embeddings, two per segment and one more per extra slip noise -family, so drawing a fault system of up to 64 segments again reuses every one. -``sampler.cache_clear()`` releases them. +command line runs that file. :func:`embed` keeps the last 128 noise embeddings, two +or three per segment, so drawing a fault system again reuses every one. +``embed.cache_clear()`` releases them. """ from importlib.metadata import version @@ -45,11 +44,13 @@ ramp_field, ) from rupture_generator.rupture.models import ( + STATIONARY, CorrelatedPerturbation, CorrelatedRise, Drawn, DrawnSlip, GaussianRake, + NoiseFamily, Perturbation, PerturbationModel, PulseSettings, @@ -61,12 +62,14 @@ SlipField, SlipModel, SpeedModel, + embed, gaussian_anamorphosis, ) from rupture_generator.rupture.propagator import JumpModel, likeliest_path, sample_path from rupture_generator.rupture.realisation import Hypocentre, Realisation from rupture_generator.rupture.source import ( MAI_BEROZA, + SUZUKI_2022, CorrelationRelation, SegmentSource, magnitude_from_moment, @@ -80,6 +83,8 @@ __all__ = [ "MAI_BEROZA", + "STATIONARY", + "SUZUKI_2022", "CorrelatedPerturbation", "CorrelatedRise", "CorrelationRelation", @@ -93,6 +98,7 @@ "JumpSeed", "Layers", "Medium", + "NoiseFamily", "Perturbation", "PerturbationModel", "PulseSettings", @@ -111,6 +117,7 @@ "SpatialField", "SpeedModel", "constant_field", + "embed", "first_arrival", "gaussian_anamorphosis", "generate", diff --git a/rupture_generator/config.py b/rupture_generator/config.py index 0436462..bd4e996 100644 --- a/rupture_generator/config.py +++ b/rupture_generator/config.py @@ -12,10 +12,9 @@ library. A ``[segments.]`` table overrides the global settings on one segment. Its -``slip``, ``rise``, ``rake``, ``speed`` and ``perturbation`` tables, and its -``correlation`` table for the source's, each start from the global table and replace -the keys they give. A depth profile, one key of its table, replaces the global profile -whole. :func:`load` does this merge. A loaded config stores each segment's tables in +``slip``, ``rise``, ``rake``, ``speed`` and ``perturbation`` tables each start from the +global table and replace the keys they give. A depth profile or a noise list, one key +of its table, replaces the global one whole. :func:`load` does this merge. A loaded config stores each segment's tables in full. :func:`load` reads and checks a file, and :func:`build` turns it into a @@ -74,7 +73,6 @@ from rupture_generator.rupture.source import ( CorrelationRelation, SegmentSource, - magnitude_from_moment, moment_from_magnitude, segment_source, split_moment, @@ -90,7 +88,7 @@ class _Strict(BaseConfig): @dataclasses.dataclass(frozen=True) class Correlation(CorrelationRelation): - """A source's ``correlation`` table, Mai & Beroza (2002) when absent.""" + """A noise family's ``relation`` table, Mai & Beroza (2002) when absent.""" Config = _Strict @@ -180,8 +178,8 @@ class Profile(_DepthProfile): """ type: Literal["profile"] = "profile" - depth_km: list[float] - values: list[float] + depth_km: tuple[float, ...] + values: tuple[float, ...] Config = _Strict def field(self) -> SpatialField: @@ -243,34 +241,37 @@ def field(self) -> SpatialField: @dataclasses.dataclass(frozen=True) -class Noise(NoiseFamily): - """``[[slip.noise]]``: :class:`~rupture_generator.rupture.models.NoiseFamily`. +class Family(NoiseFamily): + """``[[.noise]]``: :class:`~rupture_generator.rupture.models.NoiseFamily`. Parameters ---------- + relation : Correlation + ``relation``, the family's correlation lengths. weight : DepthProfile - ``[slip.noise.weight]``, how much of slip's latent the family is. - relation : Correlation, optional - ``[slip.noise.relation]``, the family's correlation lengths. + ``weight``, how much of the field the family is. """ - weight: DepthProfile - relation: Correlation | None = None + relation: Correlation = dataclasses.field(default_factory=Correlation) + weight: DepthProfile = UNMODIFIED Config = _Strict +STATIONARY = (Family(),) +"""The noise a model draws from when its table gives none.""" + + @dataclasses.dataclass(frozen=True) class Slip(DrawnSlip): """``[slip]``: :class:`~rupture_generator.rupture.models.DrawnSlip`. Parameters ---------- - noise : tuple of Noise - ``[[slip.noise]]``, the families slip's latent blends. A segment's list - replaces the global one whole. + noise : tuple of Family + ``[[slip.noise]]``, slip's noise. """ - noise: tuple[Noise, ...] = () + noise: tuple[Family, ...] = STATIONARY Config = _Strict @@ -294,8 +295,15 @@ class Rise(CorrelatedRise): @dataclasses.dataclass(frozen=True) class Rake(GaussianRake): - """``[rake]``: :class:`~rupture_generator.rupture.models.GaussianRake`.""" + """``[rake]``: :class:`~rupture_generator.rupture.models.GaussianRake`. + Parameters + ---------- + noise : tuple of Family + ``[[rake.noise]]``, rake's noise. + """ + + noise: tuple[Family, ...] = STATIONARY Config = _Strict @@ -315,8 +323,15 @@ class Speed(ShearFractionSpeed): @dataclasses.dataclass(frozen=True) class Perturbation(CorrelatedPerturbation): - """``[perturbation]``: :class:`~rupture_generator.rupture.models.CorrelatedPerturbation`.""" + """``[perturbation]``: :class:`~rupture_generator.rupture.models.CorrelatedPerturbation`. + Parameters + ---------- + noise : tuple of Family, optional + ``[[perturbation.noise]]``, the displacement's own noise. Slip's when absent. + """ + + noise: tuple[Family, ...] | None = None Config = _Strict @@ -460,14 +475,11 @@ class PerFault: Each segment's moment magnitude, by name. rakes : dict of str to float Each segment's mean rake in degrees, by name. - correlation : Correlation - How the slip correlation lengths follow each segment's magnitude. """ type: Literal["per_fault"] = "per_fault" magnitudes: dict[str, float] rakes: dict[str, float] - correlation: Correlation = dataclasses.field(default_factory=Correlation) Config = _Strict @@ -486,14 +498,11 @@ class Finite: The event's moment magnitude. rake_deg : float The mean rake in degrees, on every segment. - correlation : Correlation - How the slip correlation lengths follow each segment's magnitude. """ type: Literal["finite"] = "finite" magnitude: float rake_deg: float - correlation: Correlation = dataclasses.field(default_factory=Correlation) Config = _Strict @@ -582,8 +591,6 @@ class SegmentConfig: ``[segments..speed]``. perturbation : Perturbation, optional ``[segments..perturbation]``. - correlation : Correlation, optional - ``[segments..correlation]``, over the source's. """ slip: Slip | None = None @@ -591,7 +598,6 @@ class SegmentConfig: rake: Rake | None = None speed: Speed | None = None perturbation: Perturbation | None = None - correlation: Correlation | None = None Config = _Strict @@ -706,7 +712,6 @@ def _merged(document: dict[str, Any], overrides: object) -> dict[str, Any]: if not isinstance(overrides, dict): raise RuptureGeneratorError(f"{overrides!r} is not a table") shared = {key: _table(document.get(key)) for key in overrides} - shared["correlation"] = _table(_table(document.get("source")).get("correlation")) return { key: {**shared[key], **table} if isinstance(table, dict) else table for key, table in overrides.items() @@ -826,14 +831,8 @@ def _sources( source: PerFault | Finite, realisation: Realisation[Geometry], medium: Medium, - segments: dict[str, SegmentConfig], ) -> dict[str, SegmentSource]: names = set(realisation) - correlation = {name: source.correlation for name in realisation} | { - name: segment.correlation - for name, segment in segments.items() - if segment.correlation is not None - } match source: case PerFault(): for what, given in ( @@ -846,9 +845,7 @@ def _sources( f"{sorted(names - set(given))}, unknown {sorted(set(given) - names)}" ) return { - name: segment_source( - source.magnitudes[name], source.rakes[name], correlation[name] - ) + name: segment_source(source.magnitudes[name], source.rakes[name]) for name in realisation } case Finite(): @@ -856,11 +853,7 @@ def _sources( moment_from_magnitude(source.magnitude), realisation, medium ) return { - name: segment_source( - magnitude_from_moment(moment_nm), - source.rake_deg, - correlation[name], - ) + name: SegmentSource(moment_nm, source.rake_deg) for name, moment_nm in shares.items() } @@ -973,7 +966,7 @@ def build(config: RuptureConfig) -> Scenario: config.propagation, realisation, hypocentre, config.seed ), medium=medium, - sources=_sources(config.source, realisation, medium, config.segments), + sources=_sources(config.source, realisation, medium), model=_models(config, realisation), jump_model=(Jump() if config.propagation is None else config.propagation.model), seed=config.seed, diff --git a/rupture_generator/rupture/models.py b/rupture_generator/rupture/models.py index 74a947c..8e37649 100644 --- a/rupture_generator/rupture/models.py +++ b/rupture_generator/rupture/models.py @@ -12,7 +12,12 @@ Slip's latent is the one dependency between models. A model that correlates with slip mixes its own draw into :attr:`SlipField.latent`. Its draw comes from -:attr:`SlipField.latent_covariance`. The mixed field then has the latent's covariance. +:attr:`SlipField.sampler`. The mixed field then has the latent's covariance. + +Every model that draws noise names it: a :data:`Noise`, one or more +:class:`NoiseFamily` crossfaded by position, which :func:`embed` embeds on a chart. +Graves & Pitarka's hybrid slip is two families, long correlation lengths near the +surface and short ones at depth. Correlations with slip are between the latents: Gaussian-copula correlations, which mean the same whatever slip's marginal is, a drawn slip's or a given one's. :func:`gaussian_anamorphosis` gives any slip such a latent. @@ -22,8 +27,9 @@ """ import dataclasses +import functools from collections.abc import Callable, Mapping -from typing import Literal, Protocol +from typing import Protocol import numpy as np import scipy as sp @@ -32,19 +38,22 @@ from rupture_generator.geometry import CellArray, CellMask, Geometry from rupture_generator.rupture.medium import SpatialField, constant_field, rigidity_pa from rupture_generator.rupture.source import ( + MAI_BEROZA, CorrelationRelation, SegmentSource, - magnitude_from_moment, ) from rupture_generator.sampling import ( - Blend, + HURST, + NORMAL, Covariance, Grid, Marginal, MarginalFamily, PreCorrected, + Sampler, VonKarman, mix, + sampler, standardise, ) @@ -85,7 +94,7 @@ class SegmentContext: geometry : Geometry The segment's chart. source : SegmentSource - The segment's moment, rake and slip covariance. + The segment's moment and rake. shear_speed_km_s : CellArray The medium's shear speed at the cell centres, in kilometres per second. density_g_cm3 : CellArray @@ -119,12 +128,6 @@ def moment_nm(self, slip_m: CellArray) -> float: np.sum(self.rigidity_pa * self.geometry.areas_km2 * M2_PER_KM2 * slip_m) ) - @property - def grid(self) -> Grid: - """Grid: The chart's cells as a sampling grid, spacing down dip first.""" - strike_km, dip_km = self.geometry.spacing_km - return Grid(self.geometry.cells, (dip_km, strike_km)) - @dataclasses.dataclass(frozen=True, eq=False) class SlipField: @@ -136,9 +139,8 @@ class SlipField: Slip, in metres. The generator zeroes it off the fault. latent : CellArray A standard normal field that slip is monotone in. - latent_covariance : Covariance or Blend - The covariance the model drew ``latent`` from, or the nearest one for a slip it - didn't draw. + sampler : Sampler + What the model drew ``latent`` from, or the nearest for a slip it didn't draw. A model that mixes its own draw into ``latent`` draws from this. statistics : Mapping of str to float Numbers the model reports about its draw. @@ -146,7 +148,7 @@ class SlipField: slip_m: CellArray latent: CellArray - latent_covariance: Covariance | Blend + sampler: Sampler statistics: Mapping[str, float] = dataclasses.field(default_factory=dict) @@ -398,6 +400,22 @@ def read_at_centres( return check_cells(geometry, name, field(geometry.centres), accept, wanted) +def unit_interval(values: CellArray) -> np.ndarray: + """Whether each value lies in ``[0, 1]``. + + Parameters + ---------- + values : CellArray + The values. + + Returns + ------- + numpy.ndarray + One boolean per value. + """ + return (values >= 0.0) & (values <= 1.0) + + def positive(values: CellArray) -> np.ndarray: """Whether each value is finite and greater than zero. @@ -513,12 +531,83 @@ def _positive_marginal( return Marginal(family, spread, maximum) -def _independent( - context: SegmentContext, slip: SlipField, rng: np.random.Generator -) -> CellArray: - """A standard normal field independent of slip's latent, with its covariance.""" - independent, _ = slip.latent_covariance.sampler(context.grid).draw(rng) - return independent +# ----------------------------------------------------------------------- noise + + +@dataclasses.dataclass(frozen=True) +class NoiseFamily: + """One stationary von Karman noise, and how much of a field it is where. + + Attributes + ---------- + relation : CorrelationRelation + The correlation lengths at the segment's magnitude. + hurst : float + The Hurst exponent. + weight : SpatialField + How much of the field this family is at each position, in ``[0, 1]``. + """ + + relation: CorrelationRelation = MAI_BEROZA + hurst: float = HURST + weight: SpatialField = UNMODIFIED + + +type Noise = tuple[NoiseFamily, ...] +"""Noise families crossfaded by position, each sharing one white noise.""" + +STATIONARY: Noise = (NoiseFamily(),) +"""Mai & Beroza (2002)'s lengths and Hurst exponent everywhere.""" + + +@functools.lru_cache(maxsize=128) +def embed( + noise: Noise, geometry: Geometry, magnitude: float, marginal: Marginal = NORMAL +) -> Sampler: + """Embed a noise on a chart, or return the embedding already made. + + Every family is pre-corrected for ``marginal``, so the field it gives after the + marginal has the family's correlation. The cache keeps the last 128, so drawing the + same fault system again embeds nothing; ``embed.cache_clear()`` releases them. + + Parameters + ---------- + noise : Noise + The families. + geometry : Geometry + The chart. + magnitude : float + The segment's moment magnitude, which sets each family's lengths. + marginal : Marginal, optional + The marginal the drawn field takes. Normal by default. + + Returns + ------- + Sampler + The embedding, ready to draw from. + + Raises + ------ + RuptureGeneratorError + If a weight reads out of ``[0, 1]``, or a family doesn't embed on the chart. + """ + strike_km, dip_km = geometry.spacing_km + return sampler( + Grid(geometry.cells, (dip_km, strike_km)), + [ + Covariance( + PreCorrected(VonKarman(family.hurst), marginal), + family.relation.lengths_km(magnitude), + ) + for family in noise + ], + [ + read_at_centres( + geometry, "a noise weight", family.weight, unit_interval, "[0, 1]" + ) + for family in noise + ], + ) # ------------------------------------------------------------------------ slip @@ -537,58 +626,12 @@ def _reach(mask: CellMask, axis: int, *, reverse: bool) -> np.ndarray: return np.flip(counted, axis) if reverse else counted -@dataclasses.dataclass(frozen=True) -class NoiseFamily: - """One stationary noise slip draws from, and where on the fault it applies. - - Several families make Graves & Pitarka's hybrid slip: long correlation lengths near - the surface, short ones at depth, and a crossfade between them. - - Attributes - ---------- - weight : SpatialField - How much of slip's latent this family is at each position, in ``[0, 1]``. - relation : CorrelationRelation, optional - The family's correlation lengths at the segment's magnitude. The source's - lengths when absent. - hurst : float, optional - The family's von Karman Hurst exponent. The source's correlation function when - absent. - """ - - weight: SpatialField - relation: CorrelationRelation | None = None - hurst: float | None = None - - def covariance(self, source: SegmentSource) -> Covariance: - """The family's covariance on a segment. - - Parameters - ---------- - source : SegmentSource - The segment's source, whose moment sets the magnitude. - - Returns - ------- - Covariance - The family's correlation function and lengths. - """ - lengths = source.covariance.lengths_km - if self.relation is not None: - lengths = self.relation.lengths_km(magnitude_from_moment(source.moment_nm)) - correlation = source.covariance.correlation - if self.hurst is not None: - correlation = VonKarman(self.hurst) - return Covariance(correlation, lengths) - - @dataclasses.dataclass(frozen=True) class DrawnSlip: - """Slip drawn from the source's covariance, with an exact marginal and tapered edges. + """Slip drawn from a noise, with an exact marginal and tapered edges. - The latent field's covariance is pre-corrected for the marginal, so the slip - itself has the source's correlation lengths. The model scales slip to the source's - moment. + The noise is pre-corrected for the marginal, so the slip itself has the noise's + correlation lengths. The model scales slip to the source's moment. Attributes ---------- @@ -607,9 +650,8 @@ class DrawnSlip: bottom_taper : float The fraction of the fault's width from the bottom edge over which slip ramps to zero. - noise : tuple of NoiseFamily - The families slip's latent blends, each pre-corrected for the marginal. The - source's covariance alone when empty. + noise : Noise + The noise slip's latent comes from, and rise time's independent part. """ family: MarginalFamily = "truncated_exponential" @@ -618,7 +660,7 @@ class DrawnSlip: side_taper: float = 0.02 top_taper: float = 0.0 bottom_taper: float = 0.0 - noise: tuple[NoiseFamily, ...] = () + noise: Noise = STATIONARY def __post_init__(self) -> None: """Refuse a taper outside ``[0, 0.5]``, or a spread outside the family's range.""" @@ -659,33 +701,6 @@ def _taper(self, field: CellArray, occupied: CellMask) -> CellArray: ramp *= np.minimum(_reach(occupied, axis, reverse=reverse) / cells, 1.0) return field * ramp * occupied - def _latent_covariance(self, context: SegmentContext) -> Covariance | Blend: - """The covariance to draw the latent from, pre-corrected for the marginal.""" - - def corrected(covariance: Covariance) -> Covariance: - return Covariance( - PreCorrected(covariance.correlation, self.marginal), - covariance.lengths_km, - ) - - if not self.noise: - return corrected(context.source.covariance) - return Blend( - tuple( - corrected(family.covariance(context.source)) for family in self.noise - ), - tuple( - read_at_centres( - context.geometry, - "a slip noise family's weight", - family.weight, - lambda weight: (weight >= 0.0) & (weight <= 1.0), - "weights in [0, 1]", - ) - for family in self.noise - ), - ) - def __call__(self, context: SegmentContext, rng: np.random.Generator) -> SlipField: """Draw slip on a segment. @@ -704,12 +719,12 @@ def __call__(self, context: SegmentContext, rng: np.random.Generator) -> SlipFie Raises ------ RuptureGeneratorError - If the covariance doesn't embed on this chart, a noise family's weight - reads out of range, or the taper leaves no slip anywhere. + If the noise doesn't embed on this chart, or the taper leaves no slip + anywhere. """ geometry = context.geometry - latent_covariance = self._latent_covariance(context) - latent, _ = latent_covariance.sampler(context.grid).draw(rng) + noise = embed(self.noise, geometry, context.source.magnitude, self.marginal) + latent, _ = noise.draw(rng) pattern = self._taper(self.marginal.apply(latent), geometry.occupied) total = context.moment_nm(pattern) if not total > 0.0: @@ -720,7 +735,7 @@ def __call__(self, context: SegmentContext, rng: np.random.Generator) -> SlipFie return SlipField( slip_m=(context.source.moment_nm / total) * pattern, latent=latent, - latent_covariance=latent_covariance, + sampler=noise, ) @@ -827,7 +842,7 @@ def __call__( geometry, "the rise-time slip weight", self.slip_weight, - lambda weight: (weight >= 0.0) & (weight <= 1.0), + unit_interval, "weights in [0, 1]", ) correction = geometric_correction(geometry, context.source.rake_deg) @@ -835,7 +850,7 @@ def __call__( rho = self.correlation latent = mix( (slip_weight * rho + (1.0 - slip_weight), slip.latent), - (slip_weight * np.sqrt(1.0 - rho * rho), _independent(context, slip, rng)), + (slip_weight * np.sqrt(1.0 - rho * rho), slip.sampler.draw(rng)[0]), ) pattern = self.marginal.apply(latent) ** self.slip_exponent * factor return Drawn( @@ -852,16 +867,19 @@ class GaussianRake: """Rake that strays from the source's by ``sigma * Z`` degrees. By design, rake doesn't depend on slip. A patch with more slip has no reason to - move in a different direction. The marginal is normal, so rake keeps the source's + move in a different direction. The marginal is normal, so rake keeps the noise's correlation lengths exactly, with no pre-correction. Attributes ---------- sigma_deg : float The standard deviation of the rake about the segment's, in degrees. + noise : Noise + The noise rake strays by. Stationary by default, as in genslip. """ sigma_deg: float = 15.0 + noise: Noise = STATIONARY def __post_init__(self) -> None: """Refuse a negative spread.""" @@ -885,7 +903,9 @@ def __call__(self, context: SegmentContext, rng: np.random.Generator) -> Drawn: Drawn Rake in degrees. """ - latent, _ = context.source.covariance.sampler(context.grid).draw(rng) + latent, _ = embed(self.noise, context.geometry, context.source.magnitude).draw( + rng + ) return Drawn(context.source.rake_deg + self.sigma_deg * standardise(latent)) @@ -987,16 +1007,16 @@ class CorrelatedPerturbation: blend_sigma : float The width of the zone over which the displacement grows in from the seed, in units of the spread. - independent : {"slip", "source"} - The covariance of the displacement's own draw: slip's latent's, which follows - any slip noise families, or the source's, which genslip v5.6.3 keeps. + noise : Noise, optional + The noise of the displacement's own draw. Slip's when absent; genslip v5.6.3 + draws it from :data:`STATIONARY` while slip's is hybrid. """ offset_s: float = 0.1 coefficient: float = 0.5 correlation: float = 0.8 blend_sigma: float = 4.0 - independent: Literal["slip", "source"] = "slip" + noise: Noise | None = None def __post_init__(self) -> None: """Refuse a negative spread, or a blend of no width.""" @@ -1039,10 +1059,10 @@ def __call__( # genslip's tsfac is negative and this package reads it as a magnitude, so the # displacement subtracts: with a positive correlation, high slip ruptures early. rho = self.correlation - if self.independent == "slip": - independent = _independent(context, slip, rng) - else: - independent, _ = context.source.covariance.sampler(context.grid).draw(rng) + noise = slip.sampler + if self.noise is not None: + noise = embed(self.noise, context.geometry, context.source.magnitude) + independent, _ = noise.draw(rng) displacement_s = -scale_s * standardise( mix((rho, slip.latent), (np.sqrt(1.0 - rho * rho), independent)) ) diff --git a/rupture_generator/rupture/source.py b/rupture_generator/rupture/source.py index 6ed7905..8232dee 100644 --- a/rupture_generator/rupture/source.py +++ b/rupture_generator/rupture/source.py @@ -1,9 +1,10 @@ -"""Each segment's moment, rake and slip-patch size, derived from its magnitude. +"""Each segment's moment and rake, and the correlation lengths its magnitude implies. A :class:`SegmentSource` is the part of a rupture that differs from segment to segment. A :class:`~rupture_generator.rupture.models.RuptureModel` draws the fields from -it. The functions here -derive sources from magnitudes, which is how a hazard model states them. +it. The functions here derive sources from magnitudes, which is how a hazard model +states them. A :class:`CorrelationRelation` gives a model's noise its lengths at a +segment's magnitude. References ---------- @@ -26,7 +27,6 @@ from rupture_generator.geometry import Geometry from rupture_generator.rupture.medium import Medium from rupture_generator.rupture.realisation import Realisation -from rupture_generator.sampling import Correlation, Covariance, VonKarman DYNE_CM_PER_NM = 1.0e7 @@ -110,34 +110,27 @@ def lengths_km(self, magnitude: float) -> tuple[float, float]: tuple of float ``(down dip, along strike)`` in kilometres, the chart's axis order. """ + dip_magnitude = min(magnitude, self.dip_saturation) + strike_magnitude = min(magnitude, self.strike_saturation) return ( - 10.0 - ** ( - self.dip_exponent * min(magnitude, self.dip_saturation) - - self.dip_offset - ), - 10.0 - ** ( - self.strike_exponent * min(magnitude, self.strike_saturation) - - self.strike_offset - ), + 10.0 ** (self.dip_exponent * dip_magnitude - self.dip_offset), + 10.0 ** (self.strike_exponent * strike_magnitude - self.strike_offset), ) -VON_KARMAN = VonKarman() -"""Mai & Beroza (2002)'s correlation function, at their median Hurst exponent.""" - MAI_BEROZA = CorrelationRelation() """Mai & Beroza (2002): ``10^(Mw/2 - 2.5)`` along strike, ``10^(Mw/3 - 1.5)`` down dip.""" -SUZUKI_2022 = CorrelationRelation(0.5, 1.67, 0.5, 1.69, dip_saturation=6.3) +SUZUKI_2022 = CorrelationRelation( + strike_offset=1.67, dip_exponent=0.5, dip_offset=1.69, dip_saturation=6.3 +) """Suzuki et al. (2022): ``10^(Mw/2 - 1.67)`` along strike, ``10^(min(Mw, 6.3)/2 - 1.69)`` down dip. Smoother than Mai & Beroza, the shallow slip of Graves & Pitarka's hybrid.""" @dataclasses.dataclass(frozen=True) class SegmentSource: - """One segment's moment, mean rake, and slip covariance. + """One segment's moment and mean rake. Parameters ---------- @@ -145,14 +138,10 @@ class SegmentSource: The segment's seismic moment in newton-metres. Positive. rake_deg : float The segment's mean rake in degrees. - covariance : Covariance - The correlation of the slip *pattern*, with its lengths in the chart's axis - order, ``(down dip, along strike)``. """ moment_nm: float rake_deg: float - covariance: Covariance def __post_init__(self) -> None: """Refuse a moment of zero or less.""" @@ -161,13 +150,13 @@ def __post_init__(self) -> None: f"a segment's moment must be positive, got {self.moment_nm}" ) + @property + def magnitude(self) -> float: + """float: The segment's moment magnitude.""" + return magnitude_from_moment(self.moment_nm) + -def segment_source( - magnitude: float, - rake_deg: float, - relation: CorrelationRelation = MAI_BEROZA, - correlation: Correlation = VON_KARMAN, -) -> SegmentSource: +def segment_source(magnitude: float, rake_deg: float) -> SegmentSource: """Build a segment's source from its magnitude. Parameters @@ -176,23 +165,13 @@ def segment_source( The segment's moment magnitude. rake_deg : float The segment's mean rake in degrees. - relation : CorrelationRelation, optional - How the correlation lengths follow the magnitude. Mai & Beroza (2002) by - default. - correlation : Correlation, optional - The correlation function. Von Karman at Mai & Beroza's median Hurst exponent - by default. Returns ------- SegmentSource - The moment the magnitude gives, and the correlation lengths it implies. + The moment the magnitude gives, and the rake. """ - return SegmentSource( - moment_nm=moment_from_magnitude(magnitude), - rake_deg=rake_deg, - covariance=Covariance(correlation, relation.lengths_km(magnitude)), - ) + return SegmentSource(moment_from_magnitude(magnitude), rake_deg) def split_moment( @@ -233,7 +212,6 @@ def split_moment( __all__ = [ "MAI_BEROZA", "SUZUKI_2022", - "VON_KARMAN", "CorrelationRelation", "SegmentSource", "magnitude_from_moment", diff --git a/rupture_generator/sampling/__init__.py b/rupture_generator/sampling/__init__.py index 7bd418f..39f0d81 100644 --- a/rupture_generator/sampling/__init__.py +++ b/rupture_generator/sampling/__init__.py @@ -1,8 +1,6 @@ """Random fields on the fault: Gaussian machinery, correlation functions, NORTA.""" from rupture_generator.sampling.field import ( - Blend, - BlendedSampler, Correlation, Covariance, FieldArray, @@ -24,8 +22,6 @@ __all__ = [ "HURST", "NORMAL", - "Blend", - "BlendedSampler", "Correlation", "Covariance", "FieldArray", diff --git a/rupture_generator/sampling/field.py b/rupture_generator/sampling/field.py index 051e0b1..add0dde 100644 --- a/rupture_generator/sampling/field.py +++ b/rupture_generator/sampling/field.py @@ -2,16 +2,15 @@ This module embeds any correlation function and doesn't assume anything about what the caller does with the field. A :class:`Covariance` is a :class:`Correlation` and the -lengths that scale its argument. :func:`sampler` embeds one on a :class:`Grid` once -and draws from it as often as asked. A :class:`Blend` crossfades several covariances -cell by cell, for a field whose correlation lengths change across the grid. Every tuple -is ``(axis 0, axis 1)``. +lengths that scale its argument. :func:`sampler` embeds one or more on a :class:`Grid` +once and draws from them as often as asked. Several covariances crossfade cell by cell, +for a field whose correlation lengths change across the grid. Every tuple is +``(axis 0, axis 1)``. """ import dataclasses -import functools import math -from collections.abc import Callable +from collections.abc import Callable, Sequence import numpy as np import scipy as sp @@ -44,11 +43,6 @@ MAXIMUM_EMBEDDING_CELLS = 1 << 27 """The largest padded grid to transform. """ -SAMPLER_CACHE_SIZE = 128 -"""How many embeddings :func:`sampler` keeps. Each segment takes two, and one more for -each family a :class:`Blend` adds, so a fault system of up to 64 segments, or 42 with -two families each, draws again without embedding anything.""" - @dataclasses.dataclass(frozen=True) class Grid: @@ -85,21 +79,6 @@ class Covariance: correlation: Correlation lengths_km: tuple[float, float] - def sampler(self, grid: Grid) -> "Sampler": - """Embed the covariance on ``grid``, through :func:`sampler`'s cache. - - Parameters - ---------- - grid : Grid - The grid to draw on. - - Returns - ------- - Sampler - The embedding, ready to draw from. - """ - return sampler(grid, self) - def _decay_length(correlation: Correlation, tolerance: float) -> float: """The lag, in correlation lengths, past which the correlation is below ``tolerance``.""" @@ -204,28 +183,36 @@ def _embed( @dataclasses.dataclass(frozen=True, eq=False) class Sampler: - """One covariance embedded on one grid, the expensive part, done once. + """Covariances embedded on one grid, the expensive part, done once. + + Every covariance is embedded at one padded shape and draws from one white noise. + Each cell loads them at its own weights, so the field's correlation lengths can + change across the grid. Attributes ---------- shape : tuple of int The grid's cell counts, the shape every draw comes back in. - amplitudes : FieldArray - The square-rooted eigenvalues of the embedding, one quadrant of the padded - grid. padded : tuple of int The padded grid's shape. + amplitudes : numpy.ndarray + The square-rooted eigenvalues of each embedding, one quadrant of the padded + grid each, shaped ``(covariances, *quadrant)``. + loadings : numpy.ndarray + Each covariance's loading at each cell, shaped ``(covariances, *shape)``. """ shape: tuple[int, int] - amplitudes: FieldArray padded: tuple[int, int] + amplitudes: np.ndarray + loadings: np.ndarray def draw(self, rng: np.random.Generator) -> tuple[FieldArray, FieldArray]: """Draw a pair of independent fields on the grid, standard normal at every cell. - One complex transform gives both. Its real and imaginary parts are independent - fields with the embedded covariance (Dietrich & Newsam 1997). + One complex transform per covariance gives both. Its real and imaginary parts + are independent fields (Dietrich & Newsam 1997), and every covariance filters + the same noise, from one seed. Parameters ---------- @@ -237,150 +224,78 @@ def draw(self, rng: np.random.Generator) -> tuple[FieldArray, FieldArray]: tuple of FieldArray The real and imaginary fields, each shaped like the grid. """ - return circulant_draw(self.amplitudes, self.padded, self.shape, _seed(rng)) - - -def _seed(rng: np.random.Generator) -> int: - """Take one draw's seed from ``rng``.""" - return int(rng.integers(0, np.iinfo(np.uint64).max, dtype=np.uint64)) + seed = int(rng.integers(0, np.iinfo(np.uint64).max, dtype=np.uint64)) + pair = np.zeros((2, *self.shape)) + for amplitudes, loading in zip(self.amplitudes, self.loadings, strict=True): + pair += loading * np.asarray( + circulant_draw(amplitudes, self.padded, self.shape, seed) + ) + return pair[0], pair[1] -@functools.lru_cache(maxsize=SAMPLER_CACHE_SIZE) def sampler( - grid: Grid, covariance: Covariance, floor: tuple[int, int] = (0, 0) + grid: Grid, + covariances: Sequence[Covariance], + weights: Sequence[FieldArray] | None = None, ) -> Sampler: - """Embed ``covariance`` on ``grid``, or return the embedding already made. + """Embed covariances on ``grid`` at one padded shape, and weigh them cell by cell. + + A cell loads each covariance at ``sin(pi w / 2)`` of its weight ``w``, scaled to + unit variance from the covariances' exact cross-correlations, so the field is + standard normal at every cell however the weights overlap. Between two covariances + whose weights sum to one, that is an equal-power crossfade. Parameters ---------- grid : Grid The grid to draw on. - covariance : Covariance - The covariance to embed. - floor : tuple of int, optional - The least padded shape, so covariances that share noise share a padded grid. + covariances : sequence of Covariance + The covariances to embed. + weights : sequence of FieldArray, optional + Each covariance's weight at each cell, in ``[0, 1]``. Equal everywhere when + absent. Returns ------- Sampler - The embedding, ready to draw from. + The embeddings, ready to draw from. Raises ------ RuptureGeneratorError - If the covariance doesn't embed within the cell cap, or isn't positive - definite at any margin tried. - """ - return Sampler(grid.shape, *_embed(grid, covariance, floor)) - - -@dataclasses.dataclass(frozen=True, eq=False) -class Blend: - """Stationary covariances that filter one white noise, crossfaded cell by cell. - - It is the covariance of a field whose correlation lengths change across the grid, - such as slip with long lengths near the surface and short ones at depth. Every - component filters the same noise, so a feature changes shape through a crossfade - rather than fading out. A cell's loadings are ``sin(pi w / 2)``, scaled to unit - variance from the components' exact cross-correlations, so the field is standard - normal at every cell however the weights overlap. - - Attributes - ---------- - components : tuple of Covariance - The stationary covariances. - weights : tuple of FieldArray - Each component's weight at each cell, in ``[0, 1]``. Every cell needs a - positive weight. + If a weight is outside ``[0, 1]``, a cell has no positive weight, or a + covariance doesn't embed within the cell cap or isn't positive definite at any + margin tried. """ - - components: tuple[Covariance, ...] - weights: tuple[FieldArray, ...] - - def sampler(self, grid: Grid) -> "BlendedSampler": - """Embed every component on ``grid`` at one padded shape, and load them. - - Parameters - ---------- - grid : Grid - The grid to draw on, shaped like the weights. - - Returns - ------- - BlendedSampler - The embeddings and each one's loading, ready to draw from. - - Raises - ------ - RuptureGeneratorError - If a cell has no positive weight, or a component doesn't embed. - """ - gains = np.sin(0.5 * np.pi * np.asarray(self.weights)) - power = np.sum(np.square(gains), axis=0) - if not np.all(power > 0.0): - raise RuptureGeneratorError( - f"{int(np.sum(~(power > 0.0)))} cells have no positive weight in the blend" - ) - floor = (0, 0) - while True: - samplers = [ - sampler(grid, component, floor) for component in self.components - ] - widest = tuple(max(axis) for axis in zip(*(s.padded for s in samplers))) - if all(s.padded == widest for s in samplers): - break - floor = widest - weights = np.outer(*(_mirror_weights(extent) for extent in widest)) - amplitudes = np.array([s.amplitudes for s in samplers]) - cross = np.einsum("ikl,jkl,kl->ij", amplitudes, amplitudes, weights) - cross /= widest[0] * widest[1] - gains /= np.sqrt(power) - variance = np.einsum("ikl,ij,jkl->kl", gains, cross, gains) - return BlendedSampler(tuple(samplers), gains / np.sqrt(variance)) - - -@dataclasses.dataclass(frozen=True, eq=False) -class BlendedSampler: - """A :class:`Blend`'s components, embedded on one padded grid, and their loadings. - - Attributes - ---------- - samplers : tuple of Sampler - The components' embeddings, all at one padded shape. - loadings : numpy.ndarray - Each component's loading at each cell, shaped ``(components, *grid)``. - """ - - samplers: tuple[Sampler, ...] - loadings: np.ndarray - - def draw(self, rng: np.random.Generator) -> tuple[FieldArray, FieldArray]: - """Draw a pair of independent blended fields, standard normal at every cell. - - Every component draws from the one seed, so they filter the same noise. The - imaginary parts are independent of every real part, because the amplitudes - are even. - - Parameters - ---------- - rng : np.random.Generator - The generator the draw's seed comes from. - - Returns - ------- - tuple of FieldArray - The real and imaginary fields, each shaped like the grid. - """ - seed = _seed(rng) - real = np.zeros(self.loadings.shape[1:]) - imaginary = np.zeros(self.loadings.shape[1:]) - for loading, embedding in zip(self.loadings, self.samplers, strict=True): - parts = circulant_draw( - embedding.amplitudes, embedding.padded, embedding.shape, seed + shares = np.ones((len(covariances), *grid.shape)) + if weights is not None: + shares = np.asarray(weights, dtype=np.float64) + if not np.all((shares >= 0.0) & (shares <= 1.0)): + raise RuptureGeneratorError("a covariance's weight lies outside [0, 1]") + gains = np.sin(0.5 * np.pi * shares) + power = np.sum(np.square(gains), axis=0) + if not np.all(power > 0.0): + raise RuptureGeneratorError( + f"{int(np.sum(~(power > 0.0)))} cells weigh every covariance at zero" + ) + embedded = [_embed(grid, covariance) for covariance in covariances] + padded = tuple(np.max([shape for _, shape in embedded], axis=0).tolist()) + while any(shape != padded for _, shape in embedded): + embedded = [ + (amplitudes, shape) if shape == padded else _embed(grid, covariance, padded) + for (amplitudes, shape), covariance in zip( + embedded, covariances, strict=True ) - real += loading * parts[0] - imaginary += loading * parts[1] - return real, imaginary + ] + padded = tuple(np.max([shape for _, shape in embedded], axis=0).tolist()) + amplitudes = np.array([amplitudes for amplitudes, _ in embedded]) + amplitudes.setflags(write=False) + mirror = np.outer(*(_mirror_weights(extent) for extent in padded)) + cross = np.einsum("ikl,jkl,kl->ij", amplitudes, amplitudes, mirror) + cross /= padded[0] * padded[1] + gains /= np.sqrt(power) + variance = np.einsum("ikl,ij,jkl->kl", gains, cross, gains) + return Sampler(grid.shape, padded, amplitudes, gains / np.sqrt(variance)) def mix(*terms: tuple[float | FieldArray, FieldArray]) -> FieldArray: diff --git a/tests/test_blend.py b/tests/test_blend.py index a867f7f..56db45f 100644 --- a/tests/test_blend.py +++ b/tests/test_blend.py @@ -1,50 +1,32 @@ -"""Slip noise families: the blended sampler, the models that use it, the config.""" +"""Noise families: blended sampling, the models that draw noise, and the config.""" import dataclasses import numpy as np import pytest -from rupture_generator import RuptureGeneratorError, generate -from rupture_generator import config as rupture_config -from rupture_generator.rupture.models import ( +from rupture_generator import ( + STATIONARY, + SUZUKI_2022, CorrelatedPerturbation, DrawnSlip, + GaussianRake, NoiseFamily, + RuptureGeneratorError, RuptureModel, - SegmentContext, + generate, ) -from rupture_generator.rupture.source import MAI_BEROZA, SUZUKI_2022 -from rupture_generator.sampling import Blend, Covariance, Grid, VonKarman, sampler -from rupture_generator.sampling.field import _mirror_weights +from rupture_generator import config as rupture_config +from rupture_generator.rupture.source import MAI_BEROZA +from rupture_generator.sampling import Covariance, Grid, VonKarman, sampler from .conftest import EXAMPLE GRID = Grid((24, 60), (0.5, 0.5)) SHORT = Covariance(VonKarman(), (2.0, 3.0)) LONG = Covariance(VonKarman(), (4.0, 12.0)) - - -def _ramp(rows: int = GRID.shape[0]) -> np.ndarray: - """Weight of the deep family: 0 above row 8, 1 below row 16, linear between.""" - column = np.clip((np.arange(rows) - 8.0) / 8.0, 0.0, 1.0) - return np.repeat(column[:, None], GRID.shape[1], axis=1) - - -def _blend() -> Blend: - deep = _ramp() - return Blend((LONG, SHORT), (1.0 - deep, deep)) - - -def _cross(samplers) -> np.ndarray: - padded = samplers[0].padded - weights = np.outer(*(_mirror_weights(extent) for extent in padded)) - return np.array( - [ - [np.sum(weights * a.amplitudes * b.amplitudes) for b in samplers] - for a in samplers - ] - ) / (padded[0] * padded[1]) +DEEP = np.repeat(np.clip((np.arange(24) - 8.0) / 8.0, 0.0, 1.0)[:, None], 60, axis=1) +"""The deep family's weight: none above row 8, all below row 16, linear between.""" def test_saturation_stops_one_axis_growing(): @@ -56,52 +38,28 @@ def test_saturation_stops_one_axis_growing(): ) -def test_a_floor_only_ever_widens_the_embedding(): - alone = sampler(GRID, SHORT) - widened = sampler(GRID, SHORT, (alone.padded[0] * 2, alone.padded[1])) - assert widened.padded == (alone.padded[0] * 2, alone.padded[1]) - assert sampler(GRID, SHORT, (2, 2)).padded == alone.padded - - -def test_the_components_share_the_widest_padded_shape(): +def test_covariances_share_the_widest_padded_shape(): across = Covariance(VonKarman(), (12.0, 2.0)) - blended = Blend((across, SHORT, LONG), (np.full(GRID.shape, 1.0),) * 3).sampler( - GRID - ) - own = [sampler(GRID, c).padded for c in (across, SHORT, LONG)] - widest = tuple(max(axis) for axis in zip(*own)) - assert all(s.padded == widest for s in blended.samplers) - - -def test_every_cell_has_unit_variance_exactly(): - blended = _blend().sampler(GRID) - variance = np.einsum( - "ikl,ij,jkl->kl", blended.loadings, _cross(blended.samplers), blended.loadings - ) - np.testing.assert_allclose(variance, 1.0, rtol=1e-12) + own = [sampler(GRID, [c]).padded for c in (across, SHORT, LONG)] + widest = (max(p[0] for p in own), max(p[1] for p in own)) + assert widest not in own + assert sampler(GRID, [across, SHORT, LONG]).padded == widest def test_a_cell_of_one_family_is_that_familys_draw(): - blended = _blend().sampler(GRID) - real, imaginary = blended.draw(np.random.default_rng(3)) - # The same seed draws the same noise for every component. - seed = int( - np.random.default_rng(3).integers(0, np.iinfo(np.uint64).max, dtype=np.uint64) - ) - from rupture_generator._kernels import circulant_draw - - long_draw = circulant_draw( - blended.samplers[0].amplitudes, blended.samplers[0].padded, GRID.shape, seed - ) - scale = np.sqrt(_cross(blended.samplers)[0, 0]) - np.testing.assert_allclose(real[:8], long_draw[0][:8] / scale, rtol=1e-12) - np.testing.assert_allclose(imaginary[:8], long_draw[1][:8] / scale, rtol=1e-12) + blend = sampler(GRID, [LONG, SHORT], [1.0 - DEEP, DEEP]) + alone = sampler(GRID, [LONG]) + assert blend.padded == alone.padded + mixed, alone_draw = (s.draw(np.random.default_rng(3)) for s in (blend, alone)) + for ours, theirs in zip(mixed, alone_draw, strict=True): + np.testing.assert_allclose(ours[:8], theirs[:8], rtol=1e-12) + assert not np.allclose(ours[16:], theirs[16:]) def test_draws_are_standard_normal_and_the_pair_independent(): - blended = _blend().sampler(GRID) + blend = sampler(GRID, [LONG, SHORT], [1.0 - DEEP, DEEP]) rng = np.random.default_rng(11) - pairs = np.array([blended.draw(rng) for _ in range(600)]) + pairs = np.array([blend.draw(rng) for _ in range(600)]) real, imaginary = pairs[:, 0], pairs[:, 1] # Through the crossfade, rows 8 to 16, as well as either side of it. assert np.abs(real.var(axis=0) - 1.0).max() < 0.25 @@ -109,67 +67,66 @@ def test_draws_are_standard_normal_and_the_pair_independent(): assert abs(float(np.mean(real * imaginary))) < 0.02 -def test_a_cell_with_no_weight_is_refused(): - deep = _ramp() - with pytest.raises(RuptureGeneratorError, match="no positive weight"): - Blend((LONG, SHORT), (np.zeros(GRID.shape), deep)).sampler(GRID) +@pytest.mark.parametrize( + ("weights", "message"), + [ + ([np.zeros(GRID.shape), DEEP], "weigh every covariance at zero"), + ([1.0 - DEEP, 2.0 * DEEP], r"outside \[0, 1\]"), + ], +) +def test_bad_weights_are_refused(weights, message): + with pytest.raises(RuptureGeneratorError, match=message): + sampler(GRID, [LONG, SHORT], weights) # ----------------------------------------------------------------- the models -def _hybrid_slip() -> DrawnSlip: - def deep(positions): - return np.clip((positions[..., 2] - 4.0) / 4.0, 0.0, 1.0) +def _deep(positions): + return np.clip((positions[..., 2] - 4.0) / 4.0, 0.0, 1.0) - return DrawnSlip( - noise=( - NoiseFamily(lambda p: 1.0 - deep(p), SUZUKI_2022, hurst=0.6), - NoiseFamily(deep), - ) - ) +def _shallow(positions): + return 1.0 - _deep(positions) -def _context(scenario, name="bent"): - chart = scenario.realisation[name] - return SegmentContext( - chart, - scenario.sources[name], - scenario.medium.shear_speed_km_s(chart.centres), - scenario.medium.density_g_cm3(chart.centres), - ) + +HYBRID = (NoiseFamily(SUZUKI_2022, 0.6, _shallow), NoiseFamily(weight=_deep)) -def test_slip_draws_a_blend_its_partners_follow(scenario): - context = _context(scenario) - slip = _hybrid_slip()(context, np.random.default_rng(5)) - blend = slip.latent_covariance - assert isinstance(blend, Blend) - long, short = blend.components - assert long.correlation.target == VonKarman(0.6) - assert short.lengths_km == scenario.sources["bent"].covariance.lengths_km - assert context.moment_nm(slip.slip_m) == pytest.approx(context.source.moment_nm) - assert np.all(np.isfinite(slip.latent)) +def _draw(scenario, model): + return generate( + scenario.realisation, scenario.medium, scenario.sources, model, seed=9 + ) -def test_a_weight_out_of_range_is_refused(scenario): - slip = DrawnSlip(noise=(NoiseFamily(lambda p: np.full(p.shape[:-1], 2.0)),)) - with pytest.raises(RuptureGeneratorError, match=r"weights in \[0, 1\]"): - slip(_context(scenario), np.random.default_rng(0)) +def test_slip_draws_its_noise_and_rake_keeps_its_own(scenario): + plain = _draw(scenario, RuptureModel()) + hybrid = _draw(scenario, RuptureModel(slip=DrawnSlip(noise=HYBRID))) + rake = _draw(scenario, RuptureModel(rake=GaussianRake(noise=HYBRID))) + for name in scenario.realisation: + assert not np.allclose(hybrid[name].slip_m, plain[name].slip_m) + np.testing.assert_array_equal(hybrid[name].rake_deg, plain[name].rake_deg) + np.testing.assert_array_equal(rake[name].slip_m, plain[name].slip_m) + assert not np.allclose(rake[name].rake_deg, plain[name].rake_deg) -def test_the_perturbation_can_keep_the_sources_noise(scenario): - def onset(independent): - model = RuptureModel( - slip=_hybrid_slip(), - perturbation=CorrelatedPerturbation(independent=independent), +def test_the_perturbation_can_draw_its_own_noise(scenario): + slip = DrawnSlip(noise=HYBRID) + own, slips = ( + _draw( + scenario, + RuptureModel(slip=slip, perturbation=CorrelatedPerturbation(noise=noise)), ) - drawn = generate( - scenario.realisation, scenario.medium, scenario.sources, model, seed=9 - ) - return drawn["bent"].onset_s + for noise in (STATIONARY, None) + ) + for name in scenario.realisation: + assert not np.allclose(own[name].onset_s, slips[name].onset_s) - assert not np.allclose(onset("slip"), onset("source")) + +def test_a_weight_out_of_range_is_refused(scenario): + noise = (NoiseFamily(weight=lambda positions: np.full(positions.shape[:-1], 2.0)),) + with pytest.raises(RuptureGeneratorError, match=r"\[0, 1\]"): + _draw(scenario, RuptureModel(slip=DrawnSlip(noise=noise))) # ----------------------------------------------------------------- the config @@ -177,10 +134,13 @@ def onset(independent): NOISE = """ [[slip.noise]] relation = { strike_offset = 1.67, dip_exponent = 0.5, dip_offset = 1.69, dip_saturation = 6.3 } +hurst = 0.6 weight = { type = "ramp", centre_km = 6.5, half_width_km = 1.5, shallow = 1.0, deep = 0.0 } [[slip.noise]] weight = { type = "ramp", centre_km = 6.5, half_width_km = 1.5, shallow = 0.0, deep = 1.0 } + +[[perturbation.noise]] """ @@ -196,22 +156,25 @@ def _hybrid_file(tmp_path, extra=""): return path -def test_noise_families_load_and_round_trip(tmp_path): +def test_noise_loads_round_trips_and_draws(tmp_path): loaded = rupture_config.load(_hybrid_file(tmp_path)) shallow, deep = loaded.slip.noise assert shallow.relation == rupture_config.Correlation( strike_offset=1.67, dip_exponent=0.5, dip_offset=1.69, dip_saturation=6.3 ) - assert deep.relation is None and deep.hurst is None + assert (shallow.hurst, deep.relation) == (0.6, rupture_config.Correlation()) + assert loaded.perturbation.noise == rupture_config.STATIONARY + assert loaded.rake.noise == rupture_config.STATIONARY again = tmp_path / "again.toml" again.write_text(rupture_config.dump(loaded)) assert rupture_config.load(again) == loaded + scenario = rupture_config.build(loaded) + _draw(scenario, scenario.model) def test_a_segments_noise_replaces_the_global_list(tmp_path): - extra = '\n[segments.next.slip]\nnoise = [{ weight = { type = "constant", value = 1.0 } }]\n' + extra = "\n[segments.next.slip]\nnoise = [{ hurst = 0.5 }]\n" loaded = rupture_config.load(_hybrid_file(tmp_path, extra)) - assert len(loaded.slip.noise) == 2 segment = loaded.segments["next"].slip assert segment is not None and len(segment.noise) == 1 assert dataclasses.replace(segment, noise=loaded.slip.noise) == loaded.slip diff --git a/tests/test_config.py b/tests/test_config.py index bb3fb7d..92add1e 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -39,9 +39,16 @@ def test_defaults_are_the_librarys(): (loaded.perturbation, CorrelatedPerturbation()), ): for field in dataclasses.fields(theirs): - if field.name not in {"factor", "slip_weight"}: + if field.name not in {"factor", "slip_weight", "noise"}: assert getattr(ours, field.name) == getattr(theirs, field.name) assert loaded.slip.side_taper == DrawnSlip().side_taper + for ours, theirs in zip( + (loaded.slip.noise, loaded.rake.noise), + (DrawnSlip().noise, GaussianRake().noise), + ): + assert [(dataclasses.astuple(f.relation), f.hurst) for f in ours] == [ + (dataclasses.astuple(f.relation), f.hurst) for f in theirs + ] positions = np.zeros((2, 3)) for profile in (rupture_config.Rise().factor, rupture_config.Speed().factor): np.testing.assert_array_equal( @@ -193,11 +200,10 @@ def test_a_segment_table_overrides_the_global_one(): assert following.rise == bent.rise assert following.perturbation == bent.perturbation - default = rupture_config.build(rupture_config.load(EXAMPLE)) - assert scenario.sources["bent"] == default.sources["bent"] + (longer,) = following.slip.noise assert ( - scenario.sources["next"].covariance.lengths_km[1] - > default.sources["next"].covariance.lengths_km[1] + longer.relation.lengths_km(7.0)[1] + > bent.slip.noise[0].relation.lengths_km(7.0)[1] ) diff --git a/tests/test_models.py b/tests/test_models.py index 82c0dd5..dd170b8 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -12,9 +12,11 @@ from hypothesis.extra import numpy as hnp from rupture_generator import ( + STATIONARY, Drawn, RuptureGeneratorError, SlipField, + embed, gaussian_anamorphosis, generate, ) @@ -150,7 +152,8 @@ def test_the_generator_refuses_a_slip_that_does_not_fit( slip_m = change(ruptures[name].slip_m) def bad(context, rng): - return SlipField(slip_m, np.zeros(slip_m.shape), context.source.covariance) + noise = embed(STATIONARY, context.geometry, context.source.magnitude) + return SlipField(slip_m, np.zeros(slip_m.shape), noise) with pytest.raises(RuptureGeneratorError, match=message): _generate(scenario, dataclasses.replace(scenario.model, slip=bad)) diff --git a/tests/test_numerics.py b/tests/test_numerics.py index 7fff30d..fc01d17 100644 --- a/tests/test_numerics.py +++ b/tests/test_numerics.py @@ -8,7 +8,7 @@ import scipy as sp from rupture_generator import geometry_from_geojson, moment_from_magnitude -from rupture_generator.rupture.source import segment_source +from rupture_generator.rupture.source import MAI_BEROZA from rupture_generator.sampling import ( Covariance, Grid, @@ -130,7 +130,7 @@ def test_the_embedding_delivers_the_covariance_within_its_tolerances(shape, magn # The covariance the draws have is exactly the inverse transform of the clipped # spectrum, read at a grid lag h through the wrap as min(h, P - h), with no draws. grid = Grid(shape, (0.1, 0.1)) - covariance = segment_source(magnitude, 90.0).covariance + covariance = Covariance(VonKarman(), MAI_BEROZA.lengths_km(magnitude)) amplitudes, padded = _embed(grid, covariance) quadrant = sp.fft.idctn(amplitudes**2, type=1) wrapped = [ diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py index 3132e3b..1c110b2 100644 --- a/tests/test_pipeline.py +++ b/tests/test_pipeline.py @@ -20,9 +20,8 @@ _travel_times, first_arrival, ) +from rupture_generator.rupture.models import embed from rupture_generator.rupture.propagator import JumpModel -from rupture_generator.sampling import sampler -from rupture_generator.sampling.field import SAMPLER_CACHE_SIZE def test_every_segment_carries_its_moment(scenario, ruptures): @@ -307,7 +306,7 @@ def test_a_jump_draws_apart_from_any_segments_fields(scenario): def test_a_second_draw_reuses_every_embedding(scenario): - sampler.cache_clear() + embed.cache_clear() misses = [] for _ in range(2): generate( @@ -317,8 +316,7 @@ def test_a_second_draw_reuses_every_embedding(scenario): scenario.model, seed=scenario.seed, ) - misses.append(sampler.cache_info().misses) + misses.append(embed.cache_info().misses) # Slip's latent and rake each embed once per segment; rise time and the onset # displacement draw from slip's embedding. assert misses == [2 * len(scenario.realisation)] * 2 - assert SAMPLER_CACHE_SIZE >= 2 * 20 From 3736a8cdda2af2b9b770c0cc263d577e0f75372e Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 25 Oct 2026 01:16:58 +0000 Subject: [PATCH 3/3] Check blend weights' shape once, in the sampler; fix a stale example comment From /code-review: sampler() refuses weights that don't match the covariances and the grid with a RuptureGeneratorError, and is the one place that checks their range; embed() only asks the weights to be finite. colombia.toml no longer claims to spell out the relation. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01KrujMjfczwcCciJwtErEbs --- examples/colombia.toml | 2 +- rupture_generator/rupture/models.py | 2 +- rupture_generator/sampling/field.py | 8 +++++++- tests/test_blend.py | 1 + 4 files changed, 10 insertions(+), 3 deletions(-) diff --git a/examples/colombia.toml b/examples/colombia.toml index 8860d39..a61d0de 100644 --- a/examples/colombia.toml +++ b/examples/colombia.toml @@ -20,7 +20,7 @@ # # The parameters are not tuned to match the observations of this earthquake. # -# The correlation relation is Mai & Beroza's, spelled out. The depth profiles' shallow +# The noise is the default, Mai & Beroza's lengths. The depth profiles' shallow # ramps end at 8 km and their deep ramps start at 197.5 km, so they are flat across the # fault. # diff --git a/rupture_generator/rupture/models.py b/rupture_generator/rupture/models.py index 8e37649..934345e 100644 --- a/rupture_generator/rupture/models.py +++ b/rupture_generator/rupture/models.py @@ -603,7 +603,7 @@ def embed( ], [ read_at_centres( - geometry, "a noise weight", family.weight, unit_interval, "[0, 1]" + geometry, "a noise weight", family.weight, np.isfinite, "finite values" ) for family in noise ], diff --git a/rupture_generator/sampling/field.py b/rupture_generator/sampling/field.py index add0dde..08a0238 100644 --- a/rupture_generator/sampling/field.py +++ b/rupture_generator/sampling/field.py @@ -263,13 +263,19 @@ def sampler( Raises ------ RuptureGeneratorError - If a weight is outside ``[0, 1]``, a cell has no positive weight, or a + If the weights don't match the covariances and the grid, a weight is outside + ``[0, 1]``, a cell has no positive weight, or a covariance doesn't embed within the cell cap or isn't positive definite at any margin tried. """ shares = np.ones((len(covariances), *grid.shape)) if weights is not None: shares = np.asarray(weights, dtype=np.float64) + if shares.shape != (len(covariances), *grid.shape): + raise RuptureGeneratorError( + f"{len(covariances)} covariances on a {grid.shape} grid need weights " + f"shaped {(len(covariances), *grid.shape)}, not {shares.shape}" + ) if not np.all((shares >= 0.0) & (shares <= 1.0)): raise RuptureGeneratorError("a covariance's weight lies outside [0, 1]") gains = np.sin(0.5 * np.pi * shares) diff --git a/tests/test_blend.py b/tests/test_blend.py index 56db45f..137bffb 100644 --- a/tests/test_blend.py +++ b/tests/test_blend.py @@ -72,6 +72,7 @@ def test_draws_are_standard_normal_and_the_pair_independent(): [ ([np.zeros(GRID.shape), DEEP], "weigh every covariance at zero"), ([1.0 - DEEP, 2.0 * DEEP], r"outside \[0, 1\]"), + ([DEEP], "need weights shaped"), ], ) def test_bad_weights_are_refused(weights, message):