feat: add MPS module with DirectSampling simulation

[n0228a]

Add Multiple-Point Statistics via Direct Sampling (gstools.mps)

This PR introduces a new gstools.mps submodule implementing Multiple-Point Statistics (MPS) simulation through the Direct Sampling algorithm (Mariethoz & Caers 2010; Juda 2022). MPS is a fundamentally different paradigm from two-point geostatistics: instead of parametric covariance models, spatial structure is learned directly from a training image (TI).

What was added

gstools/mps/training_image.py — TrainingImage

The MPS analogue of CovModel. Wraps an n-dimensional NumPy array and encapsulates the distance function used to compare data events. Supports:

• Categorical variables (e.g. facies): weighted mismatch distance (Mariethoz 2010, Eq. 3)
• Continuous variables: L1, L2, Lp (arbitrary Minkowski exponent), and variation distance with automatic mean-shift correction (Mariethoz 2010, Eq. 9)
• Spatial-decay weighting of neighbours (distance power δ, Mariethoz 2010, Eq. 5)
• Upweighting of hard conditioning nodes (Mariethoz 2010, §3 ¶26)

gstools/mps/distance.py — pure distance functions

Stateless helper functions (categorical_dist, l1_dist, l2_dist, lp_dist, variation_dist, compute_node_weights) separated from class state for clarity and reusability.

gstools/mps/direct_sampling.py — DirectSampling

Subclasses gstools.field.base.Field and follows the same call interface as SRF. Key features:

• Works on n-dimensional structured grids
• Randomized simulation path (Mariethoz 2010, §3 ¶18)
• Neighbour search sorted by Euclidean distance with configurable max_offset
• Configurable scan_fraction caps TI traversal per node
• threshold=0.0 activates DSBC mode (best-candidate, full scan)
• set_condition() for hard conditioning data, with smart nearest-node snapping and collision resolution (Mariethoz 2010, §3)

examples/13_mps/channel_demo.py

An end-to-end demo using the classic Strebelle (2002) channelized fluvial training image, demonstrating conditional simulation with 100 hard-data points.

Usage

import numpy as np
from gstools import mps

ti = mps.TrainingImage(ti_array, categorical=True)
ds = mps.DirectSampling(ti, n_neighbors=30, scan_fraction=1.0, threshold=0.01)
ds.set_condition(cond_pos, cond_val)

x = y = np.arange(100, dtype=float)
field = ds([x, y], seed=42)

References

• Mariethoz, G. & Caers, J. (2010): Direct sampling method to perform multiple‑point geostatistical simulations. Water Resources Research.
• Juda, P. et al. (2022): A comparison of three recent discrete exact simulation algorithms for Gaussian random fields. Statistics and Computing.