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='independent', training_tolerance=None, propagation='joint')

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='independent', training_tolerance=None, propagation='joint')

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 (bool) – 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. “independent” (the default since 0.9.0) lets each expert speak for its own set alone – O(K) in the expert count – so duplicated points in overlapping sets may disagree, and the data-side weighting arbitrates. “consensus” (the default before, and what a model saved before 0.9.0 keeps) predicts every expert’s set from every other and combines by precision weighting – O(K^2). 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. An operation over GP nodes – Add, LinearCombination – still makes one layer: the answer and the time are the same under both rules until a GP node sits on top. The expected kernel (propagation=”joint”) takes the independent rule only.

• 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.

• propagation (str) – How a GP node reads an input that is another node’s uncertain output. “joint” (the default since 0.9.0) averages the node’s kernel over its input’s distribution – the expected kernel – each node handing its children the covariance between locations along with the variances, so two locations whose inputs move together stay correlated. “marginal” (what a model saved before 0.9.0 keeps) hands on each location’s variance alone and widens the range by half of it, which understates a deep model’s uncertainty by orders of magnitude and lets training buy a small noise with it. Only a GP node above another GP node, or above an uncertain input, is affected. Under “joint” the expert propagation must be “independent”, and a GP node reading an uncertain input takes the Gaussian, exponential, Matérn or rational quadratic kernel.

jit_predict = False

qmc_simulations = False

expert_propagation = 'consensus'

training_tolerance = None

propagation = 'marginal'

class geoml.models.VGPNetwork(data, variables, likelihoods=None, latent_network=None, 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] | dict) – The name of a variable in data to model, a list of names, or a mapping from each name to its likelihood – in which case likelihoods is left out.

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

• latent_network (_latent.network._LatentVariable | Sequence[_latent.network._LatentVariable] | None) – The tree’s leaves, from geoml.latent: a list of nodes, one per likelihood and in the same order, each sized as its likelihood; or a single node serving every likelihood, sized as their sizes summed and split among them in order. Several leaves may share parents, or sit on separate trees with roots of their own.

• 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.

leaves

The tree’s output nodes, one per likelihood or one for all.

Type:

list

latent_network

The single leaf, where there is one. A model with several leaves raises here and points at leaves.

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.viz.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.

property latent_network

The tree’s single leaf, where there is one.

Every model built before a list of leaves was accepted has exactly one, and everything that read it keeps working. A model with several leaves has no single output node to hand back, and says so rather than returning sometimes a node and sometimes a list.

property converged: bool

Whether the bound has settled in the current phase of training.

True once options.training_tolerance has stopped a call, and until the optimizer is replaced or reset or other parameters are trained. A training call made while it holds returns without a step.

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.

Consecutive calls with one optimizer on the same trained parameters continue each other: train_full(100) twice takes the steps of train_full(200) and stops where it would. set_learning_rate starts afresh.

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; a call made once it has settled (converged) takes none.

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 bound of each batch is appended to training_log.

Consecutive calls with one optimizer on the same trained parameters continue each other, the batch order included: train_svi(5) twice takes the steps of train_svi(10) and stops where it would. set_learning_rate starts afresh.

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; a call made once it has settled (converged) takes none. 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=None, 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 | None) – Number of realizations to draw per location. None takes the number newdata already holds, or 20 where it holds none.

• 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, the name of a boolean metadata column holding the same, or None for all of them. Locations left out keep whatever they hold, including their simulations, and a location never visited stays missing.

Raises:

ValueError – If where names some locations of a variable that already holds simulations, and n_sim asks for a different number of them.

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.

Where newdata holds measurements – the training data, a validation set – two things are written beside the prediction as metadata, one column per variable component, for the figures that compare a model with its data to read without the model: pit_<variable>[_<component>], where each measurement falls in the predictive distribution of a measurement (0 to 1; uniform on data the model did not see, if it is well calibrated), and warped_<variable>_<i>, the measurements through the likelihood’s warping, as the model sees them. A block model, whose locations are not points, gets neither.

predict_node(node, newdata, n_sim=None, name=None, labels=None, where=None)

Predict what a node inside the tree says, into a container.

The same prediction predict() makes, stopped at node instead of carried on to the likelihoods: the node’s mean and variance at every location, and its realizations, written as a LatentVariable with one part per output of the node. Where a GPWalk moved the coordinates, what a shared parent says before two leaves diverge, what a Linear trend adds.

Realization s of the node is the one realization s of the first leaf above it was built from, wherever operation nodes (Add, Linear, Scale, …) are all that stand between them; a GP node above reads the node’s mean and variance, never its realizations.

Parameters:

• node (_LatentVariable) – A node of this model’s tree, above its input.

• newdata (_SpatialData) – The locations to predict, of the same dimension as the training data: points, a grid, a section or a mesh. Modified in place.

• n_sim (int | None) – Number of realizations per location. None takes the number the variable already holds, or 20.

• name (str | None) – The variable to write, the node’s own name by default. A LatentVariable of that name is written into; any other variable of that name is refused.

• labels (Sequence[str] | None) – One name per output of the node, “0”, “1”, … by default.

• where (ndarray | Sequence[int] | str | None) – One boolean per location, the indices of the locations to visit, the name of a boolean metadata column holding the same, or None for all of them. Locations left out keep what they hold.

Returns:

LatentVariable – The variable written, as held by newdata.

Raises:

ValueError – If node is an input or not in this model’s tree, if newdata is a block model, if labels is not one per output, if name is held by another kind of variable or one with other labels, or if where names some locations of a variable holding a different number of realizations.

Return type:

LatentVariable

See also

predict

the prediction carried on to the likelihoods.

Notes

Point support only. A block’s value is the mean over its sub-blocks taken by the likelihood, and a node has none; predict on a Grid3D over the same ground instead.

expert_weights(container=None)

Each location’s weight for each expert, as the model blends them.

One sweep, an expert at a time, so that memory holds one expert’s computation whatever their number: at every location the expert’s explained variance gives its raw weight, the raw weights are normalized over the experts as the prediction normalizes them, and averaged over the outputs of every GP node on the expert’s input. A node that reads another reads, in the sweep, what the one expert alone gives below it – which is the blend wherever that expert carries weight.

Parameters:

container (_SpatialData | None) – The locations, the training data by default.

Returns:

ndarray – Of shape (n, n_experts), rows summing to one – on a network of several inputs, every input’s experts side by side in the network’s order, each input’s summing to one.

Return type:

NDArray[float64]

train_by_expert(epochs=10, batch_size=None, coverage=0.99, global_update='epoch', weights_every=1, visits=1, slots=True, decay='steps', sampling='assignment', quotas='equal', stepping='active')

Train an expert at a time, so memory does not grow with their number.

Each epoch every data row draws one expert from its own weights, so that a row lands in an expert’s batch with the probability its weight for the expert gives and is read once an epoch. An expert’s rows are split into batches of about N / (J x visits) rows, so a crowded expert takes more of them, and the batches come in a random order. The experts active on a batch – on each input, those holding coverage of the weight the batch carries, its own expert always – are the only ones computed, and each takes a step on the batch’s bound: its data term as it stands, less each active expert’s share of its KL divergence, shared among the batches the expert is active on by the weight it carries in each, and a share of the priors. An epoch’s batches add up to the bound. A network on several inputs assigns each row to one expert among every input’s.

Under sampling=”replacement” an epoch is visits rounds instead, each visiting every expert once in an order of its own: for expert j a batch of batch_size rows is drawn with probability proportional to their weight for j, and the batch’s data term is W_j times its mean log-likelihood, W_j the expert’s total weight, over the number of inputs, so that the experts’ terms add up to the whole data term over an epoch in expectation. Under sampling=”partition” the experts, in a random order, each draw their share of the rows still unused – by their weight, without replacement – so that a row a dense expert’s quota leaves behind falls to a later expert’s batch. Under either split of the rows, stepping=”own” steps the batch’s own expert alone, its KL counted whole in its own batches.

Under options.training_tolerance training stops once the epochs’ bound settles, as train_svi’s does, and a call made after that takes no step.

Parameters:

• epochs (int) – Number of epochs.

• batch_size (int | None) – Rows per batch under sampling=”replacement”, options.training_batch_size by default; the other samplings size their batches from visits.

• coverage (float) – The share of a batch’s expert weight its active experts hold; the rest are left out of its blend.

• global_update (str) – How the parameters every expert shares step: “batch” on each batch’s gradient; “round” once a round, on the gradients added up over it; “epoch” once an epoch.

• weights_every (int) – Epochs between sweeps of the weight table, which sets the batch distributions, the KL shares and, in slots, the active sets; without slots the sets are formed from the first sweep and kept.

• visits (int) – About how many batches an expert takes an epoch: under an assignment or a partition the batches hold N / (J x visits) rows; under “replacement” the rounds per epoch, which with batch_size divided by as much read as many rows in visits times the steps.

• slots (bool) – Compute the active experts in as many slots as the largest set holds, so that one traced step serves every set; False traces a step per set.

• decay (str) – How the learning rates decay, 0.999 a step: “steps” counts each expert’s own steps and the shared parameters’ own, as an optimizer of their own would; “epochs” counts, for both, the steps train_svi would have taken by the same point, batches of options.training_batch_size over the whole data an epoch, so that the two decay together as they do there. Each expert steps on every batch whose set it is in and the shared parameters as global_update says, so on their own counts the first falls faster than train_svi’s rate and the second slower; on train_svi’s count both have decayed before training by expert, which needs more epochs, has converged. “epochs” needs slots.

• sampling (str) – “assignment” has each row choose its expert, an epoch as many batches as the experts’ rows make, and no rounds for global_update=”round”; “replacement” draws each expert’s batch from every row; “partition” has the experts split the rows among themselves, as above.

• quotas (str) – Under “partition”, how many rows each expert draws: “equal” shares, or by “weight”, each expert’s share of the weight the rows carry, so that a dense expert leaves fewer rows behind.

• stepping (str) – Under “partition” or “assignment”, which experts step on a batch: its “own” expert alone, its KL counted whole in its own batches, or every “active” one, each expert’s KL shared among the batches it is active on by the weight it carries there.

Returns:

dict – The run’s record: the bound per epoch, seconds per epoch, the active sets (on a network of several inputs, one per input for each expert), the weight they drop, and the traces the steps took; under “partition” and “assignment” also, per epoch, the batches, the fewest and most an expert stepped on, the active sets’ sizes, the fewest and most rows a batch drew, the share of each row’s weight its batch’s sets hold, and the weight the first and the last quarter of the batches have for their own expert.

Return type:

dict[str, Any]

predict_by_expert(newdata, n_sim=None, coverage=0.99, neighbours=8, include_noise=True, grouping='home', where=None, slots=True, pack=True)

Predict with only the experts active at each location.

The training data’s expert weights are carried to the locations by an inverse-distance average of the neighbours nearest data, in each input’s transformed space. A location’s own experts, on each input, are those holding coverage of its weight, a block’s the union over its sub-blocks. Each group of locations is predicted with only its experts, the rest left out of the blend.

Under slots every group is computed in as many slots, on each input, as the largest active set training forms from the data, so that one traced refresh and one traced prediction serve them all and memory stays where training’s was. A location takes its own experts, or its leading ones where it needs more than there are slots – the weight that leaves out is in left_out. With pack, groups are merged wherever the experts of both fit in the slots, which costs nothing a padded slot would not; a location then blends experts beyond its own, so its answer depends, within the weight coverage leaves out, on what is predicted beside it. Without, it depends on the location alone.

Without slots each group is a trace of its own, and grouping keeps them few: under “exact” the locations are grouped by their own experts; under “home” a location takes, on each input, the active set training forms for the first of its leading experts whose set holds coverage of its weight, then the union of its leading experts’ sets, and its own only where neither does.

Parameters:

• newdata (_SpatialData) – The locations to predict. Modified in place, as by predict.

• n_sim (int | None) – Number of realizations, as for predict.

• coverage (float) – The share of a location’s expert weight its experts must hold.

• neighbours (int) – Data averaged into a location’s weights.

• include_noise (bool) – As for predict.

• grouping (str) – “home” or “exact”, without slots.

• where (ndarray | Sequence[int] | str | None) – The locations to predict, as for predict; the rest keep what they hold.

• slots (bool | int) – Compute the groups in a fixed number of slots – True for as many as training’s largest active set on each input, or a number for every input, memory growing with it times options.prediction_batch_size; False traces one prediction per group.

• pack (bool) – Under slots, merge groups whose experts fit the slots together.

Returns:

dict – The groups (their active experts – on a network of several inputs, one set per input – and sizes), the number of slots, and the largest share of weight left out at each location, NaN where none was predicted.

Return type:

dict[str, Any]

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. The nodes are rotated at random per location and realization from the model’s seed, so the sample is unbiased in every moment, reaches the tails, and is independent between locations.

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.

Raises:

MemoryError – If the answer would not fit. It is held whole, by design, and costs n_sim * n_nodes * 8 bytes a row for each column of each variable – 5 KB a row at the defaults, and twice that at the peak, while the batches and the assembled whole are both live. The ceiling is models.MEASUREMENT_LIMIT, and the message names the ways under it.

Return type:

dict[str, ndarray]

See also

predict

the ground, with the noise integrated out.

measurement_batches(newdata, n_sim=20, n_nodes=32, where=None)

The measurement samples, a batch of locations at a time.

What predict_measurements() returns whole, yielded in the pieces it assembles it from. The samples are the largest thing this model produces – n_sim * n_nodes values a row for each column, 5 KB a row at the defaults – and every statistic taken of them (coverage, CRPS, the point errors, a PIT) reduces the sample axis one row at a time, so a caller that accumulates as it goes never holds more than a batch. That is what cross_validate() does, and it is why this door carries no size ceiling where the other one must.

Parameters:

• newdata (_SpatialData) – Locations to ask about, as for predict_measurements().

• n_sim (int) – Latent realizations per location.

• n_nodes (int) – Equal-share noise values per realization.

• where (ndarray | Sequence[int] | str | None) – One boolean per location, or the indices of the locations to ask about; all of them by default. A location’s sample does not depend on which others are asked about with it.

Yields:

• rows (ndarray) – Indices into newdata of the locations in this batch.

• samples (dict of str to ndarray) – One (len(rows), size, n_sim * n_nodes) array per variable whose likelihood carries a warping, in the variable’s own units.

See also

predict_measurements

the same samples, assembled and returned.

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=None, split_on=None, tolerance=None, include_noise=True, where=None, meshes=None, verbose=False, by_expert=False, expert_options=None)

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: the prediction at its sub-blocks falls on both sides of a cut-off or a category boundary (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 | None) – Realizations to draw at every pass. None takes the number blocks already holds, or 20 where it holds none.

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

• tolerance (float | None) – Deprecated, and without effect: a block is divided or it is not. It was the share of realizations that had to find a block divided, and will be removed.

• 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.

• by_expert (bool) – Predict each pass with VGPNetwork.predict_by_expert(), so that memory does not grow with the number of experts.

• expert_options (dict[str, Any] | None) – Keywords for VGPNetwork.predict_by_expert() under by_expert – coverage, neighbours, grouping, slots, pack.

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, method='full', epochs=50, n_sim=20, n_nodes=32, path=None, expert_options=None)

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

The trained model is saved once and one fold model is rebuilt from the file, around the data with the first fold removed; every later fold swaps its own training rows into that same model, restores the file’s parameters and zeroes the optimizer’s memory, all in place, so each fold starts exactly where a reloaded model would while the graphs traced for the first serve them all (a model rebuilt per fold leaves its graph machinery resident for the life of the process – see Notes). Under refit=”variational” the 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, re-initialized on every node; “leaves” to re-initialize and refit it on the terminal GP nodes only – the GP nearest each likelihood, found from the leaf down through operation nodes – keeping the interior (a GPWalk’s field, a shared parent) as all the data taught it – measured to score 20% past the honest reference on Jura, the interior remembering the held-out rows, so a diagnostic of that memory rather than a score (E2 in docs/cross-validation.md); or “all” to warm-start every trainable parameter from its trained value and continue on the reduced data.

• iterations (int) – Training iterations per fold, under method=”full”. Ignored under method=”svi”, which counts in epochs.

• method (str) – “full” to refit each fold on all its data at every iteration, or “svi” to refit in minibatches of options.training_batch_size, or “by_expert” to refit with VGPNetwork.train_by_expert() and predict the held-out rows with VGPNetwork.predict_by_expert() – the same choice as VGPNetwork.train_full() against VGPNetwork.train_svi(), and worth making for the same reason: a fold refit costs the whole reduced data set per iteration whatever is frozen, so a model too large to train full-batch is too large to cross-validate that way.

• epochs (int) – Passes over each fold’s data, under method=”svi” or method=”by_expert”. A separate argument from iterations because the two count different things: an epoch is one visit to the data in batches, so it is many gradient steps, and the numbers that make sense for one are wrong for the other.

• 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.

• expert_options (dict[str, Any] | None) – Keywords for VGPNetwork.train_by_expert() under method=”by_expert”.

Returns:

• oof (container) – A copy of the training data – every variable and metadata column, the fold labels included – 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. It saves with to_zarr and opens with open like any container.

• scores (pandas.DataFrame) – One row per scored component and fold, in the folds’ sorted order, then one per component pooled over every fold, whose fold is “all”. The columns: n, the held-out measurements scored; rmse, mae and bias of the mean of the measurement samples against them, the bias being that mean minus the measurement; crps, from the samples; goodness, of the samples’ interval coverage; variable; component, the variable’s own name for a scalar one; and fold, the fold’s label.

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.

Under method=”svi” the convergence rule, if options.training_tolerance is set, reads one value an epoch – the mean bound over its batches – rather than one an iteration, so it needs a few epochs before it can fire at all. On a short refit that is worth knowing: too few epochs and the rule never speaks; the cap does the stopping.

Memory is flat across folds by construction. TensorFlow keeps the graph machinery of a differentiated function resident after the function dies, so a fold model built and dropped per fold cost 2.6 GB a fold on a 5000-row copy of a real model and took a five-fold run on the full data past a 62 GB machine; one model with its rows swapped costs one model, and the same run measured 3.2 GB in total and 3.8x faster, the rebuild, the retrace and any XLA compilation being paid once.

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]]