Repository navigation
Conversation
LSchueler
reviewed
Oct 7, 2026
LSchueler
left a comment
LSchueler
left a comment
Member
There was a problem hiding this comment.
Really impressive work, thank you so much!
When the issues below together with the other ones I marked directly in the code are fixed, I think this is ready to be merged.
-
Every postprocessing-pass rebuilds the full dependency graph, even when the visit order is the same each pass ("same", "sequential", and explicit arrays). I'm not sure how slow this is and if it's worth optimising. What do you say?
-
Don't use ipynb files for the mps_overview example.
-
I don't think
TrainingImageshould be exported to the top level, I would keep it next toVariablee.g. or does it make the mps-usage more ergonomic? -
Remove references to your internal docs from the code, a few examples:
"DS_Feature_Checklist §3.3", "spec §3", "bug B5" -
I found a few edge cases, I think it would be the most effective when I just paste code snippets for you to reproduce them.
univariate TI built from a named Variable crashes:
import numpy as np, gstools as gs
ti = gs.TrainingImage(gs.mps.Variable("f", np.zeros((9, 9))))
gs.DirectSampling(gs.MPSModel(ti))([range(3)] * 2)default max_radius=None fails on a 100^3 grid:
import numpy as np, gstools as gs
ti = gs.TrainingImage(np.zeros((9, 9, 9)))
gs.DirectSampling(gs.MPSModel(ti))([range(100)] * 3)inconsistent rotation:
import gstools as gs
from gstools.mps.nonstationary import resolve_spec
print(resolve_spec(0.5, (1, 1, 1), 3, None, "rotation")[0][0, 0, 0])
print(gs.tools.set_angles(3, 0.5))a nan-threshold is accepted:
import numpy as np, gstools as gs
print(gs.MPSModel(gs.TrainingImage(np.zeros((9, 9))), threshold=np.nan).threshold)This exception should be caught earlier not only during the simulation:
import numpy as np, gstools as gs
from gstools.mps import Variable
a = Variable("a", np.zeros((9, 9)))
ti = gs.TrainingImage([a, Variable("b", np.zeros((9, 9)))])
m = gs.MPSModel(ti, post_processing=1, post_processing_variables=[a.replace(weight=2.0)])
gs.DirectSampling(m)([range(3)] * 2)a categorical override on a cont. var. is accepted without a warning
import numpy as np, gstools as gs
from gstools.mps import Variable
ti = gs.TrainingImage([Variable("c", np.random.rand(9, 9), categorical=False)])
m = gs.MPSModel(ti, post_processing=1,
post_processing_variables=[Variable("c", np.random.rand(9, 9), categorical=True)])
print(gs.DirectSampling(m, seed=1)([range(3)] * 2)["c"])the new penalty_matrix is not checked against the conditioning values, the exception should be thrown much earlier
import numpy as np, gstools as gs
ti = gs.TrainingImage(np.eye(9))
pv = ti.variable().replace(penalty_matrix=[[0, 1], [1, 0]])
ds = gs.DirectSampling(gs.MPSModel(ti, post_processing=1, post_processing_variables=[pv]))
ds.set_condition([[0], [0]], [2])
ds([range(3)] * 2)the pp-pass applies the primary TI's mean-shift instead of the new ones distance (output should be 1, not 0)
import numpy as np, gstools as gs
data = np.random.default_rng(0).random((20, 20))
ti = gs.TrainingImage(data, categorical=False, distance="variation")
pv = ti.variable().replace(distance="l1")
f = gs.DirectSampling(gs.MPSModel(ti, post_processing=1, post_processing_variables=[pv]), seed=1)([range(10)] * 2)
print(np.isin(f, data).mean())The pp-pass drops variation to 1 neighbour, which makes every candidate distance 0
import numpy as np, gstools as gs
from gstools.mps.distance import vec_variation_dist
ti = gs.TrainingImage(np.random.rand(9, 9), categorical=False, distance="variation", n_neighbors=2)
gs.MPSModel(ti, post_processing=1, post_processing_factor=2) # post-pass: round(2/2) = 1 neighbour
print(vec_variation_dist(np.array([5.0]), np.array([[0.0], [0.7]]), np.ones(1), 1.0))
Comment on lines
+239
to
+240
| __all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"] | ||
| __all__ += ["transform", "normalizer", "config", "mps"] |
Member
There was a problem hiding this comment.
delete these lines and only add "mps" to the list on line 250
Comment on lines
+410
to
+414
| def _snap_to_nearest(self, axes, k): | ||
| """Return (grid_index_tuple, dist_sq) for the k-th conditioning point.""" | ||
| idx = tuple( | ||
| int(np.argmin(np.abs(axes[d] - self._cond_pos[d][k]))) | ||
| for d in range(self.dim) |
Member
There was a problem hiding this comment.
I think if the out-of-grid conditioning points snap to the edge of the grid, at least a warning should be printed.
| 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 |
Member
There was a problem hiding this comment.
the threshold part seems outdated, it's threshold=None, threshold=0 gets rejected
| ---------- | ||
| model : MPSModel | ||
| Algorithm configuration: training image, scan fraction, threshold, | ||
| conditioning weight, and boundary strategy. Neighbour count and |
Member
There was a problem hiding this comment.
conditioning weight lives in Variable
| patterns (the analogue of a :any:`CovModel`). Multivariate training images | ||
| are built from a list of :any:`Variable` objects, one per property. | ||
| * :any:`MPSModel` — the algorithm configuration: a training image together | ||
| with the scan fraction, threshold, conditioning weight, and the optional |
Member
There was a problem hiding this comment.
conditioning weight lives in Variable
n0228a
added a commit
to n0228a/GSTools
that referenced
this pull request
Oct 14, 2026
- R1: NumPy CSR dependency graph; post passes reuse neighbour lists - R2/R3/R4/R14: drop notebook, TrainingImage only in gs.mps, remove internal refs, fix __all__ - R5-R8, R11, R15: named univariate variables, per-axis offset bounds, set_angles rotation, non-finite/penalty-code validation, out-of-grid warning - R9/R10/R12/R13: coherent post-processing overrides (unset fields inherit, per-pass kernel and mean-shift, variation < 2 post neighbours raises) - NaN/inf TI handling, categorical-distance and max_radius checks, explicit path/node seeds no longer collide, set_condition keeps state on error - examples 05/06/08 run again; README/docstring fixes; regression tests Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
n0228a
force-pushed
the
mps-direct-sampling
branch
from
October 14, 2026 13:18
41472d4 to
c6f7ee2
Compare
- R1: NumPy CSR dependency graph; post passes reuse neighbour lists - R2/R3/R4/R14: drop notebook, TrainingImage only in gs.mps, remove internal refs, fix __all__ - R5-R8, R11, R15: named univariate variables, per-axis offset bounds, set_angles rotation, non-finite/penalty-code validation, out-of-grid warning - R9/R10/R12/R13: coherent post-processing overrides (unset fields inherit, per-pass kernel and mean-shift, variation < 2 post neighbours raises) - NaN/inf TI handling, categorical-distance and max_radius checks, explicit path/node seeds no longer collide, set_condition keeps state on error - examples 05/06/08 run again; README/docstring fixes; regression tests
# Conflicts: # CHANGELOG.md
Contributor
Author
|
All points are addressed, plus several issues found while fixing them: Review points
Additional fixes
|
OpenMP Benchmark ResultsBase: Full ASV comparison |
…eometric unchanged
jeilealr
reviewed
Oct 19, 2026
| array is constant or entirely undefined (avoids a zero/NaN normalizer). | ||
| """ | ||
| finite = arr[np.isfinite(arr)] | ||
| dmax = float(finite.max() - finite.min()) if finite.size else 0.0 |
Contributor
There was a problem hiding this comment.
Potential integer overflow in _data_range():
When the training image uses a signed integer dtype, the subtraction is performed before the conversion to float, so it can overflow.
For example:
import numpy as np
a = np.array([-128, 127], dtype=np.int8)
print(float(a.max() - a.min())) # -1.0 (overflow)
print(float(a.max()) - float(a.min())) # 255.0 (correct)one solution could be:
dmax = float(finite.max()) - float(finite.min())You could also change the type in src/gstools/mps/training_image.py, around lines 272–273 by adding if not categorical: and inside the if data_arr = data_arr.astype(np.float64), but then that could lead to memory allocations and maybe something else. so I prefer the small change in the solution I gave
jeilealr
reviewed
Oct 19, 2026
| "DirectSampling: cond_pos and cond_val length mismatch." | ||
| ) | ||
| _check_finite_positions(cond_pos_arr) | ||
| self._cond_pos = cond_pos_arr |
Contributor
There was a problem hiding this comment.
set_condition() modifies the object's stored conditioning positions before it has finished validating and converting the new conditioning values
self._cond_pos is immediately replaced with the new positions. Then, Python attempts to convert all the new values to float64; If one variable contains an invalid value, np.asarray(..., dtype=np.double) raises an exception. However, the update to self._cond_val never completes, so it still contains the old values.
Simply change the order:
new_cond_val = {
v: np.asarray(a, dtype=np.double) for v, a in cond_val.items()
}
self._cond_pos = cond_pos_arr
self._cond_val = new_cond_val
jeilealr
reviewed
Oct 19, 2026
| float64, matching the naive expression's own promotion; only an | ||
| integer-coded TI needs the explicit widen. | ||
| """ | ||
| diff = np.subtract(data_event_sim, all_de_ti) |
Contributor
There was a problem hiding this comment.
The issue is the same ordering problem I commented before for the direct_sampling.py file: integer subtraction happens before the conversion to float64.
I noticed that In training_image.py lines 1002–1007 already converts continuous data events to float64 before calling the distance functions, so the regular simulation path should be safe.
Still, would it make sense to move the integer-to-float conversion before np.subtract() here as well, to make the low-level function robust when called directly or reused elsewhere?
jeilealr
reviewed
Oct 19, 2026
| cond_pos_arr = np.asarray(cond_pos, dtype=np.double).reshape( | ||
| self.dim, -1 | ||
| ) | ||
| n_cond = len(next(iter(cond_val.values()))) |
Contributor
There was a problem hiding this comment.
len() only checks the size of the first dimension. should each variable have one scalar conditioning value per position?
If this makes sense to you then you can do something like
n_cond = cond_pos_arr.shape[1]
for v, arr in cond_val.items():
arr = np.asarray(arr, dtype=np.double)
if arr.shape != (n_cond,):
raise ValueError(
f"Conditioning values for {v!r} must have shape "
f"({n_cond},), got {arr.shape}."
)
jeilealr
reviewed
Oct 19, 2026
| "Either specify weight on every Variable or on none." | ||
| ) | ||
| else: | ||
| total = sum(raw) |
Contributor
There was a problem hiding this comment.
I think here normalising very large variables weights can lead to numerical issues. A numerically safer approach is to scale the weights by their maximum before summing
jeilealr
reviewed
Oct 19, 2026
| def _snap_to_nearest(self, axes, k): | ||
| """Return (grid_index_tuple, dist_sq) for the k-th conditioning point.""" | ||
| idx = tuple( | ||
| int(np.argmin(np.abs(axes[d] - self._cond_pos[d][k]))) |
Contributor
There was a problem hiding this comment.
Here, np.argmin() can select an invalid grid coordinate if an axis contains NaN.
For example, with axes[d] = [0.0, np.nan, 2.0] and a conditioning position of 0.0, np.argmin(np.abs(axes[d] - 0.0)) returns index 1 instead of 0.
Although set_condition() already validates that conditioning positions are finite, I couldn't find an equivalent check for the simulation-grid axes. This could result in hard-conditioning values being assigned to the wrong grid nodes.
Would it make sense to validate that all simulation-grid coordinates are finite at the beginning of DirectSampling.__call__(), before snapping the conditioning points?
Avoids 16-bit seed collisions in path/node streams. Default RNG unchanged.
Widen ints before subtracting, validate cond_val shape and assign atomically, reject non-finite grid axes, rescale overflowing weights.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds a new
gstools.mpssubpackage 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 formps; the per-node TI scan isvectorised 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.DirectSamplingandMPSModelarealso available at the top-level
gstoolsnamespace, likeSRFandCovModel;TrainingImage,Variable, andZonelive ingstools.mps.VariableHolds one variable's data and its per-variable configuration:
data— n-D NumPy array; NaN marks undefined cells, ±inf is rejectedcategorical— whether to use weighted mismatch distance (Mariethoz2010 Eq. 3);a continuous
distancerequirescategorical=Falsedistance— continuous metric:"l1"(Juda2022 Eq. 7, default),"l2","l<p>"(arbitrary Minkowski),"variation"or"variation<p>"(Mariethoz2010Eq. 9, with automatic mean-shift correction)
n_neighbors— data-event size (maximum neighbours, default 32)max_radius— optional Euclidean cap on neighbour selection, in grid cells (>= 1)weight— relative weight in the multivariate joint distance;Nonefor equalshare; weights are automatically normalized and do not need to sum to 1
penalty_matrix— optional(C, C)per-category mismatch costT[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. Validatedto 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— theseconstraints are what keep the distance in
[0, 1](there is nod_maxnormalization on this path)
distance_power— spatial-decay exponent δ for this variable's neighbourweighting (Mariethoz2010 Eq. 5);
0.0(default) gives uniform weightscond_weight— weight multiplier for this variable's conditioning nodes inthe distance (Meerschman2013 §6); default
1.0replace(**overrides)returns a newVariablewith the given search-configfields overridden, sharing the original's data buffer — no copy. This is the
cheap way to sweep search parameters, or to build a
post_processing_variablesoverride (below) without re-copying a large training image.
TrainingImageThe MPS analogue of
CovModel. Three construction modes:window(slices)returns a sub-region view of a TI, carrying over every searchparameter 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.
MPSModelConfiguration object bundling a
TrainingImagewith the search parameters andevery optional feature. Only
tiis positional; all options are keyword-only:scan_fraction,threshold,boundaryscan_path"sequential"(default) or"random"rotation,scalezonesZoneobjects for zonated simulationpost_processing,post_processing_factor,post_processing_path,post_processing_variables,post_processing_maskGeometric 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. Rotation angles follow
gstools.tools.set_angles. A stationary TI can therefore producepatterns that rotate and stretch across the grid without modifying the TI itself.
rotationis rejected on 1-D grids (no rotation planes exist);scalestaysmeaningful.
Post-processing (
post_processing=p, Meerschman et al. 2013 §4) runsppasses 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=1for categorical simulations.post_processing_factor(p_f > 0) dividesscan_fractionandn_neighborsduring the passes.p_f > 1trades fidelity forspeed — Me13's documented intent, which finds little effect and recommends
p_f=1.p_f < 1does the reverse, spending more search effort per post-pass than the mainpass; this direction has no literature precedent and is offered for experimentation.
Passes are Gauss–Seidel, but that ordering is carried by the same dependency DAG the
main path uses, so they honour
num_threadswhile staying bit-identical to serialexecution at any thread count. Each pass re-visits essentially the whole grid, so
p > 0dominates 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_pathcontrols the post-pass visit order and is animplementation extension — Me13 §4 does not prescribe one.
None(default)inherits the main-path mode;
"random","sequential","same"(reuse the mainpass's order), or an explicit
(N, dim)array are accepted.post_processing_variables(a list ofVariableoverrides, matched by nameagainst the primary TI's variables) gives post-processing passes their own search
settings — build one via
variable.replace(...). An override keeps the primaryvariable's data and categorical kind; its
distance,penalty_matrix,distance_power,cond_weight,n_neighbors,max_radius, andweightall applyin the passes, and fields not set on the override inherit the primary's value.
Weights must be set on every variable or on none. A
distance="variation"variableneeds at least 2 post-pass neighbours, checked at construction. Variables without an
override reuse the primary's search config; default
[]is bit-identical to omittingthe parameter.
post_processing_mask(a boolean, simulation-grid-shaped array) restricts whichnodes are eligible for re-simulation during post-processing;
None(default)re-simulates every non-conditioning node. Conditioning nodes are never
re-simulated regardless of the mask.
ZoneZone(ti, where)binds a training image to a region of the simulation grid forzonated (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.
whereis a boolean array of the grid shape, or a callablef(x, [y, z, ...]) -> bool mask. Non-boolean arrays raise; there is no silent!= 0coercion, so an integer zonation raster is writtenZone(ti_k, where=(labels == k)).d_max).Every search hyper-parameter and the distance kernel —
n_neighbors,max_radius,distance,weight,penalty_matrix— always come from theprimary TI, so setting them on a zone TI has no effect (during post-processing
passes, the kernel instead comes from any
post_processing_variablesoverride — main-pass zone behaviour is unaffected).
deterministic and thread-count invariant.
DirectSamplingSubclasses
gstools.field.base.Fieldand follows the same call interface asSRF:ds(pos, seed, *, ...), with every further option keyword-only. Takes anMPSModel.All search parameters (
scan_fraction,threshold,boundary,n_neighbors,max_radius,cond_weight) are read-only onDirectSampling— theyreflect the underlying
MPSModelandVariableobjects.n_neighbors,max_radius, andcond_weightare per-Variableand collapse to a scalarwhen all variables share the same value, otherwise return a
{name: value}dict. To change them, mutate the
MPSModelorVariabledirectly.Key capabilities:
"sequential"or an explicit(N, dim)visit-order array via thepathargument
scan_fractioncaps the per-node TI scan (Juda2022 §2)threshold=None— DSBC mode: full scan, best-candidate selectionboundary="strict"(default) or"partial"— partial drops lags that exceedthe TI extent and searches with the reduced template (Mariethoz2010 §6.2)
set_condition(cond_pos, cond_val)— hard-data conditioning with nearest-nodesnapping (points outside the grid warn) and per-variable collision resolution;
cond_valaccepts a{variable: array}dict for multivariate TIs, withnp.nanto leave avariable unconditioned at a point; conditioning weight is set per
Variable(
cond_weight=)set_mv_transforms()— per-variable mean / normalizer / trend post-processing.Note this is the
Fieldpipeline and is unrelated toMPSModel's DSpost_processingpassesnum_threads(orgstools.config.NUM_THREADS); outputis identical for any thread count, speed-ups are modest because the engine is pure Python
path_seed,node_seed) for independent control of visitorder vs. TI scan randomness; distinct explicit seeds give distinct streams
{variable: ndarray}dict; each variable alsostored as a named field accessible via
ds["facies"]/ds.all_fieldsUsage
Multivariate TIs are built from a list of
Variableobjects and return a{variable: ndarray}dict;rotation/scale,zones, andpost_processingareall
MPSModelarguments layered on the same call. See the examples below for each.Module layout
training_image.pyVariable,TrainingImage— data and distance functionsdata_event.pyDataEvent— internal value object bundling per-variable data-event arraysdistance.pymodel.pyMPSModel— validated search-parameter configzone.pyZone— a training image bound to a region of the gridneighbors.pyscan.pyrunner.pynonstationary.pysimulate.py_DirectSamplingEngine+ds_simulateentry pointdirect_sampling.pyDirectSampling— publicFieldsubclassExamples (
examples/13_mps/)examples/13_mps/is wired intodocs/source/tutorials.rstandsphinx_gallery_conf, so these run at documentation build time. Several fetchtheir training images from public sources at first run and cache them next to the
script.
00_simple_unconditional.py01_conditional.pyset_condition()02_continuous.py03_channel_strebelle.py04_channel_strebelle_nonstat.pyMPSModel(rotation=..., scale=...)05_channel_strebelle_radial.py06_channel_multivariate_fig8.py07_variation_distance_fig6.py08_bivariate_joint_fig4.py09_continuous_stone_fig3.pyReferences
to perform multiple-point geostatistical simulations. Water Resources Research.
https://doi.org/10.1029/2008WR007621
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
the Direct Sampling algorithm for multiple-point statistical simulations.
Applied Computing and Geosciences.
https://doi.org/10.1016/j.acags.2022.100091
using multiple-point statistics. Mathematical Geology.