# 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 variable family: `_Variable` and the concrete kinds a container holds
(continuous, vector, compositional, categorical, binary), each declaring
its own columns for the tree machinery in `base` to fold over. Constructed
by a container's `add_*_variable` methods, never directly.
"""
import collections as _col
from typing import Any as _Any
import numpy as _np
import pandas as _pd
import tensorflow as _tf
import sklearn.metrics as _skmetrics
import geoml._types as _types
import geoml.metrics as _gmlmetrics
import geoml.storage as _storage
from geoml.data.base import *
from geoml.data.base import (
_Attribute, _TreeNode, _carry_rows, _code_dtype, _copy_for_subset,
_encode,
_path_key, _subset_simulations)
# `_Variable._Attribute` is the leaf class under another name (bound at
# the end of this module). Annotations use `_Attr` so that a return type
# inside the class is not read as that alias.
_Attr = _Attribute
# What to divide a column by to turn its unit into a fraction of the whole.
# Lives here rather than in `drillhole`, where it started, because a
# composition's parts carry their unit as a fact of the variable and the
# drillhole is only one of the doors data comes in by; `drillhole` re-exports
# it, so `geoml.data.drillhole.UNITS` still resolves.
UNITS = {"fraction": 1.0, "ratio": 1.0, "1": 1.0,
"%": 100.0, "pct": 100.0, "percent": 100.0, "wt%": 100.0,
"ppm": 1e6, "g/t": 1e6, "mg/kg": 1e6,
"ppb": 1e9, "ug/kg": 1e9, "mg/t": 1e9}
def _divisor(unit):
"""What to divide a column by to turn its unit into a fraction.
A number is taken as the divisor itself; a string is looked up in
`UNITS`. `None` means undeclared, which is a fraction: that is what
makes an undeclared composition behave exactly as it did before units
existed.
"""
if unit is None:
return 1.0
if isinstance(unit, (int, float)) and not isinstance(unit, bool):
if unit <= 0:
raise ValueError(f"a unit divisor must be positive, got {unit}")
return float(unit)
key = str(unit).strip().lower()
if key not in UNITS:
raise ValueError(
f"unknown unit {unit!r}; expected one of {sorted(UNITS)}, or a "
f"number to divide the column by")
return UNITS[key]
def _units_per_label(units, labels):
"""`{label: unit}` from a mapping, a sequence, or nothing.
A variable's parts are named, so a mapping is the spelling that cannot
be got wrong; a sequence is accepted in the labels' own order, and a
label the mapping does not name is undeclared rather than an error --
a composition's `rest` is exactly that case.
"""
if units is None:
return {label: None for label in labels}
if isinstance(units, dict):
unknown = [key for key in units if key not in labels]
if unknown:
raise ValueError(
"units given for %s, which %s not among the labels %s"
% (", ".join(repr(u) for u in unknown),
"is" if len(unknown) == 1 else "are",
", ".join(repr(str(lb)) for lb in labels)))
return {label: units.get(label) for label in labels}
units = list(units)
if len(units) != len(labels):
raise ValueError(
"%d unit(s) given for %d label(s); pass a mapping from label to "
"unit where only some are declared"
% (len(units), len(labels)))
return dict(zip(labels, units))
def _store_bytes(store):
"""What `store` takes in memory read whole: float32 is read as float64."""
itemsize = _np.dtype(store.dtype).itemsize
if _np.dtype(store.dtype) == _np.float32:
itemsize = 8
return int(_np.prod(store.shape)) * itemsize
def _missing_value(attribute):
"""What an empty row of `attribute` holds: the missing code of a coded
column, False in a flag, the empty string in text, NaN in a number."""
kind = _np.dtype(attribute.values.dtype).kind
if attribute.labels is not None or kind in "iu":
return -1
if kind == "b":
return False
if kind == "O":
return ""
return _np.nan
def _measured_in_order(attribute_class, coordinates, values, labels):
"""A categorical measurement column coded against the variable's
`labels`, in their order, with anything else it holds appended after
them in the order pandas sorts it into -- so a class drawn by position
takes one colour in `predicted`, `measurements_a` and `measurements_b`
alike, and nothing measured is lost."""
# a label given twice is one category, as it always decoded
labels = list(dict.fromkeys(labels))
extra = [c for c in _pd.Categorical(_np.asarray(values)).categories
if c not in set(labels)]
full = labels + extra
codes = _pd.Categorical(_np.asarray(values), categories=full).codes
return attribute_class.encoded(
coordinates, codes.astype(_code_dtype(len(full))), labels=full)
def _refuse_past_threshold(nbytes, name):
"""Refuses to hold `nbytes` of simulations whole past the size the store
class spills to disk at -- what `get_simulations` on a block model would
ask for, hundreds of gigabytes that take the session with them. The
message names the two reads that work at any size."""
if nbytes > _storage.DEFAULT_THRESHOLD:
raise MemoryError(
"%s holds %.0f MB of simulations, past the %.0f MB the package "
"keeps in RAM; read one realization with `simulation(i)` or the "
"store in bands with `simulations.row_bands()`"
% (name, nbytes / 1024 ** 2,
_storage.DEFAULT_THRESHOLD / 1024 ** 2))
def _blank_like(store, coordinates):
"""An empty realization store of `store`'s kind on `coordinates`:
missing (NaN, or -1 for integer labels) until written."""
dtype = _np.dtype(store.dtype)
fill = _np.nan if _np.issubdtype(dtype, _np.floating) else -1
return _storage.ArrayStore.allocate(
(coordinates.n_data, store.shape[1]), dtype=dtype, fill_value=fill,
owner=coordinates)
def _write_mixture(variable, idx, kwargs):
"""What a mixture of likelihoods adds to a prediction, on the node it
belongs to: each realization's population, and the expected share of
each population, filed as the responsibilities -- the share before a
measurement, which `VGPNetwork.responsibilities` updates where one was
read."""
if "population" in kwargs.keys():
labels = _np.asarray(kwargs["population"])
if variable.population is None:
variable.population = _storage.ArrayStore.allocate(
(variable.coordinates.n_data, labels.shape[1]),
dtype=_np.int8, fill_value=-1, owner=variable.coordinates)
variable.population[idx, :] = labels
if "shares" in kwargs.keys():
shares = _np.asarray(kwargs["shares"])
for k in range(shares.shape[1]):
if k not in variable.responsibilities:
variable.responsibilities[k] = variable._Attribute(
variable.coordinates)
variable.responsibilities[k].values[idx] = shares[:, k]
class _Variable(_TreeNode):
"""Representation of a dependent random variable."""
# Subclasses that can be simulated replace this with an (n_data, n_sim)
# store; declared here so the accessors below work on every variable.
simulations: "_storage.ArrayStore | None" = None
name: str
metrics: "_pd.DataFrame | None"
# Bound at the end of the module (`_Variable._Attribute = _Attribute`),
# which is what keeps the historical `self._Attribute(...)` call sites
# working; declared here so that reading one is not an archaeology.
_Attribute: "type[_Attr]"
# not every kind has components; declaring is not assigning, so the
# `getattr` tests that ask still answer as before
components: "dict[str, _Any]"
def __init__(self, name: str, coordinates):
self.name = name
self.coordinates = coordinates
self._length = 1
self.metrics = None
# What this variable's ``labels`` are called when it prints itself.
_LABEL_KIND = "labels"
def __repr__(self):
s = "%s(%r, n_data=%s" % (
self.__class__.__name__, self.name, self.coordinates.n_data)
if self.n_sim > 0:
s += ", n_sim=%d" % self.n_sim
return s + ")"
def __str__(self):
"""
The variable, what it is made of, and what can be read off it.
``repr`` says which variable this is; this says what is on it, so that
the columns do not have to be looked up in the source -- a variable
carries different ones depending on what it models.
Only the names are listed, not whether anything has been written into
them: knowing that means reading each column, and a column on a large
model lives on disk, which is too much to do for a print.
"""
lines = [repr(self)]
labels = getattr(self, "labels", None)
if labels is not None:
lines.append(" %s: %s" % (
self._LABEL_KIND, ", ".join(str(label) for label in labels)))
columns = [name for name, value in vars(self).items()
if isinstance(value, _Attribute)]
if len(columns) > 0:
lines.append(" columns: " + ", ".join(columns))
for name in ("quantiles", "probabilities"):
keys = getattr(self, name, None)
if keys:
lines.append(" %s: %s"
% (name, ", ".join(str(key) for key in keys)))
if self.n_sim > 0:
lines.append(" simulations: %d" % self.n_sim)
return "\n".join(lines)
@property
def length(self):
return self._length
@classmethod
def from_variable(cls, coordinates, variable):
raise NotImplementedError
def split_shares(self):
"""Which decisions this variable carries, and how often each cuts a
block in two.
`{name: (n_data,) array}`, the name saying which decision *within*
this variable -- a cut-off for a grade, a boundary for a category --
and empty where the variable declares none.
Asked of the variable rather than worked out from outside, which used
to be a necessity (`divided` was an `_Attribute` on a category and an
`OrderedDict` on a grade) and is now only about the names: the
storage is one shape, but a grade's decision renders `@ 1.5` because
someone declared that number, while a category's stays bare -- its
zero is an artefact of the log-odds, and `granite @ 0` would read as
noise. `_share_label` is that one difference.
"""
shares = {}
for cutoff, attribute in (getattr(self, "divided", None) or {}).items():
shares[self._share_label(cutoff)] = attribute.values.to_numpy()
for label, component in (
getattr(self, "components", None) or {}).items():
for key, values in component.split_shares().items():
shares[("%s %s" % (label, key)).strip()] = values
return shares
def _share_label(self, cutoff):
return "@ %g" % cutoff
def _fill_pyvista(self, write, simulations=False, prefix=None,
include="**"):
"""`write(label, attribute)` for every filled column, named by path.
The same enumeration the data frame reads (`_export_leaves`), handed
to a writer instead of a frame -- in place of the twenty-one
per-class methods that each named their columns again, no two
agreeing on the spelling. A label is `render(path, "pretty")`:
`assay - Zn - noise_variance`, the same segments every other export
uses.
"""
root = None if prefix is None \
else VariablePath(prefix) / self._node_name
for path, attribute in self._export_leaves(include, simulations, root):
write(render(path, "pretty"), attribute)
def fill_pyvista_cube(self, cube, prefix=None, sigma=None,
simulations=False, include="**"):
self._fill_pyvista(
lambda label, attribute: attribute.fill_pyvista_cube(
cube, label, sigma=sigma),
simulations, prefix, include)
def fill_pyvista_points(self, points, prefix=None, simulations=False,
include="**"):
self._fill_pyvista(
lambda label, attribute: attribute.fill_pyvista_points(
points, label),
simulations, prefix, include)
def fill_pyvista_blocks(self, cube, prefix=None, sigma=None,
simulations=False, include="**"):
self._fill_pyvista(
lambda label, attribute: attribute.fill_pyvista_blocks(
cube, label, sigma=sigma),
simulations, prefix, include)
def fill_pyvista_cells(self, mesh, prefix=None, simulations=False,
include="**"):
self._fill_pyvista(
lambda label, attribute: attribute.fill_pyvista_cells(mesh, label),
simulations, prefix, include)
def carry_to(self, coordinates, keep, n_new):
"""This variable on a longer set of locations, keeping what still fits.
The locations of `coordinates` are the `keep` ones of this variable, in
their old order, followed by `n_new` that did not exist before. The
first take their values across; the rest come back missing, which is
what marks them as still to be predicted.
Used when a block set is refined: a block that was not split is the
same block on the same support, and its value is still the right
answer for it, so re-predicting it would spend the work to arrive at
the number already there.
"""
new = self.__class__.from_variable(coordinates, self)
self._copy_attrs_into(new)
self._carry_into(new, _np.asarray(keep, dtype=bool))
return new
def _carry_into(self, new, keep):
"""Fill an already-built variable with the `keep` rows of this one.
Apart from `carry_to` so that a variable holding components fills the
ones its own `from_variable` just created, rather than building a
second set that would go nowhere.
"""
n_kept = int(_np.count_nonzero(keep))
for role in self._ZARR_ATTRS:
old = getattr(self, role, None)
fresh = getattr(new, role, None)
if old is None or fresh is None or not old._has_content():
continue
fresh.labels = old.labels
fresh.values[:n_kept] = _np.asarray(old.values)[keep]
for name, store in self.realization_stores():
fresh = _blank_like(store, new.coordinates)
setattr(new, name, fresh)
_carry_rows(store, fresh, keep)
for family in self._DICT_FAMILIES:
target = getattr(new, family)
for key, attr in (getattr(self, family, None) or {}).items():
fresh = self._Attribute(new.coordinates)
fresh.values[:n_kept] = _np.asarray(attr.values)[keep]
target[key] = fresh
for name, component in (getattr(self, "components", None) or {}).items():
component._carry_into(new.components[name], keep)
# The columns that are a mean over a block's sub-blocks, and the families
# of them: a coarser block's is then the volume-weighted mean of its
# parts', exactly. What `_coarsen_into` averages; a column left out of
# these is left missing wherever a block is gathered, never wrong.
_BLOCK_MEANS: "tuple[str, ...]" = ()
_BLOCK_MEAN_FAMILIES: "tuple[str, ...]" = ()
def _coarsen_into(self, new, grouping, valid=None):
"""Fill an already-built variable on coarser blocks from this one.
`grouping` (built by `BlockSet3D.as_blocks3d`) says which coarse
block each of this variable's blocks falls in and what share of it
each holds. A coarse block that is one of them whole keeps every
column as it stands. One gathered from parts takes the
volume-weighted mean of the columns `_BLOCK_MEANS` and
`_BLOCK_MEAN_FAMILIES` declare, and of every realization index by
index; it leaves everything else missing, and a subclass recomputes
what those averages settle. `valid` marks the parts that were
predicted, for a kind whose columns cannot say so themselves.
Apart from `as_blocks3d` for the same reason `_carry_into` is apart
from `carry_to`: a variable holding components fills the ones its
own `from_variable` just built.
"""
for role in self._ZARR_ATTRS:
old = getattr(self, role, None)
if old is None or getattr(new, role, None) is None:
continue
if role in self._BLOCK_MEANS:
values = grouping.mean(old.values.to_numpy(), valid)
else:
values = grouping.kept(old.values.to_numpy(),
_missing_value(old))
fresh = self._Attribute(new.coordinates, values,
dtype=values.dtype)
fresh.labels = None if old.labels is None else list(old.labels)
setattr(new, role, fresh)
self._coarsen_realizations(new, grouping, valid)
for family in self._DICT_FAMILIES:
target = getattr(new, family)
for key, old in (getattr(self, family, None) or {}).items():
if family in self._BLOCK_MEAN_FAMILIES:
values = grouping.mean(old.values.to_numpy(), valid)
else:
values = grouping.kept(old.values.to_numpy(), _np.nan)
target[key] = self._Attribute(new.coordinates, values)
for name, component in (getattr(self, "components", None) or {}).items():
component._coarsen_into(new.components[name], grouping, valid)
def _coarsen_realizations(self, new, grouping, valid):
"""The realizations on the coarser blocks, averaged index by index."""
if self._ZARR_HAS_SIMS and self.simulations is not None:
new.allocate_simulations(self.simulations.shape[1])
grouping.realizations(self.simulations, new._sim_store(), valid)
def _subset_into(self, new, item):
"""Fill a copy with the `item` rows of this node, column by column.
Driven by the declarations rather than written out per class: a
column named in `_ZARR_ATTRS` and forgotten in a hand-written subset
survives at its old length, and only whatever reads it afterwards
finds out. That is how a `BinaryVariable` shipped a subset whose
`probability` was never cut -- the old override wrote it into a dead
`average` attribute instead.
"""
for role in self._ZARR_ATTRS:
column = getattr(self, role, None)
if column is not None:
setattr(new, role, column[item])
for family in self._DICT_FAMILIES:
target = getattr(new, family)
for key, attribute in (getattr(self, family, None) or {}).items():
target[key] = attribute[item]
for name, store in self.realization_stores():
setattr(new, name, _subset_simulations(store, item))
theirs = new.child_nodes()
for label, child in self.child_nodes().items():
child._subset_into(theirs[label], item)
def __getitem__(self, item):
new_obj = _copy_for_subset(self)
self._subset_into(new_obj, item)
return new_obj
def set_responsibilities(self, values: _types.ArrayLike) -> "_Variable":
"""File which noise component each measurement is likely to be from.
One column per component of a `Mixture` likelihood, keyed by its
position, so `assay/responsibilities/1` is how often the second
(wider) component would explain the row. One answer per location:
the mixture is over the row, not over the columns of a vector
variable.
Written by `models.VGPNetwork.responsibilities`, never by a
prediction -- it takes measurements, which a grid does not have.
Parameters
----------
values
Of shape `(n_data, n_components)`, rows summing to one, or
missing where the location carries no measurement.
"""
values = _np.asarray(values, dtype=float)
if values.ndim != 2:
raise ValueError("responsibilities come one row per location and "
"one column per component")
self.responsibilities = _col.OrderedDict()
for k in range(values.shape[1]):
attribute = self._Attribute(self.coordinates)
attribute.values[:] = values[:, k]
self.responsibilities[k] = attribute
return self
def _sim_store(self) -> "_storage.ArrayStore":
"""The simulations store, or a legible error where there is none.
Every writer goes through this rather than indexing the attribute:
a variable that was never allocated has no store, and saying so is
better than a `NoneType` error from inside a batch loop.
"""
if self.simulations is None:
raise NoDataError(
"variable %r has no simulations; call "
"`allocate_simulations` first" % str(self.name))
return self.simulations
def from_model_units(self, values):
"""Values as the model produced them, in this variable's own units.
The identity for everything the model reads directly -- a variable's
unit is a label there, and the numbers never left the units they
arrived in. A composition overrides it: its parts reach the model as
fractions of the whole, so what comes back has to be put into the
unit each part was assayed in. `values` carries the component axis
second, as everything the model hands back does.
"""
return values
def get_measurements(self):
raise NotImplementedError
def get_simulations(self):
raise NotImplementedError
@property
def n_sim(self):
"""Number of simulations available, or 0 if none were drawn."""
return 0 if self.simulations is None else self.simulations.shape[1]
def simulation(self, index: int) -> _Attr:
"""
A single simulation, in the form of an `_Attribute`.
Simulations are kept in one `(n_data, n_sim)` array, so a single one is
not an attribute in its own right. This wraps the corresponding column,
giving it the usual helpers (`as_image()`, `as_cube()`, `smooth()`,
`get_contour()`, `draw_*()`, ...).
Parameters
----------
index
Position of the simulation, from 0 to `n_sim - 1`.
Returns
-------
attribute : _Attribute
A copy of the simulation. Modifying it (with `smooth()`, for
instance) does not affect the stored simulations; assign it back
with `variable.simulations[:, index] = attribute.values` if that is
the intention.
Notes
-----
Only this column is ever held in memory, however large the store --
which is what makes processing realizations sequentially viable on a
block model. The store is chunked by location, though, so extracting
a column still visits every chunk on disk: walking all realizations
this way costs one full pass over the store *per realization*. A
computation that decomposes over locations is cheaper the other way
around -- read row bands (`ArrayStore.row_bands`) and take every
realization of each band at once.
"""
if self.simulations is None:
raise NoDataError(
f'No simulations available for variable {self.name}.')
return self._Attribute(self.coordinates, self.simulations[:, index])
def get_predictions(self):
raise NotImplementedError
def prediction_input(self):
return {}
def training_input(self, idx=None):
return {}
def copy_to(self, coordinates):
new = self.__class__.from_variable(coordinates, self)
self._copy_attrs_into(new)
coordinates.variables[self.name] = new
def _adopt_cutoffs(self, source):
"""Take the cut-offs of `source` wherever this variable's differ.
A model computes its shares against the cut-offs of the variable it
was trained on, and `update` files them by position under this
variable's own. Where the two name different ground the answer would
land under the wrong value, so this variable takes the model's,
converted to its own unit, and `set_cutoffs` drops the shares of the
ones it gives up.
Returns
-------
list of tuple
`(path, old, new)` for every node whose cut-offs changed.
"""
theirs = dict(source.walk())
changed = []
for path, node in self.walk():
other = theirs.get(path)
# the cut-offs are a grade's, a part's or a component's; a
# category's is always zero
if not isinstance(node, ContinuousVariable) \
or not isinstance(other, ContinuousVariable):
continue
if node._model_cutoffs() == other._model_cutoffs():
continue
# verbatim while the unit is the training node's, so that the
# number comes back as declared rather than rounded through the
# model's unit and back
new = other.cutoffs if node.unit == other.unit \
else node._from_model_cutoffs(other._model_cutoffs())
changed.append((path, node.cutoffs, new))
node.set_cutoffs(new)
return changed
def update(self, idx, **kwargs):
raise NotImplementedError
def allocate_simulations(self, n_sim):
raise NotImplementedError
def compute_metrics(self, **kwargs):
raise NotImplementedError
# ------------------------------------------------------------------ #
# what a prediction leaves behind
# ------------------------------------------------------------------ #
# The column `predict` always fills, whose missing values name the
# locations it never reached. Declared per class because not every kind
# of variable has a `prediction`: a rock type has an entropy, a vector
# variable an uncertainty, and reading the wrong column would call a
# predicted location unpredicted for ever.
_PREDICTED_MARKER: "str | None" = None
def unpredicted(self):
"""One boolean per location: True where no prediction reached it.
What a cancelled or partial `predict` left to do. Reading the
missing values rather than remembering a call means the answer
stays true however the container was arrived at -- reopened from a
store, subsetted, or split.
"""
if self._PREDICTED_MARKER is None:
raise ValueError(
"%s says nothing about where it was predicted; it declares "
"no marker column" % type(self).__name__)
column = getattr(self, self._PREDICTED_MARKER, None)
if column is None:
return _np.ones(self.coordinates.n_data, dtype=bool)
return _np.isnan(_np.asarray(column.values, dtype=float))
# ------------------------------------------------------------------ #
# Zarr persistence (see _SpatialData.to_zarr / _SpatialData.open)
# ------------------------------------------------------------------ #
# Scalar ``_Attribute`` roles to persist; overridden per subclass.
_ZARR_ATTRS = ()
_ZARR_HAS_SIMS = False # has a (n_data, n_sim) simulations store
# What each column holds, for a reader that colours or scales it without
# knowing the class -- `geoml.catalogue` publishes it. `(role, scale)`
# for every name in `_ZARR_ATTRS` and `_DICT_FAMILIES`, and for
# `simulations` where there are any; `test_catalogue.py` fails on a
# column left out. The role is `measurement` (the data), `value` (an
# estimate), `uncertainty` or `weight`; the scale is `unit` (the
# variable's own), `unit_squared`, `latent` (the model's standardized
# space), `latent_squared`, `unit_interval` (0 to 1), `log_odds`,
# `classes` (codes into the labels), `flag`, `ordinal` or
# `dimensionless`.
_COLUMNS: "dict[str, tuple[str, str]]" = {}
# what a family's columns are keyed by: `cutoff`, `probability` or
# `component`
_FAMILY_KEYS: "dict[str, str]" = {}
# the class of the parts a vector or categorical variable holds, whose
# columns are theirs rather than the variable's
_PART: "type | None" = None
def _save_attr(self, group, prefix, role):
"""Write one ``_Attribute``'s store into ``group``; None-valued -> skip.
String/object attributes are stored as fixed-length unicode; everything
else is streamed as a numeric/bool Zarr array.
"""
attr = getattr(self, role, None)
if attr is None:
return None
key = prefix + "/" + role
store = attr.values
if _np.dtype(store.dtype) == object:
unicode = _np.asarray(store).astype(str)
if unicode.dtype.itemsize == 0:
unicode = unicode.astype("<U1")
target = group.create_array(
name=key, shape=unicode.shape, chunks=unicode.shape,
dtype=unicode.dtype)
target[:] = unicode
return {"key": key, "encoding": "str"}
store.write_into(group, key)
info = {"key": key, "encoding": "array"}
if attr.labels is not None:
# a coded attribute's categories are its own — a measurement
# column's are not the variable's. Stringified, as the variable's
# own labels already are: this goes into JSON, and categories can
# come out of pandas as NumPy integers.
info["labels"] = [str(label) for label in attr.labels]
return info
def _load_attr(self, group, info):
z = group[info["key"]]
if info["encoding"] == "str":
return _storage.ArrayStore.from_numpy(
_np.asarray(z[:]).astype(object))
return _storage.ArrayStore.wrap_zarr(z)
def _zarr_save(self, group, prefix):
# str(name): a component is named after its category, which pandas can
# hand over as a NumPy integer, and this is JSON
meta = {"class": type(self).__name__, "name": str(self.name),
"attrs": {}}
own_labels = getattr(self, "labels", None)
if own_labels is not None:
meta["labels"] = [str(x) for x in own_labels]
for role in self._ZARR_ATTRS:
info = self._save_attr(group, prefix, role)
if info is not None:
meta["attrs"][role] = info
for name, store in self.realization_stores():
key = prefix + "/" + name
store.write_into(group, key)
if name == "simulations":
# the key every store written before the others existed has
meta["simulations"] = key
else:
meta.setdefault("stores", {})[name] = key
# the dict families and the node's own facts, off the declarations,
# with the Zarr key spelling the same string `get` takes:
# `assay/Zn/quantiles/0.5`
for family in self._DICT_FAMILIES:
entries = []
for at, attr in (getattr(self, family, None) or {}).items():
key = prefix + "/" + family + "/" + _path_key(at)
attr.values.write_into(group, key)
entries.append({"key": key, "at": float(at)})
if entries:
meta.setdefault("dicts", {})[family] = entries
facts = {name: value for name, value in self.node_attrs().items()
if value is not None}
if facts:
meta["node_attrs"] = facts
if getattr(self, "components", None):
meta["components"] = {}
for cname, comp in self.components.items():
# str: this is a JSON key, and a category can be a NumPy
# integer. The labels above are stringified for the same
# reason, so the rebuilt variable's components match.
meta["components"][str(cname)] = comp._zarr_save(
group, prefix + "/" + str(cname))
return meta
def _zarr_load(self, group, prefix, meta):
for role, info in meta.get("attrs", {}).items():
attribute = getattr(self, role)
store = self._load_attr(group, info)
if attribute.labels is not None and info["encoding"] == "str":
# written before the categorical attributes held codes; a
# value outside the variable's labels is lost here, there
# being nothing else to read the categories from
store = _storage.ArrayStore.from_numpy(
_encode(_np.asarray(store), attribute.labels))
if "labels" in info:
attribute.labels = list(info["labels"])
attribute.values = store
if meta.get("simulations") is not None:
self.simulations = _storage.ArrayStore.wrap_zarr(
group[meta["simulations"]])
for name, key in meta.get("stores", {}).items():
setattr(self, name, _storage.ArrayStore.wrap_zarr(group[key]))
for family, entries in meta.get("dicts", {}).items():
target = getattr(self, family)
for info in entries:
attr = self._Attribute(self.coordinates)
attr.values = _storage.ArrayStore.wrap_zarr(group[info["key"]])
target[info["at"]] = attr
for name, value in meta.get("node_attrs", {}).items():
setattr(self, name, value)
for cname, cmeta in meta.get("components", {}).items():
self.components[cname]._zarr_load(
group, prefix + "/" + str(cname), cmeta)
# `_Attribute` used to be defined inside `_Variable`; the alias keeps its
# `self._Attribute(...)` call sites working.
_Variable._Attribute = _Attribute
[docs]
class ContinuousVariable(_Variable):
"""
Representation of a continuous random variable.
Attributes
----------
measurements : _Attribute
The raw measurements.
latent_mean : _Attribute
The mean of the latent Gaussian representation.
latent_variance : _Attribute
The variance of the latent Gaussian representation.
dispersion : _Attribute
How much the locations inside each block differ among themselves --
the variance over a block's sub-blocks, averaged over the
realizations, in the variable's own units rather than the latent
ones. A different question from `latent_variance`, which is how sure
the model is *of* the block: a well-known block can still be
heterogeneous, and that is what decides whether cutting it finer
would tell anyone anything. Filled only where the container
discretizes; elsewhere a location has no interior and this stays
missing rather than zero.
noise_variance : _Attribute
How far a fresh *measurement* here would fall from the value above --
the likelihood noise carried into the variable's own units, averaged
over the realizations. The third of three variances and the third
question: `latent_variance` is how sure the model is of the value,
`dispersion` is how much the ground varies inside a block, and this is
how much a sample of it would scatter. A prediction reports the ground,
with the noise integrated out, so this is what has to be added back to
compare against an assay. Missing where the prediction was made with
`include_noise=False`, there being no integration to read it from.
simulations : ArrayStore
Draws from the variable's posterior distribution, in a single
`(n_data, n_sim)` array. Use `simulation()` to get one of them as an
`_Attribute`.
quantiles : dict
The variables quantiles, indexed by the corresponding percentile.
probabilities : dict
Cumulative distribution probabilities, indexed by the corresponding
quantile.
responsibilities : dict
Under a `Mixture` likelihood, how likely each measurement is to have
come from each of its noise components, indexed by the component's
position. Empty otherwise; written by `set_responsibilities`.
unit : str, float or None
What the values are measured in -- `"%"`, `"ppm"`, `"g/t"`, or a
number. On a variable a model reads directly this is a **label**:
the values go to the likelihood as they stand, and the unit travels
so that a figure can say what an axis is in and an export can record
it. On a part of a `CompositionalVariable` it is also the divisor
that turns the value into a fraction of the whole, since parts in
different units cannot be added up. `None` is undeclared.
"""
_ZARR_ATTRS = ("measurements", "latent_mean", "latent_variance",
"prediction", "dispersion", "noise_variance")
_PREDICTED_MARKER = "prediction"
_ZARR_HAS_SIMS = True
# a mixture of likelihoods' population per realization, beside them
_REALIZATION_STORES = ("simulations", "population")
_DICT_FAMILIES = ("quantiles", "probabilities", "proportions", "divided",
"responsibilities", "population_prediction")
_COLUMNS = {
"measurements": ("measurement", "unit"),
"latent_mean": ("value", "latent"),
"latent_variance": ("uncertainty", "latent_squared"),
"prediction": ("value", "unit"),
"dispersion": ("uncertainty", "unit_squared"),
"noise_variance": ("uncertainty", "unit_squared"),
"quantiles": ("value", "unit"),
# the share of the realizations at or below a cut-off
"probabilities": ("value", "unit_interval"),
# the share of a block at or below a cut-off
"proportions": ("value", "unit_interval"),
# whether the prediction's sub-blocks straddle a cut-off
"divided": ("value", "flag"),
"responsibilities": ("value", "unit_interval"),
"simulations": ("value", "unit"),
# which population of a mixture of likelihoods each realization
# takes here, and each population's own prediction
"population": ("value", "classes"),
"population_prediction": ("value", "unit"),
}
_FAMILY_KEYS = {"quantiles": "probability", "probabilities": "cutoff",
"proportions": "cutoff", "divided": "cutoff",
"responsibilities": "component",
"population_prediction": "population"}
_NODE_ATTRS = ("cutoffs", "unit")
# `dispersion` and the quantile families are read again off the
# realizations instead (`_coarsen_realizations`, `_coarsen_into`)
_BLOCK_MEANS = ("latent_mean", "latent_variance", "prediction",
"noise_variance")
_BLOCK_MEAN_FAMILIES = ("proportions",)
measurements: _Attribute
latent_mean: _Attribute
latent_variance: _Attribute
prediction: _Attribute
dispersion: _Attribute
noise_variance: _Attribute
simulations: _storage.ArrayStore | None
cutoffs: list[float] | None
unit: "str | float | None"
quantiles: dict[float, _Attribute]
probabilities: dict[float, _Attribute]
proportions: dict[float, _Attribute]
divided: dict[float, _Attribute]
responsibilities: dict[int, _Attribute]
def __init__(self, name, coordinates, measurements=None, unit=None):
super().__init__(name, coordinates)
if measurements is None:
self.measurements = self._Attribute(coordinates)
else:
self.measurements = self._Attribute(coordinates, measurements)
# What the values are measured in. A fact of the variable, so it
# rides `tree()`, Zarr and every rebuild off `_NODE_ATTRS` rather
# than being carried by hand in each of them.
self.unit = None
if unit is not None:
self.set_unit(unit)
self.latent_mean = self._Attribute(coordinates)
self.latent_variance = self._Attribute(coordinates)
self.prediction = self._Attribute(coordinates)
# How much the locations *inside* each block differ among themselves,
# which is a different question from how sure the model is of the block
# (`latent_variance`). Filled only where the container discretizes;
# everywhere else a location has no interior and this stays zero.
self.dispersion = self._Attribute(coordinates)
# What a sample taken here would read, as against what the ground
# holds: the likelihood noise in this variable's units. A prediction
# integrates that noise out, so this is the piece to add back before
# comparing with an assay.
self.noise_variance = self._Attribute(coordinates)
# A single (n_data, n_sim) store (NumPy or Zarr by size); None until
# ``allocate_simulations`` is called.
self.simulations = None
self.quantiles = _col.OrderedDict()
self.probabilities = _col.OrderedDict()
# The grades a decision turns on -- a mining cut-off, a contaminant
# limit. Declared on the data and carried to whatever is predicted
# from it; `None` means this variable takes no part in any decision,
# which is the answer for the rest component of a composition.
self.cutoffs = None
# For each of them, how much of each block sits at or below it, over
# the sub-blocks and the realizations both -- the complement of the
# recoverable share, and what a partial-block report wants.
self.proportions = _col.OrderedDict()
# And for each, how often the cut-off passes *through* the block:
# the share of realizations whose sub-blocks fall on both sides. A
# different question, and the one that says whether cutting the block
# finer would settle anything. See `likelihood._divided`.
self.divided = _col.OrderedDict()
# Which noise component each measurement came from, under a mixture
# likelihood; empty under every other one. See `set_responsibilities`.
self.responsibilities = _col.OrderedDict()
# Under a mixture of likelihoods: each realization's population at
# each location, `(n_data, n_sim)` small integers, and each
# population's own prediction. Absent under every other likelihood.
self.population = None
self.population_prediction = _col.OrderedDict()
[docs]
def set_unit(self, unit: "_types.Unit | None") -> "ContinuousVariable":
"""What the values are measured in.
A label here: the values reach the likelihood as they stand, and the
unit travels with the variable so that a figure can say what an axis
is in and an export can record it. Anything is accepted, `UNITS`
holding only the ones that can also be *divided* by -- which is what
a part of a composition needs, and what `_Component.set_unit`
insists on.
"""
self.unit = unit
return self
[docs]
def divisor(self) -> float:
"""What to divide this variable's values by to make them fractions."""
return _divisor(self.unit)
[docs]
def set_cutoffs(self, cutoffs: _types.Cutoffs) -> "ContinuousVariable":
"""The grades this variable is judged against.
They travel with the variable, so a model trained on data that
declares them hands them to every block model predicted from it, and
`refine` knows what the blocks have to be resolved against without
being told a second time. The `proportions` and `divided` columns of
a cut-off no longer declared are dropped.
"""
self.cutoffs = None if cutoffs is None else \
[float(c) for c in _np.atleast_1d(cutoffs)]
# a share filed under a cut-off nobody declares any more would keep
# voting on where `refine` splits
kept = set(self.cutoffs or [])
for family in (self.proportions, self.divided):
for cutoff in [c for c in family if c not in kept]:
del family[cutoff]
return self
def _model_cutoffs(self):
"""The cut-offs as the model reads them.
The same numbers here -- a variable a model reads directly is in
whatever units it was measured in, and its unit is a label. A part
of a composition overrides this: its cut-off is declared in its own
unit and the model works in fractions.
"""
return list(self.cutoffs or [])
def _from_model_cutoffs(self, cutoffs):
"""`_model_cutoffs` undone: cut-offs as the model reads them, back
in this variable's own unit."""
return list(cutoffs)
[docs]
def get_measurements(self):
values = self.measurements.values.copy()[:, None]
has_value = (~ _np.isnan(values)) * 1.0
values[_np.isnan(values)] = 0
return values, has_value
[docs]
def get_simulations(self):
# The one deliberate materializer of the (n_data, n_sim) store, for
# data small enough to hold whole; past the size the store spills to
# disk at it refuses, naming `simulation(i)` and the row bands
if self.simulations is not None:
_refuse_past_threshold(_store_bytes(self.simulations), self.name)
return _np.asarray(self.simulations)
[docs]
def get_predictions(self):
return self.prediction.values.to_numpy()
[docs]
def reset_quantiles(
self, probabilities: _types.ArrayLike | None = None) -> None:
"""
Resets the variable's quantiles.
Parameters
----------
probabilities
Probabilities between 0 and 1, exclusive, at which to take
the quantiles.
"""
if self.simulations is None:
raise NoDataError(f'No simulations available for variable {self.name}.')
self.quantiles = _col.OrderedDict()
if probabilities is not None:
probabilities = _np.atleast_1d(
_np.asarray(probabilities, dtype=float))
# All quantiles are computed lazily in a single chunk-by-chunk
# pass over the simulations; the (n_data, n_sim) array is never
# fully materialized.
columns = self.simulations.row_quantiles(probabilities)
targets = []
for p in probabilities:
attr = self._Attribute(self.coordinates)
self.quantiles[p] = attr
targets.append(attr.values)
_storage.store_columns(columns, targets)
[docs]
def reset_probabilities(
self, quantiles: _types.ArrayLike | None = None) -> None:
"""
Resets the variable's probabilities.
Parameters
----------
quantiles
Values in the variable's own units, at which to take the
cumulative probabilities.
"""
if self.simulations is None:
raise NoDataError(f'No simulations available for variable {self.name}.')
self.probabilities = _col.OrderedDict()
if quantiles is not None:
quantiles = _np.atleast_1d(
_np.asarray(quantiles, dtype=float))
# Empirical CDF, the inverse of reset_quantiles: for each cutoff,
# the fraction of simulations at or below it, in (0, 1). Computed
# lazily in a single chunk-by-chunk pass. (The previous
# implementation misused np.percentile, treating the cutoff as a
# percent rank.)
columns = self.simulations.row_cdf(quantiles)
targets = []
for q in quantiles:
attr = self._Attribute(self.coordinates)
self.probabilities[q] = attr
targets.append(attr.values)
_storage.store_columns(columns, targets)
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
# the facts the variable carries -- its cut-offs -- follow separately,
# by `_copy_attrs_into` in `copy_to`/`carry_to`, off the `_NODE_ATTRS`
# declaration rather than named here again
return cls(variable.name, coordinates)
[docs]
def update(self, idx, **kwargs):
# The likelihood speaks in `(rows, components, ...)` whatever the
# number of components; a scalar variable takes its single column.
# Flat arrays still arrive from the legacy closed-form model and from
# a vector variable distributing columns to its parts.
def column(key):
arr = kwargs[key].numpy()
return arr[:, 0] if arr.ndim > 1 else arr
self.prediction.values[idx] = column("average_sim")
if "mean" in kwargs.keys():
self.latent_mean.values[idx] = column("mean")
self.latent_variance.values[idx] = column("variance")
if "dispersion" in kwargs.keys():
self.dispersion.values[idx] = column("dispersion")
if "noise_variance" in kwargs.keys():
self.noise_variance.values[idx] = column("noise_variance")
for key, target in (("proportions", self.proportions),
("divided", self.divided)):
if key not in kwargs.keys():
continue
values = kwargs[key].numpy()
if values.ndim == 3:
values = values[:, 0, :]
cutoffs = self.cutoffs or []
if values.shape[1] == 0:
# a component of a vector variable that declared none, whose
# columns its parent has already trimmed away
continue
# the model asked the *training* variable what the cut-offs were,
# and the answer is being filed against this one; `predict` has
# made this one declare the same (`_adopt_cutoffs`)
if len(cutoffs) != values.shape[1]:
raise ValueError(
"%r was predicted against %d cut-off(s) but declares %s; "
"set them on the data the model was trained from, and let "
"`copy_to` carry them"
% (self.name, values.shape[1], self.cutoffs))
for i, cutoff in enumerate(cutoffs):
if cutoff not in target:
target[cutoff] = self._Attribute(self.coordinates)
target[cutoff].values[idx] = values[:, i]
if "simulations" in kwargs.keys():
# Whole (batch, n_sim) block written as one region into the store.
sims = kwargs["simulations"].numpy()
if sims.ndim == 3:
sims = sims[:, 0, :]
self._sim_store()[idx, :] = sims
_write_mixture(self, idx, kwargs)
if "population_prediction" in kwargs.keys():
values = _np.asarray(kwargs["population_prediction"])
if values.ndim == 3:
values = values[:, 0, :]
for k in range(values.shape[1]):
if k not in self.population_prediction:
self.population_prediction[k] = self._Attribute(
self.coordinates)
self.population_prediction[k].values[idx] = values[:, k]
[docs]
def allocate_simulations(self, n_sim):
self.simulations = _storage.ArrayStore.allocate(
(self.coordinates.n_data, n_sim),
dtype=_storage.realization_dtype(), fill_value=_np.nan,
owner=self.coordinates)
def _coarsen_realizations(self, new, grouping, valid):
# The spread between the parts comes out of the same pass as their
# mean: a block's dispersion is its parts' mean dispersion plus how
# far they sit from one another, realization by realization -- the
# variance of a mixture, and the "larger by exactly what the grouping
# absorbed" of `group`. About the block's prediction, which the
# scalar columns have already averaged.
if self.simulations is None:
return
new.allocate_simulations(self.simulations.shape[1])
centre = new.prediction.values.to_numpy()
between = grouping.realizations(
self.simulations, new._sim_store(), valid,
shift=_np.where(_np.isfinite(centre), centre, 0.0))
within = grouping.mean(self.dispersion.values.to_numpy(), valid)
new.dispersion = self._Attribute(new.coordinates, within + between)
def _coarsen_into(self, new, grouping, valid=None):
super()._coarsen_into(new, grouping, valid)
if self.simulations is None:
return
# read again off the averaged realizations, at the same keys; a block
# kept whole keeps its own, which the same reading would reproduce
for family, reset in (("quantiles", new.reset_quantiles),
("probabilities", new.reset_probabilities)):
old = getattr(self, family)
if not old:
continue
reset(list(old))
fresh = getattr(new, family)
for key in fresh:
kept = grouping.kept(old[key].values.to_numpy(), _np.nan)
fresh[key] = self._Attribute(new.coordinates, _np.where(
grouping.whole, kept, fresh[key].values.to_numpy()))
[docs]
def compute_metrics(self, alpha=0.05):
"""
Scores this variable's prediction against its own measurements.
Parameters
----------
alpha
Significance level for the interval-based scores.
Returns
-------
dict
One entry per score, named.
Notes
-----
The spread-based scores here -- goodness, coverage, CRPS, the
interval score -- are of the **ground**: a container's simulations
have the likelihood's noise integrated out, so they describe a
quantity no sample observes, while the measurements they are
compared against carry it. Those scores therefore read pessimistic
on held-out data, by the share of the variance the model calls
noise, and increasingly so for a model with more capacity, which
calls less of it noise. Measured on Jura, a nominal 90% band read
0.59 here against 0.94 through the measurement distribution.
For calibration on data the model has not seen, use
:func:`geoml.models.cross_validate`, whose scores come from
:meth:`geoml.models.VGPNetwork.predict_measurements`, or the
`accuracy` figure, which asks the model for the same thing. The
location-wise scores (rmse, mae, bias) are unaffected: integrating
the noise out changes the spread, not the value.
See Also
--------
geoml.models.VGPNetwork.predict_measurements : the distribution an
assay is drawn from.
geoml.models.cross_validate : out-of-fold scores, of measurements.
"""
y_true, has_value = self.get_measurements()
if _np.sum(has_value) == 0:
raise ValueError('No measurements available')
y_pred = self.prediction.values.to_numpy()
has_value = has_value[:, 0]
y_true = y_true[has_value == 1]
y_pred = y_pred[has_value == 1]
# Only the measured rows, indexed out of the chunks: materializing the
# store to cut a sliver from it is what kills a session on a large
# model (same reasoning as `_subset_simulations`).
sims = _np.asarray(self._sim_store().as_dask()[has_value == 1])
metrics = {
'Root Mean Square Error (prediction)': _skmetrics.root_mean_squared_error(y_true, y_pred),
'Mean Absolute Error (prediction)': _skmetrics.mean_absolute_error(y_true, y_pred),
'Median Absolute Error (prediction)': _skmetrics.median_absolute_error(y_true, y_pred),
'Bias (prediction)': _gmlmetrics.bias(y_true, y_pred),
'Root Mean Square Error (simulations)': _skmetrics.root_mean_squared_error(
_np.broadcast_to(y_true, sims.shape), sims),
'Mean Absolute Error (simulations)': _skmetrics.mean_absolute_error(
_np.broadcast_to(y_true, sims.shape), sims),
'Median Absolute Error (simulations)': _skmetrics.median_absolute_error(
_np.broadcast_to(y_true, sims.shape), sims),
'CRPS (simulations)': _gmlmetrics.crps(y_true, sims),
# declustered: the pairs come from wherever the drilling went, and
# unweighted they would describe the sampling as much as the
# field. The locations are indexed out of the chunks like the
# simulations above, rather than materialized and then cut down
'Variogram score (simulations)': _gmlmetrics.variogram_score(
y_true, sims,
coordinates=_np.asarray(
self.coordinates.coordinates.as_dask()[has_value == 1],
dtype=float)),
}
bias_2, variance = _gmlmetrics.bias_variance_decomposition(y_true, sims)
metrics['Bias squared (simulations)'] = bias_2
metrics['Variance (simulations)'] = variance
nominal, observed = _gmlmetrics.coverage(y_true, sims)
metrics['Goodness (simulations)'] = _gmlmetrics.goodness(
nominal, observed)
if not isinstance(alpha, (list, tuple)):
alpha = [alpha]
for a in alpha:
metrics[f'Interval score ({a})'] = _gmlmetrics.interval_score(y_true, sims, a)
self.metrics = _pd.Series(metrics, name=self.name)
return self.metrics
[docs]
class DerivedVariable(ContinuousVariable):
"""
A variable computed from others, realization by realization.
The middle ground between metadata (a constant the models never see) and
a modelled variable (measured, likelihooded, written by a model): it
carries a full set of simulations and everything built on them --
quantiles, cut-offs, contours, grade-tonnage -- but every bit of its
uncertainty is inherited from the variables it was derived from. Built by
`derive` on the container, never fed to a model. Applying the function to
each realization and summarizing afterwards is what keeps a nonlinear
function honest: `f(E[grades])` is not `E[f(grades)]`, and the second is
the answer.
The recipe -- the function itself -- lives in the script that ran
`derive`, not here: functions do not survive a Zarr store honestly. A
reloaded container has the values, fully usable; re-deriving is running
the script again. `parents` records which paths it came from.
"""
_NODE_ATTRS = ContinuousVariable._NODE_ATTRS + ("parents",)
def __init__(self, name, coordinates, parents=None, unit=None):
super().__init__(name, coordinates, unit=unit)
self.parents = list(parents) if parents is not None else None
def _from(self):
return ("derived from %s" % ", ".join(map(repr, self.parents))
if self.parents else "a derived variable")
[docs]
def get_measurements(self):
raise TypeError("%r is %s; it holds no measurements"
% (self.name, self._from()))
[docs]
def update(self, idx, **kwargs):
raise TypeError("%r is %s; a model cannot predict into it -- "
"derive it again after predicting its parents"
% (self.name, self._from()))
class _LatentPart(_Variable):
"""One output of a latent node: its moments and its realizations."""
latent_mean: _Attribute
latent_variance: _Attribute
_ZARR_ATTRS = ("latent_mean", "latent_variance")
_PREDICTED_MARKER = "latent_mean"
_ZARR_HAS_SIMS = True
_COLUMNS = {"latent_mean": ("value", "latent"),
"latent_variance": ("uncertainty", "latent_squared"),
"simulations": ("value", "latent")}
def __init__(self, name, coordinates):
super().__init__(name, coordinates)
self.latent_mean = self._Attribute(coordinates)
self.latent_variance = self._Attribute(coordinates)
self.simulations = None
@classmethod
def from_variable(cls, coordinates, variable):
return cls(variable.name, coordinates)
def get_simulations(self):
if self.simulations is not None:
_refuse_past_threshold(_store_bytes(self.simulations), self.name)
return _np.asarray(self.simulations)
def get_predictions(self):
return self.latent_mean.values.to_numpy()
def allocate_simulations(self, n_sim):
self.simulations = _storage.ArrayStore.allocate(
(self.coordinates.n_data, n_sim),
dtype=_storage.realization_dtype(), fill_value=_np.nan,
owner=self.coordinates)
[docs]
class LatentVariable(_Variable):
"""
What a node inside a model's tree says, at every location.
Written by :meth:`geoml.models.VGPNetwork.predict_node`: one part per
output of the node, each holding the node's mean and variance there and
its realizations. Everything is on the latent scale -- a node has no
likelihood, so nothing is back-transformed, and there are no
measurements, units or cut-offs.
The moments are the node's own, carried up the tree from the inputs.
The realizations carry only the variance a GP node's inducing points
explain, so their spread can fall short of `latent_variance`; above a
nonlinear node (`Exponentiation`, `Multiply`, `GaussianMixture`) the
moments are an approximation and the realizations are the reference.
A model does not train on it: it holds a prediction, not measurements.
Attributes
----------
components : dict
One part per output, keyed by label, each with `latent_mean`,
`latent_variance` and `simulations`.
"""
components: "dict[str, _LatentPart]"
_PART = _LatentPart
_LABEL_KIND = "components"
def __init__(self, name, coordinates, labels):
super().__init__(name, coordinates)
self.labels = [str(label) for label in labels]
self._length = len(self.labels)
self.components = {label: _LatentPart(label, coordinates)
for label in self.labels}
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
return cls(variable.name, coordinates, variable.labels)
@property
def n_sim(self):
return self.components[self.labels[0]].n_sim
[docs]
def get_measurements(self):
raise TypeError("%r is a latent node's prediction; it holds no "
"measurements" % self.name)
[docs]
def get_simulations(self):
stores = [self.components[label].simulations for label in self.labels]
if all(store is not None for store in stores):
_refuse_past_threshold(sum(_store_bytes(s) for s in stores),
self.name)
return _np.stack([self.components[label].get_simulations()
for label in self.labels], axis=2)
[docs]
def get_predictions(self):
return _np.stack([self.components[label].get_predictions()
for label in self.labels], axis=1)
[docs]
def allocate_simulations(self, n_sim):
for part in self.components.values():
part.allocate_simulations(n_sim)
[docs]
def update(self, idx, **kwargs):
"""Writes one batch: `latent_mean` and `latent_variance` as
`(rows, size)`, `simulations` as `(rows, size, n_sim)`."""
if "latent_mean" not in kwargs:
raise TypeError("%r is a latent node's prediction; a model's "
"prediction cannot be written into it"
% self.name)
mean = _np.asarray(kwargs["latent_mean"])
variance = _np.asarray(kwargs["latent_variance"])
sims = _np.asarray(kwargs["simulations"])
for i, label in enumerate(self.labels):
part = self.components[label]
part.latent_mean.values[idx] = mean[:, i]
part.latent_variance.values[idx] = variance[:, i]
part._sim_store()[idx, :] = sims[:, i, :]
[docs]
def unpredicted(self):
"""One boolean per location: True where any part is missing its
mean."""
return _np.any([part.unpredicted()
for part in self.components.values()], axis=0)
[docs]
class VectorVariable(_Variable):
uncertainty: _Attribute
components: "dict[str, ContinuousVariable]"
responsibilities: "dict[int, _Attribute]"
_ZARR_ATTRS = ("uncertainty",)
_PREDICTED_MARKER = "uncertainty"
# the mixture is over the row, so the responsibilities -- and a mixture
# of likelihoods' population per realization -- belong to the variable
# rather than to its components: one answer per location
_DICT_FAMILIES = ("responsibilities",)
_REALIZATION_STORES = ("simulations", "population")
# the mean of the components' latent variances
_COLUMNS = {"uncertainty": ("uncertainty", "latent_squared"),
"responsibilities": ("value", "unit_interval"),
"population": ("value", "classes")}
_FAMILY_KEYS = {"responsibilities": "component"}
_BLOCK_MEANS = ("uncertainty",)
_LABEL_KIND = "components"
def __init__(self, name, coordinates, labels, measurements=None,
units=None):
super().__init__(name, coordinates)
if measurements is not None \
and isinstance(measurements, _pd.DataFrame):
measurements = measurements.values
self.labels = labels
self._length = len(labels)
units = _units_per_label(units, labels)
self.components = {}
for i, label in enumerate(labels):
self.components[label] = ContinuousVariable(
label,
coordinates,
measurements[:, i] if measurements is not None else None,
unit=units[label],
)
self.uncertainty = self._Attribute(coordinates)
# Which noise component each measurement came from, under a mixture
# likelihood; empty under every other one.
self.responsibilities = _col.OrderedDict()
self.population = None
[docs]
def get_measurements(self):
# not allowing partial missing data
out = [self.components[label].measurements.values.to_numpy()
for label in self.labels]
out = _np.stack(out, axis=1)
has_value = _np.all(~ _np.isnan(out), axis=1, keepdims=True) * 1.0
has_value = _np.tile(has_value, [1, self.length])
out = _np.where(has_value == 0.0, 1.0, out)
return out, has_value
[docs]
def get_simulations(self):
# Materializes every component at once -- n_data x n_sim x n_comp in
# RAM -- so the refusal is on the total, not on one component
stores = [self.components[v].simulations for v in self.labels]
if all(store is not None for store in stores):
_refuse_past_threshold(sum(_store_bytes(s) for s in stores),
self.name)
sims = _np.stack([self.components[v].get_simulations() for v in self.labels], axis=2)
return sims
[docs]
def get_predictions(self):
pred = _np.stack([self.components[v].get_predictions() for v in self.labels], axis=1)
return pred
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
# the components are built fresh by `__init__`; what they *know* --
# their cut-offs -- follows by `_copy_attrs_into`, which walks the
# two trees in step so nothing is named here again
return cls(variable.name, coordinates, variable.labels)
[docs]
@classmethod
def from_data_frame(cls, name, coordinates, df, columns=None,
units=None, *args, **kwargs):
new_var = cls(
name,
coordinates,
labels=columns,
measurements=df.loc[:, columns].values,
units=units,
)
return new_var
[docs]
def update(self, idx, elementwise=False, **kwargs):
prediction = _tf.unstack(kwargs["average_sim"], axis=1)
simulations = _tf.unstack(kwargs["simulations"], axis=1)
# each component is dispersed inside a block, and measured, on its own
# account
dispersion = _tf.unstack(kwargs["dispersion"], axis=1)
noise = _tf.unstack(kwargs["noise_variance"], axis=1)
blank = [None] * len(self.labels)
shares = {
key: (_tf.unstack(kwargs[key], axis=1)
if key in kwargs.keys() else blank)
for key in ("proportions", "divided")}
# The latent field's columns are the components' own only under an
# elementwise warping, which is what the model says with
# `elementwise`; a rotation or a projection leaves no column that
# is any one component's, and then the components' latent moments
# stay empty rather than carry a mixture under one label.
moments = {
key: (_tf.unstack(kwargs[key], axis=1)
if elementwise and key in kwargs.keys() else blank)
for key in ("mean", "variance")}
for i, (lb, p, s, d, nv) in enumerate(zip(self.labels, prediction,
simulations, dispersion,
noise)):
values = {"average_sim": p, "simulations": s, "dispersion": d,
"noise_variance": nv}
for key, unstacked in moments.items():
if unstacked[i] is not None:
values[key] = unstacked[i]
for key, unstacked in shares.items():
column = unstacked[i]
if column is not None:
# the matrix was padded out to the widest component, so
# this one takes only the cut-offs it declared
declared = len(self.components[lb].cutoffs or [])
values[key] = column[:, :declared]
if "population_prediction" in kwargs.keys():
values["population_prediction"] = _np.asarray(
kwargs["population_prediction"])[:, i, :]
self.components[lb].update(idx, **values)
_write_mixture(self, idx, kwargs)
self.uncertainty.values[idx] = kwargs["uncertainty"].numpy()
[docs]
def allocate_simulations(self, n_sim):
for comp in self.labels:
self.components[comp].allocate_simulations(n_sim)
[docs]
def reset_quantiles(self, probabilities=None):
for el in self.labels:
self.components[el].reset_quantiles(probabilities)
[docs]
def reset_probabilities(self, quantiles=None):
for el in self.labels:
self.components[el].reset_probabilities(quantiles)
[docs]
def compute_metrics(self, alpha=0.05):
metrics = [self.components[comp].compute_metrics(alpha) for comp in self.labels]
metrics = _pd.concat(metrics, axis=1)
metrics.columns = self.labels
self.metrics = metrics
return self.metrics
class _Component(ContinuousVariable):
# a component reads the composition's latent field rather than one of
# its own, so these two stay empty where its parent fills them
latent_mean: "_Attr | None"
latent_variance: "_Attr | None"
def __init__(self, name, coordinates, measurements=None, unit=None):
super().__init__(name, coordinates, measurements, unit=unit)
self.latent_mean = None
self.latent_variance = None
def set_unit(self, unit):
"""The part's unit, which here is also a divisor.
Parts in different units cannot be added up, so a composition's
unit has to be one the package can convert: a name from `UNITS` or
a number. Checked when it is declared rather than at the door,
where the message would arrive a training run late.
"""
_divisor(unit)
self.unit = unit
return self
def _model_cutoffs(self):
# declared in the part's own unit; the model works in fractions
divisor = self.divisor()
return [c / divisor for c in (self.cutoffs or [])]
def _from_model_cutoffs(self, cutoffs):
divisor = self.divisor()
return [c * divisor for c in cutoffs]
def update(self, idx, **kwargs):
# The model speaks in fractions of the whole -- it has to, since
# the parts are added up -- and this part is stored, reported and
# contoured in the unit it was assayed in. So the crossing happens
# here, once, and everything downstream of it (quantiles, cut-off
# shares, grade-tonnage, the plots) is already in the right unit
# without knowing that units exist. A variance carries the square.
scale = self.divisor()
self.prediction.values[idx] = kwargs["prediction"].numpy() * scale
self._sim_store()[idx, :] = kwargs["simulations"].numpy() * scale
if "dispersion" in kwargs.keys():
self.dispersion.values[idx] = \
kwargs["dispersion"].numpy() * scale ** 2
if "noise_variance" in kwargs.keys():
self.noise_variance.values[idx] = \
kwargs["noise_variance"].numpy() * scale ** 2
for key, target in (("proportions", self.proportions),
("divided", self.divided)):
# a share is a share whatever the unit; the cut-offs they are
# indexed by are this part's own, as declared
if key not in kwargs.keys():
continue
values = kwargs[key].numpy()
for i, cutoff in enumerate(self.cutoffs or []):
if cutoff not in target:
target[cutoff] = self._Attribute(self.coordinates)
target[cutoff].values[idx] = values[:, i]
def allocate_simulations(self, n_sim):
self.simulations = _storage.ArrayStore.allocate(
(self.coordinates.n_data, n_sim),
dtype=_storage.realization_dtype(), fill_value=_np.nan,
owner=self.coordinates)
def get_simulations(self):
return _np.asarray(self.simulations)
[docs]
class CompositionalVariable(VectorVariable):
"""A composition, each part in its own unit.
The parts are stored, reported and simulated in the units they were
measured in -- percent, ppm, g/t -- and turned into fractions of the
whole only where the model reads them, since parts in different units
cannot be added up. `add_compositional_variable` is the door that
prepares them; the two crossings are `get_measurements` here and
`_Component.update` on the way back.
"""
def __init__(self, name, coordinates, labels, measurements=None,
units=None):
super().__init__(name, coordinates, labels, measurements, units)
units = _units_per_label(units, labels)
for i, label in enumerate(labels):
self.components[label] = _Component(
label,
coordinates,
measurements[:, i] if measurements is not None else None,
unit=units[label])
[docs]
def divisors(self):
"""What each part is divided by to become a fraction, in order."""
return _np.array([self.components[label].divisor()
for label in self.labels], dtype=float)
[docs]
def from_model_units(self, values):
"""Model-space values (fractions) in the parts' own units.
`values` has the component axis second, as everything the model
hands back does.
"""
scale = self.divisors()
values = _np.asarray(values, dtype=float)
shape = [1] * values.ndim
shape[1] = scale.size
return values * scale.reshape(shape)
[docs]
def get_measurements(self):
# not allowing partial missing data
out = [self.components[label].measurements.values.to_numpy()
for label in self.labels]
# each part divided by its own unit, so that they can be added up:
# 1 % and 10 000 ppm are the same fraction, and only in fractions
# does a row sum to one
out = _np.stack(out, axis=1) / self.divisors()[None, :]
total = _np.sum(out, axis=1, keepdims=True)
total = _np.where(_np.abs(total - 1) < 1e-10, 1.0, _np.nan)
out = out * total
has_value = _np.all(~ _np.isnan(out), axis=1, keepdims=True) * 1.0
has_value = _np.tile(has_value, [1, self.length]).astype(float)
# out[_np.isnan(out)] = 1
out = _np.where(has_value == 0.0, 1.0, out).astype(float)
return out, has_value
# allowing partial missing data
# out = [self.components[label].measurements.values
# for label in self.labels]
# out = _np.stack(out, axis=1)
# has_value = ~ _np.isnan(out)
# out[_np.isnan(out)] = 1
# return out, has_value
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
new_var = cls(variable.name, coordinates, variable.labels)
return new_var
[docs]
@classmethod
def from_data_frame(cls, name, coordinates, df, columns=None,
units=None, *args, **kwargs):
new_var = cls(
name,
coordinates,
labels=columns,
measurements=df.loc[:, columns].values,
units=units)
return new_var
[docs]
def update(self, idx, **kwargs):
prediction = _tf.unstack(kwargs["average_sim"], axis=1)
simulations = _tf.unstack(kwargs["simulations"], axis=1)
# a part varies inside a block, and is assayed, on its own account
dispersion = _tf.unstack(kwargs["dispersion"], axis=1)
noise = _tf.unstack(kwargs["noise_variance"], axis=1)
blank = [None] * len(self.labels)
shares = {
key: (_tf.unstack(kwargs[key], axis=1)
if key in kwargs.keys() else blank)
for key in ("proportions", "divided")}
for i, (lb, p, s, d, nv) in enumerate(zip(
self.labels,
prediction, simulations, dispersion, noise)):
values = {
"prediction": p,
"simulations": s,
"dispersion": d,
"noise_variance": nv
}
for key, unstacked in shares.items():
column = unstacked[i]
if column is not None:
# padded out to the widest part, as `prediction_input`
# built it; this one takes the cut-offs it declared
declared = len(self.components[lb].cutoffs or [])
values[key] = column[:, :declared]
self.components[lb].update(idx, **values)
self.uncertainty.values[idx] = kwargs["uncertainty"].numpy()
[docs]
def compute_metrics(self, alpha=0.05):
metrics = [self.components[comp].compute_metrics(alpha) for comp in self.labels]
metrics = _pd.concat(metrics, axis=1)
metrics.columns = self.labels
comp_true, has_value = self.get_measurements()
# `get_measurements` closes the parts into fractions, and the
# Aitchison distance is not invariant to scaling one part alone, so
# the predictions are put on the same footing before comparing
comp_pred = _np.stack([self.components[c].prediction.values.to_numpy()
for c in self.labels], axis=1)
comp_pred = comp_pred / self.divisors()[None, :]
comp_true = comp_true[has_value[:, 0] == 1]
comp_pred = comp_pred[has_value[:, 0] == 1]
ad = _gmlmetrics.aitchison_distance(comp_true, comp_pred)
metrics.loc['Aitchison distance', :] = ad
self.metrics = metrics
return self.metrics
class _Category(_Variable):
probability: _Attribute
indicator: _Attribute
indicator_mean: _Attribute
indicator_variance: _Attribute
indicator_predicted: _Attribute
proportions: "dict[float, _Attribute]"
divided: "dict[float, _Attribute]"
_ZARR_ATTRS = ("probability", "indicator", "indicator_mean",
"indicator_variance", "indicator_predicted")
_ZARR_HAS_SIMS = True
_DICT_FAMILIES = ("proportions", "divided")
_COLUMNS = {
"probability": ("value", "unit_interval"),
# 1 where measured in the category, 0 where not, a half at a contact
"indicator": ("measurement", "unit_interval"),
"indicator_mean": ("value", "latent"),
"indicator_variance": ("uncertainty", "latent_squared"),
# the category's log-odds against its best rival
"indicator_predicted": ("value", "log_odds"),
"proportions": ("value", "unit_interval"),
"divided": ("value", "flag"),
"simulations": ("value", "latent"),
}
_FAMILY_KEYS = {"proportions": "cutoff", "divided": "cutoff"}
_BLOCK_MEANS = ("probability", "indicator_mean", "indicator_variance",
"indicator_predicted")
_BLOCK_MEAN_FAMILIES = ("proportions",)
def __init__(self, name, coordinates, indicator):
super().__init__(name, coordinates)
n_data = coordinates.n_data
self.probability = self._Attribute(coordinates, _np.zeros(n_data))
self.indicator = self._Attribute(coordinates, indicator)
self.indicator_mean = self._Attribute(coordinates)
self.indicator_variance = self._Attribute(coordinates)
self.indicator_predicted = self._Attribute(coordinates)
# How much of each block this category holds, and whether the block
# is cut in two by this category's boundary -- see
# `likelihood._divided`. The same dicts a grade keeps, keyed by the
# one cut-off a category has: zero on `ind_skew`, its log-odds
# against its best rival, which is not a number anyone declared but
# the level set the contact *is*. One shape for both kinds, so
# nothing downstream asks which it is holding. `proportions` is a
# different thing from `probability`, which is how sure the model is
# that the block as a whole belongs here.
self.proportions = _col.OrderedDict()
self.divided = _col.OrderedDict()
self.simulations = None
def _share_label(self, cutoff):
# the zero is an artefact of the log-odds, not a number anyone
# declared, so the share keeps a bare name where a grade's says
# `@ 1.5`
return ""
def update(self, idx, **kwargs):
self.indicator_predicted.values[idx] = kwargs["indicator"].numpy()
self.indicator_mean.values[idx] = kwargs["mean"].numpy()
self.indicator_variance.values[idx] = kwargs["variance"].numpy()
self.probability.values[idx] = kwargs["probability"].numpy()
for family in ("proportions", "divided"):
if family in kwargs.keys():
target = getattr(self, family)
if 0.0 not in target:
target[0.0] = self._Attribute(self.coordinates)
target[0.0].values[idx] = kwargs[family].numpy()
self._sim_store()[idx, :] = kwargs["simulations"].numpy()
def allocate_simulations(self, n_sim):
self.simulations = _storage.ArrayStore.allocate(
(self.coordinates.n_data, n_sim),
dtype=_storage.realization_dtype(), fill_value=_np.nan,
owner=self.coordinates)
[docs]
class RockTypeVariable(_Variable):
predicted: _Attribute
entropy: _Attribute
uncertainty: _Attribute
measurements_a: _Attribute
measurements_b: _Attribute
boundary: _Attribute
components: "dict[str, _Category]"
_ZARR_ATTRS = ("predicted", "entropy", "uncertainty",
"measurements_a", "measurements_b", "boundary")
_COLUMNS = {
"predicted": ("value", "classes"),
# divided by the log of the number of categories
"entropy": ("uncertainty", "unit_interval"),
# the root of the latent variance times the entropy
"uncertainty": ("uncertainty", "latent"),
"measurements_a": ("measurement", "classes"),
"measurements_b": ("measurement", "classes"),
"boundary": ("measurement", "flag"),
}
_PREDICTED_MARKER = "entropy"
# taken per sub-block and then averaged (`_resolve`); `predicted` is read
# again off the averaged probabilities (`_coarsen_into`)
_BLOCK_MEANS = ("entropy", "uncertainty")
_LABEL_KIND = "categories"
def __init__(self, name, coordinates, labels=None, measurements_a=None,
measurements_b=None):
# Coerce to numpy arrays: with plain Python lists the element-wise
# comparisons below (indicator building, boundary detection) would
# silently collapse to scalars.
if measurements_a is not None:
measurements_a = _np.asarray(measurements_a)
if measurements_b is not None:
measurements_b = _np.asarray(measurements_b)
if measurements_b is None:
measurements_b = measurements_a
if labels is None:
if measurements_a is None:
raise Exception("either the labels or measurements"
"must be provided")
cat_a = _pd.Categorical(measurements_a)
cat_b = _pd.Categorical(measurements_b)
labels = _pd.api.types.union_categoricals([cat_a, cat_b])
labels = labels.categories.values
n_cat = len(labels)
n_data = coordinates.n_data
avg_vals = None
if measurements_a is not None:
vals_a = _np.zeros([n_data, n_cat])
vals_b = vals_a.copy()
for i, label in enumerate(labels):
vals_a[measurements_a == label, i] = 1
vals_b[measurements_b == label, i] = 1
avg_vals = 0.5 * (vals_a + vals_b)
super().__init__(name, coordinates)
self.labels = labels
self._length = len(labels)
self.components = {}
for i, label in enumerate(labels):
self.components[label] = _Category(
label,
coordinates,
avg_vals[:, i] if avg_vals is not None else None,
)
# the categories are known here, so these three hold codes into
# `labels` rather than one string object per data location
self.predicted = self._Attribute.encoded(coordinates, labels=labels)
self.entropy = self._Attribute(coordinates)
self.uncertainty = self._Attribute(coordinates)
if measurements_a is None:
self.measurements_a = self._Attribute.encoded(
coordinates, labels=labels)
self.measurements_b = self._Attribute.encoded(
coordinates, labels=labels)
self.boundary = self._Attribute(
coordinates, [False]*n_data, dtype=bool)
else:
# the variable's categories first, in its order, so that a class
# has one code in every column; then whatever else was measured,
# since a measurement outside `labels` is still worth keeping as
# what it says
self.measurements_a = _measured_in_order(
self._Attribute, coordinates, measurements_a, labels)
self.measurements_b = _measured_in_order(
self._Attribute, coordinates, measurements_b, labels)
self.boundary = self._Attribute(
coordinates, measurements_a != measurements_b, dtype=bool)
[docs]
def get_measurements(self):
# not allowing partial missing data
out = [self.components[label].indicator.values.to_numpy()
for label in self.labels]
out = _np.stack(out, axis=1)
total = _np.sum(out, axis=1, keepdims=True)
total = _np.where(_np.abs(total - 1) < 1e-10, 1.0, _np.nan)
out = out * total
has_value = _np.all(~ _np.isnan(out), axis=1, keepdims=True) * 1.0
has_value = _np.tile(has_value, [1, self.length]).astype(float)
out = _np.where(has_value == 0.0, 1.0, out).astype(float)
return out, has_value
[docs]
def allocate_simulations(self, n_sim):
for comp in self.labels:
self.components[comp].allocate_simulations(n_sim)
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
new_var = cls(variable.name, coordinates, variable.labels)
return new_var
[docs]
@classmethod
def from_data_frame(cls, name, coordinates, df, col_a=None, col_b=None,
*args, **kwargs):
labels = _pd.api.types.union_categoricals(
[_pd.Categorical(df[col_a]),
_pd.Categorical(df[col_b])])
labels = labels.categories.values
new_var = cls(name, coordinates, labels,
measurements_a=df[col_a].values,
measurements_b=df[col_b].values)
return new_var
def __getitem__(self, item):
new_obj = _copy_for_subset(self)
new_obj.entropy = self.entropy[item]
new_obj.uncertainty = self.uncertainty[item]
new_obj.boundary = self.boundary[item]
new_obj.measurements_a = self.measurements_a[item]
new_obj.measurements_b = self.measurements_b[item]
# Only the categories still present: slicing is a pre-processing step
# before training, and a category with no data left in it has nothing
# to teach the model. The parent's order is kept, since it is the order
# of the components and of the codes below.
present = set(new_obj.measurements_a.to_numpy()) \
| set(new_obj.measurements_b.to_numpy())
labels = [label for label in self.labels if label in present]
if len(labels) == 0:
# nothing is measured here at all — a prediction target, say, whose
# categories come from the model rather than from the data
labels = list(self.labels)
new_obj.labels = labels
new_obj._length = len(labels)
# the components and the prediction follow the labels: `update` writes
# the winning label's position, which the dropped ones would shift
new_obj.components = {label: self.components[label][item]
for label in labels}
predicted = self.predicted[item]
predicted.values = _encode(predicted.to_numpy(), labels)
predicted.labels = list(labels)
new_obj.predicted = predicted
return new_obj
[docs]
def update(self, idx, **kwargs):
self.entropy.values[idx] = kwargs["entropy"].numpy()
self.uncertainty.values[idx] = kwargs["uncertainty"].numpy()
# the winning category's position is the code
self.predicted.values[idx] = _np.argmax(
kwargs["probability"].numpy(), axis=1)
mean = _tf.unstack(kwargs["mean"], axis=1)
variance = _tf.unstack(kwargs["variance"], axis=1)
indicators = _tf.unstack(kwargs["indicators"], axis=1)
probability = _tf.unstack(kwargs["probability"], axis=1)
simulations = _tf.unstack(kwargs["simulations"], axis=1)
# the share of each block this category holds, one column per label
blank = [None] * len(self.labels)
proportions = _tf.unstack(kwargs["proportions"], axis=1) \
if "proportions" in kwargs.keys() else blank
divided = _tf.unstack(kwargs["divided"], axis=1) \
if "divided" in kwargs.keys() else blank
for lb, m, v, i, p, s, share, cut in zip(
self.labels, mean, variance, indicators,
probability, simulations, proportions, divided):
values = {"mean": m, "variance": v, "indicator": i,
"probability": p, "simulations": s}
if share is not None:
values["proportions"] = share
values["divided"] = cut
self.components[lb].update(idx, **values)
def _coarsen_into(self, new, grouping, valid=None):
# A category's probability reads 0 where nothing was predicted, not
# missing, so the parts that were are told apart by their label
predicted = self.predicted.values.to_numpy() >= 0
valid = predicted if valid is None else valid & predicted
super()._coarsen_into(new, grouping, valid)
# the winner of the averaged probabilities, as `update` names it
probability = _np.stack(
[new.components[label].probability.values.to_numpy()
for label in self.labels], axis=1)
gathered = ~grouping.whole & _np.all(_np.isfinite(probability),
axis=1)
codes = new.predicted.values.to_numpy().copy()
codes[gathered] = _np.argmax(probability[gathered], axis=1)
new.predicted.values[:] = codes
[docs]
def compute_metrics(self, decluster: bool = False) -> "_pd.DataFrame":
"""
Scores this variable's prediction against its own measurements.
Every score is of one category against the rest, so the table has
a column per category. The locations scored are the ones the
confusion matrix counts: measured, predicted, and off the contacts,
where a location holds two measurements and no one truth.
- Balanced accuracy, Jaccard, Matthews, Cohen's kappa, precision,
recall and F1 score read the predicted category. Recall is the
share of the locations measured as the category that the model
calls it; precision is the share of the locations the model
calls it that were measured as it.
- Quantity and allocation disagreement split the category's errors
into the part a wrong proportion explains and the part a wrong
place explains (Pontius and Millones, 2011). Summed over the
categories and halved, they add up to one minus the accuracy.
- The Brier score and the log score read the predicted
probability, as a forecast of whether a location is the
category. Both are proper: neither hedging nor overconfidence
improves them. Lower is better.
Parameters
----------
decluster
Weight each location by the container's `"declustering"`
metadata column, which the container's `decluster` method
writes, so that densely drilled ground does not dominate.
Returns
-------
pandas.DataFrame
One row per score, one column per category.
Raises
------
ValueError
If no location holds both a measurement and a prediction, or if
`decluster` is asked for and the container keeps no weights.
Notes
-----
A score with no value is NaN: the precision of a category the model
never calls, the recall of one never measured. At the locations a
model was trained on, every score flatters it; the out-of-fold
container :func:`geoml.models.cross_validate` fills is the honest
input.
References
----------
Brier, G. W. (1950). Verification of forecasts expressed in terms of
probability. *Monthly Weather Review*, 78(1), 1-3.
Cohen, J. (1960). A coefficient of agreement for nominal scales.
*Educational and Psychological Measurement*, 20(1), 37-46.
Pontius, R. G., & Millones, M. (2011). Death to Kappa: birth of
quantity disagreement and allocation disagreement for accuracy
assessment. *International Journal of Remote Sensing*, 32(15),
4407-4429.
"""
y_pred = self.predicted.to_numpy()
y_true_a = self.measurements_a.to_numpy()
y_true_b = self.measurements_b.to_numpy()
# the missing code decodes to the empty string on both sides
valid = (y_true_a == y_true_b) & (y_true_a != "") & (y_pred != "")
if not valid.any():
raise ValueError(
"no location holds both a measured category and a "
"prediction of %r; predict on the data first" % self.name)
weights = None
if decluster:
column = self.coordinates.metadata.get("declustering")
if column is None:
raise ValueError(
"the container keeps no 'declustering' column to weight "
"the locations by; run its decluster() method first")
weights = _np.asarray(column.values, dtype=float)[valid]
share = _np.ones(valid.sum()) if weights is None else weights
share = share / share.sum()
y_pred = y_pred[valid]
y_true = y_true_a[valid]
series = []
for lab in self.labels:
yp = _np.where(y_pred == lab, 1, 0)
yt = _np.where(y_true == lab, 1, 0)
claim = self.components[lab].probability.values.to_numpy()[valid]
# the category's hits and its two errors, as shares of the
# locations: called it where another was measured (commission),
# measured as it where another was called (omission)
hits = share @ (yt * yp)
commission = share @ (yp * (1 - yt))
omission = share @ (yt * (1 - yp))
with _np.errstate(invalid="ignore", divide="ignore"):
# 0/0 where there is nothing to score: the precision of a
# category never called, the recall of one never measured,
# the kappa of one neither
precision = hits / (hits + commission)
recall = hits / (hits + omission)
f1 = 2 * hits / (2 * hits + commission + omission)
kappa = _skmetrics.cohen_kappa_score(
yt, yp, labels=[0, 1], sample_weight=weights)
# a claim of exactly zero for the category measured would score
# infinity; clipped at machine precision, as scikit-learn clips
eps = _np.finfo(float).eps
clipped = _np.clip(claim, eps, 1 - eps)
d = {
'Balanced accuracy': _skmetrics.balanced_accuracy_score(
yt, yp, sample_weight=weights),
'Jaccard': _skmetrics.jaccard_score(
yt, yp, sample_weight=weights),
'Matthews': _skmetrics.matthews_corrcoef(
yt, yp, sample_weight=weights),
"Cohen's kappa": kappa,
'Precision': precision,
'Recall': recall,
'F1 score': f1,
'Quantity disagreement': abs(commission - omission),
'Allocation disagreement': 2 * min(commission, omission),
'Brier score': share @ (claim - yt) ** 2,
'Log score': -share @ (yt * _np.log(clipped)
+ (1 - yt) * _np.log1p(-clipped)),
}
series.append(_pd.Series(d, name=lab))
self.metrics = _pd.concat(series, axis=1)
return self.metrics
[docs]
class CategoricalVariable(RockTypeVariable):
def __init__(self, name, coordinates, labels=None, measurements=None):
super().__init__(name, coordinates, labels=labels, measurements_a=measurements)
[docs]
@classmethod
def from_data_frame(cls, name, coordinates, df, measurements_col=None,
*args, **kwargs):
labels = _pd.Categorical(df[measurements_col])
labels = labels.categories.values
new_var = cls(name, coordinates, labels,
measurements=df[measurements_col].values)
return new_var
[docs]
class OrderedRockType(RockTypeVariable):
_ZARR_ATTRS = RockTypeVariable._ZARR_ATTRS + ("implicit_values",)
_COLUMNS = dict(RockTypeVariable._COLUMNS,
implicit_values=("measurement", "ordinal"))
def __init__(self, name, coordinates, labels=None, measurements_a=None,
measurements_b=None):
super().__init__(name, coordinates, labels, measurements_a,
measurements_b)
self._length = 1
if measurements_a is None:
# Labels-only construction (prediction targets, reloading from
# disk): no measured contacts, all implicit values missing.
self.implicit_values = self._Attribute(coordinates)
return
measurements_a = _np.asarray(measurements_a)
if measurements_b is None:
measurements_b = measurements_a
else:
measurements_b = _np.asarray(measurements_b)
# Labels may have been derived from the measurements by the parent.
self.implicit_values = self._Attribute(
coordinates,
self._implicit_values(measurements_a, measurements_b, self.labels))
@staticmethod
def _implicit_values(measurements_a, measurements_b, labels):
"""Where each pair of measurements sits in the label sequence.
A function of the labels, so it has to be recomputed whenever they
change — slicing away a category shifts every position after it.
"""
# Built as a float array from the start; the previous
# ``-0.5 * _np.ones_like(measurements_a)`` crashed on string inputs.
implicit_values = _np.full(len(measurements_a), -0.5)
for i in range(len(labels[:-1])):
implicit_values = _np.where(
(measurements_a == labels[i])
& (measurements_b == labels[i + 1]),
i,
implicit_values
)
implicit_values = _np.where(
(measurements_a == labels[i + 1])
& (measurements_b == labels[i]),
i,
implicit_values
)
implicit_values = _np.where(
(measurements_a == labels[i])
& (measurements_b == labels[i]),
i - 0.5,
implicit_values
)
implicit_values = _np.where(
(measurements_a == labels[i + 1])
& (measurements_b == labels[i + 1]),
i + 0.5,
implicit_values
)
# implicit_values[(measurements_a == labels[i])
# & (measurements_b == labels[i + 1])] = i
# implicit_values[(measurements_a == labels[i + 1])
# & (measurements_b == labels[i])] = i
# implicit_values[(measurements_a == labels[i])
# & (measurements_b == labels[i])] = i - 0.5
# implicit_values[(measurements_a == labels[i + 1])
# & (measurements_b == labels[i + 1])] = i + 0.5
return implicit_values
[docs]
def get_measurements(self):
values = self.implicit_values.values.copy()[:, None]
has_value = (~ _np.isnan(values)) * 1.0
values[_np.isnan(values)] = 0
return values, has_value
def __getitem__(self, item):
new_obj = super().__getitem__(item)
implicit = self.implicit_values[item]
if new_obj.measurements_a._has_content():
# positions in the label sequence, and the slice may have dropped
# labels, so these are recomputed rather than carried over
implicit.values = self._implicit_values(
new_obj.measurements_a.to_numpy(),
new_obj.measurements_b.to_numpy(),
new_obj.labels)
new_obj.implicit_values = implicit
return new_obj
[docs]
class BinaryVariable(_Variable):
indicator: _Attribute
measurements: _Attribute
weights: _Attribute
predicted: _Attribute
probability: _Attribute
entropy: _Attribute
uncertainty: _Attribute
latent_mean: _Attribute
latent_variance: _Attribute
_ZARR_ATTRS = ("indicator", "measurements", "weights", "predicted",
"probability", "entropy", "uncertainty",
"latent_mean", "latent_variance")
_COLUMNS = {
"indicator": ("measurement", "unit_interval"),
"measurements": ("measurement", "classes"),
"weights": ("weight", "dimensionless"),
"predicted": ("value", "classes"),
"probability": ("value", "unit_interval"),
# divided by the log of two
"entropy": ("uncertainty", "unit_interval"),
"uncertainty": ("uncertainty", "latent"),
"latent_mean": ("value", "latent"),
"latent_variance": ("uncertainty", "latent_squared"),
# the probability each realization gives
"simulations": ("value", "unit_interval"),
}
_PREDICTED_MARKER = "entropy"
# `entropy` and `uncertainty` are functions of the block's own
# probability, not means over it, and stay out; `predicted` is read
# again off the averaged probability (`_coarsen_into`)
_BLOCK_MEANS = ("probability", "latent_mean", "latent_variance")
_LABEL_KIND = "categories"
_ZARR_HAS_SIMS = True
def __init__(self, name, coordinates, labels=None, measurements=None):
super().__init__(name, coordinates)
n_data = coordinates.n_data
# Coerce to a numpy array: with a plain Python list the element-wise
# comparisons below (indicator and weights) silently collapse to
# scalars and the indicators stay NaN.
if measurements is not None:
measurements = _np.asarray(measurements)
if labels is None:
if measurements is None:
raise Exception("either the labels or measurements"
"must be provided")
cat = _pd.Categorical(measurements)
labels = cat.categories.values
self.labels = labels
self._length = 1
if len(labels) != 2:
raise ValueError(f"There must be exactly 2 labels - found {len(labels)}.")
self.indicator = self._Attribute(
coordinates, _np.array([_np.nan]*n_data))
if measurements is None:
self.measurements = self._Attribute.encoded(
coordinates, labels=labels)
self.weights = self._Attribute(coordinates)
else:
# their own categories: an `AnomalyVariable` labels everything that
# is not the anomaly `_dummy`, but the measurements say what it was
self.measurements = self._Attribute.encoded(
coordinates, measurements)
self.weights = self._Attribute(coordinates, _np.ones(n_data))
self.indicator.values[measurements == labels[0]] = 1
self.indicator.values[measurements == labels[1]] = 0
self.predicted = self._Attribute.encoded(coordinates, labels=labels)
self.probability = self._Attribute(coordinates, _np.zeros(n_data))
self.entropy = self._Attribute(coordinates)
self.uncertainty = self._Attribute(coordinates)
self.latent_mean = self._Attribute(coordinates)
self.latent_variance = self._Attribute(coordinates)
self.simulations = None
if measurements is not None:
for label in labels:
idx = measurements == label
n_in_label = _np.sum(idx)
if n_in_label > 0:
self.weights.values[idx] = \
n_data / (n_in_label * len(labels))
[docs]
def get_measurements(self):
values = self.indicator.values.copy()[:, None]
has_value = (~ _np.isnan(values)) * 1.0
values[_np.isnan(values)] = 0
return values, has_value
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
new_var = cls(variable.name, coordinates, variable.labels)
return new_var
[docs]
def update(self, idx, **kwargs):
prob = kwargs["probability"].numpy()
mean = kwargs["mean"].numpy()
var = kwargs["variance"].numpy()
entropy = kwargs["entropy"].numpy()
uncertainty = kwargs["uncertainty"].numpy()
sims = kwargs["simulations"].numpy()
if len(prob.shape) > 1:
prob = prob[:, 0]
mean = mean[:, 0]
var = var[:, 0]
# entropy = entropy[:, 0]
# uncertainty = uncertainty[:, 0]
sims = sims[:, 0, :]
label_idx = _np.zeros(prob.shape, dtype=int) # positive class
label_idx[prob < 0.5] = 1 # negative class
# the label's position is the code
self.predicted.values[idx] = label_idx
self.latent_mean.values[idx] = mean
self.latent_variance.values[idx] = var
self.entropy.values[idx] = entropy
self.uncertainty.values[idx] = uncertainty
self.probability.values[idx] = prob
self._sim_store()[idx, :] = sims
def _coarsen_into(self, new, grouping, valid=None):
# the probability reads 0 where nothing was predicted, not missing,
# so the parts that were are told apart by their label
predicted = self.predicted.values.to_numpy() >= 0
valid = predicted if valid is None else valid & predicted
super()._coarsen_into(new, grouping, valid)
# the positive class at one half or more, as `update` has it
probability = new.probability.values.to_numpy()
gathered = ~grouping.whole & _np.isfinite(probability)
codes = new.predicted.values.to_numpy().copy()
codes[gathered] = _np.where(probability[gathered] < 0.5, 1, 0)
new.predicted.values[:] = codes
[docs]
def allocate_simulations(self, n_sim):
self.simulations = _storage.ArrayStore.allocate(
(self.coordinates.n_data, n_sim),
dtype=_storage.realization_dtype(), fill_value=_np.nan,
owner=self.coordinates)
[docs]
@classmethod
def from_data_frame(cls, name, coordinates, df, col, positive_class):
labels = _pd.Categorical(df[col])
labels = labels.categories.values.tolist()
pos = None
for i, label in enumerate(labels):
if label == positive_class:
pos = i
labels.pop(pos)
labels.append(positive_class)
labels = labels[::-1]
new_var = cls(name, coordinates, labels,
measurements=df[col].values)
return new_var
[docs]
class AnomalyVariable(BinaryVariable):
def __init__(self, name, coordinates, label, measurements=None):
labels = [label, "_dummy"]
super().__init__(name, coordinates, labels, measurements)
[docs]
@classmethod
def from_data_frame(cls, name, coordinates, df, col, positive_class):
new_var = cls(name, coordinates, positive_class,
measurements=df[col].values)
return new_var
[docs]
@classmethod
def from_variable(cls, coordinates, variable):
new_var = cls(variable.name, coordinates, variable.labels[0])
return new_var
# the parts' own columns: a vector variable's components and a categorical
# one's categories, declared once where their classes are
VectorVariable._PART = ContinuousVariable
CompositionalVariable._PART = _Component
RockTypeVariable._PART = _Category