GeoStat-Framework · n0228a · May 12, 2026 · May 12, 2026 · May 12, 2026 · May 12, 2026
diff --git a/.claude/agents/code-reviewer.md b/.claude/agents/code-reviewer.md
@@ -0,0 +1,40 @@
+---
+name: code-reviewer
+description: "Read-only code quality and domain correctness reviewer. Use on any new implementation or PR before merging. Reviews for geostatistical correctness, API consistency, and docstring completeness."
+tools: [Read, Glob, Grep]
+model: sonnet
+---
+
+You are a scientific Python developer with expertise in geostatistics. You review code for correctness and quality. You cannot modify files.
+
+## What you review
+
+**Domain correctness**
+- Does the algorithm match its mathematical definition? Check against the docstring references (Mariethoz2010, Juda2022, etc.).
+- Are edge cases handled correctly or explicitly documented as unsupported?
+- Are physical constraints enforced (non-negative variance, valid distance range [0,1], etc.)?
+
+**API consistency**
+- Does the public API follow the GSTools pattern? Compare against `field/srf.py` and `krige/base.py` as references.
+- Are parameter names consistent with the rest of GSTools (`pos`, `seed`, `mesh_type`, not `positions`, `random_seed`, `grid_type`)?
+- Does `__call__` accept `pos` tuples in the standard GSTools convention?
+
+**Docstrings**
+- NumPy style with explicit parameter shapes.
+- `Returns` section describes the output array shape.
+- `Examples` section contains runnable code.
+- Mathematical references cited where the formula is non-obvious.
+
+**Implementation**
+- No `np.random` calls — uses `gstools.random.RNG`.
+- No hardcoded 2D/3D assumptions — uses `len(shape)` for dimensionality.
+- No Python loops over individual array elements where a vectorized operation exists.
+
+## Output format
+
+Report issues grouped by severity:
+- **Critical**: wrong result, silent data corruption, incorrect algorithm
+- **Major**: API inconsistency, missing edge case handling, broken docstring
+- **Minor**: style, naming, missing example
+
+For each issue: file, line range, what is wrong, what it should be.
diff --git a/.git_archival.txt b/.git_archival.txt
diff --git a/.gitattributes b/.gitattributes
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
diff --git a/.gitignore b/.gitignore
diff --git a/.readthedocs.yml b/.readthedocs.yml
diff --git a/AUTHORS.md b/AUTHORS.md
diff --git a/CITATION.bib b/CITATION.bib
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
diff --git a/LICENSE b/LICENSE
diff --git a/README.md b/README.md
diff --git a/docs/Makefile b/docs/Makefile
diff --git a/docs/source/_static/custom.css b/docs/source/_static/custom.css
diff --git a/docs/source/_templates/autosummary/class.rst b/docs/source/_templates/autosummary/class.rst
diff --git a/docs/source/_templates/autosummary/module.rst b/docs/source/_templates/autosummary/module.rst
diff --git a/docs/source/_templates/layout.html b/docs/source/_templates/layout.html
diff --git a/docs/source/api.rst b/docs/source/api.rst
diff --git a/docs/source/authors.rst b/docs/source/authors.rst
diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
diff --git a/docs/source/conf.py b/docs/source/conf.py
diff --git a/docs/source/contents.rst b/docs/source/contents.rst
diff --git a/docs/source/index.rst b/docs/source/index.rst
diff --git a/docs/source/pics/05_ordinary.png b/docs/source/pics/05_ordinary.png
diff --git a/docs/source/pics/05_simple.png b/docs/source/pics/05_simple.png
diff --git a/docs/source/pics/06_ensemble.png b/docs/source/pics/06_ensemble.png
diff --git a/docs/source/pics/07_00_std.png b/docs/source/pics/07_00_std.png
diff --git a/docs/source/pics/07_01_lognormal.png b/docs/source/pics/07_01_lognormal.png
diff --git a/docs/source/pics/07_02_binary.png b/docs/source/pics/07_02_binary.png
diff --git a/docs/source/pics/07_03_zinnharvey.png b/docs/source/pics/07_03_zinnharvey.png
diff --git a/docs/source/pics/07_04_arcsin.png b/docs/source/pics/07_04_arcsin.png
diff --git a/docs/source/pics/07_05_combine.png b/docs/source/pics/07_05_combine.png
diff --git a/docs/source/pics/09_cond_ens.png b/docs/source/pics/09_cond_ens.png
diff --git a/docs/source/pics/20_gstools.png b/docs/source/pics/20_gstools.png
diff --git a/docs/source/pics/20_pykrige.png b/docs/source/pics/20_pykrige.png
diff --git a/docs/source/pics/2d_multi_tree_pgs.png b/docs/source/pics/2d_multi_tree_pgs.png
diff --git a/docs/source/pics/2d_pgs.png b/docs/source/pics/2d_pgs.png
diff --git a/docs/source/pics/2d_tree_pgs.png b/docs/source/pics/2d_tree_pgs.png
diff --git a/docs/source/pics/3d_gau_field.png b/docs/source/pics/3d_gau_field.png
diff --git a/docs/source/pics/3d_pgs.png b/docs/source/pics/3d_pgs.png
diff --git a/docs/source/pics/3d_tree_pgs.png b/docs/source/pics/3d_tree_pgs.png
diff --git a/docs/source/pics/GS_3d_vector_field.png b/docs/source/pics/GS_3d_vector_field.png
diff --git a/docs/source/pics/GS_pyvista_cut.png b/docs/source/pics/GS_pyvista_cut.png
diff --git a/docs/source/pics/cond_ens.png b/docs/source/pics/cond_ens.png
diff --git a/docs/source/pics/cov_model_vario.png b/docs/source/pics/cov_model_vario.png
diff --git a/docs/source/pics/demonstrator.png b/docs/source/pics/demonstrator.png
diff --git a/docs/source/pics/exp_vario_fit.png b/docs/source/pics/exp_vario_fit.png
diff --git a/docs/source/pics/gau_field.png b/docs/source/pics/gau_field.png
diff --git a/docs/source/pics/gstools.ico b/docs/source/pics/gstools.ico
diff --git a/docs/source/pics/gstools.png b/docs/source/pics/gstools.png
diff --git a/docs/source/pics/gstools_150.png b/docs/source/pics/gstools_150.png
diff --git a/docs/source/pics/paraview.png b/docs/source/pics/paraview.png
diff --git a/docs/source/pics/pyvista_export.png b/docs/source/pics/pyvista_export.png
diff --git a/docs/source/pics/srf_tut_exp_ani_rot.png b/docs/source/pics/srf_tut_exp_ani_rot.png
diff --git a/docs/source/pics/srf_tut_gau_field.png b/docs/source/pics/srf_tut_gau_field.png
diff --git a/docs/source/pics/srf_tut_gau_field_ens.png b/docs/source/pics/srf_tut_gau_field_ens.png
diff --git a/docs/source/pics/srf_tut_merge.png b/docs/source/pics/srf_tut_merge.png
diff --git a/docs/source/pics/srf_tut_unstr.png b/docs/source/pics/srf_tut_unstr.png
diff --git a/docs/source/pics/srf_vector_field.png b/docs/source/pics/srf_vector_field.png
diff --git a/docs/source/pics/stab_vario_fit.png b/docs/source/pics/stab_vario_fit.png
diff --git a/docs/source/pics/tplstable_field.png b/docs/source/pics/tplstable_field.png
diff --git a/docs/source/pics/vario_tut_aniso_fit_exp.png b/docs/source/pics/vario_tut_aniso_fit_exp.png
diff --git a/docs/source/pics/vario_tut_fit_exp.png b/docs/source/pics/vario_tut_fit_exp.png
diff --git a/docs/source/pics/vario_tut_herten.png b/docs/source/pics/vario_tut_herten.png
diff --git a/docs/source/pics/vario_tut_new_herten.png b/docs/source/pics/vario_tut_new_herten.png
diff --git a/docs/source/pics/vec_srf_tut_exp.png b/docs/source/pics/vec_srf_tut_exp.png
diff --git a/docs/source/pics/vec_srf_tut_gau.png b/docs/source/pics/vec_srf_tut_gau.png
diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst
diff --git a/examples/00_misc/00_tpl_stable.py b/examples/00_misc/00_tpl_stable.py
diff --git a/examples/00_misc/01_export.py b/examples/00_misc/01_export.py
diff --git a/examples/00_misc/02_check_rand_meth_sampling.py b/examples/00_misc/02_check_rand_meth_sampling.py
diff --git a/examples/00_misc/04_herten.py b/examples/00_misc/04_herten.py
diff --git a/examples/00_misc/05_standalone_field.py b/examples/00_misc/05_standalone_field.py
diff --git a/examples/00_misc/README.rst b/examples/00_misc/README.rst
diff --git a/examples/00_misc/grid_dim_origin_spacing.txt b/examples/00_misc/grid_dim_origin_spacing.txt
diff --git a/examples/00_misc/herten_transmissivity.gz b/examples/00_misc/herten_transmissivity.gz
diff --git a/examples/01_random_field/00_gaussian.py b/examples/01_random_field/00_gaussian.py
diff --git a/examples/01_random_field/01_srf_ensemble.py b/examples/01_random_field/01_srf_ensemble.py
diff --git a/examples/01_random_field/02_fancier.py b/examples/01_random_field/02_fancier.py
diff --git a/examples/01_random_field/03_unstr_srf_export.py b/examples/01_random_field/03_unstr_srf_export.py
diff --git a/examples/01_random_field/04_srf_merge.py b/examples/01_random_field/04_srf_merge.py
diff --git a/examples/01_random_field/06_pyvista_support.py b/examples/01_random_field/06_pyvista_support.py
diff --git a/examples/01_random_field/08_fourier.py b/examples/01_random_field/08_fourier.py
diff --git a/examples/01_random_field/09_fourier_trans.py b/examples/01_random_field/09_fourier_trans.py
diff --git a/examples/01_random_field/README.rst b/examples/01_random_field/README.rst
diff --git a/examples/01_random_field/hexagon_cells.npy b/examples/01_random_field/hexagon_cells.npy
diff --git a/examples/01_random_field/hexagon_points.npy b/examples/01_random_field/hexagon_points.npy
diff --git a/examples/02_cov_model/00_intro.py b/examples/02_cov_model/00_intro.py
diff --git a/examples/02_cov_model/07_roughness.py b/examples/02_cov_model/07_roughness.py
diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst
diff --git a/examples/03_variogram/00_fit_variogram.py b/examples/03_variogram/00_fit_variogram.py
diff --git a/examples/03_variogram/05_auto_fit_variogram.py b/examples/03_variogram/05_auto_fit_variogram.py
diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py
diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst
diff --git a/examples/04_vector_field/00_2d_vector_field.py b/examples/04_vector_field/00_2d_vector_field.py
diff --git a/examples/04_vector_field/README.rst b/examples/04_vector_field/README.rst
diff --git a/examples/05_kriging/01_ordinary_kriging.py b/examples/05_kriging/01_ordinary_kriging.py
diff --git a/examples/05_kriging/README.rst b/examples/05_kriging/README.rst
diff --git a/examples/06_conditioned_fields/00_condition_ensemble.py b/examples/06_conditioned_fields/00_condition_ensemble.py
diff --git a/examples/06_conditioned_fields/01_2D_condition_ensemble.py b/examples/06_conditioned_fields/01_2D_condition_ensemble.py
diff --git a/examples/06_conditioned_fields/README.rst b/examples/06_conditioned_fields/README.rst
diff --git a/examples/07_transformations/README.rst b/examples/07_transformations/README.rst
diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst
diff --git a/examples/08_geo_coordinates/de_borders.txt b/examples/08_geo_coordinates/de_borders.txt
diff --git a/examples/08_geo_coordinates/temp_obs.txt b/examples/08_geo_coordinates/temp_obs.txt
diff --git a/examples/09_spatio_temporal/01_precip_1d.py b/examples/09_spatio_temporal/01_precip_1d.py
diff --git a/examples/09_spatio_temporal/02_precip_2d.py b/examples/09_spatio_temporal/02_precip_2d.py
diff --git a/examples/09_spatio_temporal/03_geographic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py
diff --git a/examples/09_spatio_temporal/README.rst b/examples/09_spatio_temporal/README.rst
diff --git a/examples/10_normalizer/00_lognormal_kriging.py b/examples/10_normalizer/00_lognormal_kriging.py
diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py
diff --git a/examples/10_normalizer/02_compare.py b/examples/10_normalizer/02_compare.py
diff --git a/examples/10_normalizer/README.rst b/examples/10_normalizer/README.rst
diff --git a/examples/11_plurigaussian/00_simple.py b/examples/11_plurigaussian/00_simple.py
diff --git a/examples/11_plurigaussian/01_pgs.py b/examples/11_plurigaussian/01_pgs.py
diff --git a/examples/11_plurigaussian/02_spatial_relations.py b/examples/11_plurigaussian/02_spatial_relations.py
diff --git a/examples/11_plurigaussian/03_correlations.py b/examples/11_plurigaussian/03_correlations.py
diff --git a/examples/11_plurigaussian/04_3d_pgs.py b/examples/11_plurigaussian/04_3d_pgs.py
diff --git a/examples/11_plurigaussian/05_conditioned.py b/examples/11_plurigaussian/05_conditioned.py
diff --git a/examples/11_plurigaussian/06_periodic.py b/examples/11_plurigaussian/06_periodic.py
diff --git a/examples/11_plurigaussian/07_tree_pgs.py b/examples/11_plurigaussian/07_tree_pgs.py
diff --git a/examples/11_plurigaussian/08_multi_field_pgs.py b/examples/11_plurigaussian/08_multi_field_pgs.py
diff --git a/examples/11_plurigaussian/09_3d_tree_pgs.py b/examples/11_plurigaussian/09_3d_tree_pgs.py
diff --git a/examples/11_plurigaussian/README.rst b/examples/11_plurigaussian/README.rst
diff --git a/examples/11_plurigaussian/conditional_values.npz b/examples/11_plurigaussian/conditional_values.npz
diff --git a/examples/12_sum_model/00_simple_sum_model.py b/examples/12_sum_model/00_simple_sum_model.py
diff --git a/examples/12_sum_model/01_fitting_sum_model.py b/examples/12_sum_model/01_fitting_sum_model.py
diff --git a/examples/12_sum_model/README.rst b/examples/12_sum_model/README.rst
diff --git a/examples/13_mps/00_mps_overview.ipynb b/examples/13_mps/00_mps_overview.ipynb
diff --git a/examples/13_mps/00_simple_unconditional.py b/examples/13_mps/00_simple_unconditional.py
@@ -0,0 +1,56 @@
+r"""
+A first Direct Sampling simulation
+----------------------------------
+
+This is the minimal Multiple Point Statistics example: build a training image,
+wrap it in a :any:`TrainingImage`, and generate one unconditional realization
+with :any:`DirectSampling`.
+
+We use a small, synthetic *channelized* training image generated with NumPy, so
+the example is fast and needs no downloads.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+###############################################################################
+# Create a synthetic binary training image with curvilinear "channels".
+# The two facies (0 and 1) form connected, meandering bands — exactly the kind
+# of structure two-point statistics struggles to reproduce.
+
+gx, gy = np.meshgrid(np.arange(60), np.arange(60), indexing="ij")
+ti_data = ((np.sin(gx / 5.0) + np.sin((gx + gy) / 8.0)) > 0).astype(float)
+
+###############################################################################
+# Wrap the array in a :any:`TrainingImage`. For a categorical variable (facies
+# codes) the distance is the fraction of mismatching neighbours, so the
+# ``distance`` argument is ignored here.
+
+ti = gs.TrainingImage(ti_data, categorical=True)
+print(ti)
+
+###############################################################################
+# Create the :any:`DirectSampling` generator and simulate on a 40x40 grid.
+#
+# * ``n_neighbors`` — how many already-known cells define each data event.
+# * ``scan_fraction`` — fraction of the training image scanned per cell
+#   (smaller is faster, slightly noisier).
+# * ``threshold=0.0`` — the recommended DSBC mode: always take the best match.
+
+ds = gs.DirectSampling(
+    gs.MPSModel(ti, n_neighbors=12, scan_fraction=0.3, threshold=0.0)
+)
+field = ds([np.arange(40, dtype=float)] * 2, seed=20250616)
+
+###############################################################################
+# Plot the training image next to the realization. The realization is not a
+# copy of the TI, but it reproduces the same channel patterns.
+
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))
+ax0.imshow(ti_data, cmap="cividis", origin="lower")
+ax0.set_title("Training image")
+ax1.imshow(field, cmap="cividis", origin="lower")
+ax1.set_title("DS realization")
+fig.tight_layout()
diff --git a/examples/13_mps/01_conditional.py b/examples/13_mps/01_conditional.py
@@ -0,0 +1,58 @@
+r"""
+Conditioning to hard data
+-------------------------
+
+Real simulations must honour measured data (boreholes, samples). Direct
+Sampling supports this through :meth:`~gstools.DirectSampling.set_condition`:
+conditioning values are pinned into the grid before simulation and preserved
+exactly, while the rest of the field is filled in around them.
+
+We reuse the synthetic channel training image from the first example.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+###############################################################################
+# Same synthetic channelized training image as before.
+
+gx, gy = np.meshgrid(np.arange(60), np.arange(60), indexing="ij")
+ti_data = ((np.sin(gx / 5.0) + np.sin((gx + gy) / 8.0)) > 0).astype(float)
+ti = gs.TrainingImage(ti_data, categorical=True)
+
+###############################################################################
+# Draw 40 random "hard data" points and read their facies from the TI. In a
+# real study these would be field measurements; here we sample the TI so the
+# conditioning data are consistent with the patterns.
+
+rng = np.random.default_rng(0)
+cond_x = rng.integers(0, 40, 40).astype(float)
+cond_y = rng.integers(0, 40, 40).astype(float)
+cond_val = ti_data[cond_x.astype(int), cond_y.astype(int)]
+
+###############################################################################
+# Set the conditioning data and simulate. ``set_condition`` snaps each point to
+# its nearest grid node, so the values are honoured exactly at those cells.
+
+ds = gs.DirectSampling(
+    gs.MPSModel(ti, n_neighbors=12, scan_fraction=0.3, threshold=0.0)
+)
+ds.set_condition([cond_x, cond_y], cond_val)
+field = ds([np.arange(40, dtype=float)] * 2, seed=7)
+
+honored = int(
+    (field[cond_x.astype(int), cond_y.astype(int)] == cond_val).sum()
+)
+print(f"conditioning honoured: {honored}/{cond_val.size}")
+
+###############################################################################
+# Plot the realization with the conditioning points overlaid. Every marker sits
+# on a cell whose simulated facies matches the datum.
+
+fig, ax = plt.subplots(figsize=(6, 5))
+ax.imshow(field, cmap="cividis", origin="lower")
+ax.scatter(cond_y, cond_x, c=cond_val, cmap="cividis", edgecolors="red", s=30)
+ax.set_title("Conditional DS realization (red-edged = hard data)")
+fig.tight_layout()
diff --git a/examples/13_mps/02_continuous.py b/examples/13_mps/02_continuous.py
@@ -0,0 +1,62 @@
+r"""
+Continuous variables and distance metrics
+-----------------------------------------
+
+Direct Sampling is not limited to categorical facies: with ``categorical=False``
+it simulates continuous variables (permeability, porosity, elevation, ...). The
+``distance`` argument then selects how two patterns are compared:
+
+* ``"l1"`` / ``"l2"`` — Manhattan / Euclidean distance on the raw values
+  (Mariethoz et al., 2010, Eq. 6 / Eq. 4).
+* ``"variation"`` — compares only the *relative* variations within a pattern
+  (Eq. 9), tolerating a locally varying mean. Useful for non-stationary data.
+
+Here we use a smooth continuous training image and a small acceptance
+``threshold`` (standard DS mode) to allow approximate matches.
+
+.. note::
+
+    The acceptance ``threshold`` is **not comparable across metrics**: the
+    ``"variation"`` distance is normalised by ``2·d_max``, so the same value is
+    stricter than it would be for ``"l1"``/``"l2"``. Re-tune it when you switch
+    the distance metric.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+###############################################################################
+# A smooth, continuous synthetic training image.
+
+gx, gy = np.meshgrid(np.arange(60), np.arange(60), indexing="ij")
+ti_data = np.sin(gx / 6.0) * np.cos(gy / 8.0)
+
+###############################################################################
+# Build a continuous training image with the Euclidean (``"l2"``) distance.
+
+ti = gs.TrainingImage(ti_data, categorical=False, distance="l2")
+print(ti)
+
+###############################################################################
+# Simulate. ``threshold=0.03`` accepts the first pattern within that distance
+# (standard DS), which is faster than the exhaustive best-candidate search for
+# continuous variables.
+
+ds = gs.DirectSampling(
+    gs.MPSModel(ti, n_neighbors=12, scan_fraction=0.3, threshold=0.03)
+)
+field = ds([np.arange(32, dtype=float)] * 2, seed=3)
+
+###############################################################################
+# Plot. The realization reproduces the smooth wavy structure of the TI without
+# copying it. A shared colour scale makes the comparison fair.
+
+vmin, vmax = ti_data.min(), ti_data.max()
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))
+im = ax0.imshow(ti_data, cmap="RdBu_r", origin="lower", vmin=vmin, vmax=vmax)
+ax0.set_title("Training image (continuous)")
+ax1.imshow(field, cmap="RdBu_r", origin="lower", vmin=vmin, vmax=vmax)
+ax1.set_title("DS realization")
+fig.colorbar(im, ax=(ax0, ax1), shrink=0.7)
diff --git a/examples/13_mps/03_channel_strebelle.py b/examples/13_mps/03_channel_strebelle.py
@@ -0,0 +1,85 @@
+r"""
+A real training image: the Strebelle channels
+----------------------------------------------
+
+The previous examples used tiny synthetic training images. Here we use the
+classic **Strebelle (2002) channelized fluvial training image**, the de-facto
+benchmark for MPS, and condition the simulation on random hard data.
+
+.. note::
+
+    **Data source / license.** The training image is downloaded from the
+    `GeoDataSets <https://github.com/GeostatsGuy/GeoDataSets>`_ repository by
+    Michael Pyrcz (GeostatsGuy), which is distributed under the **MIT license**
+    (redistribution permitted with attribution). The underlying channel TI is
+    due to Strebelle, S. (2002), *Conditional simulation of complex geological
+    structures using multiple-point statistics*, Mathematical Geology, 34(1),
+    1-21. If the download is unavailable, the example falls back to a synthetic
+    training image so it still runs offline.
+"""
+
+import os
+import urllib.request
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.colors import ListedColormap
+
+import gstools as gs
+
+###############################################################################
+# Load the Strebelle training image, with a synthetic fallback for offline use.
+
+TI_URL = (
+    "https://raw.githubusercontent.com/GeostatsGuy/"
+    "GeoDataSets/master/MPS_Training_image_and_Realizations_500.npz"
+)
+CACHE = "mps_strebelle.npz"
+try:
+    if not os.path.exists(CACHE):
+        urllib.request.urlretrieve(TI_URL, CACHE)
+    ti_arr = np.load(CACHE)["array1"].astype(float)
+    source = "Strebelle (2002) via GeoDataSets (MIT)"
+except Exception as err:  # pragma: no cover - network fallback
+    print(f"download failed ({err}); using a synthetic channel TI instead")
+    gx, gy = np.meshgrid(np.arange(150), np.arange(150), indexing="ij")
+    ti_arr = ((np.sin(gx / 6.0) + np.sin((gx + gy) / 10.0)) > 0).astype(float)
+    source = "synthetic fallback"
+
+ti = gs.TrainingImage(ti_arr, categorical=True)
+print(f"TI {ti.shape} ({source}), sand fraction = {ti_arr.mean():.3f}")
+
+###############################################################################
+# Take 80 random conditioning points from the training image patterns.
+
+sg_size = 80
+rng = np.random.default_rng(0)
+cond_x = rng.integers(0, sg_size, 80).astype(float)
+cond_y = rng.integers(0, sg_size, 80).astype(float)
+cond_val = ti_arr[cond_x.astype(int), cond_y.astype(int)]
+
+###############################################################################
+# Simulate with DSBC-style parameters (best-candidate + partial scan).
+
+ds = gs.DirectSampling(
+    gs.MPSModel(ti, n_neighbors=30, scan_fraction=0.2, threshold=0.0)
+)
+ds.set_condition([cond_x, cond_y], cond_val)
+field = ds([np.arange(sg_size, dtype=float)] * 2, seed=42)
+
+honored = int(
+    (field[cond_x.astype(int), cond_y.astype(int)] == cond_val).sum()
+)
+print(f"conditioning honoured: {honored}/{cond_val.size}")
+
+###############################################################################
+# Plot the training image crop next to the conditional realization.
+
+cmap = ListedColormap(["#c9a96e", "#2b6cb0"])  # shale / sand
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(11, 5.5))
+ax0.imshow(ti_arr[:sg_size, :sg_size], cmap=cmap, origin="lower")
+ax0.set_title("Training image (crop)")
+ax1.imshow(field, cmap=cmap, origin="lower")
+ax1.scatter(cond_y, cond_x, c=cond_val, cmap=cmap, edgecolors="k", s=18)
+ax1.set_title("Conditional DS realization")
+fig.tight_layout()
diff --git a/examples/13_mps/README.rst b/examples/13_mps/README.rst
@@ -0,0 +1,38 @@
+Multiple Point Statistics
+=========================
+
+Two-point geostatistics (covariance models, kriging, SRFs) describes spatial
+structure through pairs of points and a variogram. This is powerful, but it
+cannot reproduce *curvilinear* or *connected* features such as meandering
+channels, fractures, or other patterns that depend on the joint configuration
+of many points at once.
+
+**Multiple Point Statistics (MPS)** addresses this by learning patterns
+directly from a **training image (TI)** — an example image deemed
+representative of the spatial structure to simulate. Instead of fitting a
+variogram, MPS borrows whole patterns from the TI.
+
+GSTools provides the **Direct Sampling (DS)** algorithm
+(`Mariethoz et al., 2010 <https://doi.org/10.1029/2008WR007621>`_), together
+with the **Direct Sampling Best Candidate (DSBC)** parametrization
+(`Juda et al., 2022 <https://doi.org/10.1016/j.acags.2022.100091>`_), through
+two classes:
+
+* :any:`TrainingImage` — the MPS model: the training image plus the distance
+  used to compare patterns (the analogue of a :any:`CovModel`).
+* :any:`DirectSampling` — the generator that produces realizations on a
+  structured grid (the analogue of :any:`SRF`).
+
+The core idea: to fill each grid cell, DS looks at the values already present
+around it (its *data event*), scans the training image for a location whose
+surroundings look similar enough, and copies that cell's value over.
+"Similar enough" is decided by a **distance** between the two surroundings,
+controlled by three parameters — the number of neighbours ``n``, the scan
+fraction ``f``, and the acceptance threshold ``t`` (with ``t = 0`` giving the
+recommended DSBC mode).
+
+The following tutorials build up from a minimal unconditional simulation to
+conditioning and continuous variables.
+
+Examples
+--------
diff --git a/pyproject.toml b/pyproject.toml