feat: add gstools.mps — Multiple Point Statistics via Direct Sampling

[n0228a]

Summary

Adds a new gstools.mps subpackage implementing Multiple Point Statistics (MPS)
simulation via the Direct Sampling (DS) algorithm (Mariethoz et al. 2010) with
the DSBC parametrization (Juda et al. 2022). MPS is a fundamentally different
paradigm from two-point geostatistics: instead of a parametric covariance model,
spatial structure is learned directly from a training image (TI), enabling
reproduction of curvilinear, connected patterns — channels, fractures, categorical
facies — that variogram-based methods cannot capture.

The implementation is pure Python + NumPy. Unlike field/krige/variogram,
there is no gstools_core (Rust) or Cython path for mps; the per-node TI scan is
vectorised over candidate blocks and parallelised across nodes with threads.
Node-level parallelism is a dependency DAG rather than a barrier — a node is
dispatched as soon as the neighbours it reads are settled — and the same machinery
drives the main path and every post-processing pass, so output is bit-identical for
any num_threads.

Supersedes #409.

---

Public API

Five classes exported from gstools.mps. DirectSampling, MPSModel,
TrainingImage, and Zone are also available at the top-level gstools
namespace; Variable is available as gstools.mps.Variable.

Variable

Holds one variable's data and its per-variable configuration:

• data — n-D NumPy array
• categorical — whether to use weighted mismatch distance (Mariethoz2010 Eq. 3)
• distance — continuous metric: "l1" (Juda2022 Eq. 7, default), "l2",
  "l<p>" (arbitrary Minkowski), "variation" or "variation<p>" (Mariethoz2010
  Eq. 9, with automatic mean-shift correction)
• n_neighbors — data-event size (maximum neighbours, default 32)
• max_radius — optional Euclidean cap on neighbour selection
• weight — relative weight in the multivariate joint distance; None for equal
  share; weights are automatically normalized and do not need to sum to 1
• penalty_matrix — optional (C, C) per-category mismatch cost T[u, v]
  replacing the flat binary mismatch for categorical variables, so that
  "sand mistaken for silt" can cost less than "sand mistaken for shale":
  Σ w_i · T[Z(x_i), Z(y_i)] / Σ w_i. Asymmetric matrices are allowed. Validated
  to be 2-D square with an all-zero diagonal and off-diagonal entries in (0, 1],
  and the variable's data must hold integer-valued category codes < C — these
  constraints are what keep the distance in [0, 1] (there is no d_max
  normalization on this path)

TrainingImage

The MPS analogue of CovModel. Three construction modes:

# Univariate — plain array + keyword args
ti = gs.TrainingImage(ti_array, categorical=True, n_neighbors=30)

# Univariate — pre-configured Variable
ti = gs.TrainingImage(gs.mps.Variable("z", ti_array, categorical=True, n_neighbors=30))

# Multivariate — list of Variable objects
ti = gs.TrainingImage([
    gs.mps.Variable("facies",   facies_arr, categorical=True,  weight=0.6, n_neighbors=30),
    gs.mps.Variable("porosity", por_arr,    categorical=False, weight=0.4,
                    distance="l2", n_neighbors=15),
])

Shared option: distance_power — spatial-decay exponent δ for neighbour
weighting (Mariethoz2010 Eq. 5); 0.0 (default) gives uniform weights.

window(slices) returns a sub-region view of a TI, carrying over every search
parameter and the distance kernel — useful for cutting several zone TIs out of one
large image.

NaN masking — undefined TI cells (NaN) are excluded from distance
comparisons and never pasted into the simulation grid.

MPSModel

Configuration object bundling a TrainingImage with the search parameters and
every optional feature:

model = gs.MPSModel(ti, scan_fraction=0.2, threshold=0.0)

Argument	Purpose
scan_fraction, threshold, cond_weight, boundary	core DS search parameters
rotation, scale	geometric non-stationarity (see below)
zones	list of Zone objects for zonated simulation
post_processing, post_processing_factor, post_processing_path	re-simulation passes (see below)

Geometric non-stationarity — rotation= / scale= accept a scalar, a vector,
a full per-node array, or a callable evaluated on the grid coordinates, and apply
the lag transform of Mariethoz2010 §6.2. A stationary TI can therefore produce
patterns that rotate and stretch across the grid without modifying the TI itself.
rotation is rejected on 1-D grids (no rotation planes exist); scale stays
meaningful.

Post-processing (post_processing=p, Meerschman et al. 2013 §4) runs p
passes after the main path, each re-simulating every non-conditioning node against
a fully informed neighbourhood to remove simulation noise. Me13 §7 advises at least
p=1 for categorical simulations. post_processing_factor (p_f ≥ 1) divides
scan_fraction and n_neighbors during the passes to trade fidelity for speed;
Me13 finds little effect and recommends p_f=1. Passes are Gauss–Seidel, but that
ordering is carried by the same dependency DAG the main path uses, so they honour
num_threads while staying bit-identical to serial execution at any thread count.
Each pass re-visits essentially the whole grid, so p > 0 dominates the cost of a
run — budget accordingly. post_processing=0 (default) is bit-identical to
omitting the parameter and consumes no extra RNG draws.

post_processing_path controls the post-pass visit order and is an
implementation extension — Me13 §4 does not prescribe one. None (default)
inherits the main-path mode; "random", "sequential", "same" (reuse the main
pass's order), or an explicit (N, dim) array are accepted.

Zone

Zone(ti, where) binds a training image to a region of the simulation grid for
zonated (non-stationary) simulation, Mariethoz2010 para [40]. Nodes inside the
region scan this zone's TI instead of the primary one, while the data event is
still built from simulation-grid neighbours regardless of their zone — which keeps
zone boundaries coherent instead of cutting the pattern at the seam.

• where is a boolean array of the grid shape, or a callable
  f(x, [y, z, ...]) -> bool mask. Non-boolean arrays raise; there is no silent
  != 0 coercion, so an integer zonation raster is written
  Zone(ti_k, where=(labels == k)).
• A zone contributes only its data, shape, and continuous range (d_max).
  Every search hyper-parameter and the distance kernel — n_neighbors,
  max_radius, distance, weight, penalty_matrix — always come from the
  primary TI, so setting them on a zone TI has no effect.
• Zoning does not perturb RNG consumption or the DAG, so output stays
  deterministic and thread-count invariant.

DirectSampling

Subclasses gstools.field.base.Field and follows the same call interface as
SRF. Takes an MPSModel.

All search parameters (scan_fraction, threshold, cond_weight, boundary,
n_neighbors, max_radius) are read-only on DirectSampling — they
reflect the underlying MPSModel and Variable objects. n_neighbors and
max_radius collapse to a scalar when all variables share the same value,
otherwise return a {name: value} dict. To change them, mutate the MPSModel
or Variable directly.

Key capabilities:

• n-dimensional structured grids — fully dimension-agnostic
• Randomised simulation path (Mariethoz2010 §3 ¶13); also accepts
  "sequential" or an explicit (N, dim) visit-order array via the path
  argument
• scan_fraction caps the per-node TI scan (Juda2022 §2)
• threshold=0.0 — DSBC mode: full scan, best-candidate selection
• boundary="strict" (default) or "partial" — partial drops lags that exceed
  the TI extent and searches with the reduced template (Mariethoz2010 §6.2)
• set_condition(cond_pos, cond_val) — hard-data conditioning with nearest-node
  snapping and per-variable collision resolution; cond_val accepts a
  {variable: array} dict for multivariate TIs, with np.nan to leave a
  variable unconditioned at a point; conditioning weight is set via MPSModel
• set_mv_transforms() — per-variable mean / normalizer / trend post-processing.
  Note this is the Field pipeline and is unrelated to MPSModel's DS
  post_processing passes
• DAG-based parallelism via num_threads (or gstools.config.NUM_THREADS)
• Dual RNG seeds (path_seed, node_seed) for independent control of visit
  order vs. TI scan randomness
• Multivariate return value: {variable: ndarray} dict; each variable also
  stored as a named field accessible via ds["facies"] / ds.all_fields

---

Usage

import numpy as np
import gstools as gs

ti = gs.TrainingImage(ti_array, categorical=True, n_neighbors=30)
model = gs.MPSModel(ti, scan_fraction=0.2, threshold=0.0)
ds = gs.DirectSampling(model)
ds.set_condition([cond_x, cond_y], cond_val)
field = ds([np.arange(100, dtype=float)] * 2, seed=42)

Multivariate TIs are built from a list of Variable objects and return a
{variable: ndarray} dict; rotation/scale, zones, and post_processing are
all MPSModel arguments layered on the same call. See the examples below for each.

---

Module layout

File	Responsibility
training_image.py	Variable, TrainingImage — data and distance functions
data_event.py	DataEvent — internal value object bundling per-variable data-event arrays
distance.py	Stateless vectorised distance helpers
model.py	MPSModel — validated search-parameter config
zone.py	Zone — a training image bound to a region of the grid
neighbors.py	Neighbour selection, geometric lag transform, search-window bounds
scan.py	TI scanning loop
runner.py	Simulation path and DAG thread executor
nonstationary.py	Rotation / scale map resolution
simulate.py	_DirectSamplingEngine + ds_simulate entry point
direct_sampling.py	DirectSampling — public Field subclass

---

Examples (examples/13_mps/)

examples/13_mps/ is now wired into docs/source/tutorials.rst and
sphinx_gallery_conf, so these run at documentation build time. Several fetch
their training images from public sources at first run and cache them next to the
script.

File	What it shows
00_simple_unconditional.py	Minimal unconditional simulation on a synthetic channel TI
01_conditional.py	Hard-data conditioning with set_condition()
02_continuous.py	Continuous variable simulation
03_channel_strebelle.py	Strebelle (2002) fluvial TI, conditional
04_channel_strebelle_nonstat.py	Non-stationary simulation via MPSModel(rotation=...)
05_channel_strebelle_radial.py	Radial non-stationarity with custom spiral path; reproduces Mariethoz2010 Fig. 7
06_channel_multivariate_fig8.py	Multivariate co-simulation; reproduces Mariethoz2010 Fig. 8
07_variation_distance_fig6.py	Variation-distance metric; reproduces Mariethoz2010 Fig. 6
08_bivariate_joint_fig4.py	Bivariate joint simulation; reproduces Mariethoz2010 Fig. 4
09_continuous_stone_fig3.py	Continuous simulation on the GAIA stone TI; reproduces Mariethoz2010 Fig. 3
10_zonated_facies.py	Zonation with Zone — a separate TI per region, via window() views
00_mps_overview.ipynb	End-to-end notebook covering all modes

---

Changes outside gstools.mps

MPS's per-node rotation/scale non-stationarity (Mariethoz2010 §6.2) needs a
per-axis scale with no pinned axis — including uniform dilation (M10's
affinity r), which the existing anis convention (set_anis, axis 0 fixed
to ratio 1, every other axis a ratio to it) cannot express. gstools.tools
gains a second, more general primitive family for this:

• set_scale — validates a scalar-or-dim-vector scale, no padding, no pin
• matrix_scale, matrix_transform, matrix_detransform — diagonal scale
  matrix, rotate∘scale, and its exact inverse
• great_circle_to_chordal added to tools/geometric.py's module __all__,
  alongside chordal_to_great_circle (already present, just not exported)

gstools.mps.neighbors uses matrix_detransform directly for the per-node
SG→TI lag pullback. Separately, the existing anis-pinned helpers —
matrix_isotropify, matrix_anisotropify, matrix_isometrize, and
matrix_anisometrize — are now all expressed as the axis-0-pinned special
case of the same general primitives, instead of each duplicating its own
diag/matmul construction. matrix_isometrize and matrix_anisometrize are
CovModel's actual call path (covmodel/base.py::isometrize/anisometrize),
so this is verified behaviour-preserving two ways: tests/test_geometric.py::TestWrapperBitIdentical
checks exact (assert_array_equal, not allclose) equality against a frozen
copy of every pre-refactor implementation over 150 randomized cases across
dims 2–4, and the full test_covmodel.py suite passes unmodified.

---

References

• Mariethoz, G., Renard, P., Straubhaar, J. (2010). The Direct Sampling method
  to perform multiple-point geostatistical simulations. Water Resources Research.
  https://doi.org/10.1029/2008WR007621
• Meerschman, E., Pirot, G., Mariethoz, G., Straubhaar, J., Van Meirvenne, M.,
  Renard, P. (2013). A practical guide to performing multiple-point statistical
  simulations with the Direct Sampling algorithm. Computers & Geosciences.
  https://doi.org/10.1016/j.cageo.2012.09.019
• Juda, P., Renard, P., Straubhaar, J. (2022). A parsimonious parametrization of
  the Direct Sampling algorithm for multiple-point statistical simulations.
  Applied Computing and Geosciences.
  https://doi.org/10.1016/j.acags.2022.100091
• Strebelle, S. (2002). Conditional simulation of complex geological structures
  using multiple-point statistics. Mathematical Geology.