mixle.ppl.field module¶
Latent fields observed through many proxies, fit jointly, with a posterior readable off any node.
This is the dedicated builder behind the “Cox foundation” extension: a shared latent field (a vector over an index grid) carrying a Gaussian-process / Gaussian-Markov-random-field prior, observed through an arbitrary list of heterogeneous proxy likelihoods (each its own forward model + noise), all coupled to the one field and fit in a single joint optimization. The Laplace posterior is then read off any node – the field itself or any proxy’s latent parameters – because information is additive: the joint posterior precision is the prior precision plus every proxy’s Fisher information, evaluated at the joint MAP.
The motivating instance is the earth-field engine: a latent temperature curve T(t) inferred jointly
from a benthic delta18O forward model (Gaussian) and foram thermal niches (logistic occupancy), where
combining proxies sharpens the posterior on the shared climate field. A log-Gaussian Cox process (latent
log-intensity field -> Poisson counts) is the same pattern with a Poisson proxy.
Example:
field = GaussianField(index=np.arange(50), kernel=RandomWalk(scale=0.3, ridge=3.0), name="T")
post = fit_field(field, [
GaussianProxy(d18O, index=obs_idx, slope=free, intercept=free, scale=0.15),
LogisticNicheProxy(presence), # presence[taxon, bin]
], how="laplace")
mean, sd = post.posterior("T") # the field marginal
post.summary() # mean/sd for every node
- class FieldKernel[source]
Bases:
objectA Gaussian field prior over an index grid, expressed as a precision matrix.
- covariance(index)[source]
The prior covariance matrix, when available cheaply (a kernel defined by its covariance). Used by the low-rank Gauss-Newton posterior to get field marginals without a dense precision inverse;
None(the default) falls back to the dense path.
- class RandomWalk(scale=1.0, order=1, ridge=None)[source]
Bases:
FieldKernelA Gaussian-Markov-random-field smoothness prior via finite differences.
order=1penalizessum (f[i+1]-f[i])^2 / scale^2(a random-walk / integrated-noise prior);order=2penalizes the discrete curvature (an integrated-Wiener / thin-plate prior).ridgeaddsI / ridge^2so the otherwise-improper prior is proper (it anchors the level/trend).- scale: float = 1.0
- order: int = 1
- class RBF(lengthscale=1.0, amplitude=1.0, jitter=1e-06)[source]
Bases:
FieldKernelA squared-exponential (RBF) GP prior:
K[i,j] = amplitude^2 exp(-0.5 (d_ij/lengthscale)^2).The precision is
inv(K + jitter*I).indexmay be 1-D (a grid) or 2-D (coordinates, one row per node) – distances are Euclidean, so this is the spatial-field prior as well.- lengthscale: float = 1.0
- amplitude: float = 1.0
- jitter: float = 1e-06
- covariance(index)[source]
The prior covariance matrix, when available cheaply (a kernel defined by its covariance). Used by the low-rank Gauss-Newton posterior to get field marginals without a dense precision inverse;
None(the default) falls back to the dense path.
- class AnisotropicRBF(ranges=(1.0, 1.0), angle=0.0, amplitude=1.0, jitter=1e-06, metric=None)[source]
Bases:
FieldKernelA geometrically-anisotropic squared-exponential GP prior – correlation that is longer along one direction than another (geological layering / bedding, faulted fabric, flow channels).
The squared distance is the Mahalanobis form
(x-y)^T M (x-y)withMbuilt from per-axis correlationrangesafter rotating the coordinates: in 2-D,angle(radians, counter-clockwise from the x-axis) sets the principal direction andranges=(major, minor)the correlation length along/across it. For >2-D, supplyrangesper axis (axis-aligned) or a fullmetricmatrixMdirectly. Provably positive-definite – it is an ordinary RBF on linearly-transformed coordinates.- Parameters:
- angle: float = 0.0
- amplitude: float = 1.0
- jitter: float = 1e-06
- covariance(index)[source]
The prior covariance matrix, when available cheaply (a kernel defined by its covariance). Used by the low-rank Gauss-Newton posterior to get field marginals without a dense precision inverse;
None(the default) falls back to the dense path.
- class GreatCircleRBF(lengthscale=1.0, amplitude=1.0, jitter=1e-06, radius=1.0)[source]
Bases:
FieldKernelSquared-exponential GP prior on the sphere for a
(lat, lon)index (degrees).K[i,j] = amplitude^2 exp(-0.5 chord_ij^2 / lengthscale^2)wherechordis the straight-line distance through the R^3 embedding – so the kernel is positive-definite on the sphere at every lengthscale (a geodesic-distance RBF is not). For nearby points the chord matches the great-circle distance, solengthscalereads as a spatial correlation length in the same units asradius(radius=6371.0088-> kilometres on Earth; default unit sphere -> radians of arc).- lengthscale: float = 1.0
- amplitude: float = 1.0
- jitter: float = 1e-06
- radius: float = 1.0
- covariance(index)[source]
The prior covariance matrix, when available cheaply (a kernel defined by its covariance). Used by the low-rank Gauss-Newton posterior to get field marginals without a dense precision inverse;
None(the default) falls back to the dense path.
- class GreatCircleMatern(lengthscale=1.0, amplitude=1.0, nu=1.5, jitter=1e-06, radius=1.0)[source]
Bases:
FieldKernelMatern GP prior on the sphere for a
(lat, lon)index (degrees) – the geostatistics default.A Matern covariance of the R^3 chordal distance (so PD on the sphere at every lengthscale), giving rougher, more realistic spatial fields than the infinitely-smooth RBF.
nuin {0.5, 1.5, 2.5} (0.5 = exponential / Ornstein-Uhlenbeck, 1.5 and 2.5 are the common once/twice-differentiable choices).lengthscaleandradiuscarry the same spatial-units meaning as inGreatCircleRBF.- lengthscale: float = 1.0
- amplitude: float = 1.0
- nu: float = 1.5
- jitter: float = 1e-06
- radius: float = 1.0
- covariance(index)[source]
The prior covariance matrix, when available cheaply (a kernel defined by its covariance). Used by the low-rank Gauss-Newton posterior to get field marginals without a dense precision inverse;
None(the default) falls back to the dense path.
- great_circle_distance(latlon_a, latlon_b=None, *, radius=1.0)[source]
Great-circle (geodesic) distance between
(lat, lon)points in degrees.Returns
radius * thetawherethetais the central angle (radians). Withradius=1the result is the angle itself;radius=6371.0088gives kilometres on Earth.latlon_b=Nonereturns the full pairwise(n, n)matrix forlatlon_a; otherwise the(na, nb)cross matrix (a float when both are single points). Uses the haversine formula for accuracy at small angles.
- class GaussianField(index, kernel, name='field')[source]
Bases:
objectA latent field: an index grid plus a Gaussian (GP/GMRF) prior over its node values.
- index: ndarray
- kernel: FieldKernel
- name: str = 'field'
- class FieldSystem(fields, coregion=None)[source]
Bases:
objectSeveral named latent fields fit jointly – the spine for coupled multiphysics / multivariate earth-science models (e.g. an ore-grade field constrained by several geophysical surveys, or coupled paleo-environment fields like temperature + salinity + pCO2 read through many proxies).
By default the prior is block-diagonal over the fields (each keeps its own kernel precision) and cross-field dependence is expressed in the proxy likelihoods: a proxy
loglik(field_t, params, torch)receives the fullparamsdict, so it can read any field by name and couple them (a coupling PDE residual, a grade that depends on density + susceptibility, …). Attach a proxy to a particular field withproxy.on('name'); an unattached proxy defaults to the first field.Pass
coregion=B(aK x Kbetween-field covariance for theKfields) for prior-level coregionalization – the intrinsic coregionalization model, joint prior precisionB^-1 (x) Lambda, where the fields share one spatial structureLambda(the first field’s kernel) and are correlated a priori throughB(e.g. temperature and salinity covary before any data). Requires all fields to share one index/dim.Bmust be symmetric positive-definite.- fields: Sequence[GaussianField]
- class Proxy[source]
Bases:
objectA likelihood hung off the shared field. Subclasses declare params and score the data in torch.
- prefix: str = 'proxy'
- class GaussianProxy(y, index=None, slope=1.0, intercept=0.0, scale=1.0, prefix='gauss')[source]
Bases:
ProxyA linear-Gaussian forward model:
y_j ~ N(intercept + slope * field[index_j], scale).slope,interceptandscaleare each a fixed float or thefreetoken (estimated). With everything fixed this is the exact linear-Gaussian observation that makes the joint posterior exactly Gaussian – the closed-form check on the Laplace covariance.- y: ndarray
- slope: Any = 1.0
- intercept: Any = 0.0
- scale: Any = 1.0
- prefix: str = 'gauss'
- loglik(field_t, params, torch)[source]
- residual(field_t, params, torch)[source]
The standardized Gaussian residual
(y - prediction) / scale(1-D), orNoneif this proxy is not a Gaussian-misfit observation. Used by the Gauss-Newton posterior (how='gauss_newton').
- class LogisticNicheProxy(presence, mu_scale=2.0, prefix='niche')[source]
Bases:
ProxyThermal-niche occupancy:
presence[i,t] ~ Bernoulli(sigmoid(b - 0.5 kappa_i (field[t] - mu_i)^2)).Each row of
presenceis one taxon’s presence/absence across the field’s index bins; the latent niche locationmu_iand precisionkappa_iare co-estimated, andbis a shared baseline. The assemblage acts as a community thermometer – a unimodal response peaked at the niche optimum.- presence: ndarray
- mu_scale: float = 2.0
- prefix: str = 'niche'
- loglik(field_t, params, torch)[source]
- class PoissonProxy(counts, index=None, offset=0.0, prefix='cox')[source]
Bases:
ProxyA log-Gaussian Cox process:
counts[j] ~ Poisson(exp(offset_j + field[index_j])).The field is the latent log-intensity; counts are the point/aggregate observations. This is the canonical Cox-process proxy – the namesake of the foundation.
- counts: ndarray
- offset: Any = 0.0
- prefix: str = 'cox'
- loglik(field_t, params, torch)[source]
- class CustomProxy(loglik_fn, param_specs=(), prefix='custom')[source]
Bases:
ProxyAn arbitrary proxy: supply a torch log-likelihood
loglik_fn(field_t, params, torch)and the parameter specs it reads fromparams(each(name, shape, support, init)).- loglik_fn: Callable
- prefix: str = 'custom'
- loglik(field_t, params, torch)[source]
- fit_field(field, proxies, *, how='laplace', max_iter=500, lr=0.4, init=None, vi_steps=400, vi_lr=0.05, vi_samples=4)[source]
Fit a latent field jointly to a list of proxy likelihoods.
field=Noneruns pure-parameter inference (the proxies carry all the latents, e.g. ODE/PDE coefficients) with no shared field.how='map'returns the joint MAP (no covariance);how='laplace'adds the Gaussian posterior (the inverse-Hessian covariance, exact when every factor is Gaussian; needs a twice-differentiable forward);how='gauss_newton'builds the posterior fromJ^T J + priorwithJthe Jacobian of the standardized residual – first-order only, so it is the posterior for the sparse adjoint solve (wherehow='laplace'cannot run) and is exact for a linear forward;how='vi'fits a mean-field Gaussian variational posterior (reparameterized ADVI on the unconstrained vector), the calibrated approximation for genuinely non-Gaussian posteriors (e.g. with total-variation / Potts priors), and it too is first-order so it works through the sparse solve. Posteriors over any node are read off the returnedFieldPosterior.
- class FieldPosterior(map_values, _cov, _layout, _field_name, _hessian, objective, _field_prior=<factory>, _proxy_info=<factory>, _supports=<factory>, _marg_var=<factory>)[source]
Bases:
objectThe joint posterior.
posterior(node)returns(mean, sd)for the field or any proxy param.The Laplace covariance is the inverse of the joint negative-log-posterior Hessian at the MAP – the prior precision plus every proxy’s Fisher information.
posterior(field, coupling=False)instead inverts only the field’s own Hessian block (the posterior conditional on the other nodes at their MAP), which is the per-proxy additive-information picture.- Parameters:
- map_values: dict
- objective: float
- cov(node)[source]
Posterior covariance of
nodein its natural space (delta method applied for positive nodes).
- posterior(node, *, coupling=True)[source]
(mean, sd)fornodein its natural space.coupling=True(default) marginalizes the other nodes (the honest marginal);coupling=Falsefixes them at the MAP (additive information).
- field_posterior(include=None)[source]
The field posterior
(mean, sd)under a subset of proxies, evaluated at the one joint MAP.Because information is additive, the field posterior precision under any subset of proxies is the prior precision plus those proxies’ Fisher-information blocks.
include=Noneuses every proxy (the joint posterior);include=["gauss"]uses only that proxy – so you can read how much each proxy sharpens the shared field without re-fitting (which would be ill-posed for a single proxy with a free forward-model gain). The mean is the joint MAP for every subset.
- sample(size=1, rng=None, *, nodes=None, given=None)[source]
Draw joint samples from the Gaussian posterior, returned in each node’s natural space.
The Laplace / Gauss-Newton posterior is the Gaussian
N(map, _cov)in the unconstrained parameter space; this draws via a Cholesky factor of that joint covariance (so cross-node and within-field correlations are preserved) and maps each draw back through the node’s support transform (expfor positive nodes). That back-transform is exact – unlikecov(), which linearizes it with the delta method – so a positive node’s draws are properly lognormal. A mean-fieldhow='vi'(or low-rank Woodbury) fit exposes only per-node marginal variances, so there the draws are independent across nodes (which is exactly the mean-field assumption).given(a{node: value}dict) conditions the draw on fixed node values – the closed-form Gaussian conditional of the remaining coordinates given the observed ones (e.g. draw the field consistent with a pinned measurement or a fixed proxy parameter). Conditioning couples nodes through the joint covariance, so it requires the full-covariance fit (how='laplace'/'gauss_newton'); the fixed nodes come back at their given values.- Parameters:
size (int) – number of joint draws.
rng – a
numpy.random.RandomState, an integer seed, orNone.nodes (Sequence[str] | None) – which nodes to return (default: all). The draw is always joint over the full vector.
given (dict | None) – optional
{node: value}to condition on (values in each node’s natural space).
- Returns:
{node: ndarray}with shape(size,)for scalar nodes and(size, dim)otherwise.- Return type:
- class GP(name, index, kernel)[source]
Bases:
objectA latent field node for the equation-style surface:
T = GP("T", index=grid, kernel=...).Supports affine algebra (
c0 - c1*T,T + b) so a linear forward model reads as math; the result carries the field, the gain and the offset tojoint().- Parameters:
name (str)
index (np.ndarray)
kernel (FieldKernel)
- Gaussian(y, *, mean, sd)[source]
A linear-Gaussian observation
y ~ N(gain*field + offset, sd);meanis an affine in a GP.- Return type:
- Niche(presence, *, over, mu_scale=2.0)[source]
Logistic thermal-niche occupancy of
presenceover the fieldover(a community thermometer).
- Cox(counts, *, log_intensity, offset=0.0)[source]
A log-Gaussian Cox process
counts ~ Poisson(exp(offset + field));log_intensityis a GP.- Return type:
- joint(observations)[source]
Assemble a latent-field model from equation-style observations; call
.fit(how=...)to fit.Each item is the
(field, proxy)pair returned byGaussian(),Niche(),Cox()ormixle.ppl.Differential(). Observations sharing a field must name the same one (this surface targets one shared field); field-free observations (a pure-parameter ODE inverse problem) are allowed.
- class FieldModel(field, proxies)[source]
Bases:
objectA latent-field model built from equation-style observations; fit with
.fit(how=...).The one fit verb, matching the rest of
mixle.ppl:joint([...]).fit(how='map'|'laplace')returns aFieldPosteriorwith a posterior over any node. Delegates tofit_field().- Parameters:
field (GaussianField | None)
proxies (list)
- field: GaussianField | None
- proxies: list
- multistart(model, inits, *, how='map', **kw)[source]
Fit
modelfrom several initializations and keep the best (lowest-objective) fit.For the multimodal posteriors of nonlinear inverse problems (e.g. scattering / FWI cycle-skipping), a single optimization can land in a poor local mode;
initsis a list of{node: value}dicts and the fit with the smallestobjectiveis returned. (Frequency continuation – fit a coarse/low- frequency model, then seed a finer one viafit(init=...)– is the complementary strategy.)