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

examples/13_mps/00_simple_unconditional.py

@@ -0,0 +1,57 @@
+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, n_neighbors=12)
+print(ti)
+
+###############################################################################
+# Create the :any:`DirectSampling` generator and simulate on a 40x40 grid.
+#
+# * ``n_neighbors`` — how many already-known cells define each data event
+#   (set on the :any:`TrainingImage`).
+# * ``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, 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()

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, n_neighbors=12)
+
+###############################################################################
+# 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, 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()

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", n_neighbors=12)
+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, 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)

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, n_neighbors=30)
+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, 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()

examples/13_mps/04_channel_strebelle_nonstat.py

@@ -0,0 +1,71 @@
+r"""
+A non-stationary real training image
+------------------------------------
+
+Here we take the Strebelle channelized fluvial training image and apply a
+spatially varying rotation and scaling factor to show how Direct Sampling
+can generate non-stationary features from a stationary training image.
+"""
+
+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
+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:
+    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, n_neighbors=30)
+print(f"TI {ti.shape} ({source}), sand fraction = {ti_arr.mean():.3f}")
+
+# Grid size for simulation
+sg_size = 80
+xs = np.arange(sg_size, dtype=float)
+ys = np.arange(sg_size, dtype=float)
+gx, gy = np.meshgrid(xs, ys, indexing="ij")
+
+# Create non-stationary rotation and anisotropy (scaling) maps
+# Rotation varies smoothly along the Y-axis from 0 to 45 degrees
+rotation = (gy / sg_size) * (np.pi / 4.0)
+
+# Anisotropy (channel width) varies along the X-axis
+# > 1 means wider channels, < 1 means narrower
+anis = 0.8 + (gx / sg_size) * 0.4
+
+ds = gs.DirectSampling(
+    gs.MPSModel(ti, scan_fraction=0.1, threshold=0.0)
+)
+ds.set_nonstationary(rotation=rotation, anis=anis)
+
+print("Simulating non-stationary field...")
+# Generate the field
+field = ds([xs, ys], seed=42)
+
+# Plot the result
+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.set_title("Non-stationary DS realization\n(varying rotation & scale)")
+fig.tight_layout()
+plt.savefig("nonstat_strebelle.png")
+print("Saved nonstat_strebelle.png")