# geoML - machine learning models for geospatial data
# Copyright (C) 2026 Ítalo Gomes Gonçalves
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR a PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
The numbers behind the figures.
Nothing here draws anything, and nothing here imports matplotlib. A figure is
one way of looking at these arrays; a dashboard is another, and it reads the
same functions rather than working the values out a second time. It also means
the arithmetic -- which components a fraction of variance asks for, what a
composition looks like once it is opened up -- can be tested against numbers
instead of against pictures.
"""
from collections.abc import Sequence
import numpy as _np
import scipy.sparse as _sparse
import scipy.spatial as _spatial
import scipy.stats as _stats
import geoml._types as _types
import geoml.data as _data
import geoml.likelihood as _lk
import geoml.math.geometry as _geom
import geoml.metrics as _gmet
import geoml.storage as _storage
[docs]
def variable(container: "_data._SpatialData", name: str):
"""The variable called `name`, or a message saying what there is."""
try:
return container.variables[name]
except KeyError:
raise KeyError(
"no variable named %r; found %s"
% (name, ", ".join(sorted(container.variables)) or "none"))
[docs]
def variable_or_component(container: "_data._SpatialData", name: str):
"""
The variable called `name`, or the component of a vector variable.
A composition is held as one variable, so its parts are not in
`container.variables` and asking for `"Zn"` would otherwise mean reaching
into `Elements` by hand. The search is the container's own -- this module
kept a near-copy of it for years, which is the duplication the path work
was started to remove.
"""
try:
named, _ = container._variable_or_component(str(name))
except ValueError as err:
raise KeyError(str(err))
return named
[docs]
def numeric_values(var):
"""
A continuous or vector variable as a matrix, and the rows that hold one.
Returns
-------
values : array
`(n_data, n_columns)`, whatever the variable's length.
measured : array
Boolean, `(n_data,)`. A vector variable is measured at a location only
where every component is.
labels : list
A name per column.
"""
_, has_value = var.get_measurements()
measured = _np.all(has_value == 1.0, axis=1)
# the stored columns rather than what `get_measurements` hands the
# model: a composition's parts reach it as fractions of the whole, and
# a figure shows the units they were measured in -- the same units every
# other figure reads `prediction` and `noise_variance` in
components = getattr(var, "components", None)
parts = [var] if components is None \
else [components[label] for label in var.labels]
values = _np.stack([part.measurements.values.to_numpy() for part in parts],
axis=1)
labels = [axis_label(part) for part in parts]
return _np.asarray(values, dtype=float), measured, labels
[docs]
def axis_label(part):
"""A continuous variable's name, with its unit where it declares one."""
unit = getattr(part, "unit", None)
return str(part.name) if unit is None else "%s (%s)" % (part.name, unit)
[docs]
def category_values(var):
"""
A categorical variable as one label per location.
Returns
-------
values : array
The label at each location; the empty string where nothing was
measured, which is how a missing code reads.
measured : array
Boolean, `(n_data,)`.
labels : list
The categories actually present, in the variable's own order.
"""
column = getattr(var, "measurements_a", None)
if column is None:
column = getattr(var, "measurements", None)
if column is None:
raise TypeError(
"%s holds no measurements to group by" % type(var).__name__)
values = _np.asarray(column.to_numpy(), dtype=object)
measured = values != ""
present = set(values[measured])
labels = [str(label) for label in var.labels if label in present]
return values, measured, labels
[docs]
def groups(values, measured, labels):
"""A boolean mask per category, in the order the labels come in."""
return [(label, measured & (values == label)) for label in labels]
[docs]
def centred_log_ratio(values: _types.ArrayLike) -> _types.FloatArray:
"""
A composition opened up into ordinary numbers.
Proportions carry a constant sum, so their covariance is negative by
construction and a PCA of them describes the constraint as much as the
data. The log-ratios have no such constraint. This is `warping`'s
`CenteredLogRatio` in NumPy, for exploring rather than for modelling.
"""
logged = _np.log(values)
return logged - _np.mean(logged, axis=1, keepdims=True)
[docs]
def logarithm(values: _types.ArrayLike,
compositional: bool = False) -> _types.FloatArray:
"""
The data on a log scale, by the road that suits it.
A composition takes the centred log-ratio: its parts carry a constant sum,
and logging them one at a time leaves that constraint in place, whereas
dividing each by the row's geometric mean first removes it. Anything else
takes an ordinary logarithm.
Non-positive values stop this rather than quietly becoming infinities:
which columns they are in is worth knowing, since a zero in an assay
usually means below detection rather than absent, and what to put in its
place is a decision about the data, not about the figure.
"""
values = _np.asarray(values, dtype=float)
if _np.any(values <= 0):
bad = _np.where(_np.any(values <= 0, axis=0))[0]
raise ValueError(
"a logarithm needs positive values; columns %s hold zeros or "
"negatives. Replace them first -- half the smallest positive "
"value is the usual choice for a detection limit"
% ", ".join(str(int(column)) for column in bad))
return centred_log_ratio(values) if compositional else _np.log(values)
[docs]
def principal_components(values: _types.ArrayLike,
explained: float = 0.9):
"""
The principal components of `values`, down to a share of the variance.
Follows `warping.PCA`: an eigendecomposition of the covariance, taken
largest first. The scores are the plain projection, which is what a biplot
puts on its axes; the loadings are the eigenvectors, for the arrows.
Parameters
----------
values : array
`(n_data, n_columns)`.
explained : float
The share of the total variance to reach. The number of components is
the fewest that reach it -- never fewer than two, so that there is a
plot to draw, and never more than there are columns.
Returns
-------
dict
`scores` `(n_data, n_components)`, `loadings`
`(n_columns, n_components)`, `ratio` (the share each component carries,
all of them), `n_components`, and -- for drawing the components back
onto the data they came from -- `mean`, the centre they turn about,
and `eigenvalues`, the variance along each, whose square root is a
length in the data's own units.
"""
values = _np.asarray(values, dtype=float)
mean = _np.mean(values, axis=0, keepdims=True)
centred = values - mean
covariance = centred.T @ centred / centred.shape[0]
eigenvalues, eigenvectors = _np.linalg.eigh(covariance)
order = _np.argsort(eigenvalues)[::-1]
eigenvalues = _np.maximum(eigenvalues[order], 0.0)
eigenvectors = eigenvectors[:, order]
total = _np.sum(eigenvalues)
ratio = eigenvalues / total if total > 0 else _np.zeros_like(eigenvalues)
reached = int(_np.searchsorted(_np.cumsum(ratio), explained) + 1)
n_components = int(_np.clip(reached, 2, values.shape[1]))
return {"scores": centred @ eigenvectors[:, :n_components],
"loadings": eigenvectors[:, :n_components],
"ratio": ratio,
"n_components": n_components,
"mean": mean.ravel(),
"eigenvalues": eigenvalues}
[docs]
def component_analysis(var, explained: float = 0.9, log: bool = False):
"""
`principal_components` for a variable, on a log scale if asked.
A composition is opened up whether or not `log` is set, and that is not an
oversight: its parts carry a constant sum, so their covariance is singular
and the first component of the raw proportions would spend itself
describing the closure rather than the data. There is no useful PCA of a
composition to refuse it. Everything else is left as measured unless `log`
says otherwise.
"""
values, measured, labels = numeric_values(var)
values = values[measured]
compositional = isinstance(var, _data.CompositionalVariable)
if compositional or log:
values = logarithm(values, compositional)
labels = ["%s(%s)" % ("clr" if compositional else "log", label)
for label in labels]
analysis = principal_components(values, explained)
analysis["labels"] = labels
analysis["measured"] = measured
return analysis
[docs]
def warped_values(model, name: str):
"""
The measurements as the model sees them, after its warping.
A model does not work on the data as measured: the likelihood's warping
takes it somewhere it can be treated as Gaussian and, for several
variables at once, uncorrelated. Whether it succeeded is a thing to look
at rather than assume, and looking at it means passing the data through
the same warping the model trained with -- already fitted, since
`VGPNetwork` initializes each likelihood from the measurements.
Only measured rows go through. `get_measurements` fills the rest with 1.0
to keep the array rectangular, and a warping has no reason to be kind to a
value that was never there.
Returns
-------
values : array
`(n_measured, warping.size_out)`.
measured : array
Boolean, `(n_data,)`, saying which rows those are.
labels : list
A name per column.
"""
names = list(model.variables)
if name not in names:
raise KeyError("the model does not hold %r; it holds %s"
% (name, ", ".join(names)))
likelihood = model.likelihoods[names.index(name)]
if isinstance(likelihood, _lk.LikelihoodMixture):
raise TypeError(
"%s warps each population its own way, so there is no one "
"warped space to show" % type(likelihood).__name__)
warping = getattr(likelihood, "warping", None)
if warping is None:
raise TypeError(
"%s carries no warping, so there is nothing to transform"
% type(likelihood).__name__)
# What the model fed the warping -- `get_measurements`, a composition's
# parts as fractions of the whole, the rows closed -- rather than the
# stored columns `numeric_values` shows in their measured units. The
# warping was initialized on the former; sending ppm and percent
# through it put the log of each part's divisor on every warped column
# as an offset, and the figure showed components centred at 6 and -5.
values, has_value = variable(model.data, name).get_measurements()
values = _np.asarray(values, dtype=float)
measured = _np.all(_np.asarray(has_value) == 1.0, axis=1)
warped, _ = warping.forward(values[measured])
warped = _np.asarray(warped, dtype=float)
# The columns are numbered rather than named after what was measured. A
# warping may rotate the data, or bend it, or take seven columns to five,
# and then a column is a mixture of the measurements rather than one of
# them -- calling it Zn would be wrong in a way nobody would catch.
labels = ["Variable %d" % (i + 1) for i in range(warped.shape[1])]
return warped, measured, labels
def _strided(store, stride):
"""Every `stride`-th value of a store, flattened, read a band at a time.
The same values as `asarray(store).reshape(-1)[::stride]`, arrived at
without the materialization. Each band picks up where the last one left
off, so which values are taken depends on the stride alone and not on where
the chunk boundaries happen to fall -- two stores of the same shape are
thinned to the same positions, which is what lets the columns be paired.
"""
n_columns = store.shape[1] if len(store.shape) > 1 else 1
pieces = []
for band in store.row_bands():
values = _np.asarray(store[band], dtype=float).reshape(-1)
pieces.append(values[(-band.start * n_columns) % stride::stride])
return _np.concatenate(pieces) if len(pieces) > 0 else _np.zeros(0)
[docs]
def simulation_sample(var, most: int = 100000) -> _types.FloatArray:
"""
The simulated values, thinned to a number that can be drawn.
A simulated variable holds `n_data x n_sim` values per component, which on
a block model runs to hundreds of millions -- more than a figure can show
and, spread across every pair of a matrix, more than memory should be
asked to hold. What a distribution looks like is settled by far fewer, so
a stride is taken through the block -- and taken while reading, a band of
locations at a time, so the values that are thrown away are never all in
memory together either.
The stride is the same for every component, and the mask that drops
non-finite values is applied to whole rows. Thinning the components
separately would pair the value simulated at one location with the value
simulated at another, which is not a pair the model ever produced -- the
joint shape, which is the thing being looked at, would be an artefact.
Returns
-------
values : array
`(n_sample, n_columns)`, columns in the variable's own order.
"""
parts = continuous_parts(var)
stride, columns = None, []
for part in parts:
# asked for simulations it does not have, a variable hands back
# `asarray(None)` rather than raising, which would go on to fail
# somewhere less informative
store = getattr(part, "simulations", None)
if store is None:
raise ValueError(
"%r carries no simulations to compare against; predict with "
"n_sim greater than zero first" % var.name)
if stride is None:
stride = max(1, int(_np.prod(store.shape)) // int(most))
columns.append(_strided(store, stride))
values = _np.column_stack(columns)
return values[_np.all(_np.isfinite(values), axis=1)]
[docs]
def padded_range(values: _types.ArrayLike, margin: float = 0.1):
"""
The span of each column, opened up by `margin` at both ends.
Returns a `(low, high)` pair per column. A column of one repeated value
has no span to open up, so it is given a unit of room rather than a window
of zero width.
"""
values = _np.asarray(values, dtype=float)
low = _np.min(values, axis=0)
high = _np.max(values, axis=0)
room = _np.where(high > low, (high - low) * margin, 1.0)
return list(zip(low - room, high + room))
[docs]
def color_limits(values: _types.ArrayLike,
clip: "_types.ArrayLike | None" = None):
"""
The ends of a colour scale, set by where the data mostly is.
One value far from the rest takes the whole of a colour scale with it and
leaves everything else in a single shade -- the usual fate of a
geochemical assay, whose tail is long and whose interest is not in it.
Naming a pair of quantiles ends the scale where the data mostly ends
instead: `[0, 0.99]` for a variable skewed to the right, `[0.01, 0.99]` to
take both tails.
Nothing is dropped and nothing is altered. The values keep their own
numbers and a hover still reports what was measured; what is bounded is
the scale, so the few beyond it take the end colour rather than setting
it.
Returns `(low, high)`, or None when `clip` is None -- and None again when
the two ends come out equal, since a scale of zero width is no scale.
"""
if clip is None:
return None
clip = _np.asarray(clip, dtype=float).ravel()
if len(clip) != 2 or not _np.all((clip >= 0) & (clip <= 1)) \
or clip[0] >= clip[1]:
raise ValueError(
"clip takes two quantiles between 0 and 1, the lower first: "
"[0, 0.99] keeps the long right tail off the scale, [0.01, 0.99] "
"takes both ends. Got %r" % (clip.tolist(),))
values = _np.asarray(values, dtype=float).ravel()
values = values[_np.isfinite(values)]
if len(values) == 0:
return None
low, high = _np.quantile(values, clip)
return None if low >= high else (float(low), float(high))
[docs]
def color_choices(container: "_data._SpatialData",
names: "Sequence[str]"
) -> "tuple[_types.FloatArray, _types.IndexArray, list]":
"""
Several variables' values over the locations that carry all of them.
For a figure that keeps one set of points and swaps the values over them.
The locations are those measured in *every* one of the variables named,
which is the point: a cloud that gained and lost points as the choice
changed would be two things changing at once, and the comparison -- this
variable against that one, here -- is exactly what would be lost.
A name may be a variable, one component of a vector variable, or a vector
variable itself, which stands for all of its components in order. So
`["Elements"]` names seven grades and `["Cd", "Zn"]` names two.
Returns
-------
values : array
`(n_kept, n_columns)`, a column per choice in the order asked for.
rows : array
Which locations those are, as indices into the container.
labels : list
A name per column.
"""
columns, labels = [], []
measured = _np.ones(container.n_data, dtype=bool)
for name in names:
var = variable_or_component(container, name)
if isinstance(var, (_data.RockTypeVariable, _data.BinaryVariable)):
raise TypeError(
"%r is categorical, and a menu of choices colours by a scale "
"rather than by a legend; name continuous variables here and "
"draw the categorical one on its own" % name)
values, has_value, names_here = numeric_values(var)
for i, label in enumerate(names_here):
columns.append(values[:, i])
labels.append(label)
measured &= has_value
if not _np.any(measured):
raise ValueError(
"no location carries all of %s at once, and a menu holds one set "
"of points for every choice on it" % ", ".join(labels))
rows = _np.where(measured)[0]
return _np.column_stack(columns)[rows], rows, labels
[docs]
def cells(n_points: int, most: int):
"""
How many cells to count `n_points` into, at most `most` of them.
The two halves of a comparison rarely hold comparable numbers -- a few
hundred measurements against a hundred thousand simulated values -- and
binning both the same way leaves the sparse half as scattered single counts
that read as noise. Roughly the square root of the count keeps several
points in a typical cell either way.
"""
return int(_np.clip(_np.sqrt(n_points), 10, most))
[docs]
def counts_2d(x: _types.ArrayLike, y: _types.ArrayLike,
bins: int = 60, log: bool = False):
"""
Points counted into cells, with the empty ones left empty.
A cell holding nothing comes back as NaN rather than as zero, so that
whatever draws it can leave it unpainted: painted, it takes the bottom of
the colour scale and fills the panel with a background that looks like
data.
Returns
-------
dict
`x` and `y`, the cell centres; `z`, what the colour is to be taken
from, which is the count or its base-ten logarithm; and `count`, the
count itself, so that a hover can say how many points a cell holds
whichever of the two is being coloured. `z` and `count` are shaped
`(n_y, n_x)`, the way an image is indexed.
"""
counted, x_edges, y_edges = _np.histogram2d(
_np.asarray(x, dtype=float), _np.asarray(y, dtype=float), bins=bins)
# a heatmap is indexed by row, and a row runs along y
counted = counted.T
empty = counted == 0
values = _np.where(empty, _np.nan, counted)
return {"x": 0.5 * (x_edges[:-1] + x_edges[1:]),
"y": 0.5 * (y_edges[:-1] + y_edges[1:]),
"z": _np.log10(values) if log else values,
"count": counted.astype(int)}
[docs]
def density_grid(x: _types.ArrayLike, y: _types.ArrayLike,
grid: int = 60, most: int = 4000):
"""
Smoothed density over a pair of columns, on a regular mesh.
The estimate costs one pass over the data for every mesh cell, so a long
column is thinned first: a density is a shape, and a few thousand points
settle it as well as a few million do.
Returns `(x_axis, y_axis, density)` with `density` shaped
`(grid, grid)` -- rows along `y_axis`, as an image is -- or None when
there is nothing to smooth: fewer than three points, or a column that
never varies, which has no width for a kernel to sit in.
"""
x = _np.asarray(x, dtype=float)
y = _np.asarray(y, dtype=float)
if len(x) > most:
step = len(x) // most
x, y = x[::step], y[::step]
if len(x) < 3 or _np.ptp(x) == 0 or _np.ptp(y) == 0:
return None
x_axis = _np.linspace(_np.min(x), _np.max(x), grid)
y_axis = _np.linspace(_np.min(y), _np.max(y), grid)
mesh_x, mesh_y = _np.meshgrid(x_axis, y_axis)
kernel = _stats.gaussian_kde(_np.vstack([x, y]))
density = kernel(_np.vstack([mesh_x.ravel(), mesh_y.ravel()]))
return x_axis, y_axis, density.reshape(grid, grid)
[docs]
def density_curve(values: _types.ArrayLike, limits, points: int = 200,
most: int = 5000):
"""
A smoothed distribution as a line, across `limits`.
Thinned as `density_grid` is, and for the same reason. Returns
`(grid, density)`, or None when there is nothing to smooth.
"""
values = _np.asarray(values, dtype=float)
if len(values) > most:
values = values[::len(values) // most]
if len(values) < 3 or _np.ptp(values) == 0:
return None
grid = _np.linspace(limits[0], limits[1], points)
return grid, _stats.gaussian_kde(values)(grid)
[docs]
def normal_curve(values: _types.ArrayLike, low: float, high: float,
points: int = 200):
"""
The normal of the same mean and spread as `values`, across a window.
Fitted rather than standard, so that the curve asks about the shape alone.
A warping's parameters are trained along with everything else and need not
leave the data at unit variance -- the GP's amplitude absorbs a scale
factor -- so a standard normal would call a perfectly symmetric result
skewed. Returns `(x, density)`, or None for a column that never varies.
"""
values = _np.asarray(values, dtype=float)
mean, deviation = _np.mean(values), _np.std(values)
if deviation <= 0:
return None
x = _np.linspace(low, high, points)
density = _np.exp(-0.5 * ((x - mean) / deviation) ** 2) \
/ (deviation * _np.sqrt(2 * _np.pi))
return x, density
[docs]
def summary_statistics(values: _types.ArrayLike) -> "dict[str, float]":
"""
The numbers a distribution is summed up by.
The moments are those of the values themselves, in their population
form, and the kurtosis is the excess over a normal's, so a normal reads
0 on the skewness and the kurtosis alike. The coefficient of variation
is the standard deviation over the mean, and missing where the mean is
not positive, which leaves it meaning nothing. Values that are not
finite are left out.
Returns
-------
dict
`n`, `mean`, `std`, `cv`, `skewness`, `kurtosis`, `min`, `q1`,
`median`, `q3` and `max`.
"""
values = _np.asarray(values, dtype=float).ravel()
values = values[_np.isfinite(values)]
if values.size == 0:
raise ValueError("there are no finite values to sum up")
mean = float(values.mean())
deviation = float(values.std())
centred = values - mean
if deviation > 0:
skewness = float(_np.mean(centred ** 3)) / deviation ** 3
kurtosis = float(_np.mean(centred ** 4)) / deviation ** 4 - 3.0
else:
skewness = kurtosis = _np.nan
q1, median, q3 = (float(q) for q in _np.quantile(values,
[0.25, 0.5, 0.75]))
return {"n": int(values.size), "mean": mean, "std": deviation,
"cv": deviation / mean if mean > 0 else _np.nan,
"skewness": skewness, "kurtosis": kurtosis,
"min": float(values.min()), "q1": q1, "median": median,
"q3": q3, "max": float(values.max())}
[docs]
def statistics_lines(statistics: "dict[str, float]") -> "list[str]":
"""
Summary statistics as lines of text for a box in a monospace font.
Two columns, aligned: the count and the moments on the left, the order
statistics on the right. A value in the variable's units keeps four
significant figures and a ratio -- the coefficient of variation, the
skewness, the kurtosis -- two decimals. A number that is missing reads
as a dash.
"""
def number(value, pattern="%.4g"):
if not _np.isfinite(value):
return "-"
text = pattern % value
# a value that rounds to nothing carries no sign
return text[1:] if text.startswith("-") and float(text) == 0 else text
left = [("n", "%d" % statistics["n"]),
("mean", number(statistics["mean"])),
("sd", number(statistics["std"])),
("CV", number(statistics["cv"], "%.2f")),
("skew", number(statistics["skewness"], "%.2f")),
("kurt", number(statistics["kurtosis"], "%.2f"))]
right = [("min", number(statistics["min"])),
("Q1", number(statistics["q1"])),
("median", number(statistics["median"])),
("Q3", number(statistics["q3"])),
("max", number(statistics["max"])),
("", "")]
def width(pairs, item):
return max(len(pair[item]) for pair in pairs)
lines = []
for (name, value), (other, figure) in zip(left, right):
line = "%-*s %*s" % (width(left, 0), name, width(left, 1), value)
if other:
line += " %-*s %*s" % (width(right, 0), other,
width(right, 1), figure)
lines.append(line)
return lines
[docs]
def statistics_side(counts: _types.ArrayLike) -> str:
"""
The upper corner a box of statistics takes over a histogram.
The one over whichever half of the bins has the lower bars: the right,
for the long tail of an assay.
Parameters
----------
counts
The height of the tallest bar in each bin, left to right.
Returns
-------
str
`"left"` or `"right"`.
"""
counts = _np.asarray(counts, dtype=float)
half = len(counts) // 2
left = float(counts[:half].max()) if half else 0.0
right = float(counts[half:].max()) if len(counts) > half else 0.0
return "right" if right <= left else "left"
[docs]
def statistics_top(counts: _types.ArrayLike, edges: _types.ArrayLike,
low: float, high: float, height: float) -> float:
"""
How tall a histogram's axis has to be for its bars to clear a box.
Every bar whose bin reaches under the box has to end below it, with a
little room to spare; the bars beside it may run as high as they like.
Parameters
----------
counts
The height of the tallest bar in each bin, left to right.
edges
The bins' edges.
low, high
Where the box starts and ends along the axis, in the axis's units.
height
How much of the axis's height the box covers, from the top.
Returns
-------
float
The top the axis needs.
"""
counts = _np.asarray(counts, dtype=float)
edges = _np.asarray(edges, dtype=float)
under = (edges[1:] > low) & (edges[:-1] < high)
tallest = float(counts[under].max()) if _np.any(under) else 0.0
return max(1.05 * float(counts.max()),
1.02 * tallest / max(1.0 - height, 0.1))
[docs]
def continuous_parts(var) -> "list[_data.ContinuousVariable]":
"""The continuous columns a figure draws, one per component.
A scalar variable is its own single part; a vector or compositional one
comes back as its components, in label order. Everything that reads
`measurements`, `prediction`, `noise_variance` or the simulations asks
through this, so a figure handed a categorical variable is refused by
name rather than failing later at whichever column it reached for first.
Parameters
----------
var
A variable, as :func:`variable` returns it.
Returns
-------
list of geoml.data.ContinuousVariable
One part for a scalar variable, one per component otherwise.
Raises
------
ValueError
If the variable carries no continuous columns of its own.
"""
components = getattr(var, "components", None)
if components is None:
parts = [var]
else:
parts = [components[label] for label in var.labels]
for part in parts:
if not isinstance(part, _data.ContinuousVariable):
raise ValueError(
"%r is a %s, which carries no continuous values to draw; "
"this figure is of measured and predicted numbers"
% (str(var.name), type(var).__name__))
return parts
[docs]
def component_names(var) -> str:
"""The components a caller can name, ready to drop into a message.
Every figure that needs one grade out of several says so the same way,
and each of them used to reach for `labels` -- an attribute only some
variables have. Asking through `continuous_parts` keeps that reach in
one place.
"""
return ", ".join(str(part.name) for part in continuous_parts(var))
[docs]
def prediction_values(container: "_data._SpatialData", name: str):
"""
What was measured against what was predicted, component by component.
Returns
-------
measured_values, predicted_values : array
`(n_compared, n_columns)`, holding only the locations that carry both.
labels : list
A name per column.
rows : array
Which locations those are, as indices into the container. A figure
that can be brushed needs to say which row each of its points came
from, and only the rows carrying both a measurement and a prediction
are drawn.
"""
parts = continuous_parts(variable(container, name))
measured_values, predicted_values, labels = [], [], []
for part in parts:
if getattr(part, "prediction", None) is None:
raise ValueError(
"%r carries no prediction; run the model over this data first"
% name)
measured_values.append(part.measurements.values.to_numpy())
predicted_values.append(part.prediction.values.to_numpy())
labels.append(axis_label(part))
measured_values = _np.stack(measured_values, axis=1).astype(float)
predicted_values = _np.stack(predicted_values, axis=1).astype(float)
both = _np.all(~_np.isnan(measured_values), axis=1) \
& _np.all(~_np.isnan(predicted_values), axis=1)
if not _np.any(both):
# a variable is built with a prediction column already in place, empty
# until a model fills it, so this is what never having run one looks
# like -- as well as what a genuine mismatch looks like
raise ValueError(
"%r has no location carrying both a measurement and a prediction; "
"run the model over this data first" % name)
return (measured_values[both], predicted_values[both], labels,
_np.where(both)[0])
[docs]
def inside_trim(measured: _types.ArrayLike, predicted: _types.ArrayLike,
trim: "_types.ArrayLike | None" = None):
"""
Which locations a predicted-against-measured panel keeps under a trim.
A few values far from the rest set the limits of such a panel and
squeeze everything else into a corner of it. Naming a pair of quantiles
sets a window instead, running from the lower quantile of the measured
or the predicted values, whichever is lower, to the upper quantile of
whichever is higher; a location outside it on either axis is left out.
`[0, 0.99]` takes off the long right tail of an assay, `[0.01, 0.99]`
both ends.
Returns a boolean mask over the locations, every one kept when `trim`
is None.
"""
measured = _np.asarray(measured, dtype=float)
predicted = _np.asarray(predicted, dtype=float)
if trim is None:
return _np.ones(len(measured), dtype=bool)
levels = _np.asarray(trim, dtype=float).ravel()
if len(levels) != 2 or not _np.all((levels >= 0) & (levels <= 1)) \
or levels[0] >= levels[1]:
raise ValueError(
"trim takes two quantiles between 0 and 1, the lower first: "
"[0, 0.99] leaves out the long right tail, [0.01, 0.99] both "
"ends. Got %r" % (levels.tolist(),))
# each end from whichever side reaches further, so that only what would
# stretch the window is left out: trimming each axis by its own
# quantiles would also drop the highest predictions, which a smoothing
# model packs well inside the measured range -- points that stretch
# nothing, and the very ones that show the smoothing
low = min(_np.quantile(measured, levels[0]),
_np.quantile(predicted, levels[0]))
high = max(_np.quantile(measured, levels[1]),
_np.quantile(predicted, levels[1]))
return ((measured >= low) & (measured <= high)
& (predicted >= low) & (predicted <= high))
def _bin_edges(values, bins):
"""Where to cut, from a count or from the positions themselves."""
if _np.ndim(bins) > 0:
edges = _np.unique(_np.asarray(bins, dtype=float))
else:
bins = int(bins)
if bins < 1:
raise ValueError("bins must be at least 1, got %r" % bins)
# equal *count*, not equal width: a predicted grade is skewed, and
# equal-width bins would put nine tenths of the data in the first one
# and a single sample in the last. Pass the positions for equal width.
edges = _np.unique(_np.quantile(values, _np.linspace(0, 1, bins + 1)))
if len(edges) < 2:
raise ValueError(
"there is nothing to bin: every predicted value is %g" % edges[0])
return edges
[docs]
def step_path(lo, hi, values):
"""A per-bin value as a polyline that steps at the edges.
One number per bin is not a curve, and drawing it through the bin centres
says it is. Stepping at the edges shows where the bins are without a
second thing on the figure to say so. A gap between bins comes back as a
break rather than a line across it.
"""
lo, hi, values = (_np.asarray(a, dtype=float) for a in (lo, hi, values))
x, y = [], []
for i in range(len(values)):
if i > 0 and lo[i] > hi[i - 1]:
x.append(_np.nan)
y.append(_np.nan)
x.extend([lo[i], hi[i]])
y.extend([values[i], values[i]])
return _np.array(x), _np.array(y)
def _measured_categories(var, name):
"""The measured label per location, and the rows a score may use.
The truth side shared by `confusion_matrix` and `reliability`: labels
as strings — how the figures name categories everywhere, with the
missing code decoding to the empty string — and the rows narrowed to
the unambiguous ones: a contact carries two measurements, and where
they disagree (`boundary`) there is no one truth for a prediction to
be right about.
"""
truth = getattr(var, "measurements_a", None)
if truth is None:
truth = getattr(var, "measurements", None)
if truth is None or getattr(truth, "labels", None) is None:
raise TypeError(
"%s %r holds no measured categories to score against"
% (type(var).__name__, str(name)))
measured = _np.asarray(truth.to_numpy()).astype(str)
keep = measured != ""
boundary = getattr(var, "boundary", None)
if boundary is not None:
keep &= ~_np.asarray(boundary.values, dtype=bool)
return measured, keep
[docs]
def reliability(container: "_data._SpatialData", name: str,
bins: _types.Bins = 10) -> list[dict]:
"""
Whether a claimed probability is the frequency it claims.
The categorical half of what the accuracy figure asks of a continuous
variable: for each category, the locations are binned by the
probability the model assigned to it, and each bin's mean claim is
set against the share of its locations actually measured as that
category. On the diagonal a 70% claim is that category 70% of the
time; below it the model is overconfident, above it hedging. The
expected calibration error summarizes a curve as the count-weighted
mean distance from the diagonal.
The same locations count as in `confusion_matrix`: a contact carries
two measurements and no one truth, and locations missing either the
measurement or the prediction are left out.
Only honest on data the model has not seen: at a training location
the claim was fitted to its own outcome. The out-of-fold container
`models.cross_validate` fills is the honest input.
Parameters
----------
container
Any container from the `data` module.
name
The name of a categorical variable holding measurements and
predicted probabilities.
bins
How many bins, or where their edges are. A count gives
**equal-count** bins over the claimed probabilities — for a rare
category the claims pile up near zero and for a dominant one near
one, and equal width would leave most bins holding nothing — so
pass explicit edges to ask for equal width instead.
Returns
-------
panels : list of dict
One per category, in the variable's own order: `label`;
`claimed` and `observed`, one pair per non-empty bin; `count`,
the locations in each; `ece`, the count-weighted mean
`|observed - claimed|`.
Raises
------
TypeError
If the variable holds no measured categories or no components.
ValueError
If no location carries both a measurement and a probability.
"""
var = variable(container, name)
measured, keep = _measured_categories(var, name)
components = getattr(var, "components", None)
called = getattr(var, "predicted", None)
if components is None or called is None:
raise TypeError(
"%s %r holds no categories with predicted probabilities"
% (type(var).__name__, str(name)))
# the predicted label is what says a location was predicted at all --
# a category's `probability` initializes to zero, not to absence, so
# an unpredicted container would otherwise read as a model claiming
# zero everywhere
keep = keep & (_np.asarray(called.to_numpy()).astype(str) != "")
panels = []
for label in getattr(var, "labels", []):
claimed = _np.asarray(
components[label].probability.values, dtype=float).ravel()
rows = keep & _np.isfinite(claimed)
if not _np.any(rows):
raise ValueError(
"no location carries both a measured category and a "
"predicted probability of %r; predict on the data first -- "
"and note this figure is only honest out of fold (see "
"`models.cross_validate`)" % str(name))
hit = (measured[rows] == str(label)).astype(float)
claimed = claimed[rows]
if _np.all(claimed == claimed[0]):
# one constant claim is one point, not nothing to draw
panels.append({"label": str(label),
"claimed": claimed[:1].copy(),
"observed": _np.array([float(hit.mean())]),
"count": _np.array([len(hit)]),
"ece": float(abs(hit.mean() - claimed[0]))})
continue
edges = _bin_edges(claimed, bins)
which = _np.clip(
_np.searchsorted(edges, claimed, side="right") - 1,
0, len(edges) - 2)
count = _np.bincount(which, minlength=len(edges) - 1)
filled = count > 0
safe = _np.where(filled, count, 1)
claim = _np.bincount(which, weights=claimed,
minlength=len(edges) - 1) / safe
freq = _np.bincount(which, weights=hit,
minlength=len(edges) - 1) / safe
ece = float(_np.sum(count[filled]
* _np.abs(freq - claim)[filled])
/ count.sum())
panels.append({"label": str(label), "claimed": claim[filled],
"observed": freq[filled], "count": count[filled],
"ece": ece})
return panels
[docs]
def confusion_matrix(container: "_data._SpatialData", name: str) -> dict:
"""
Measured categories against predicted ones, counted.
Rows are what was measured, columns what the model called there. Only
the unambiguous locations count: a rock type variable carries two
measurements per location, and where they disagree the point lies on a
contact and there is no one truth for the prediction to be right
about; locations missing either the measurement or the prediction are
left out likewise.
Only honest on data the model has not seen: at a training location the
prediction interpolates its own measurement. The out-of-fold container
`models.cross_validate` fills is the honest input.
Parameters
----------
container
Any container from the `data` module.
name
The name of a categorical variable holding measurements and
predictions.
Returns
-------
table : dict
`counts` -- an integer matrix, measured by predicted; `share` --
each cell as a share of its measured row, zero where a row holds
nothing; `labels` -- one name per row and column, the categories
present in either role in the variable's own order, with any
measured category the variable does not model appended after;
`agreement` -- the diagonal's share of the total.
Raises
------
TypeError
If the variable holds no measured and predicted categories.
ValueError
If no location carries both.
"""
var = variable(container, name)
measured, keep = _measured_categories(var, name)
called = getattr(var, "predicted", None)
if called is None:
raise TypeError(
"%s %r holds no predicted categories to compare"
% (type(var).__name__, str(name)))
predicted = _np.asarray(called.to_numpy()).astype(str)
keep = keep & (predicted != "")
if not _np.any(keep):
raise ValueError(
"no location carries both a measured and a predicted category "
"of %r; predict on the data first -- and note this table is "
"only honest out of fold (see `models.cross_validate`)"
% str(name))
measured, predicted = measured[keep], predicted[keep]
present = set(measured) | set(predicted)
labels = [str(label) for label in getattr(var, "labels", [])
if str(label) in present]
labels += sorted(present - set(labels))
index = {label: i for i, label in enumerate(labels)}
k = len(labels)
values, codes = _np.unique(
_np.concatenate([measured, predicted]), return_inverse=True)
codes = _np.array([index[value] for value in values])[codes]
counts = _np.bincount(
codes[:len(measured)] * k + codes[len(measured):],
minlength=k * k).reshape(k, k)
totals = counts.sum(axis=1, keepdims=True)
return {"counts": counts,
"share": counts / _np.where(totals == 0, 1, totals),
"labels": labels,
"agreement": float(_np.trace(counts)) / float(counts.sum())}
[docs]
def spread_check(container: "_data._SpatialData", name: str,
bins: _types.Bins = 8) -> list[dict]:
"""
What a model claims a value's spread is, against what it turned out to be.
A residual holds two things at once -- how wrong the model was about the
ground, and how far the assay fell from the ground -- so it can only be
read against the two together. This lays all three out along the predicted
value: the noise the model fitted, the whole spread it claims, and the
spread the errors actually had.
Reading it: the observed points on the claimed line means calibrated,
below it means hedging, above it means over-confident. The level axis is
what says *which* term is at fault. A warping bends, so the noise grows
with the value while the model's own uncertainty does not, and a shortfall
that widens with the grade is the noise where a flat one is the posterior.
Observed points sitting inside the noise band alone are the plainest case
of all: the fitted noise over-explains the errors by itself.
Only honest on data the model has not seen. At a training location the
model interpolates its own measurement and the residual is not an error.
Parameters
----------
container :
Point data carrying measurements and a prediction.
name : str
The variable.
bins : int or sequence
How many bins, or where their edges are. A count gives **equal-count**
bins; positions are taken as they come, so `np.linspace(...)` is how
to ask for equal width.
Returns
-------
list of dict
One per component, with `label`, the bin bounds `lo`/`hi`, the mean
predicted value in each `centre`, the `count`, the `observed` root
mean square residual and its `observed_error`, and the claimed
`noise` and `total` spreads.
"""
parts = continuous_parts(variable(container, name))
panels = []
for part in parts:
if not part.noise_variance._has_content():
raise ValueError(
"%r carries no noise variance, so there is nothing to check "
"the residuals against; predict with `include_noise=True`, "
"which is the default" % str(part.name))
measured = part.measurements.values.to_numpy().astype(float)
predicted = part.prediction.values.to_numpy().astype(float)
noise = part.noise_variance.values.to_numpy().astype(float)
keep = ~(_np.isnan(measured) | _np.isnan(predicted) | _np.isnan(noise))
if not _np.any(keep):
raise ValueError(
"%r has no location carrying both a measurement and a "
"prediction; this figure is a comparison against what was "
"observed" % str(part.name))
# the model's own uncertainty, in the variable's units, a band of rows
# at a time -- a block model holds more simulations than memory does
store = realization_store(part)
signal = _np.empty(len(measured), dtype=float)
for band in store.row_bands():
signal[band] = _np.var(_np.asarray(store[band, :], dtype=float),
axis=1)
measured, predicted = measured[keep], predicted[keep]
noise, signal = noise[keep], signal[keep]
residual = measured - predicted
edges = _bin_edges(predicted, bins)
index = _np.clip(_np.searchsorted(edges, predicted, side="right") - 1,
0, len(edges) - 2)
panel = {"label": str(part.name), "lo": [], "hi": [], "centre": [],
"count": [], "observed": [], "observed_error": [],
"noise": [], "total": []}
for i in range(len(edges) - 1):
here = index == i
n = int(_np.count_nonzero(here))
if n == 0:
continue
rms = float(_np.sqrt(_np.mean(residual[here] ** 2)))
panel["lo"].append(float(edges[i]))
panel["hi"].append(float(edges[i + 1]))
panel["centre"].append(float(_np.mean(predicted[here])))
panel["count"].append(n)
panel["observed"].append(rms)
# the sampling error of a root mean square over n values, so that
# a gap can be read as real rather than as a short bin
panel["observed_error"].append(rms / _np.sqrt(2 * n))
panel["noise"].append(float(_np.sqrt(_np.mean(noise[here]))))
panel["total"].append(
float(_np.sqrt(_np.mean(noise[here] + signal[here]))))
panels.append({key: (val if key == "label" else _np.asarray(val))
for key, val in panel.items()})
return panels
def _axis_index(container, axis):
"""Which coordinate to slab along, by index or by label."""
labels = [str(label) for label in
(getattr(container, "coordinate_labels", None) or [])]
if isinstance(axis, str):
if axis not in labels:
raise ValueError("no coordinate named %r; found %s"
% (axis, ", ".join(labels) or "none"))
return labels.index(axis), axis
index = int(axis)
if not 0 <= index < container.n_dim:
raise ValueError("axis must be between 0 and %d, got %d"
% (container.n_dim - 1, index))
return index, labels[index] if index < len(labels) else "axis %d" % index
def _slab_edges(positions, bins):
"""Where the slabs are cut, from a count or from the positions.
A count gives equal *width*: a slab is a place, not a value, and the
question a swath asks is where along the deposit the model and the
data part company. `_bin_edges` gives equal count, which is right for
a predicted grade and wrong here.
"""
if _np.ndim(bins) > 0:
edges = _np.unique(_np.asarray(bins, dtype=float))
else:
bins = int(bins)
if bins < 1:
raise ValueError("bins must be at least 1, got %r" % bins)
edges = _np.linspace(positions.min(), positions.max(), bins + 1)
if len(edges) < 2 or edges[-1] <= edges[0]:
raise ValueError("there is nothing to slab: every position is %g"
% edges[0])
return edges
def _slab_of(positions, edges):
"""The slab each position falls in, and whether it falls in any."""
index = _np.searchsorted(edges, positions, side="right") - 1
# the last edge belongs to the last slab, so the range is closed
index = _np.where(positions == edges[-1], len(edges) - 2, index)
inside = (index >= 0) & (index < len(edges) - 1)
return _np.clip(index, 0, len(edges) - 2), inside
def _where_mask(container, where):
"""A stored or given filter as a boolean row mask."""
if where is None:
return _np.ones(container.n_data, dtype=bool)
if isinstance(where, str):
return _np.asarray(container.get_metadata(where)).astype(bool)
mask = _np.asarray(where).astype(bool)
if mask.shape != (container.n_data,):
raise ValueError("where must hold one boolean per location (%d), "
"got shape %s" % (container.n_data, mask.shape))
return mask
def _volume_rows(container):
"""One volume per location, or None where locations have no size."""
try:
volume = block_volume(container)
except TypeError:
return None
if _np.ndim(volume) == 0:
return _np.full(container.n_data, float(volume))
return _np.asarray(volume, dtype=float)
def _declustering(container, weights, values):
"""One weight per location: explicit, then the stored column, then
computed from `values` -- the precedence every declustered consumer
shares. `values` None means nothing to compute from, and the weights
come back None: an unweighted mean, said rather than guessed."""
if weights is not None:
return _np.asarray(weights, dtype=float).ravel(), None, False
column = container.metadata.get("declustering")
if column is not None:
return _np.asarray(column.values, dtype=float), None, True
if values is None:
return None, None, False
points = _np.asarray(container.coordinates, dtype=float)
finite = _np.isfinite(values)
computed = _np.ones(container.n_data)
computed[finite], cell = _geom.declustering_weights(
points[finite], values[finite])
return computed, cell, False
[docs]
def swath(data: "_data._SpatialData", predicted: "_data._SpatialData",
name: str, axis: "int | str" = 0, bins: _types.Bins = 12,
weights: "_types.ArrayLike | None" = None,
where: _types.Where = None,
quantiles=(0.05, 0.95)) -> list[dict]:
"""
The data's mean against the model's, slab by slab along one axis.
The check that *localizes* conditional bias instead of aggregating it
away: a model that runs high in the north and low in the south can
score unbiased overall, and only a mean per slab shows where it drifts.
Two corrections make the comparison fair, and without them the figure
describes the drilling rather than the deposit. The data's means are
**declustered** -- explicit `weights`, else the `"declustering"` column
`container.decluster()` keeps, else cell weights computed here -- so a
crowded patch of holes speaks once. The model's means run only over
the ground its data informs, which `where` names: the boolean column
`assign_from_data` writes, or any mask. Over a block model each block
counts at its own volume.
Where the model carries simulations, the slab mean of every realization
is taken as well and the band between two of their quantiles reported,
which a kriging swath cannot draw.
Parameters
----------
data
The samples, carrying measurements of `name`.
predicted
The grid or block model carrying the model's prediction of it.
name
The variable.
axis
Which coordinate the slabs cut across, by index or by label.
bins
How many slabs, or where their edges are. A count gives slabs of
**equal width** over the range both containers span.
weights
One declustering weight per sample, overriding the stored column.
where
Which locations of `predicted` take part: a boolean mask, or the
name of a boolean metadata column. Everything, by default.
quantiles
The two quantiles of the realizations' slab means drawn as a band.
Returns
-------
list of dict
One per component, with `label`, the slab bounds `lo`/`hi` and
`centre`, `axis`, then `data_mean`, `data_weight` (the summed
weights, an effective count) and `data_count`, `model_mean` and
`model_count`, `band_lo`/`band_hi` (None without simulations),
and `cell`, the declustering cell when the weights were computed
here.
See Also
--------
categorical_swath : the same figure for a categorical variable.
geoml.data.PointData.decluster : the stored weights this reads.
geoml.data.PointData.assign_from_data : the stored reach `where` names.
"""
axis_index, axis_label = _axis_index(data, axis)
if predicted.n_dim != data.n_dim:
raise ValueError("data and predicted must have the same dimension")
data_parts = continuous_parts(variable(data, name))
model_parts = {str(part.name): part
for part in continuous_parts(variable(predicted, name))}
data_pos = _np.asarray(data.coordinates, dtype=float)[:, axis_index]
model_pos = _np.asarray(predicted.coordinates, dtype=float)[:, axis_index]
keep = _where_mask(predicted, where)
if not keep.any():
raise ValueError("`where` leaves no location of the model to compare")
edges = _slab_edges(_np.concatenate([data_pos, model_pos[keep]]), bins)
n_slabs = len(edges) - 1
data_slab, data_in = _slab_of(data_pos, edges)
model_slab, model_in = _slab_of(model_pos, edges)
model_in &= keep
volume = _volume_rows(predicted)
first = data_parts[0].measurements.values.to_numpy().astype(float)
all_weights, cell, _ = _declustering(data, weights, first)
panels = []
for part in data_parts:
label = str(part.name)
if label not in model_parts:
raise ValueError("%r has no component %r in the predicted "
"container" % (name, label))
measured = part.measurements.values.to_numpy().astype(float)
rows = _np.isfinite(measured) & data_in
share = _np.ones(rows.sum()) if all_weights is None \
else all_weights[rows]
weight = _np.bincount(data_slab[rows], weights=share,
minlength=n_slabs)
total = _np.bincount(data_slab[rows], weights=share * measured[rows],
minlength=n_slabs)
with _np.errstate(invalid="ignore", divide="ignore"):
data_mean = _np.where(weight > 0, total / weight, _np.nan)
model_part = model_parts[label]
if getattr(model_part, "prediction", None) is None:
raise ValueError("%r carries no prediction; run the model over "
"the predicted container first" % name)
prediction = model_part.prediction.values.to_numpy().astype(float)
cells = _np.isfinite(prediction) & model_in
size = _np.ones(cells.sum()) if volume is None else volume[cells]
mass = _np.bincount(model_slab[cells], weights=size,
minlength=n_slabs)
held = _np.bincount(model_slab[cells], weights=size * prediction[cells],
minlength=n_slabs)
with _np.errstate(invalid="ignore", divide="ignore"):
model_mean = _np.where(mass > 0, held / mass, _np.nan)
band_lo = band_hi = None
store = realization_store(model_part)
if store.shape[1] > 1:
sums = _np.zeros([n_slabs, store.shape[1]])
for band in store.row_bands():
chunk = _np.asarray(store[band], dtype=float)
inside = cells[band]
size_band = 1.0 if volume is None \
else volume[band][inside][:, None]
_np.add.at(sums, model_slab[band][inside],
size_band * chunk[inside])
with _np.errstate(invalid="ignore", divide="ignore"):
means = _np.where(mass[:, None] > 0, sums / mass[:, None],
_np.nan)
band_lo, band_hi = _np.nanquantile(means, quantiles, axis=1)
panels.append({
"label": label, "axis": axis_label,
"lo": edges[:-1], "hi": edges[1:],
"centre": 0.5 * (edges[:-1] + edges[1:]),
"data_mean": data_mean, "data_weight": weight,
"data_count": _np.bincount(data_slab[rows], minlength=n_slabs),
"model_mean": model_mean,
"model_count": _np.bincount(model_slab[cells], minlength=n_slabs),
"band_lo": band_lo, "band_hi": band_hi, "cell": cell})
return panels
[docs]
def categorical_swath(data: "_data._SpatialData",
predicted: "_data._SpatialData", name: str,
axis: "int | str" = 0, bins: _types.Bins = 12,
weights: "_types.ArrayLike | None" = None,
where: _types.Where = None) -> dict:
"""
The data's category shares against the model's, slab by slab.
The categorical `swath`: per slab, the declustered share of each
category among the samples against the model's mean predicted
probability of it -- its expected share -- over the locations `where`
names, each block at its own volume. The shares of a slab sum to one on
both sides, which is what lets them stack.
The declustering weights are explicit `weights` or the stored
`"declustering"` column; a categorical variable offers no values to
compute a cell from, so without either the shares are raw, which is
said in the result.
Parameters
----------
data, predicted, name, axis, bins, weights, where
As in `swath`.
Returns
-------
dict
`labels`, `axis`, the slab bounds `lo`/`hi` and `centre`,
`data_share` and `model_share` (`(n_slabs, n_labels)`),
`data_weight`, `data_count`, `model_count`, and `declustered`.
"""
axis_index, axis_label = _axis_index(data, axis)
if predicted.n_dim != data.n_dim:
raise ValueError("data and predicted must have the same dimension")
values, measured, labels = category_values(variable(data, name))
components = variable(predicted, name).components
data_pos = _np.asarray(data.coordinates, dtype=float)[:, axis_index]
model_pos = _np.asarray(predicted.coordinates, dtype=float)[:, axis_index]
keep = _where_mask(predicted, where)
if not keep.any():
raise ValueError("`where` leaves no location of the model to compare")
edges = _slab_edges(_np.concatenate([data_pos, model_pos[keep]]), bins)
n_slabs = len(edges) - 1
data_slab, data_in = _slab_of(data_pos, edges)
model_slab, model_in = _slab_of(model_pos, edges)
model_in &= keep
volume = _volume_rows(predicted)
all_weights, _, _ = _declustering(data, weights, None)
rows = measured & data_in
share = _np.ones(rows.sum()) if all_weights is None else all_weights[rows]
weight = _np.bincount(data_slab[rows], weights=share, minlength=n_slabs)
data_share = _np.full([n_slabs, len(labels)], _np.nan)
model_share = _np.full([n_slabs, len(labels)], _np.nan)
model_count = _np.zeros(n_slabs, dtype=int)
for k, label in enumerate(labels):
hit = (values[rows] == label).astype(float)
with _np.errstate(invalid="ignore", divide="ignore"):
data_share[:, k] = _np.bincount(
data_slab[rows], weights=share * hit, minlength=n_slabs) \
/ weight
if label not in components:
raise ValueError("%r has no category %r in the predicted "
"container" % (name, label))
probability = components[label].probability.values.to_numpy() \
.astype(float)
cells = _np.isfinite(probability) & model_in
size = _np.ones(cells.sum()) if volume is None else volume[cells]
mass = _np.bincount(model_slab[cells], weights=size,
minlength=n_slabs)
with _np.errstate(invalid="ignore", divide="ignore"):
model_share[:, k] = _np.bincount(
model_slab[cells], weights=size * probability[cells],
minlength=n_slabs) / mass
model_count = _np.maximum(
model_count, _np.bincount(model_slab[cells], minlength=n_slabs))
return {"labels": labels, "axis": axis_label,
"lo": edges[:-1], "hi": edges[1:],
"centre": 0.5 * (edges[:-1] + edges[1:]),
"data_share": data_share, "model_share": model_share,
"data_weight": weight,
"data_count": _np.bincount(data_slab[rows], minlength=n_slabs),
"model_count": model_count,
"declustered": all_weights is not None}
[docs]
def proportions(data: "_data._SpatialData", predicted: "_data._SpatialData",
name: str, weights: "_types.ArrayLike | None" = None,
where: _types.Where = None) -> dict:
"""
The data's category shares against the model's, over the whole model.
`categorical_swath` without the slabs: the declustered share of each
category among the samples against the model's mean predicted
probability of it -- its expected share -- over the locations `where`
names, each block at its own volume. Both sides sum to one. It is the
global check the row-normalized confusion matrix cannot make: a model
can place every measured sample in its right category and still call
the dominant rock over ground the data never reached.
Parameters
----------
data, predicted, name, weights, where
As in `swath`.
Returns
-------
dict
`labels`, `data_share` and `model_share` (one per label),
`data_count`, `model_count`, and `declustered`.
"""
values, measured, labels = category_values(variable(data, name))
components = variable(predicted, name).components
keep = _where_mask(predicted, where)
if not keep.any():
raise ValueError("`where` leaves no location of the model to compare")
volume = _volume_rows(predicted)
all_weights, _, _ = _declustering(data, weights, None)
share = _np.ones(measured.sum()) if all_weights is None \
else all_weights[measured]
with _np.errstate(invalid="ignore", divide="ignore"):
data_share = _np.array([
share[values[measured] == label].sum() / share.sum()
for label in labels])
model_share = _np.full(len(labels), _np.nan)
model_count = 0
for k, label in enumerate(labels):
if label not in components:
raise ValueError("%r has no category %r in the predicted "
"container" % (name, label))
probability = components[label].probability.values.to_numpy() \
.astype(float)
cells = _np.isfinite(probability) & keep
size = _np.ones(cells.sum()) if volume is None else volume[cells]
with _np.errstate(invalid="ignore", divide="ignore"):
model_share[k] = (size * probability[cells]).sum() / size.sum()
model_count = max(model_count, int(cells.sum()))
return {"labels": labels, "data_share": data_share,
"model_share": model_share,
"data_count": int(measured.sum()), "model_count": model_count,
"declustered": all_weights is not None}
def _drillhole_column(container, key, what):
column = container.metadata.get(key)
if column is None:
raise ValueError("%s must carry the %r metadata column a drillhole "
"conversion writes" % (what, key))
return _np.asarray(column.values)
[docs]
def variogram(container: "_data._SpatialData", name: str,
n_lags: int = 15, max_lag: "float | None" = None,
direction: "_types.ArrayLike | None" = None,
tolerance: float = 45.0, max_pairs: int = 2_000_000,
residuals: bool = False,
decluster: "bool | float" = True) -> list[dict]:
"""
The data's spatial structure, against the fan the simulations reproduce.
The experimental semivariogram of the measured values, and one curve per
realization computed on the same pairs: a model that learned the spatial
structure scatters its realizations *around* the data's curve, a kernel
too smooth sags below it at short lags, and a nugget fitted into the
range lifts it there. Neither shows in the marginal checks
(`accuracy`, `spread_check`), which is what this figure exists for.
The two curves are put on the same footing before being compared. The
measurements carry the likelihood noise and the stored realizations do
not (`predict` integrates it out), so the fan is raised by what that
noise adds to a semivariogram. For noise independent between locations
that is exactly the pair-averaged `(var_i + var_j) / 2`, taken from the
`noise_variance` column and added bin by bin -- no draw, no seed, and
exact in expectation rather than approximated. Without the correction
the fan sits a nugget below the data at every lag and every model looks
over-smooth. A container predicted without `include_noise` carries no
such column and its fan is left where it is.
Pairs are **declustered by default**, which is the other half of making
the comparison fair and matters more than it sounds. Samples follow the
ore, so an experimental variogram computed on raw pairs describes the
sampling as much as the field: on the bundled Walker Lake set, whose
exhaustive truth is known, the raw sample curve runs 1.4 times the true
variogram at the sill and 2.3 times at the shortest lag, and a model
matching the field perfectly would look far too smooth against it. Each
pair is therefore weighted by `w_i * w_j` from
:func:`geoml.math.geometry.declustering_weights`, and the sill likewise,
with the same weights used for the fan so that both sides estimate the
same thing. Pass `decluster=False` for the raw curve, or a number to fix
the cell size rather than let it be chosen.
With `residuals=True` the variogram is of `measured - predicted` instead
and the fan is omitted: structure left in the residuals is structure the
model missed, and on cross-validated predictions (`models.cross_validate`)
it is honest. Whatever the flag, at a training location a model
interpolates its own measurement, so the fan is only worth reading
against data the model has not seen or with the data curve as the anchor.
Parameters
----------
container :
Point data carrying measurements (and simulations, for the fan).
name : str
The variable.
n_lags : int
Number of equal-width lag bins between zero and `max_lag`.
max_lag : float, optional
The longest separation considered. Half the bounding-box diagonal by
default.
direction : array-like, optional
A direction vector for a directional variogram. Omnidirectional when
absent. Call once per direction to compare them -- the anisotropy
ellipsoid's principal axes are the ones worth asking about.
tolerance : float
Angular tolerance around `direction`, in degrees.
max_pairs : int
The pair budget. Past it the locations are strided down --
deterministic, so two calls agree.
residuals : bool
Variogram of `measured - predicted` rather than of the measurements,
with no fan.
decluster : bool or float
Weight pairs by cell-declustering weights, so that the curve
estimates the field's variogram rather than the sampling's. `True`
chooses the cell size, a number fixes it, `False` leaves the pairs
raw.
Returns
-------
list of dict
One per component, with `label`, the bin centres `lag`, the pair
`count` per bin, the data curve `data`, the sample variance `sill`,
and `realizations` -- an `(n_realizations, n_lags)` array, or None
when there is nothing simulated or `residuals` was asked. `noise` is
what was added to the fan per bin to put it on the measurements'
footing, or None when the container could not say. `cell` is the
declustering cell used, or None when the pairs were left raw.
`score` is :func:`geoml.metrics.variogram_score` over the same
locations and weights, or None when there was no fan to score.
Notes
-----
The `score` is the figure's verdict as one number, and it is a ranking
rather than a measurement: unlike the curves it cannot be put on the
measurements' footing, since `|difference| ** p` is not a second moment
and has no constant to add. It never reaches zero, and only comparisons
between models on the same data mean anything. It is also taken over
every pair of the locations kept, not only the binned ones, so it does
not answer to `max_lag` or move when `direction` does.
"""
parts = continuous_parts(variable(container, name))
coords = _np.asarray(container.coordinates, dtype=float)
n = coords.shape[0]
# the pair budget, met by striding the locations down -- deterministic,
# where sampling pairs at random would put a seed inside a figure
most = int(_np.floor(_np.sqrt(2 * max_pairs))) + 1
stride = max(1, int(_np.ceil(n / most)))
rows = _np.arange(0, n, stride)
points = coords[rows]
if max_lag is None:
box = container.bounding_box
span = _np.asarray(box.max, dtype=float).ravel() \
- _np.asarray(box.min, dtype=float).ravel()
max_lag = 0.5 * float(_np.sqrt((span ** 2).sum()))
distance = _spatial.distance.pdist(points)
i_idx, j_idx = _np.triu_indices(len(points), k=1)
keep = distance <= max_lag
if direction is not None:
u = _np.asarray(direction, dtype=float).ravel()
u = u / _np.linalg.norm(u)
separation = points[j_idx] - points[i_idx]
along = _np.abs(separation @ u) / _np.maximum(distance, 1e-30)
keep &= along >= _np.cos(_np.deg2rad(tolerance))
i_idx, j_idx, distance = i_idx[keep], j_idx[keep], distance[keep]
edges = _np.linspace(0.0, max_lag, n_lags + 1)
lag_bin = _np.clip(_np.searchsorted(edges, distance, side="right") - 1,
0, n_lags - 1)
centre = 0.5 * (edges[:-1] + edges[1:])
# declustering weights, one per location, shared by every component: the
# sampling geometry is the same for all of them. The column
# `container.decluster()` keeps is preferred — one consistent set of
# weights for every declustered consumer — and only in its absence is
# the cell swept here, from the first component that has values
weights, cell, stored = None, None, False
if decluster is not False:
column = container.metadata.get("declustering") \
if decluster is True else None
if column is not None:
weights = _np.asarray(column.values, dtype=float)[rows]
stored = True
else:
first = parts[0].measurements.values.to_numpy() \
.astype(float)[rows]
weights, cell = _geom.declustering_weights(
points, first,
cell=None if decluster is True else float(decluster))
def _binned(values, per_pair):
"""One weighted average per lag bin, over the pairs that are usable.
`per_pair` builds the quantity being averaged from the pair indices,
so the same machinery serves the semivariogram, the declustered one
and the noise correction.
"""
ok = _np.isfinite(values)
pair_ok = ok[i_idx] & ok[j_idx]
pi, pj, pb = i_idx[pair_ok], j_idx[pair_ok], lag_bin[pair_ok]
share = _np.ones(pi.shape) if weights is None \
else weights[pi] * weights[pj]
count = _np.bincount(pb, minlength=n_lags)
mass = _np.bincount(pb, weights=share, minlength=n_lags)
total = _np.bincount(pb, weights=share * per_pair(pi, pj),
minlength=n_lags)
averaged = total / _np.maximum(mass, _np.finfo(float).tiny)
averaged[count == 0] = _np.nan
return averaged, count
def curve(values):
return _binned(values,
lambda pi, pj: 0.5 * (values[pi] - values[pj]) ** 2)
def noise_lift(variance):
"""What independent measurement noise adds to a semivariogram.
A measurement is the ground plus an error of its own at each
location, so a pair differs by `(g_i - g_j) + (e_i - e_j)`. The
cross term averages away and the rest is `(var_i + var_j) / 2` per
pair, which is what a realization of the ground has to be raised by
before it can be laid against the data's curve. Only the variance
enters, never the shape of the noise, which is why nothing has to be
drawn here.
"""
return _binned(variance,
lambda pi, pj: 0.5 * (variance[pi] + variance[pj]))[0]
panels = []
for part in parts:
measured = part.measurements.values.to_numpy().astype(float)[rows]
if residuals:
predicted = part.prediction.values.to_numpy().astype(float)[rows]
values = measured - predicted
else:
values = measured
if not _np.any(_np.isfinite(values)):
raise ValueError(
"%r has nothing to compute a variogram from"
% str(part.name))
gamma, count = curve(values)
fan = None
lift = None
score = None
store = getattr(part, "simulations", None)
if not residuals and store is not None \
and len(getattr(store, "shape", ())) == 2 \
and store.shape[1] > 1:
# the realizations are of the ground; the data carries the
# measurement noise, so the fan is raised onto its footing
noise = getattr(part, "noise_variance", None)
noise = None if noise is None else noise.values
if noise is not None and getattr(noise, "shape", None) is not None:
noise = noise.to_numpy().astype(float).ravel()[rows]
if _np.any(_np.isfinite(noise)):
lift = noise_lift(noise)
fan = _np.full((store.shape[1], n_lags), _np.nan)
draws = []
for r in range(store.shape[1]):
# one realization is one column, read without materializing
# the store -- the same discipline as everywhere else here
column = _np.asarray(store[:, r], dtype=float).ravel()[rows]
fan[r] = curve(column)[0]
if lift is not None:
fan[r] = fan[r] + lift
# kept for the score, which needs the realizations side by
# side rather than one at a time. Only the strided rows, so
# this is bounded by the pair budget however large the
# container is, and it costs no reading that the fan has not
# already done.
draws.append(column)
keep = _np.isfinite(values)
if keep.sum() > 1:
if stored and weights is not None:
# the same weights the curve used, exactly
score = _gmet.variogram_score(
values[keep], _np.column_stack(draws)[keep],
weights=weights[keep])
else:
score = _gmet.variogram_score(
values[keep], _np.column_stack(draws)[keep],
coordinates=points[keep],
decluster=cell if cell else False)
# the sill is weighted the same way, or the line drawn across the
# figure would belong to a different population from the curve
usable = _np.isfinite(values)
if weights is None:
sill = float(_np.var(values[usable]))
else:
share = weights[usable]
centred = values[usable] \
- (share * values[usable]).sum() / share.sum()
sill = float((share * centred ** 2).sum() / share.sum())
panels.append({
"label": str(part.name),
"lag": centre,
"count": count,
"data": gamma,
"sill": sill,
"realizations": fan,
"noise": lift,
"cell": cell,
"score": score,
})
return panels
[docs]
def moving_average(values: _types.ArrayLike, window: int
) -> "tuple[_types.IndexArray, _types.FloatArray]":
"""
The running mean of `values`, and where each point belongs.
Returns the positions as well: a mean over `window` points only exists
once there are that many, so the curve starts later than the one it
smooths and has to be drawn against its own x.
"""
values = _np.asarray(values, dtype=float)
window = int(window)
if window <= 1 or window > len(values):
return _np.arange(1, len(values) + 1), values.copy()
kernel = _np.ones(window) / window
smoothed = _np.convolve(values, kernel, mode="valid")
return _np.arange(window, len(values) + 1), smoothed
[docs]
def training_curve(model, window: "int | None" = None):
"""
The training log, with a running mean over it.
Each entry is one optimizer step -- and under `train_svi` one *batch*, so
the curve is noisy by construction: every value is the ELBO estimated from
a sample of the data and a sample of the latent variables. The mean is what
says whether it is still climbing.
Parameters
----------
window : int
Points to average over. Defaults to a fiftieth of the log, which keeps
the smoothing proportionate to however long training ran, but never
fewer than five: a mean of two or three smooths nothing, and draws a
second line on top of the first saying the same thing. A window longer
than the log leaves it alone, and then there is no second line at all.
"""
values = _np.asarray(getattr(model, "training_log", []), dtype=float)
if len(values) == 0:
raise ValueError("this model has not been trained yet")
if window is None:
window = int(max(5, len(values) // 50))
position, smoothed = moving_average(values, window)
return {"iteration": _np.arange(1, len(values) + 1), "value": values,
"smooth_iteration": position, "smooth": smoothed,
"window": window}
[docs]
def realizations(var):
"""
A variable's simulations, or its prediction as the only one there is.
Returns `(n_data, n_realizations)` either way, so that whatever reads it
does not have to care which it got.
"""
values = None
try:
values = _np.asarray(var.get_simulations(), dtype=float)
except Exception:
# a vector variable with nothing simulated fails stacking its
# components rather than handing anything back
pass
# and a continuous one hands back `asarray(None)`, a dimensionless nan,
# which reads as a value right up until it is asked for its length
if values is None or values.ndim == 0:
prediction = getattr(var, "prediction", None)
if prediction is None:
raise ValueError(
"%r carries neither simulations nor a prediction" % var.name)
values = _np.asarray(prediction.values.to_numpy(), dtype=float)
return values.reshape(len(values), -1)
[docs]
def realization_store(var):
"""
A variable's realizations as a store to read in bands, not as an array.
The same `(n_data, n_realizations)` that `realizations` hands back, without
asking memory for all of it at once. A block model's simulations are the
one thing a container holds that will not fit -- hundreds of gigabytes is
an ordinary size for them -- and every use of them here is a reduction over
locations, which never needs more than a band of rows at a time.
A variable with no simulations falls back on its prediction, which is a
single column and already in RAM, and comes back as a one-band store so
that the caller has only the one path to write.
"""
store = getattr(var, "simulations", None)
if store is not None and len(getattr(store, "shape", ())) == 2:
return store
return _storage.ArrayStore.from_numpy(realizations(var))
[docs]
def block_density(container, density, n_data, n_realizations):
"""
A density per block, to turn a volume into a tonnage.
`density` may be a number, the name of a metadata column, or the name of a
`ContinuousVariable`. Only the last of these can be uncertain, and when it
is, its realizations are matched one to one with the grade's: simulation
`i` of the density belongs with simulation `i` of the grade, and pairing
them any other way would invent a correlation that was never modelled.
Comes back as a single number where the density is one, and otherwise as a
store to read in bands alongside the grade. A simulated density is exactly
as big as the grade, so materializing it would give back everything not
materializing the grade saves.
"""
if density is None:
return 1.0
if isinstance(density, (int, float)):
return float(density)
if density in container.metadata:
values = _np.asarray(container.get_metadata(density), dtype=float)
return _storage.ArrayStore.from_numpy(values.reshape(n_data, 1))
if density in container.variables:
store = realization_store(variable(container, density))
if store.shape[1] not in (1, n_realizations):
raise ValueError(
"%r carries %d realizations and the grade carries %d; they "
"have to be matched one to one, or the density has to be the "
"same in all of them"
% (density, store.shape[1], n_realizations))
return store
raise KeyError(
"nothing named %r to take a density from; metadata holds %s and the "
"variables are %s"
% (density, ", ".join(sorted(container.metadata)) or "nothing",
", ".join(sorted(container.variables))))
[docs]
def containing_variable(container, var):
"""The variable `var` is a component of, if it is a component of one."""
for candidate in container.variables.values():
components = getattr(candidate, "components", None) or {}
if any(component is var for component in components.values()):
return candidate
return None
def _column_of(obj, name):
"""
An attribute of `obj` that is a column of values, or None.
A column that was never filled does not count as found. A vector
variable's components are ordinary continuous variables, so each carries a
`latent_variance` whether or not a model ever wrote one; stopping at the
empty one would hide the column on the parent that does hold something,
and then every block would fail the filter for being NaN.
"""
attribute = getattr(obj, name, None)
if attribute is None or not hasattr(attribute, "values"):
return None
values = _np.asarray(attribute.values.to_numpy(), dtype=float)
return None if _np.all(_np.isnan(values)) else values
[docs]
def uncertainty_values(container, name, variables=()):
"""
A number per location saying how much the model doubts itself there.
Which number that is depends on what was modelled, so it is named rather
than guessed at. `name` may be:
- **an array**, one value per location. Whatever the number is and wherever
it came from, it can always be handed over directly, which is the way out
of every case the names below do not reach.
- **a path**, naming exactly where to read it -- `"Elements/uncertainty"`
while the grade is one of its components, or a column belonging to some
other variable entirely. (The old dotted `"Variable.column"` is refused
with the replacement spelled out: a `.` inside a label would make a
wrong guess look like a working one.)
- **a bare name**, looked for on each of `variables` in turn and then among
the metadata columns.
`variables` is what a bare name is tried against: the grade, and then the
variable containing it. That second one is not a nicety. A component of a
vector variable has `latent_variance` set to None and no uncertainty of its
own -- the column that exists belongs to the parent, so grading `Zn` and
asking for `"uncertainty"` has to reach `Elements` to find anything.
"""
if not isinstance(name, str):
values = _np.asarray(name, dtype=float).ravel()
if len(values) != container.n_data:
raise ValueError(
"an uncertainty given as an array needs one value per "
"location: got %d for %d" % (len(values), container.n_data))
return values
if _data.PATH_SEP in name:
found = container.get(name) # says what is there when it is not
if not isinstance(found, _data._Attribute):
raise KeyError(
"%r is a %s, not a column of values"
% (name, type(found).__name__))
values = _np.asarray(found.values.to_numpy(), dtype=float)
if _np.all(_np.isnan(values)):
raise KeyError("%r was never filled" % name)
return values
if "." in name:
owner, _, column = name.partition(".")
raise KeyError(
"%r is no longer accepted; use the path %r"
% (name, "%s/%s" % (owner, column)))
for var in variables:
values = _column_of(var, name)
if values is not None:
return values
if name in container.metadata:
return _np.asarray(container.get_metadata(name), dtype=float)
carried = sorted({key for var in variables for key, value
in vars(var).items() if hasattr(value, "values")})
raise KeyError(
"nothing named %r to take an uncertainty from. The variables in hand "
"carry %s; the metadata holds %s. Name a column on another variable "
"by its path, 'Variable/column', or pass the values themselves"
% (name, ", ".join(carried) or "no columns",
", ".join(sorted(container.metadata)) or "nothing"))
def _grade_band(store, band, keep):
"""One band of realizations, with the blocks that were filtered out gone."""
values = _np.asarray(store[band], dtype=float)
return values if keep is None else values[keep[band]]
def _per_block(value, band, keep):
"""One quantity for the blocks of a band: a number, a column, or a store.
A number stays a number -- a model whose blocks are all the same size has
one volume, not a copy of it per block -- and anything longer is banded and
filtered like the grade beside it.
"""
if _np.ndim(value) == 0:
return float(value)
values = _np.asarray(value[band], dtype=float)
if values.ndim == 1:
values = values[:, None]
if keep is not None:
values = values[keep[band]]
return values
def _mass_band(volume, density, band, keep, shape):
"""What every block of a band weighs: its size times what fills it.
Both may be one number for the whole model or one per block, and either
way the product broadcasts over the realizations.
"""
return _np.broadcast_to(
_per_block(volume, band, keep) * _per_block(density, band, keep),
shape)
[docs]
def block_volume(container):
"""What one block is worth, as a number or as one value per block.
A regular grid has a single spacing and so a single volume; a `BlockSet3D`
carries a size per block and answers with a column. Kept apart from the
density because only one of them is a property of the container.
"""
volume = getattr(container, "block_volume", None)
if volume is not None:
return _np.asarray(volume, dtype=float)
step = getattr(container, "step_size", None)
if step is None:
raise TypeError(
"grade-tonnage needs to know how big a block is, and %s says "
"neither `block_volume` nor `step_size`"
% type(container).__name__)
return float(_np.prod(step))
def _cutoff_range(store, bands, keep, name):
"""The span of the finite realizations, in one pass over the store."""
low, high = _np.inf, -_np.inf
for band in bands:
values = _grade_band(store, band, keep)
finite = values[_np.isfinite(values)]
if finite.size > 0:
low = min(low, float(finite.min()))
high = max(high, float(finite.max()))
if not _np.isfinite(low):
raise ValueError("%r holds no values to cut" % name)
return low, high
[docs]
def grade_tonnage(container, name, density=None, cutoffs=30,
uncertainty=None, max_uncertainty=None):
"""
How much material sits above a cut-off, and how good it is.
Each block contributes its volume, or its mass where a density is given,
to every cut-off its grade clears. Simulations are carried through
separately rather than averaged first: the curve of the mean model is not
the mean of the curves, since a cut-off is a threshold and averaging either
side of it gives different answers.
The simulations are read a band of blocks at a time and never held whole:
a block model runs to hundreds of gigabytes of them, and what comes out is
one small number per cut-off per realization. Each block is placed at the
highest cut-off it clears and the curve is the running total from the top
down, so the cost is one pass over the grade rather than one per cut-off.
Giving `cutoffs` as values rather than as a count saves the pass that would
otherwise be needed to find their range.
Parameters
----------
container
A gridded container -- the volume of a block comes from its spacing.
name : str
The variable to take as the grade.
density : float or str
A number, a metadata column, or a `ContinuousVariable`. Without one the
curve is in volume.
cutoffs : int or array-like
The grades to cut at, or how many of them to spread evenly across the
range of the data.
uncertainty : str or array
Where to read how sure the model is at each block: a column name, a
a path (`"Variable/column"`) naming which variable it belongs to, or the values
themselves. See `uncertainty_values`. A bare name is looked for on the
grade, then on the variable containing it, then in the metadata.
max_uncertainty : float
Blocks doubted more than this are left out altogether -- not counted
at any cut-off, and not counted towards the grade above one. A block
the model cannot speak for is not tonnage.
Returns
-------
dict
`cutoff` `(n_cutoffs,)`; `tonnage`, `grade` and `metal`
`(n_cutoffs, n_realizations)`; `extent`, what is accumulated above
the cut-off -- a block's own extent, or `"mass"` where a density
was given; `unit`, what the grade is measured in where it says so;
and `kept` and `total`, the blocks that survived the uncertainty
filter and the blocks there were.
"""
volume = block_volume(container)
# a block of a two-dimensional grid has an area, not a volume, and saying
# otherwise on the axis of a figure someone is reading off is not harmless
extent = {1: "length", 2: "area", 3: "volume"}.get(
container.n_dim, "volume")
var = variable_or_component(container, name)
grade = realization_store(var)
n_data, n_realizations = grade.shape
mass_per_volume = block_density(
container, density, n_data, n_realizations)
total, keep = n_data, None
if max_uncertainty is not None:
if uncertainty is None:
raise ValueError(
"there is no uncertainty column to filter by; name one when "
"the Explorer is built, or pass uncertainty=")
owner = containing_variable(container, var)
doubt = uncertainty_values(
container, uncertainty,
[var] + ([owner] if owner is not None else []))
keep = _np.isfinite(doubt) & (doubt <= float(max_uncertainty))
if not _np.any(keep):
raise ValueError(
"no block is certain enough to keep: the smallest %r is %g, "
"above the %g asked for"
% (uncertainty, _np.nanmin(doubt), max_uncertainty))
kept = total if keep is None else int(_np.count_nonzero(keep))
bands = grade.row_bands()
if isinstance(cutoffs, (int, _np.integer)):
low, high = _cutoff_range(grade, bands, keep, name)
cutoffs = _np.linspace(low, high, int(cutoffs))
cutoffs = _np.asarray(cutoffs, dtype=float)
# What each band adds is the mass sitting *at* a cut-off -- above it and
# below the next one up. Summing those from the top down at the end turns
# them into the mass above each cut-off, which is the curve.
n_cutoffs = len(cutoffs)
at_cutoff = _np.zeros([n_cutoffs, n_realizations])
metal_at_cutoff = _np.zeros_like(at_cutoff)
columns = _np.arange(n_realizations)
for band in bands:
values = _grade_band(grade, band, keep)
if values.shape[0] == 0:
continue
weight = _mass_band(volume, mass_per_volume, band, keep, values.shape)
# the highest cut-off each block clears; -1 for one that clears none,
# which is where a block with no value belongs as well
index = _np.searchsorted(cutoffs, values, side="right") - 1
index[~_np.isfinite(values)] = -1
above = index >= 0
# one bin per (cut-off, realization) pair, so every realization is
# accumulated in the same pass over the band
binned = (index * n_realizations + columns)[above]
values, weight = values[above], weight[above]
at_cutoff += _np.bincount(
binned, weights=weight,
minlength=n_cutoffs * n_realizations
).reshape(n_cutoffs, n_realizations)
metal_at_cutoff += _np.bincount(
binned, weights=weight * values,
minlength=n_cutoffs * n_realizations
).reshape(n_cutoffs, n_realizations)
# the volume is already in the weights: with a block size per block it
# cannot be saved for the end the way one shared size could
tonnage = _np.cumsum(at_cutoff[::-1], axis=0)[::-1]
metal = _np.cumsum(metal_at_cutoff[::-1], axis=0)[::-1]
with _np.errstate(invalid="ignore", divide="ignore"):
mean_grade = _np.where(tonnage > 0, metal / tonnage, _np.nan)
return {"cutoff": cutoffs, "tonnage": tonnage, "grade": mean_grade,
"metal": metal,
# what is being accumulated above the cut-off -- a length, an
# area, a volume, or a mass where a density was given
"extent": extent if density is None else "mass",
# and what the grade itself is measured in, where it says so
"unit": getattr(var, "unit", None),
"kept": kept, "total": total}
[docs]
def dispersion_by_support(container: "_data.BlockSet3D", name: str,
component: "str | None" = None) -> "list[dict]":
"""
The within-block standard deviation at every block size a block set has.
Each block's `dispersion` says how much the ground varies inside it, on
the block's own support. Merging the blocks into their parents, level by
level up to the coarsest, gives the same number at every size the
lattice has, so it can be read against the block size.
A parent is never predicted; it is put together from the blocks inside
it. Its dispersion is theirs, volume-weighted, plus how much their block
values differ among themselves -- taken realization by realization, a
realization being the one thing that comes across a regrouping exactly,
and averaged over the realizations as a block's own dispersion is. Every
realization the variable holds is used, read a band of blocks at a time.
A parent missing any of its ground -- a block never predicted, such as
one a `where=` filter left out -- is left out at that size and at every
size above it, as `BlockSet3D.group` refuses a partial family.
Parameters
----------
container
The block set.
name
The continuous variable.
component
One component of a vector or compositional variable; every one by
default.
Returns
-------
list of dict
One per component, in label order, with `name`, `label` and `sizes`:
one dict per level, the finest first, holding `level`, `size` (the
block size along each axis), `deviation` (the within-block standard
deviation of every block of that size), `depth` (how many times the
refinement split inside each: 0 for a block it left whole), `rms`
(the root mean square of `deviation`, the dispersion of the ground
within blocks of that size), `count`, `share` (of the set's volume
the blocks cover) and `left_out` (blocks of that size missing some
of their ground).
Raises
------
ValueError
If the container is not a `BlockSet3D`, or the variable carries no
dispersion or no simulations.
KeyError
If `component` names no component of the variable.
"""
if not isinstance(container, _data.BlockSet3D):
raise ValueError(
"the dispersion by block size merges blocks into their parents, "
"which needs a BlockSet3D; got a %s" % type(container).__name__)
var = variable(container, name)
parts = continuous_parts(var)
if component is not None:
parts = [part for part in parts if str(part.name) == str(component)]
if not parts:
raise KeyError("no component %r in %r; found %s"
% (component, name, component_names(var)))
return [{"name": str(part.name), "label": axis_label(part),
"sizes": _dispersion_sizes(container, part)} for part in parts]
def _dispersion_sizes(blocks, part):
"""`dispersion_by_support` for one column."""
dispersion = _np.asarray(part.dispersion.values.to_numpy(), dtype=float)
leaf = _np.isfinite(dispersion)
if not _np.any(leaf):
raise ValueError(
"%r carries no dispersion; a model writes one when it predicts "
"onto blocks, and a derived variable has none" % str(part.name))
store = getattr(part, "simulations", None)
if store is None or len(getattr(store, "shape", ())) != 2:
raise ValueError(
"%r carries no simulations, which a parent's dispersion is put "
"together from" % str(part.name))
level = _np.asarray(blocks.level)
ratio = _np.array(blocks.discretization, dtype=_np.int64)
volume = _np.prod(blocks._size, axis=1) # in base cells
finest = int(level[leaf].max())
n_sim = int(store.shape[1])
def cells(coarse):
"""A block's volume at level `coarse`, in base cells."""
return int(_np.prod(blocks._coarse_size // ratio ** coarse))
# every block finer than a level, numbered by its ancestor there
parents = {}
for coarse in range(finest):
inside = leaf & (level > coarse)
_, inverse = _np.unique(blocks._ancestor(coarse)[inside],
return_inverse=True)
index = _np.full(blocks.n_data, -1, dtype=_np.int64)
index[inside] = inverse.ravel()
parents[coarse] = (index, int(index.max()) + 1)
# The realizations enter shifted by their mean: the spread between
# blocks is a difference of two sums of squares, and a grade far from
# zero would otherwise leave it at the mercy of rounding.
prediction = _np.asarray(part.prediction.values.to_numpy(), dtype=float)
finite = _np.isfinite(prediction) & leaf
shift = float(prediction[finite].mean()) if _np.any(finite) else 0.0
# A parent's realizations are sums over blocks that sit anywhere in the
# store -- a split appends its children at the end -- so they are
# gathered band by band into one row per parent and realization. Where
# that many rows would not fit, the realizations are taken a slice at a
# time, each slice another pass over the store.
rows = sum(n for _, n in parents.values())
per_pass = max(1, int(_storage.DEFAULT_THRESHOLD // (8 * max(rows, 1))))
square = _np.zeros(blocks.n_data)
between = {coarse: _np.zeros(n) for coarse, (_, n) in parents.items()}
for first in range(0, n_sim, per_pass):
columns = slice(first, min(first + per_pass, n_sim))
sums = {coarse: _np.zeros((n, columns.stop - columns.start))
for coarse, (_, n) in parents.items()}
for band in store.row_bands():
keep = _np.flatnonzero(leaf[band])
if keep.size == 0:
continue
values = _np.asarray(store[band, columns], dtype=float)
values = values[keep] - shift
keep = keep + band.start
square[keep] += _np.sum(values ** 2, axis=1)
for coarse, (index, n) in parents.items():
parent = index[keep]
used = parent >= 0
if not _np.any(used):
continue
count = int(_np.count_nonzero(used))
gather = _sparse.csr_matrix(
(volume[keep[used]] / cells(coarse),
(parent[used], _np.arange(count))), shape=(n, count))
sums[coarse] += gather @ values[used]
for coarse in parents:
between[coarse] += _np.sum(sums[coarse] ** 2, axis=1)
total = float(volume.sum())
sizes = []
for coarse in range(finest, -1, -1):
whole = leaf & (level == coarse)
variance = [dispersion[whole]]
depth = [_np.zeros(int(_np.count_nonzero(whole)), dtype=_np.int64)]
left_out = int(_np.count_nonzero(~leaf & (level == coarse)))
if coarse in parents:
index, n = parents[coarse]
inside = index >= 0
weight = volume[inside] / cells(coarse)
# the law of total variance, realization by realization: what
# varies inside the blocks, plus how their values differ
within = _np.bincount(index[inside], minlength=n,
weights=weight * dispersion[inside])
spread = _np.bincount(index[inside], minlength=n,
weights=weight * square[inside] / n_sim) \
- between[coarse] / n_sim
covered = _np.bincount(index[inside], minlength=n,
weights=volume[inside])
complete = covered == cells(coarse)
deepest = _np.zeros(n, dtype=_np.int64)
_np.maximum.at(deepest, index[inside], level[inside] - coarse)
variance.append((within + spread)[complete])
depth.append(deepest[complete])
left_out += int(_np.count_nonzero(~complete))
variance = _np.maximum(_np.concatenate(variance), 0.0)
count = len(variance)
sizes.append({
"level": coarse,
"size": blocks.base_step * (blocks._coarse_size
// ratio ** coarse),
"deviation": _np.sqrt(variance),
"depth": _np.concatenate(depth),
"rms": float(_np.sqrt(variance.mean())) if count else _np.nan,
"count": count,
"share": count * cells(coarse) / total,
"left_out": left_out})
return sizes
[docs]
def split_label(depth: int) -> str:
"""How many times the refinement split inside a block, in words."""
names = {0: "left whole", 1: "split once", 2: "split twice"}
return names.get(int(depth), "split %d times" % int(depth))
[docs]
def support_tick(entry: dict) -> "list[str]":
"""The lines naming one block size on the axis of
`dispersion_by_support`: the size, then how much of the set it is."""
lines = [" × ".join("%g" % side for side in entry["size"]),
"%d block%s, %.0f%%" % (entry["count"],
"" if entry["count"] == 1 else "s",
100 * entry["share"])]
if entry["left_out"]:
lines.append("%d left out" % entry["left_out"])
return lines
[docs]
def jitter(n: int, width: float = 0.5) -> _types.FloatArray:
"""Offsets spreading `n` points across a strip `width` wide.
Taken from the golden-ratio sequence rather than at random, so a figure
comes out the same every time and the points spread evenly whatever
their number.
"""
return ((_np.arange(n) * 0.6180339887498949 % 1.0 - 0.5)
* width).astype(_np.float64)
[docs]
def support_strip(sizes: "list[dict]", most: int = 5000) -> "list[tuple]":
"""The points of `dispersion_by_support` drawn as a strip.
One `(depth, x, y)` per depth of splitting, `x` being the position of
the block's size plus a spread across the strip. At most about `most`
blocks of each size are kept, taken by striding through them.
"""
x, y, depth = [], [], []
for position, entry in enumerate(sizes):
if entry["count"] == 0:
continue
pick = _np.arange(0, entry["count"],
max(1, int(_np.ceil(entry["count"] / most))))
x.append(position + jitter(len(pick)))
y.append(entry["deviation"][pick])
depth.append(entry["depth"][pick])
if not x:
return []
x, y, depth = (_np.concatenate(a) for a in (x, y, depth))
return [(int(d), x[depth == d], y[depth == d]) for d in _np.unique(depth)]
[docs]
def grid_shape(n_panels):
"""Rows and columns for `n_panels`, as square as they go."""
columns = int(_np.ceil(_np.sqrt(n_panels)))
rows = int(_np.ceil(n_panels / columns))
return rows, columns
# --------------------------------------------------------------------------- #
# mesh sets
# --------------------------------------------------------------------------- #
[docs]
def key_label(shells: "_data.MeshSet", key) -> str:
"""A mesh set's key as a figure writes it: a category's name, or the
cut-off as a number."""
return str(key) if shells.kind == "category" else "%g" % key
[docs]
def key_axis(shells: "_data.MeshSet") -> str:
"""What the keys of a mesh set are, for the axis they run along."""
if shells.kind == "category":
return "category"
if "probability_of" in shells.provenance:
return "probability"
return "cut-off" if shells.unit is None else "cut-off (%s)" % shells.unit
[docs]
def volume_dispersion(shells: "_data.MeshSet",
relative: bool = False) -> dict:
"""
What `volume_dispersion` draws: every realization's mesh volume at each
cut-off or category, beside the prediction's.
Read off what the set measured as each realization was made, so no mesh
is loaded.
Parameters
----------
shells
A set built with its realizations.
relative
Whether to divide every volume by the prediction's, so that the
prediction reads one.
Returns
-------
dict
`labels`, one per key; `values`, the realizations' volumes at each
key; `prediction`, the prediction's; `title`; `axis`, what the
volumes are; `keys`, what the keys are.
"""
table = shells.realization_volumes()
summary = shells.volume_dispersion()
values, prediction = [], []
for key in shells:
held = table[key].to_numpy(dtype=float)
held = held[_np.isfinite(held)]
own = float(summary.loc[key, "prediction"])
if relative:
held = held / own if own > 0 else _np.full(held.shape, _np.nan)
own = 1.0 if own > 0 else _np.nan
values.append(held)
prediction.append(own)
return {"labels": [key_label(shells, key) for key in shells],
"values": values, "prediction": _np.asarray(prediction),
"relative": relative,
"title": "%s: volume by realization" % shells.path,
"axis": "volume over the prediction's" if relative
else "volume",
"keys": key_axis(shells)}
[docs]
def connectivity(shells: "_data.MeshSet") -> dict:
"""
What `connectivity` draws: how many pieces each mesh is in, and the
largest one's share of its volume.
Returns
-------
dict
`labels`; `x`, the cut-offs or the categories' positions;
`numeric`, whether `x` is a scale; `largest` and `pieces`, the
prediction's; `band`, the realizations' P10 and P90 of the largest
share, `median` their P50 and `pieces_median` the pieces' P50 --
None without realizations; `title`; `keys`.
"""
frame = shells.connectivity()
numeric = shells.kind != "category"
if numeric:
x = _np.asarray(list(shells), dtype=float)
else:
x = _np.arange(len(shells), dtype=float)
held = "largest_p10" in frame
band = None
if held:
band = (frame["largest_p10"].to_numpy(dtype=float),
frame["largest_p90"].to_numpy(dtype=float))
return {"labels": [key_label(shells, key) for key in shells],
"x": x, "numeric": numeric,
"largest": frame["largest"].to_numpy(dtype=float),
"pieces": frame["pieces"].to_numpy(dtype=float),
"band": band,
"median": frame["largest_p50"].to_numpy(dtype=float)
if held else None,
"pieces_median": frame["pieces_p50"].to_numpy(dtype=float)
if held else None,
"title": "%s: connectivity" % shells.path,
"keys": key_axis(shells)}
[docs]
def mesh_section(shells: "_data.MeshSet", axis, value: float,
variable=None, resolution: "float | None" = None) -> dict:
"""
What `section` draws: where each mesh of a set crosses a plane, and
what the model holds on it.
Parameters
----------
shells
The set.
axis
The coordinate held fixed, by index or by label.
value
Where along it the plane sits.
variable
A continuous variable or component of the set's container, drawn
beneath the lines; nothing is drawn beneath without one.
resolution
The spacing to sample the variable at on the plane. The finest
block by default, or a grid's own spacing.
Returns
-------
dict
`lines`, for every key's label a list of `(n, 2)` polylines in the
plane; `axes`, the labels of the two in-plane coordinates; `image`,
the sampled variable or None; `title`; `keys`.
"""
data = shells.data
labels = [str(label) for label in
(getattr(data, "coordinate_labels", None) or ("X", "Y", "Z"))]
if isinstance(axis, (int, _np.integer)):
index = int(axis)
else:
lowered = [label.lower() for label in labels]
wanted = str(axis).lower()
index = lowered.index(wanted) if wanted in lowered \
else "xyz".index(wanted)
plane = [i for i in range(3) if i != index]
cut = shells.section(index, value)
lines = {key_label(shells, key): [line[:, plane] for line in cut[key]]
for key in shells}
image = None
if variable is not None and data is not None:
image = section_image(data, variable, index, value, resolution)
return {"lines": lines, "axes": (labels[plane[0]], labels[plane[1]]),
"image": image, "keys": key_axis(shells),
"title": "%s: section at %s = %g" % (shells.path, labels[index],
value)}
[docs]
def section_image(container, var, index: int, value: float,
resolution: "float | None" = None) -> "dict | None":
"""
A variable's prediction sampled on a plane across one axis.
A block model answers exactly which block holds each sample; any other
container its nearest location, left empty past half a cell's diagonal.
At most 600 samples a side.
Returns
-------
dict or None
`values` `(n_v, n_u)`, `extent` `(u_min, u_max, v_min, v_max)` and
`label`; None when the variable holds no prediction.
"""
column = getattr(var, "prediction", None)
if column is None or not column._has_content():
return None
values = _np.asarray(column.values, dtype=float).ravel()
low = _np.ravel(container.bounding_box.min)
high = _np.ravel(container.bounding_box.max)
plane = [i for i in range(3) if i != index]
if resolution is None:
step = getattr(container, "base_step", None)
if step is None:
step = getattr(container, "step_size", None)
if step is not None:
resolution = float(_np.min(_np.asarray(step, dtype=float)[plane]))
else:
resolution = float(_np.max(high - low)) / 200.0
counts = [int(min(600, max(2, _np.ceil((high[i] - low[i]) / resolution))))
for i in plane]
u = _np.linspace(low[plane[0]], high[plane[0]], counts[0] + 1)
v = _np.linspace(low[plane[1]], high[plane[1]], counts[1] + 1)
grid_u, grid_v = _np.meshgrid((u[:-1] + u[1:]) / 2, (v[:-1] + v[1:]) / 2)
points = _np.zeros((grid_u.size, 3))
points[:, plane[0]] = grid_u.ravel()
points[:, plane[1]] = grid_v.ravel()
points[:, index] = float(value)
if isinstance(container, _data.BlockSet3D):
found = container.index_data(_data.PointData.from_array(points))
sampled = _np.where(found >= 0, values[_np.maximum(found, 0)],
_np.nan)
else:
coordinates = _np.asarray(container.coordinates, dtype=float)
step = _np.asarray(getattr(container, "step_size", resolution),
dtype=float) * _np.ones(3)
distance, nearest = _spatial.KDTree(coordinates).query(points)
sampled = _np.where(distance <= 0.5 * float(_np.linalg.norm(step)),
values[nearest], _np.nan)
return {"values": sampled.reshape(grid_u.shape),
"extent": (u[0], u[-1], v[0], v[-1]), "label": var.name}