Geoml

geoML

How to use the geoML package: the object model, the workflow, the recommendations and the geostatistical intuitions its API takes for granted. It says what to build and why. It does not list arguments: the installed package answers that exactly, through its catalogue (§0). Where anything here disagrees with the code, trust the code.

Writing a notebook or example code? Read references/notebook-style.md first, for the import aliases, the workflow arc, the plotting conventions and the tree paths. Worked examples are the manual's chapters, in references/manual/, each run against the package at every release (§4).

0. Looking things up: the catalogue

Every class and function a model or a script may use is described by the installed geoML itself: arguments, defaults, bounds, how a node's size follows from its arguments, which parents it takes, and how stable it is. Ask it rather than guessing, and before writing any constructor call:

python -m geoml.catalogue --show BoxCox          # one entry, as text
python -m geoml.catalogue --show BasicGP --json  # the same, as JSON
python -m geoml.catalogue --list                 # the categories
python -m geoml.catalogue --list warping         # one category's entries

A short name that means two things (Gaussian is a kernel and a likelihood) lists the dotted candidates instead; pass the dotted one. The answer is for the version installed where the command runs, which is the version the code will run against.

Stability. Every entry is one of three:

• public: use freely.
• experimental: works and is tested, but its interface or its recommendation may still change. Say so when you use one.
• internal: kept so that old saved models load. Never build a new model with one; --list leaves them out.

Where geoML is not installed, the same descriptions are on the reference pages of the documentation site, https://italo-goncalves.github.io/geoML/reference/.

1. The package

Install:

pip install git+https://github.com/italo-goncalves/geoML

Backend: TensorFlow 2.x + TensorFlow-Probability, GPU-accelerated. All computation is in float64: geostatistical matrices are ill-conditioned and float32 breaks the Cholesky factorizations. Realizations are stored as float32 (half the disk, every read widened back to float64); geoml.set_realization_dtype("float64") keeps them wide.

License: GPL-3 (dual-licensed; see README). Version: 0.9.0.

Layout: five subpackages (data, latent, math, stats, viz) plus the older plots, around modules left flat on purpose: models, likelihood, kernels, transform, warping, parameter, datasets, metrics, persistence, storage. A saved model records the dotted path of every class in it and imports it by that path when it loads, which is why those modules never move.

1.1 The object model

• Containers hold data and predictions. PointData, the grids (Grid1D/Grid2D/Grid3D), the variable-size block model BlockSet3D, and triangulated meshes (Surface3D, Solid3D). A container holds variables (continuous, vector, compositional, categorical), and each variable holds columns: measurements, a prediction, variances, quantiles, realizations.
• Everything in a container is reached by tree path. container.values("V/prediction") is a numpy array, container.get("V/quantiles/0.5") the column itself (for as_image() or get_contour()), container.tree() the picture of what is there. _metadata/... holds per-location facts the model never reads: hole ids, depths, folds, filters. Paths are in references/notebook-style.md §0.
• A model is a latent network plus one likelihood per variable. models.VGPNetwork is the core. The network is built from latent nodes: an input node carrying the inducing points, GP nodes, and operations that combine them. Its leaves are where the likelihoods attach: pass one leaf per likelihood, or name them together with variables={"Rock": likelihood, ...}. A likelihood carries a warping, the monotone map from the data to the Gaussian scale the model works on.
• Predictions are written into the target container, not returned. model.predict(grid) fills the grid's columns. grid.unpredicted() names what no prediction reached yet, so a cancelled or partial prediction is finished with model.predict(grid, where=grid.unpredicted()).
• Long calls report and can be cancelled. Inside with geoml.progress(callback): training, prediction, refinement, cross-validation and mesh sets report each unit they finish; the callback raising cancels the call.
• Everything trainable is a Parametric holding RealParameters, constrained through their transforms and trained by Adam. Bounds a constructor clips to are in the catalogue.

1.2 Worked example (Jura, categorical)

import geoml
import geoml.latent as gl, geoml.transform as tr
import geoml.kernels as kr, geoml.likelihood as lk

geoml.set_seed(1234)                      # BEFORE constructing anything
jura_train, _ = geoml.datasets.jura()
labels = list(jura_train.get("Landuse").labels)

network_input = gl.BasicInput(
    inducing_points=jura_train, transform=tr.Isotropic(0.5))
network_output = gl.BasicGP(
    parent=network_input, size=len(labels), kernel=kr.Matern32(),
    fix_range=True)

model = geoml.models.VGPNetwork(
    data=jura_train, variables="Landuse", latent_network=network_output,
    likelihoods=lk.CategoricalGaussianIndicator(n_components=len(labels)),
    options=geoml.models.GPOptions(verbose=False))
model.train_full(5)

grid = geoml.data.Grid2D(start=[0, 0], end=[6, 6], n=[21, 21])
model.predict(grid)                       # writes into the container
image = grid.get("Landuse/predicted").as_image()

model.to_dot() renders the network as a diagram. The full modelling arcs are the case-study chapters: Walker Lake (15), Jura with a non-stationary multivariate network (16), and a folded quartz vein in 3D (17).

1.3 Where things live

python -m geoml.catalogue --list CATEGORY gives the classes of each.

1.4 Recommendations

What was measured to work, as of this version. Arguments for each class are in the catalogue.

• Warpings. A single positive grade: BoxCox then ZScore; where the data hold zeros, give BoxCox a shift of the order of the smallest positive value, not the default. A variable centred on zero: YeoJohnson then ZScore. A vector of grades: BoxCox, RobustPCA, ZScore, SinhArcsinh, ZScore, chained with ChainedWarping. These parametric links replaced Spline, which stays for saved models (with backbone="rq" if you use it anyway). One marginal transform is where the benefit stops: stacking rotation-and-spline pairs made held-out scores worse.
• Noise law. For a skewed grade under a log-like link, prefer a Gaussian likelihood: a heavy-tailed law (Laplace, EpsilonInsensitive) pushed back through the link can have an unbounded variance in data units and over-wide intervals, though its bound may look better. For data with gross errors, likelihood.Mixture (experimental) names the bad readings, with its warping led by ZScore(robust=True).
• Several populations. Where one variable is drawn from populations that differ in level or skew, ore and waste say, use likelihood.LikelihoodMixture over continuous likelihoods, each with its own warping. With shares="latent" the shares change from place to place: read them from a GP node of their own, apart from the populations', and where a domain was logged let a CategoricalGaussianIndicator with bias=True, trained on every logged interval rather than the assayed ones alone, read the same node. On one deposit it beat a single likelihood on every metal out of fold, but its intervals came out narrower than the held-out data supported; check them.
• Inducing points. The data's own locations plus a regular backbone, divided into overlapping experts: data.inducing.experts over data.inducing.combine. Where the survey does not fill its box -- a fan of drillholes -- take the backbone from data.inducing.from_hull, which keeps the lattice inside the data's hull and a margin around it. Experts come out compact and share points with every neighbour once each borrows at least as many points as it has neighbours: keep experts large enough for their overlap (an overlap of 0.1 on 25-point experts cannot reach six neighbours). How many points a model can absorb depends on the whole configuration; measure it on held-out data rather than carrying a number across problems.
• Structure. One variable influencing another: Linear off the first's field, added to the second's GP through LinearCombination, with unit_norm=False so the model can decline the influence. A field that is not stationary: GPWalk moves the coordinates and a stationary kernel reads the moved space. Build the GP reading the walk with isotropic=True: a range per direction there is a second description of the anisotropy, which belongs in the input's transform. Independent parts of a model can sit on independent trees, one leaf each. Chapter 16 builds the first two, one leaf per variable.
• Depth. Since 0.9.0 a GP node reading another node's uncertain output averages its kernel over it, the covariance between locations included (the expected kernel, GPOptions(propagation="joint"), the default); a model saved before keeps the old rule, which made deep networks overconfident (propagation="marginal"). Build a deep network with the coordinates concatenated beside the inner GP (Concatenate(root, inner)). A GP on an uncertain input takes the Gaussian, exponential, Matérn or rational quadratic kernel; spherical and cubic are refused there, and the error message says so. On a GaussianInput a BasicGP takes the mixture's moments (the variance of the mean over the input included) and the likelihood integrates what its realizations leave out; UncertainInputGP is deprecated. A location error is still better left untold: a noise term absorbs it. That variance costs time and memory in proportion to the inducing points; experts do not shorten it, train_by_expert holds a fraction of it in memory.
• Training. GPOptions(training_tolerance=0.01) stops once the bound has settled; the last few percent of the bound buys sharpness held-out data does not support.
• Many experts. When memory grows with their number, train and predict an expert at a time: model.train_by_expert, model.predict_by_expert, and refine(..., by_expert=True) for a block model. It matched or beat train_svi on CRPS and on a deposit's rock types in a third to a half of the device memory, a little behind on rmse on a smooth synthetic field. Fewer, larger experts scored better than many small ones at a fixed total of inducing points. A deep GP fed by another needs the coordinates concatenated into its input to stay local; one fed by a GP alone spreads every expert over the field. Five nodes are refused; the message names them.
• Validation. PointData.spatial_k_fold writes folds that resemble the real prediction task, models.cross_validate scores the model out of fold, and models.conformalize recalibrates the intervals. In-sample scores flatter; only out-of-fold ones speak for the ground between samples.
• Block models. A BlockSet3D refined by models.refine: predict coarse, split only the blocks the prediction's surface runs through (a block the model is merely unsure about is left whole), predict what the split made. MeshSet contours a column for the prediction and every realization at once, with volumes measured as it goes.