geoml.models

class geoml.models.GP(data, variable, covariance, warping=None, directional_data=None, interpolation=False, use_trend=False, options=None)

Bases: _GPModel

Basic Gaussian process model.

__init__(data, variable, covariance, warping=None, directional_data=None, interpolation=False, use_trend=False, options=None)

Basic Gaussian process model.

This model is based on the standard Gaussian process for a single output variable. It supports warping for non-Gaussian variables and directional data as gradients of the modelled field, but not both simultaneously.

Parameters:

• data – A PointData object from the ´data´ module.

• variable (str) – The name of the variable to be modelled. Must be a continuous variable.

• covariance – The covariance function to build the covariance matrices.

• warping – An object from the warping module. If None, the data is assumed to have zero mean and unit variance.

• directional_data – A DirectionalData object from the ´data´ module. The corresponding variable will be used as the gradient of the modelled field.

• interpolation (bool) – If True, will assume that the data is noiseless and try to honor the data points.

• use_trend (bool) – If True, will model a linear trend in the data in addition to the GP.

• options (GPOptions) – Additional configurations.

set_learning_rate(rate)

Resets the model’s optimizer with the provided learning rate. Will erase the optimizer’s memory.

Parameters:

rate (float) – The learning rate to use.

refresh(jitter=1e-09)

Updates the model’s internal state.

If called within TensorFlow’s eager mode, will allow inspection of the internal tensors.

Parameters:

jitter (float) – Small value added to the covariance matrices for numerical stability.

log_likelihood(jitter=1e-09)

Computes the model’s log-likelihood with the current parameters.

Parameters:

jitter (float) – Small value added to the covariance matrices for numerical stability.

predict_raw(x_new, jitter=1e-09, n_sim=50)

predict(newdata, n_sim=50)

Makes a prediction on the specified coordinates.

Parameters:

newdata – A reference to a spatial points object of compatible dimension. The object’s variables will be updated.

train(max_iter=1000)

Model training.

The standard GP does not support batches of data, allways using the full data instead. This is feasible for up to a few thousand data points.

Parameters:

max_iter (int) – The number of iterations to train.

class geoml.models.GPEnsemble(data, variable, covariance, warping=None, directional_data=None, use_trend=False, options=None)

Bases: _EnsembleModel

An ensemble of Gaussian processes.

__init__(data, variable, covariance, warping=None, directional_data=None, use_trend=False, options=None)

An ensemble of Gaussian processes.

This model combines independent GPs into a consolidated prediction using the Product of Experts approach. It is preferable to divide the data spatially instead of randomly, so that each expert can focus on a specific region of the space.

Parameters:

• data – A list or tuple of PointData objects.

• variable (str) – The name of the variable to be modelled. Must be present in all data objects.

• covariance – The covariance function to build the covariance matrices.

• warping – An object from the warping module. If None, the data is assumed to have zero mean and unit variance.

• directional_data – A DirectionalData object from the ´data´ module. The corresponding variable will be used as the gradient of the modelled field.

• use_trend (bool) – If True, will model a linear trend in the data in addition to the GP.

• options (GPOptions) – Additional configurations.

train(max_iter=1000)

class geoml.models.Normalizer(warping, options=None)

Bases: _GPModel

Trainable data normalizer.

__init__(warping, options=None)

Trainable data normalizer.

This model will fit a warping object to a data vector, allowing its transformation to a Gaussian distribution with zero mean and unit variance.

Parameters:

• warping – An object from the warping module.

• options (GPOptions) – Additional configurations.

normalize(x, max_iter=250)

Model training.

Parameters:

• x (array-like) – The data vector to train on.

• max_iter (int) – The number of iterations to run.

class geoml.models.StructuralField(tangents, covariance, normals=None, mean_vector=None, options=None)

Bases: _GPModel

Structural field modeling based on gradient data.

set_learning_rate(rate)

Resets the model’s optimizer with the provided learning rate. Will erase the optimizer’s memory.

Parameters:

rate (float) – The learning rate to use.

refresh(jitter=1e-09)

log_likelihood(jitter=1e-09)

train(max_iter=1000)

predict_raw(x_new, jitter=1e-09)

predict_raw_directions(x_new, x_new_dir, jitter=1e-09)

predict(newdata, variable)

Makes a prediction on the specified coordinates.

Parameters:

• newdata – A reference to a spatial points object of compatible dimension. The object’s variables are updated.

• variable (str) – Name of output variable.

class geoml.models.GPOptions(verbose=True, prediction_batch_size=20000, jitter=1e-09, training_batch_size=2000, training_samples=20, jit_predict=False, qmc_simulations=False, expert_propagation='consensus', training_tolerance=None)

Bases: _ModelOptions

__init__(verbose=True, prediction_batch_size=20000, jitter=1e-09, training_batch_size=2000, training_samples=20, jit_predict=False, qmc_simulations=False, expert_propagation='consensus', training_tolerance=None)

Configuration of Gaussian process models.

This object can be passed on to models based on the Gaussian process in order to control their behavior.

The seed that training and the simulations read is not a parameter: it is drawn from the package generator when the object is built, so geoml.set_seed before construction governs it the way it already governs parameter initialization, and a saved model keeps the number it drew.

Parameters:

• verbose – Whether to show the training process on screen.

• prediction_batch_size (int) – Batch size for prediction/inference.

• jitter (float) – Small value added to covariance matrices for numerical stability.

• training_batch_size (int) – Number of data points per batch during training.

• training_samples (int) – Number of Monte Carlo samples to be drawn when training requires it.

• jit_predict (bool) – Whether to compile the prediction graph with XLA. Worth 3-5x on a grid of any size, and more on a GPU, but XLA compiles per distinct batch shape (so a grid that prediction_batch_size does not divide pays that cost twice) and refuses to run anything it cannot compile, rather than falling back. Prediction only: the same treatment makes training both slower and unstable.

• qmc_simulations (bool) – Whether to draw the posterior simulations from a seeded-scramble Sobol sequence instead of pseudo-random normals, so the same number of them covers the predictive distribution evenly rather than by chance. Measured on the Walker Lake model at 16-256 simulations: the ensemble mean lands 7-37x closer to the exact posterior mean, proportions below a cut-off and the outer quantiles about a quarter closer – the accuracy of half again to twice the simulations – while the ensemble’s own spread and the correlation between locations gain nothing, and the cost is not measurable. Deterministic given seed, batch-invariant either way. The sequence covers at most 21201 dimensions (size times the inducing points of a node), beyond which scipy refuses.

• expert_propagation (str) – How a deep network’s experts see each other’s inducing sets, for every BasicGP in the network at once. “consensus” (the default, and the historical behavior) predicts every expert’s set from every other and combines by precision weighting – O(K^2) in the expert count. “independent” lets each expert speak for its own set alone – O(K) – so duplicated points in overlapping sets may disagree, and the data-side weighting arbitrates. Measured (Walker deep model and a 3000-point synthetic, K = 5-40): training 1.6x to 6.3x faster and prediction up to 8x as K grows; quality within a few percent of consensus and sometimes ahead, the consensus coupling appearing to slow optimization at large K. Only deep (multi-layer) networks are affected: below a terminal node the propagation never runs.

• training_tolerance (float, optional) –

  When to stop training before its iteration count runs out, as a fraction. The bound is smoothed, and training stops once its gain over the last twenty iterations falls below this share of everything gained since the call began. None (the default) trains for exactly as long as it is told, which is what every version before 0.6.5 did. 0.01 is a reasonable setting.

  The count passed to train_full/train_svi remains the cap, and the criterion only ever ends training sooner. It is deliberately unable to stop a model that begins the call already converged: with nothing gained, there is nothing to take a fraction of, and the run goes to its cap. Training in phases is the pattern this protects – a smaller learning rate makes progress the previous phase could not, and each phase is judged against its own starting point.

jit_predict = False

qmc_simulations = False

expert_propagation = 'consensus'

training_tolerance = None

class geoml.models.VGPNetwork(data, variables, likelihoods, latent_network, directional_data=None, options=None)

Bases: _GPModel

Variational Gaussian process network.

A generalization of the standard Gaussian process: variables of any kind (continuous, categorical, compositional, directional) are modelled through a network of latent Gaussian processes, fitted by maximizing the evidence lower bound on inducing points.

Parameters:

• data (_data.PointData) – The training data, a container from the geoml.data module.

• variables (str | Sequence[str]) – The name of a variable in data to model, or a list of names.

• likelihoods (_lk._Likelihood | Sequence[_lk._Likelihood]) – A likelihood from geoml.likelihood, or a list of one per variable, in the same order as variables.

• latent_network (_latent.network._LatentVariable) – The network’s terminal node, from geoml.latent. Its size must match the likelihoods’ sizes summed.

• directional_data (_data.DirectionalData | None) – Structural measurements, whose variable is taken as the gradient of the modelled field.

• options (GPOptions | None) – Training and prediction settings.

training_log

The evidence lower bound at each iteration of the last training run.

Type:

list of float

Examples

>>> import geoml
>>> geoml.set_seed(1234)
>>> walker, grid = geoml.datasets.walker()
>>> inducing = geoml.data.inducing.from_kmeans(walker, 100, seed=0)
>>> gp = geoml.latent.BasicGP(
...     geoml.latent.BasicInput(inducing), size=1)
>>> model = geoml.models.VGPNetwork(
...     walker, "V", geoml.likelihood.Gaussian(), gp)
>>> model.train_full(max_iter=100)
>>> model.predict(grid, n_sim=20)

to_dot(legend=True, rankdir='BT')

Writes the model as a Graphviz diagram: coordinates, latent network, warpings and output variables.

See geoml.graphviz.to_dot.

set_learning_rate(rate)

Resets the model’s optimizer with the provided learning rate. Will erase the optimizer’s memory.

Parameters:

rate (float) – The learning rate to use.

train_full(max_iter=1000)

Train on the whole data set at every iteration.

Feasible while the data and the latent network fit in memory together; past that, use train_svi(). The evidence lower bound of each iteration is appended to training_log.

Parameters:

max_iter (int) – Number of iterations, and a cap rather than a count when options.training_tolerance asks training to stop once the bound settles.

See also

train_svi

minibatch training, for larger data sets.

train_svi(epochs=100)

Train in minibatches, by stochastic variational inference.

Each epoch visits the data once, in batches of options.training_batch_size, drawn in an order reproducible from the model’s seed. The mean bound over an epoch’s batches is appended to training_log.

Parameters:

epochs (int) – Number of passes over the data, and a cap rather than a count when options.training_tolerance asks training to stop once the bound settles. The criterion reads one value an epoch, the mean over its batches.

See also

train_full

one gradient step per iteration on all the data.

predict_raw(*args, **kwargs)

_predict_raw in a graph, XLA-compiled when the options ask for it.

jit_compile is settled when a tf.function is built, and the simulation draws are baked into the graph, so honouring either option means holding one function per combination of settings rather than a flag on a single one. Each is traced at most once per model, and None (rather than False) leaves the uncompiled path exactly as it was. The simulation_rule/propagation_rule contexts wrap the call rather than the trace because a retrace (a new batch shape, a new n_sim) can happen on any call, and has to see the flags the cache key promised. The propagation rule joins the key for the networks whose graphs refresh internally (ProjectedVGP); on this class the graph reads snapshotted state, and the extra key is merely unused.

predict(newdata, n_sim=20, include_noise=True, where=None)

Predict at new locations, writing the answer into the container.

The variables the model was trained on are created on newdata if absent and filled in place: prediction, latent moments, simulations, and the columns each variable kind adds to those.

Parameters:

• newdata (_SpatialData) – The locations to predict, of the same dimension as the training data. Modified in place.

• n_sim (int) – Number of realizations to draw per location.

• include_noise (bool) – Whether to integrate the likelihood noise out of the answer. The prediction then reports the value the ground would show once measurement error and sub-resolution variability are averaged over – a deterministic correction, and a large one on a skewed variable. Turn it off to see the latent field alone.

• where (ndarray | Sequence[int] | str | None) – One boolean per location, the indices of the locations to visit, or None for all of them. Locations left out keep whatever they hold, including their simulations, and a location never visited stays missing.

Notes

A location’s simulated values do not depend on what else is in its batch, so predicting a subset gives the same answer as predicting everything and reading that subset back.

predict_measurements(newdata, n_sim=20, n_nodes=32)

Draw the predictive distribution of a measurement per location.

predict() reports the ground, with the likelihood noise integrated out, so its simulations describe a quantity no sample observes. This keeps the noise instead, giving n_sim * n_nodes equally likely readings per location – the distribution to compare against measured data, as an accuracy plot or a cross-validation does.

Nothing is stored: newdata is not modified.

Parameters:

• newdata (_SpatialData) – Locations to ask about, of the same dimension as the training data. Intended for the locations that carry measurements, not for a block model.

• n_sim (int) – Latent realizations per location.

• n_nodes (int) – Equal-share noise values per realization, so that the two axes pool into one sample.

Returns:

dict of str to ndarray – One (n_data, size, n_sim * n_nodes) array per variable whose likelihood carries a warping. Categorical variables are skipped: their noise lives in the probabilities, leaving no value for a measurement to scatter around.

Return type:

dict[str, ndarray]

See also

predict

the ground, with the noise integrated out.

responsibilities(newdata, store=True)

Posterior probability that each measurement came from each noise component of a mixture likelihood.

One answer per location: a Mixture is a mixture over the row, so a measurement wrong in one component of a vector variable is a wrong measurement.

Parameters:

• newdata (_SpatialData) – Point data carrying the variables’ measurements – the training data, a validation set, or the out-of-fold container cross_validate() returns. Modified in place if store.

• store (bool) – Whether to file the answer on each variable, as <variable>/responsibilities/<component>, as well as return it.

Returns:

dict of str to ndarray – One (n_data, n_components) array per variable whose likelihood is a mixture, rows summing to one, and missing at locations without a measurement.

Raises:

ValueError – If no variable has a mixture likelihood, if newdata lacks one of those variables, or if it fans each location into several rows, as a block model does.

Return type:

dict[str, ndarray]

Notes

Read out of fold. At a training location the model interpolates its own measurement, so an outlier is partly absorbed into the fit and its responsibility understates it.

geoml.models.refine(model, blocks, n_sim=20, split_on=None, tolerance=0.05, include_noise=True, where=None, meshes=None, verbose=False)

Predict on a block model, cutting finer wherever it cannot decide.

Predicts on the coarse blocks, splits the ones still in doubt, predicts only what the split created, and repeats. Three criteria mark a block for splitting: its realizations disagree about which side of a cut-off or a category boundary it falls on (needs_splitting()); a neighbour is more than one level finer than it (unbalanced()); or a mesh passes through it (crossed_by()).

The loop ends by itself. Each pass takes the blocks it splits one level finer, and no criterion marks a block already at the lattice’s max_levels, so within that many passes every block is either settled or as fine as the block model was built to go.

Parameters:

• model – A trained model with a predict method, normally a VGPNetwork.

• blocks (BlockSet3D) – The coarse block model to start from. Not modified: each pass builds a new one.

• n_sim (int) – Realizations to draw at every pass.

• split_on (str | Sequence[str] | None) – Which variables have a say in the decision. All of them by default.

• tolerance (float) – The share of realizations that must find a block divided before it is cut.

• include_noise (bool) – Passed to VGPNetwork.predict().

• where (ndarray | Sequence[int] | str | None) – One boolean per block, the indices of the blocks worth modelling, or the name of a boolean metadata column holding the same. The rest are never predicted and never cut, and keep their missing values. Given once, against the blocks as they stand: the mask is carried across each split.

• meshes (Sequence[Mesh3D] | None) – Surfaces and closed bodies whose crossings force a split. Costs a side test per sub-block of every splittable block, each pass.

• verbose (bool) – Print what each pass cut.

Returns:

BlockSet3D – The refined model, predicted throughout.

Return type:

BlockSet3D

Notes

What decides a split carries no noise: the likelihood noise is integrated out rather than drawn, so a block never straddles a cut-off on account of spread that splitting cannot resolve.

Only the blocks a pass creates are predicted. A block that was not split is the same block on the same support, and its value still stands.

Written out, the loop is three calls – predict, ask, split – which is the way to stop part way and inspect a pass. See docs/variable-block-models.md.

geoml.models.cross_validate(model, folds='fold', refit='variational', iterations=200, n_sim=20, n_nodes=32, path=None)

Score a model on folds it never saw, with one short refit per fold.

The trained model is saved once. Each fold gets a copy rebuilt around the data with that fold removed; under refit=”variational” its variational state – the part of a trained model that encodes the data – is re-initialized and every other parameter frozen, so the fold model starts ignorant of the held-out rows, and only that state is refitted. The fold model then predicts its held-out rows, and only those, into one shared copy of the training data. Folds partition the data, so every location ends up predicted by a model that never saw it.

Scores are of measurements: the held-out values are samples, so each fold model is asked through VGPNetwork.predict_measurements(). Categorical variables get no rows – subset the returned container by fold and use the variable’s own compute_metrics.

Parameters:

• model (VGPNetwork) – A trained model. It is saved and copied; the original is untouched.

• folds (str) – Name of the metadata column holding the fold labels, as spatial_k_fold() writes. Any labelling works: a hole-id column gives leave-one-hole-out.

• refit (str) – “variational” to refit the variational state alone, or “all” to warm-start every trainable parameter from its trained value and continue on the reduced data.

• iterations (int) – Training iterations per fold.

• n_sim (int) – Latent realizations, for the out-of-fold predictions and the measurement samples alike.

• n_nodes (int) – Noise values per realization in the measurement samples.

• path (str | PathLike | None) – Where to keep the saved model and its fold copies. A temporary directory, removed at the end, unless one is given.

Returns:

• oof (container) – A copy of the training data carrying out-of-fold predictions and simulations, plus one metadata column per scored component (pit_<variable>, or pit_<variable>_<component>) holding where each measurement fell inside its own predictive distribution.

• scores (pandas.DataFrame) – One row per component and fold, plus a pooled “all” row, with rmse, mae, bias, crps and goodness against the held-out measurements.

Return type:

tuple[_SpatialData, DataFrame]

See also

conformalize

calibrate interval widths on the PIT columns.

geoml.data.PointData.spatial_k_fold

folds that mimic a prediction task.

Notes

The hyperparameters and the warping were fitted on all the data, including each fold’s – the concession kriging cross-validation also makes when it keeps the variogram fixed. Design record and measurements: docs/cross-validation.md.

geoml.models.conformalize(oof, name, component=None)

Build a conformal calibration from a cross-validation.

Reads the out-of-fold PIT column cross_validate() left on its container.

Parameters:

• oof (_SpatialData) – The container cross_validate() returned.

• name (str) – The variable.

• component (str | None) – Which component, when the variable is a vector or compositional one.

Returns:

ConformalCalibration – Calibrated on that component’s out-of-fold measurements.

Return type:

ConformalCalibration

Examples

>>> oof, scores = geoml.models.cross_validate(model)
>>> calibration = geoml.models.conformalize(oof, "grade")
>>> samples = model.predict_measurements(new_points)["grade"]
>>> lower, upper = calibration.interval(samples[:, 0, :], coverage=0.9)

class geoml.models.ConformalCalibration(pit)

Bases: object

Interval coverage repaired from out-of-fold PIT values.

Split conformal prediction on the score \(|u - 1/2|\), where u is where a held-out measurement fell inside its own predictive distribution. nominal() answers: at what level must a central interval be cut so that it covers a given share of fresh measurements? For a calibrated model that is the share itself; an overconfident model is told to cut wider and a hedging one narrower. The conformal quantile carries the finite-sample guarantee – coverage at least the share asked for – under exchangeability of the calibration scores with the prediction’s.

Parameters:

pit (ArrayLike) – Probability integral transform values, one per calibration measurement, as cross_validate() stores them.

Raises:

ValueError – If no finite PIT values are given.

Notes

Spatial data is not exchangeable point by point, which is what the folds are for: built to mimic the prediction task, they draw the calibration scores from conditions like deployment’s. The intervals are of measurements – the ground is never observed.

The repair is bounded by the ensemble: an interval cut from samples cannot reach past their range, so nominal(q) == 1.0 means the model was too sure for its samples to say how much wider the interval should be. Raise n_sim or n_nodes, or reconsider the model.

References

Vovk, V., Gammerman, A. and Shafer, G. (2005) Algorithmic Learning in a Random World. Springer.

nominal(coverage)

The level to cut a central interval at, to cover coverage.

Parameters:

coverage (float) – The share of fresh measurements the interval should contain.

Returns:

float – The level to pass to interval(), between 0 and 1. A value of 1 means the calibration data cannot say how much wider the interval must be.

Return type:

float

interval(samples, coverage=0.9)

A calibrated central interval per row of measurement samples.

Parameters:

• samples (ArrayLike) – One component’s measurement samples, of shape (n_data, n_samples), as VGPNetwork.predict_measurements() returns them once sliced to the component.

• coverage (float) – The share of fresh measurements the interval should cover.

Returns:

lower, upper (ndarray) – One bound per location, of shape (n_data,).

Return type:

tuple[NDArray[float64], NDArray[float64]]