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..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. # @@ -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 028e3d0..6e2fe4d 100644 --- a/rupture_generator/__init__.py +++ b/rupture_generator/__init__.py @@ -14,9 +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 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. +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 @@ -44,11 +44,13 @@ ramp_field, ) from rupture_generator.rupture.models import ( + STATIONARY, CorrelatedPerturbation, CorrelatedRise, Drawn, DrawnSlip, GaussianRake, + NoiseFamily, Perturbation, PerturbationModel, PulseSettings, @@ -60,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, @@ -79,6 +83,8 @@ __all__ = [ "MAI_BEROZA", + "STATIONARY", + "SUZUKI_2022", "CorrelatedPerturbation", "CorrelatedRise", "CorrelationRelation", @@ -92,6 +98,7 @@ "JumpSeed", "Layers", "Medium", + "NoiseFamily", "Perturbation", "PerturbationModel", "PulseSettings", @@ -110,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 3dd6208..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 @@ -64,6 +63,7 @@ CorrelatedRise, DrawnSlip, GaussianRake, + NoiseFamily, PulseSettings, RuptureModel, ShearFractionSpeed, @@ -73,7 +73,6 @@ from rupture_generator.rupture.source import ( CorrelationRelation, SegmentSource, - magnitude_from_moment, moment_from_magnitude, segment_source, split_moment, @@ -89,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 @@ -179,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: @@ -241,10 +240,38 @@ def field(self) -> SpatialField: # ----------------------------------------------------------- the default models +@dataclasses.dataclass(frozen=True) +class Family(NoiseFamily): + """``[[.noise]]``: :class:`~rupture_generator.rupture.models.NoiseFamily`. + + Parameters + ---------- + relation : Correlation + ``relation``, the family's correlation lengths. + weight : DepthProfile + ``weight``, how much of the field the family is. + """ + + 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`.""" + """``[slip]``: :class:`~rupture_generator.rupture.models.DrawnSlip`. + + Parameters + ---------- + noise : tuple of Family + ``[[slip.noise]]``, slip's noise. + """ + noise: tuple[Family, ...] = STATIONARY Config = _Strict @@ -268,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 @@ -289,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 @@ -434,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 @@ -460,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 @@ -556,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 @@ -565,7 +598,6 @@ class SegmentConfig: rake: Rake | None = None speed: Speed | None = None perturbation: Perturbation | None = None - correlation: Correlation | None = None Config = _Strict @@ -680,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() @@ -800,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 ( @@ -820,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(): @@ -830,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() } @@ -947,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 54e322e..934345e 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,6 +27,7 @@ """ import dataclasses +import functools from collections.abc import Callable, Mapping from typing import Protocol @@ -31,13 +37,21 @@ 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 ( + MAI_BEROZA, + CorrelationRelation, + SegmentSource, +) from rupture_generator.sampling import ( + HURST, + NORMAL, Covariance, Grid, Marginal, MarginalFamily, PreCorrected, + Sampler, + VonKarman, mix, sampler, standardise, @@ -80,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 @@ -114,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: @@ -131,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 - 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. @@ -141,7 +148,7 @@ class SlipField: slip_m: CellArray latent: CellArray - latent_covariance: Covariance + sampler: Sampler statistics: Mapping[str, float] = dataclasses.field(default_factory=dict) @@ -393,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. @@ -508,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, _ = sampler(context.grid, slip.latent_covariance).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, np.isfinite, "finite values" + ) + for family in noise + ], + ) # ------------------------------------------------------------------------ slip @@ -534,11 +628,10 @@ def _reach(mask: CellMask, axis: int, *, reverse: bool) -> np.ndarray: @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 ---------- @@ -557,6 +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 : Noise + The noise slip's latent comes from, and rise time's independent part. """ family: MarginalFamily = "truncated_exponential" @@ -565,6 +660,7 @@ class DrawnSlip: side_taper: float = 0.02 top_taper: float = 0.0 bottom_taper: float = 0.0 + noise: Noise = STATIONARY def __post_init__(self) -> None: """Refuse a taper outside ``[0, 0.5]``, or a spread outside the family's range.""" @@ -623,14 +719,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 noise doesn't embed on this chart, 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 + 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: @@ -641,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, ) @@ -748,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) @@ -756,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( @@ -773,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.""" @@ -806,7 +903,9 @@ def __call__(self, context: SegmentContext, rng: np.random.Generator) -> Drawn: Drawn Rake in degrees. """ - latent, _ = sampler(context.grid, context.source.covariance).draw(rng) + latent, _ = embed(self.noise, context.geometry, context.source.magnitude).draw( + rng + ) return Drawn(context.source.rake_deg + self.sigma_deg * standardise(latent)) @@ -908,12 +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. + 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 + noise: Noise | None = None def __post_init__(self) -> None: """Refuse a negative spread, or a blend of no width.""" @@ -956,11 +1059,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 + 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(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..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 ---------- @@ -12,6 +13,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 @@ -22,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 @@ -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. @@ -98,22 +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 * magnitude - self.dip_offset), - 10.0 ** (self.strike_exponent * magnitude - 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( + 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 ---------- @@ -121,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.""" @@ -137,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 @@ -152,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( @@ -208,7 +211,7 @@ def split_moment( __all__ = [ "MAI_BEROZA", - "VON_KARMAN", + "SUZUKI_2022", "CorrelationRelation", "SegmentSource", "magnitude_from_moment", diff --git a/rupture_generator/sampling/field.py b/rupture_generator/sampling/field.py index 284a94b..08a0238 100644 --- a/rupture_generator/sampling/field.py +++ b/rupture_generator/sampling/field.py @@ -2,14 +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. 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 @@ -42,10 +43,6 @@ 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.""" - @dataclasses.dataclass(frozen=True) class Grid: @@ -126,19 +123,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) @@ -179,28 +183,36 @@ def _embed(grid: Grid, covariance: Covariance) -> tuple[FieldArray, tuple[int, i @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 ---------- @@ -213,32 +225,83 @@ def draw(self, rng: np.random.Generator) -> tuple[FieldArray, 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) + 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] + +def sampler( + grid: Grid, + covariances: Sequence[Covariance], + weights: Sequence[FieldArray] | None = None, +) -> Sampler: + """Embed covariances on ``grid`` at one padded shape, and weigh them cell by cell. -@functools.lru_cache(maxsize=SAMPLER_CACHE_SIZE) -def sampler(grid: Grid, covariance: Covariance) -> Sampler: - """Embed ``covariance`` on ``grid``, or return the embedding already made. + 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. + 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. + 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. """ - return Sampler(grid.shape, *_embed(grid, covariance)) + 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) + 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 + ) + ] + 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 new file mode 100644 index 0000000..137bffb --- /dev/null +++ b/tests/test_blend.py @@ -0,0 +1,181 @@ +"""Noise families: blended sampling, the models that draw noise, and the config.""" + +import dataclasses + +import numpy as np +import pytest + +from rupture_generator import ( + STATIONARY, + SUZUKI_2022, + CorrelatedPerturbation, + DrawnSlip, + GaussianRake, + NoiseFamily, + RuptureGeneratorError, + RuptureModel, + generate, +) +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)) +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(): + 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_covariances_share_the_widest_padded_shape(): + across = Covariance(VonKarman(), (12.0, 2.0)) + 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(): + 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(): + blend = sampler(GRID, [LONG, SHORT], [1.0 - DEEP, DEEP]) + rng = np.random.default_rng(11) + 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 + assert abs(float(real.var() - 1.0)) < 0.05 + assert abs(float(np.mean(real * imaginary))) < 0.02 + + +@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\]"), + ([DEEP], "need weights shaped"), + ], +) +def test_bad_weights_are_refused(weights, message): + with pytest.raises(RuptureGeneratorError, match=message): + sampler(GRID, [LONG, SHORT], weights) + + +# ----------------------------------------------------------------- the models + + +def _deep(positions): + return np.clip((positions[..., 2] - 4.0) / 4.0, 0.0, 1.0) + + +def _shallow(positions): + return 1.0 - _deep(positions) + + +HYBRID = (NoiseFamily(SUZUKI_2022, 0.6, _shallow), NoiseFamily(weight=_deep)) + + +def _draw(scenario, model): + return generate( + scenario.realisation, scenario.medium, scenario.sources, model, seed=9 + ) + + +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_draw_its_own_noise(scenario): + slip = DrawnSlip(noise=HYBRID) + own, slips = ( + _draw( + scenario, + RuptureModel(slip=slip, perturbation=CorrelatedPerturbation(noise=noise)), + ) + for noise in (STATIONARY, None) + ) + for name in scenario.realisation: + assert not np.allclose(own[name].onset_s, slips[name].onset_s) + + +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 + +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]] +""" + + +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_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 (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 = [{ hurst = 0.5 }]\n" + loaded = rupture_config.load(_hybrid_file(tmp_path, extra)) + 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