# metrics/distribution.py
"""Metrics on probability distributions."""
from __future__ import annotations
import inspect
import math
from collections.abc import Callable
from collections.abc import Iterable
from contextlib import ContextDecorator
from math import log2
from typing import Any
import numpy as np
from numpy.typing import ArrayLike
from scipy.spatial.distance import cdist
from scipy.special import entr
from scipy.special import rel_entr
from pyphi import numerics
from pyphi import utils
from pyphi import validate
from pyphi.cache import joblib_memory
from pyphi.conf import config
from pyphi.direction import Direction
from pyphi.distribution import flatten
from pyphi.distribution import marginal
from pyphi.distribution import marginal_zero
from pyphi.exceptions import MissingOptionalDependenciesError
from pyphi.measures.protocols import CompositeMeasure
from pyphi.measures.protocols import DistributionMeasure
from pyphi.measures.protocols import StateAwareMeasure
from pyphi.measures.protocols import StatefulDistributionMeasure
from pyphi.measures.protocols import satisfies_composite_measure
from pyphi.measures.protocols import satisfies_distribution_measure
from pyphi.measures.protocols import satisfies_state_aware_measure
from pyphi.measures.protocols import satisfies_stateful_distribution_measure
from pyphi.registry import Registry
from pyphi.types import Repertoire
from pyphi.types import State
_LN_OF_2 = np.log(2)
# ---------------------------------------------------------------------------
# Alphabet-support helpers
# ---------------------------------------------------------------------------
def _any_alphabet(_alphabet_sizes: tuple[int, ...]) -> bool:
"""Return ``True`` for any combination of node cardinalities."""
return True
[docs]
class DistanceResult(float):
"""A :class:`float` that carries auxiliary data about its computation.
A distance or φ value with arbitrary metadata attached, so the provenance of
a result (the measure used, the direction, the winning partition, and so on)
remains inspectable after the number has been produced.
Comparison and arithmetic follow exact :class:`float` semantics. Metadata is
preserved through :func:`min` / :func:`max` and copying, but is dropped by
plain arithmetic, which returns an ordinary float. Tolerant comparison of φ
values at decision points — where two values equal to within
``config.numerics.precision`` should be treated as equal — is performed
explicitly through :mod:`pyphi.numerics` and
:func:`pyphi.resolve_ties.resolve_ties`, not by the comparison operators of
this type.
Parameters
----------
value
The numeric value.
**kwargs
Arbitrary keyword arguments stored as metadata attributes.
Attributes
----------
All :class:`float` attributes are available, plus any metadata passed as
keyword arguments.
Notes
-----
To build a NumPy float64 array from a collection of ``DistanceResult``
objects, use ``DistanceResult.values_array(results)``: it extracts the
float values and drops the metadata. The metadata stays available on
the individual results.
Examples
--------
Basic usage with metadata:
>>> from pyphi.measures.distribution import DistanceResult
>>> result = DistanceResult(0.5, method='EMD', direction='CAUSE')
>>> float(result) # Extract numeric value
0.5
>>> result.method # Access metadata
'EMD'
>>> result.direction
'CAUSE'
Mathematical operations preserve the numeric value:
>>> result + 0.3
0.8
>>> result * 2
1.0
>>> result > 0.3
True
Type preservation in min/max with metadata:
>>> results = [
... DistanceResult(0.5, method='EMD'),
... DistanceResult(0.3, method='L1'),
... DistanceResult(0.7, method='GID')
... ]
>>> min_result = min(results)
>>> float(min_result)
0.3
>>> min_result.method # Metadata from the minimum value is preserved
'L1'
NumPy array creation (explicit float extraction):
>>> results = [DistanceResult(0.5), DistanceResult(0.3), DistanceResult(0.7)]
>>> arr = DistanceResult.values_array(results)
>>> arr.dtype
dtype('float64')
>>> arr
array([0.5, 0.3, 0.7])
Serialization through :mod:`pyphi.serialize` preserves the auxiliary
metadata:
>>> from pyphi import serialize
>>> result = DistanceResult(0.5, method='EMD', direction='CAUSE')
>>> restored = serialize.loads(serialize.dumps(result))
>>> restored.method
'EMD'
Typical scientific workflow:
>>> # Compute multiple phi values
>>> phi_values = [
... DistanceResult(0.5, method='EMD', system='ABC'),
... DistanceResult(0.3, method='L1', system='ABC'),
... DistanceResult(0.7, method='GID', system='DEF')
... ] # doctest: +SKIP
>>> # Find maximum
>>> max_phi = max(phi_values) # doctest: +SKIP
>>> print(f"Max φ = {max_phi:.3f} using {max_phi.method}") # doctest: +SKIP
Max φ = 0.700 using GID
>>> # Statistical analysis
>>> # Explicitly extract floats for statistics (metadata dropped):
>>> phi_array = DistanceResult.values_array(phi_values) # doctest: +SKIP
>>> np.mean(phi_array) # doctest: +SKIP
0.5
"""
def __new__(cls, value, **kwargs):
instance = super().__new__(cls, value)
for key, val in kwargs.items():
setattr(instance, key, val)
return instance
def _public_aux_data(self) -> dict:
"""Auxiliary data the user attached at construction or via setattr.
Excludes underscore-prefixed names, which are reserved for internal
use."""
return {k: v for k, v in self.__dict__.items() if not k.startswith("_")}
def __repr__(self):
aux_data = self._public_aux_data()
if aux_data:
aux_str = ", ".join(f"{k}={v!r}" for k, v in aux_data.items())
return f"DistanceResult({float(self)}, {aux_str})"
return f"DistanceResult({float(self)})"
def __str__(self):
"""Short representation for use in formatted output."""
from pyphi.models import fmt
return fmt.fmt_number(float(self))
def __copy__(self):
"""Ensure auxiliary data is preserved when copying."""
return DistanceResult(float(self), **self._public_aux_data())
def __deepcopy__(self, memo):
"""Ensure auxiliary data is preserved when deep copying."""
import copy
aux_data = {
k: copy.deepcopy(v, memo) for k, v in self._public_aux_data().items()
}
return DistanceResult(float(self), **aux_data)
[docs]
@classmethod
def values_array(
cls, results: Iterable[DistanceResult], dtype: Any = None
) -> np.ndarray:
"""Return the float values of an iterable of ``DistanceResult``\\ s
as a NumPy array.
Auxiliary metadata (``method``, ``state``, etc.) is intentionally
dropped — callers that need it should iterate the input directly.
Use this method when you want explicit control over the metadata-loss
boundary, rather than relying on implicit ``np.array(results)``
coercion, which silently drops metadata.
Parameters
----------
results
Iterable of DistanceResult objects.
dtype
Optional NumPy dtype for the array (default: float64).
Returns
-------
np.ndarray
1-D array of the float values.
Examples
--------
>>> results = [DistanceResult(0.5, method='EMD'),
... DistanceResult(0.3, method='L1')]
>>> arr = DistanceResult.values_array(results)
>>> arr.dtype
dtype('float64')
>>> arr
array([0.5, 0.3])
"""
import numpy as np
if dtype is None:
dtype = np.float64
return np.fromiter((float(r) for r in results), dtype=dtype)
[docs]
class OptionalEMD:
"""Class to handle EMD computations via POT (Python Optimal Transport).
Allows deferring import of ``ot`` in case it is not needed.
"""
def __init__(self) -> None:
self._ot = None
@property
def ot(self) -> Any:
if self._ot is None:
try:
import ot
self._ot = ot
except ModuleNotFoundError as exc:
raise ModuleNotFoundError(
MissingOptionalDependenciesError.MSG.format(dependencies="pot")
) from exc
return self._ot
def compute(self, first_histogram, second_histogram, distance_matrix) -> float:
first = np.asarray(first_histogram, dtype=float)
second = np.asarray(second_histogram, dtype=float)
# POT's network simplex normalizes by total mass and returns NaN (with a
# divide-by-zero warning) when both signatures are empty; pyemd returned
# 0, and the EMD between two zero-mass distributions is 0 by definition
# (e.g. the CES-distance between two empty constellations).
if not first.any() and not second.any():
return 0.0
# ``ot.emd2`` returns the exact EMD *cost* (the unique optimal-transport
# objective). ``check_marginals=False`` skips POT's mass-balance warning:
# callers pass equal-mass distributions (normalized repertoires, or the
# balanced phi signatures in ``ces.py``), and the check only warns -- it
# does not affect the computed cost.
return float(self.ot.emd2(first, second, distance_matrix, check_marginals=False))
# Usage
EMD = OptionalEMD()
[docs]
class DistributionMeasureRegistry(Registry):
"""Storage for ``(p, q) -> float | DistanceResult`` distance functions.
Each registered function is validated at registration time to have
exactly two required positional parameters named ``p`` and ``q``;
signature drift fails at import rather than deep in a phi
computation. The ``asymmetric`` flag attaches to the function as an
attribute so callers can filter without consulting a parallel list.
Examples
--------
>>> @distribution_measures.register('ALWAYS_ZERO') # doctest: +SKIP
... def always_zero(p, q):
... return 0
"""
# pylint: disable=arguments-differ
desc = "distribution-to-distribution distance functions"
[docs]
def register( # type: ignore[override]
self,
name: str,
asymmetric: bool = False,
supports_alphabet: Callable[[tuple[int, ...]], bool] = _any_alphabet,
) -> Callable[[Callable[..., float]], Callable[..., float]]:
"""Decorator for registering a :class:`DistributionMeasure`.
Parameters
----------
name
The name of the measure.
asymmetric
``True`` if the measure is asymmetric. Stored as
an attribute on the function.
supports_alphabet
Callable ``(alphabet_sizes) -> bool``
indicating whether the measure handles the given node
cardinalities. Defaults to :data:`_any_alphabet`.
"""
def register_func(func: Callable[..., float]) -> Callable[..., float]:
if not satisfies_distribution_measure(func):
raise TypeError(
f"Cannot register {func!r} as DistributionMeasure {name!r}: "
f"required params must be exactly (p, q); got "
f"{list(inspect.signature(func).parameters)}."
)
func.name = name # type: ignore[attr-defined]
func.asymmetric = asymmetric # type: ignore[attr-defined]
func.supports_alphabet = supports_alphabet # type: ignore[attr-defined]
self.store[name] = func
return func
return register_func
[docs]
class StateAwareMeasureRegistry(Registry):
"""Storage for ``(p, state) -> float | DistanceResult`` measures.
The function reads off a single state's value from a single
distribution.
"""
# pylint: disable=arguments-differ
desc = "pointwise state-aware measures"
[docs]
def register( # type: ignore[override]
self,
name: str,
supports_alphabet: Callable[[tuple[int, ...]], bool] = _any_alphabet,
) -> Callable[[Callable[..., float]], Callable[..., float]]:
"""Decorator for registering a :class:`StateAwareMeasure`.
Parameters
----------
name
The name of the measure.
supports_alphabet
Callable ``(alphabet_sizes) -> bool``
indicating whether the measure handles the given node
cardinalities. Defaults to :data:`_any_alphabet`.
"""
def register_func(func: Callable[..., float]) -> Callable[..., float]:
if not satisfies_state_aware_measure(func):
raise TypeError(
f"Cannot register {func!r} as StateAwareMeasure {name!r}: "
f"required params must be exactly (p, state); got "
f"{list(inspect.signature(func).parameters)}."
)
func.name = name # type: ignore[attr-defined]
func.supports_alphabet = supports_alphabet # type: ignore[attr-defined]
self.store[name] = func
return func
return register_func
[docs]
class CompositeMeasureRegistry(Registry):
"""Storage for composite measures of shape
``(forward, partitioned, selectivity, *, state) -> DistanceResult``.
Used at the system / mechanism boundary by GID, INTRINSIC_SPECIFICATION,
and INTRINSIC_INFORMATION.
"""
# pylint: disable=arguments-differ
desc = "composite measures"
[docs]
def register( # type: ignore[override]
self,
name: str,
asymmetric: bool = False,
applies_intrinsic_information_requirement: bool = False,
partition_measure: Callable[..., Any] | None = None,
supports_alphabet: Callable[[tuple[int, ...]], bool] = _any_alphabet,
) -> Callable[[Callable[..., Any]], Callable[..., Any]]:
"""Decorator for registering a :class:`CompositeMeasure`.
Parameters
----------
name
The name of the measure.
asymmetric
``True`` if the measure is asymmetric.
applies_intrinsic_information_requirement
``True`` if the measure applies the intrinsic-information
requirement (Mayner et al. 2026, Eq. 23;
``INTRINSIC_INFORMATION`` only).
partition_measure
The composite measure used to score
partitions when this measure is the system measure;
``None`` means "use self".
supports_alphabet
Callable ``(alphabet_sizes) -> bool``
indicating whether the measure handles the given node
cardinalities. Defaults to :data:`_any_alphabet`.
"""
def register_func(func: Callable[..., Any]) -> Callable[..., Any]:
if not satisfies_composite_measure(func):
raise TypeError(
f"Cannot register {func!r} as CompositeMeasure {name!r}: "
f"first three params must include 'forward', "
f"'partitioned', and 'selectivity'; got "
f"{list(inspect.signature(func).parameters)}."
)
func.name = name # type: ignore[attr-defined]
func.asymmetric = asymmetric # type: ignore[attr-defined]
func.applies_intrinsic_information_requirement = ( # type: ignore[attr-defined]
applies_intrinsic_information_requirement
)
func.partition_measure = partition_measure # type: ignore[attr-defined]
func.supports_alphabet = supports_alphabet # type: ignore[attr-defined]
self.store[name] = func
return func
return register_func
[docs]
class StatefulDistributionMeasureRegistry(Registry):
"""Storage for ``(p, q, state) -> float | DistanceResult`` measures.
Both distributions contribute to the result; the state selects an
element. Used by IIT_4.0_SMALL_PHI variants and APMI.
"""
# pylint: disable=arguments-differ
desc = "two-distribution state-aware measures"
[docs]
def register( # type: ignore[override]
self,
name: str,
asymmetric: bool = False,
supports_alphabet: Callable[[tuple[int, ...]], bool] = _any_alphabet,
) -> Callable[[Callable[..., float]], Callable[..., float]]:
"""Decorator for registering a :class:`StatefulDistributionMeasure`.
Parameters
----------
name
The name of the measure.
asymmetric
``True`` if the measure is asymmetric. Stored as
an attribute on the function.
supports_alphabet
Callable ``(alphabet_sizes) -> bool``
indicating whether the measure handles the given node
cardinalities. Defaults to :data:`_any_alphabet`.
"""
def register_func(func: Callable[..., float]) -> Callable[..., float]:
if not satisfies_stateful_distribution_measure(func):
raise TypeError(
f"Cannot register {func!r} as StatefulDistributionMeasure "
f"{name!r}: required params must be exactly (p, q, state); "
f"got {list(inspect.signature(func).parameters)}."
)
func.name = name # type: ignore[attr-defined]
func.asymmetric = asymmetric # type: ignore[attr-defined]
func.supports_alphabet = supports_alphabet # type: ignore[attr-defined]
self.store[name] = func
return func
return register_func
distribution_measures = DistributionMeasureRegistry()
state_aware_measures = StateAwareMeasureRegistry()
composite_measures = CompositeMeasureRegistry()
stateful_distribution_measures = StatefulDistributionMeasureRegistry()
[docs]
class ActualCausationMeasureRegistry(Registry):
"""Storage for distance functions used in :mod:`pyphi.actual`.
Users can define custom measures:
Examples
--------
>>> @actual_causation_measures.register('ALWAYS_ZERO') # doctest: +SKIP
... def always_zero(a, b):
... return 0
And use them by setting, *e.g.*, ``config.repertoire_distance = 'ALWAYS_ZERO'``.
"""
# pylint: disable=arguments-differ
desc = "distance functions for use in actual causation calculations"
def __init__(self) -> None:
super().__init__()
self._asymmetric: list[str] = []
[docs]
def register(
self,
name: str,
asymmetric: bool = False,
supports_alphabet: Callable[[tuple[int, ...]], bool] = _any_alphabet,
) -> Callable[[Callable[..., float]], Callable[..., float]]:
"""Decorator for registering an actual causation measure with PyPhi.
Parameters
----------
name : string
The name of the measure.
asymmetric : boolean
``True`` if the measure is asymmetric.
supports_alphabet
Callable ``(alphabet_sizes) -> bool``
indicating whether the measure handles the given node
cardinalities. Defaults to :data:`_any_alphabet`.
"""
def register_func(func: Callable[..., float]) -> Callable[..., float]:
if asymmetric:
self._asymmetric.append(name)
func.supports_alphabet = supports_alphabet # type: ignore[attr-defined]
self.store[name] = func
return func
return register_func
[docs]
def asymmetric(self) -> list[str]:
"""Return a list of asymmetric measures."""
return self._asymmetric
actual_causation_measures = ActualCausationMeasureRegistry()
[docs]
class np_suppress(np.errstate, ContextDecorator):
"""Decorator to suppress NumPy warnings about divide-by-zero and
multiplication of ``NaN``.
.. note::
This should only be used in cases where you are *sure* that these
warnings are not indicative of deeper issues in your code.
"""
def __init__(self) -> None:
super().__init__(divide="ignore", invalid="ignore")
# Load precomputed hamming matrices.
_NUM_PRECOMPUTED_HAMMING_MATRICES = 10
_hamming_matrices = utils.load_data(
"hamming_matrices", _NUM_PRECOMPUTED_HAMMING_MATRICES
)
# TODO extend to nonbinary nodes
def _hamming_matrix(N: int) -> np.ndarray:
"""Return a matrix of Hamming distances for the possible states of ``N``
binary nodes.
Parameters
----------
N : int
The number of nodes under consideration.
Returns
-------
np.ndarray
A ``2**N`` by ``2**N`` matrix whose ``(i, j)`` element is the Hamming
distance between state ``i`` and state ``j``.
Examples
--------
>>> _hamming_matrix(2)
array([[0., 1., 1., 2.],
[1., 0., 2., 1.],
[1., 2., 0., 1.],
[2., 1., 1., 0.]])
"""
if N < _NUM_PRECOMPUTED_HAMMING_MATRICES:
return _hamming_matrices[N]
return _compute_hamming_matrix(N)
@joblib_memory.cache
def _compute_hamming_matrix(N: int) -> np.ndarray:
"""Compute and store a Hamming matrix for ``N`` nodes.
Hamming matrices have the following sizes::
N MBs
== ===
9 2
10 8
11 32
12 128
13 512
Because these matrices grow large and are needed infrequently, computed
matrices are stored in the joblib filesystem cache rather than added to
the in-memory ``_hamming_matrices`` global.
This function is only called when ``N`` exceeds
``_NUM_PRECOMPUTED_HAMMING_MATRICES``. Do not call it directly; use
``_hamming_matrix`` instead.
"""
possible_states = np.array(list(utils.all_states(N)))
return cdist(possible_states, possible_states, "hamming") * N
@joblib_memory.cache
def _kary_hamming_matrix(alphabet_sizes: tuple[int, ...]) -> np.ndarray:
"""Hamming ground-distance matrix over a (possibly non-binary) state space.
Generalizes ``_hamming_matrix`` to heterogeneous alphabets: the ``(i, j)``
entry is the number of nodes whose state differs between the ``i``th and
``j``th joint states, with states enumerated in the little-endian order used
by :func:`~pyphi.distribution.flatten` (so the matrix indices align with a
flattened repertoire).
"""
states = np.array(list(utils.all_states(alphabet_sizes)))
return cdist(states, states, "hamming") * len(alphabet_sizes)
# In-memory k-ary ground metrics, mirroring the precomputed binary matrices:
# matrices for small state spaces are kept in memory (a filesystem-cache load
# per EMD call is ~500x slower than the binary path's dict lookup); larger
# ones stay in the joblib filesystem cache only.
_kary_hamming_matrices: dict[tuple[int, ...], np.ndarray] = {}
_MAX_IN_MEMORY_GROUND_METRIC_STATES = 2**_NUM_PRECOMPUTED_HAMMING_MATRICES
def _ground_metric(alphabet_sizes: tuple[int, ...]) -> np.ndarray:
"""Return the Hamming ground-distance matrix for the given alphabet sizes.
Binary substrates (every alphabet size 2) use the precomputed and cached
``_hamming_matrix`` path; non-binary substrates use
``_kary_hamming_matrix``, memoized in memory for state spaces up to
``_MAX_IN_MEMORY_GROUND_METRIC_STATES`` states.
"""
if all(k == 2 for k in alphabet_sizes):
return _hamming_matrix(len(alphabet_sizes))
if math.prod(alphabet_sizes) < _MAX_IN_MEMORY_GROUND_METRIC_STATES:
matrix = _kary_hamming_matrices.get(alphabet_sizes)
if matrix is None:
matrix = _kary_hamming_matrices[alphabet_sizes] = _kary_hamming_matrix(
alphabet_sizes
)
return matrix
return _kary_hamming_matrix(alphabet_sizes)
[docs]
def hamming_emd(p: ArrayLike, q: ArrayLike) -> float:
"""Return the Earth Mover's Distance between two distributions (indexed
by state, one dimension per node) using the Hamming distance between states
as the transportation cost function.
Supports non-binary (k-ary) substrates: the ground metric is the number of
nodes whose state differs, taken over the substrate's actual (possibly
heterogeneous) state space. Singleton dimensions are squeezed out.
"""
p = np.asarray(p)
q = np.asarray(q)
alphabet_sizes = p.squeeze().shape
p_flat = flatten(p)
q_flat = flatten(q)
assert p_flat is not None
assert q_flat is not None
if p_flat.shape != q_flat.shape:
raise ValueError(
"hamming_emd requires distributions of equal shape; got "
f"{p_flat.shape} and {q_flat.shape}"
)
return EMD.compute(p_flat, q_flat, _ground_metric(alphabet_sizes))
[docs]
def effect_emd(p: ArrayLike, q: ArrayLike) -> float:
"""Compute the EMD between two effect repertoires.
Because the nodes are independent, the EMD between effect repertoires equals
the sum of the EMDs between each node's marginal distributions (this is exact
for product distributions under the additive Hamming ground metric). For a
node, that per-node EMD is the total variation between its marginals; for a
binary node this reduces to the absolute difference in the probabilities that
the node is OFF.
Parameters
----------
p : np.ndarray
The first repertoire.
q : np.ndarray
The second repertoire.
Returns
-------
float
The EMD between ``p`` and ``q``.
"""
p = np.asarray(p)
q = np.asarray(q)
total = 0.0
for i in range(p.ndim):
if p.shape[i] == 2:
total += abs(marginal_zero(p, i) - marginal_zero(q, i))
else:
# Total variation = 1/2 * L1 between the node's k-ary marginals.
total += 0.5 * np.abs(marginal(p, i) - marginal(q, i)).sum()
return float(total)
[docs]
@distribution_measures.register("EMD")
def emd(p: ArrayLike, q: ArrayLike, direction: Direction | None = None) -> float:
"""Compute the EMD between two repertoires for a given direction.
The full EMD computation (:func:`hamming_emd`) is used for cause
repertoires. A fast analytic solution (:func:`effect_emd`) is used for
effect repertoires, exploiting the independence of the nodes in an effect
repertoire.
Parameters
----------
p : np.ndarray
The first repertoire.
q : np.ndarray
The second repertoire.
direction : Direction or None
``CAUSE`` or ``EFFECT``. If ``EFFECT``, :func:`effect_emd` is used;
otherwise :func:`hamming_emd` is used. ``None`` is treated as
``CAUSE``.
Returns
-------
DistanceResult
The EMD between ``p`` and ``q``, carrying ``method="EMD"`` and the
``direction`` as metadata. The value is not rounded.
Raises
------
ValueError
If ``direction`` is invalid.
"""
if (direction == Direction.CAUSE) or (direction is None):
func = hamming_emd
elif direction == Direction.EFFECT:
func = effect_emd
else:
# TODO: test that ValueError is raised
validate.direction(direction)
raise ValueError(f"Invalid direction: {direction}")
return DistanceResult(func(p, q), method="EMD", direction=direction)
[docs]
@distribution_measures.register("L1")
def l1(p: ArrayLike, q: ArrayLike) -> float:
"""Return the L1 distance between two distributions.
Parameters
----------
p : np.ndarray
The first probability distribution.
q : np.ndarray
The second probability distribution.
Returns
-------
float
The sum of absolute differences of ``p`` and ``q``.
"""
p = np.asarray(p)
q = np.asarray(q)
return DistanceResult(np.abs(p - q).sum(), method="L1")
[docs]
@distribution_measures.register("ENTROPY_DIFFERENCE")
def entropy_difference(p: ArrayLike, q: ArrayLike) -> float:
"""Return the difference in entropy between two distributions."""
hp = entr(p).sum() / _LN_OF_2
hq = entr(q).sum() / _LN_OF_2
return DistanceResult(abs(hp - hq), method="ENTROPY_DIFFERENCE")
[docs]
@distribution_measures.register("PSQ2")
def psq2(p: ArrayLike, q: ArrayLike) -> float:
r"""Compute the PSQ2 measure.
This is defined as :math:`\mid f(p) - f(q) \mid`, where
.. math::
f(x) = \sum_{i=0}^{N-1} p_i^2 \log_2 (p_i N)
Parameters
----------
p : np.ndarray
The first distribution.
q : np.ndarray
The second distribution.
"""
p = np.asarray(p)
q = np.asarray(q)
fp = (p * (-1.0 * entr(p))).sum() / _LN_OF_2 + (p**2 * log2(p.size)).sum()
fq = (q * (-1.0 * entr(q))).sum() / _LN_OF_2 + (q**2 * log2(q.size)).sum()
return DistanceResult(abs(fp - fq), method="PSQ2")
[docs]
@distribution_measures.register("MP2Q", asymmetric=True)
@np_suppress()
def mp2q(p: ArrayLike, q: ArrayLike) -> float:
r"""Compute the MP2Q measure.
This is defined as
.. math::
\frac{1}{N}
\sum_{i=0}^{N-1} \frac{p_i^2}{q_i} \log_2\left(\frac{p_i}{q_i}\right)
Parameters
----------
p : np.ndarray
The first distribution.
q : np.ndarray
The second distribution.
Returns
-------
float
The distance.
"""
p = np.asarray(p, dtype=float)
q = np.asarray(q, dtype=float)
# There is already a factor of p in the `information_density`, so we only
# multiply by p, not p**2. Terms with p_i = 0 contribute nothing (the
# p_i^2 factor dominates the 1/q_i pole), so they are masked to zero
# rather than left as 0/0 = NaN; p_i > 0 with q_i = 0 stays +inf.
density = information_density(p, q)
ratio = np.divide(p, q, out=np.zeros_like(p), where=p > 0)
terms = np.where(p > 0, ratio * density, 0.0)
return DistanceResult(
np.sum(terms / p.size),
method="MP2Q",
asymmetric=True,
)
[docs]
@distribution_measures.register("KLD", asymmetric=True)
def kld(p: ArrayLike, q: ArrayLike) -> float:
"""Return the Kullback-Leibler Divergence (KLD) between two distributions.
Parameters
----------
p : np.ndarray
The first probability distribution.
q : np.ndarray
The second probability distribution.
Returns
-------
float
The KLD of ``p`` from ``q``.
"""
return DistanceResult(information_density(p, q).sum(), method="KLD", asymmetric=True)
[docs]
def approximate_specified_state(
repertoire: ArrayLike, partitioned_repertoire: ArrayLike
) -> np.ndarray:
"""Estimate the purview state that maximizes the AID between the repertoires.
This returns only the state of the purview nodes (i.e., there is one element
in the state vector for each purview node, not for each node in the
substrate).
.. note::
Although deterministic, *results are only a good guess*. This function
should only be used in cases where running :func:`specified_state`
becomes unfeasible.
This algorithm runs in linear time as a function of purview size, as opposed
to the exponential (on average) exhaustive exact search. Single-node
(i.e. marginal) repertoires are considered one by one, and their state is
determined according to the following heuristics:
If the most probable state in the unpartitioned repertoire (:math:`p > 1/2`)
becomes less probable in the partitioned one (:math:`p > q`), we should pick
that state for that node. Note that there can be ties. In that case, the
state with the lowest index is arbitrarily chosen.
Now suppose that was enough to specify the state of only :math:`k` nodes,
with joint point unpartitioned probability :math:`p_k` and partitioned
probability :math:`q_k`, and suppose we add node :math:`z`. Let the node
:math:`z` have probability :math:`p_z` for the state ``0``. For the
complementary state ``1``, the probability is :math:`1 - p_z`. We want to
know which state of :math:`z` gives higher intrinsic information when it is
added to the :math:`k` nodes. In other words, we want to compare :math:`I_x`
and :math:`I_y`:
.. math::
I_x = \\left( p_k p_z \\right) \\log_2 \\left( \\frac{p_k p_z}{q_k q_z} \\right)
.. math::
I_y = \\left( p_k (1-p_z) \\right) \\log_2 \\left(
\\frac{p_k (1-p_z)}{q_k(1-q_z)} \\right)
For state ``1`` to give higher intrinsic information (i.e., :math:`I_y >
I_x`), :math:`p_z` and :math:`q_z` must satisfy two equations:
.. math::
p_z < 1/2
.. math::
\\log_2 \\left( \\frac{p_k}{q_k} \\right) < \\left( \\frac{1}{1-2p_z}
\\right) \\left( p_z \\log_2 \\left( \\frac{p_z}{q_z} \\right) - (1-p_z)
\\log_2 \\left( \\frac{1-p_z}{1-q_z} \\right) \\right)
Otherwise, we should pick the state ``0`` as the state of node :math:`z`.
Parameters
----------
repertoire : np.ndarray
The first probability distribution.
partitioned_repertoire : np.ndarray
The second probability distribution.
Returns
-------
np.ndarray
A 2D array where the single row is the approximate
:func:`specified_state`.
"""
# TODO: All the marginalization defeats the whole purpose. Config option
# must prevent calculating outer product at `system`, and pass node
# marginal repertoires instead.
def joint_to_marginals(repertoire: np.ndarray) -> np.ndarray:
"""Converts a joint repertoire in multidimensional form to a 2D array of
single-node marginal repertoires.
Parameters
----------
repertoire : np.ndarray
The joint repertoire of a purview in
multidimensional form, e.g., as obtained from
:mod:`pyphi.system`. Note that `repertoire` is assumed to be
a well-formed probability distribution whose sum over all states
equals one.
Returns
-------
np.ndarray
A 2D array with one row per node in the purview (the
marginalized repertoires) and one column per state, in the same
order as the argument.
"""
# Remove singleton dimensions.
repertoire = repertoire.squeeze()
# Map each dimension in the squeezed repertoire to a local node index.
node_indices = set(range(repertoire.ndim))
# All the sets of indices of size n - 1 (i.e. combinations(n, n - 1)).
complements = [node_indices - {n} for n in tuple(node_indices)]
# Marginalize out all the complementary dimensions for each
# node in the repertoire.
marginals = [repertoire.sum(tuple(c)) for c in complements]
return np.vstack(marginals)
P = joint_to_marginals(np.asarray(repertoire))
Q = joint_to_marginals(np.asarray(partitioned_repertoire))
# Preallocate arrays for the specified states and their corresponding point
# probabilities in P and Q.
purview_size = P.shape[0]
specified_states = np.zeros((purview_size, 1)) * np.nan
specified_P = np.ones((purview_size, 1))
specified_Q = np.ones((purview_size, 1))
# Find "fixed" nodes. A fixed node is defined as one for which its most
# selective state according to the unpartitioned repertoire (p > 0.5) is
# also informative for the node itself, as revealed by partitioning the
# mechanism (i.e., p > q).
is_selective = P >= (1 / 2)
informativeness = P / Q
is_informative = informativeness >= 1
fixed_nodes = np.where(np.sum(is_selective * is_informative, axis=1))[0]
def informative_state(node):
return np.where(informativeness[node, :] == informativeness[node, :].max())[0]
for fixed_node in fixed_nodes:
specified_state = np.where(
is_selective[fixed_node, :] * is_informative[fixed_node, :]
)[0]
# TODO: state ties.
# If P[ON] == P[OFF] == Q[ON] == Q[OFF] then |specified_state| > 1.
# Arbitrarily pick the first state.
specified_state = specified_state[0]
specified_states[fixed_node] = specified_state
specified_P[fixed_node] = P[fixed_node, specified_state]
specified_Q[fixed_node] = Q[fixed_node, informative_state(fixed_node)[0]]
if fixed_nodes.size == purview_size:
return specified_states.astype(int).T
# Estimate the state of the remaining (i.e. non-fixed) nodes, one by one,
# based on a greedy search on their impact on "temporary informativeness".
nonfixed_nodes = np.setdiff1d(np.arange(purview_size), fixed_nodes)
# First, compute discriminant values for every non-fixed node. This
# discriminant will be compared to the temporary informativeness.
p = np.array([P[n, informative_state(n)] for n in nonfixed_nodes]).flatten()
q = np.array([Q[n, informative_state(n)] for n in nonfixed_nodes]).flatten()
discriminants = (p * np.log2(p / q) - (1 - p) * np.log2((1 - p) / (1 - q))) / (
1 - 2 * p
)
# The smaller the discriminant of a purview node, the more likely its true
# specified state is to violate p > q. Thus we consider nodes in that order.
discriminant_indices = np.argsort(discriminants)
discriminants = np.sort(discriminants) # ascending
for index, discriminant in zip(discriminant_indices, discriminants, strict=False):
# The temporary-informativeness, updated as new nodes are included.
tmp_inform = np.log2(specified_P.prod()) - np.log2(specified_Q.prod())
nonfixed_node = nonfixed_nodes[index]
# TODO: nonbinary states.
# If discriminant < tmp_inform, select the state that gives p < q.
# Exact comparison is intended: this is a documented approximation, so a
# sub-precision tie degrades an explicit guess, not an exact result.
if discriminant < tmp_inform:
specified_state = int(not informative_state(nonfixed_node)[0])
else:
specified_state = informative_state(nonfixed_node)[0]
specified_states[nonfixed_node] = specified_state
specified_P[nonfixed_node] = P[nonfixed_node, specified_state]
specified_Q[nonfixed_node] = Q[nonfixed_node, specified_state]
return specified_states.astype(int).T
[docs]
@distribution_measures.register("ID", asymmetric=True)
def intrinsic_difference(p: ArrayLike, q: ArrayLike) -> float:
r"""Compute the intrinsic difference (ID) between two distributions.
This is defined as
.. math::
\max_i \left\{
p_i \log_2 \left( \frac{p_i}{q_i} \right)
\right\}
where :math:`p_i \log_2 \left( \frac{p_i}{q_i} \right)` is :math:`0`
when :math:`p_i = 0`, and :math:`+\infty` when :math:`p_i > 0` and
:math:`q_i = 0`.
See the following paper:
Barbosa LS, Marshall W, Streipert S, Albantakis L, Tononi G (2020).
A measure for intrinsic information.
*Sci Rep*, 10, 18803. https://doi.org/10.1038/s41598-020-75943-4
Parameters
----------
p : np.ndarray
The first probability distribution.
q : np.ndarray
The second probability distribution.
Returns
-------
float
The intrinsic difference.
"""
return DistanceResult(
np.max(information_density(p, q)), method="ID", asymmetric=True
)
[docs]
@distribution_measures.register("AID", asymmetric=True)
@distribution_measures.register("KLM", asymmetric=True) # Backwards-compatible alias
@distribution_measures.register("BLD", asymmetric=True) # Backwards-compatible alias
def absolute_intrinsic_difference(p: ArrayLike, q: ArrayLike) -> float:
"""Compute the absolute intrinsic difference (AID) between two
distributions.
This is the same as the ID, but with the absolute value taken before the
maximum is taken.
See documentation for :func:`intrinsic_difference` for further details
and references.
Parameters
----------
p : float
The first probability distribution.
q : float
The second probability distribution.
Returns
-------
float
The absolute intrinsic difference.
"""
return DistanceResult(
np.max(absolute_information_density(p, q)), method="AID", asymmetric=True
)
[docs]
@stateful_distribution_measures.register("IIT_4.0_SMALL_PHI", asymmetric=True)
def iit_4_small_phi(p: ArrayLike, q: ArrayLike, state: State) -> float:
r"""Absolute intrinsic difference between two distributions at one state.
Returns :math:`\left| p_i \log_2 (p_i / q_i) \right|` evaluated at
``state`` (the pointwise absolute information density; see
:func:`absolute_information_density`). Unlike
:func:`absolute_intrinsic_difference`, which maximizes over states, this
reads off the value at a single specified purview state.
Parameters
----------
p : np.ndarray
The first probability distribution.
q : np.ndarray
The second probability distribution.
state : State
Index into the squeezed pointwise array selecting the value returned.
Returns
-------
DistanceResult
The absolute information density at ``state``, with
``method="IIT_4.0_SMALL_PHI"``.
"""
return DistanceResult(
absolute_information_density(p, q).squeeze()[state],
method="IIT_4.0_SMALL_PHI",
asymmetric=True,
state=state,
)
[docs]
@stateful_distribution_measures.register(
"IIT_4.0_SMALL_PHI_NO_ABSOLUTE_VALUE", asymmetric=True
)
def iit_4_small_phi_no_absolute_value(p: ArrayLike, q: ArrayLike, state: State) -> float:
r"""Signed intrinsic difference between two distributions at one state.
Like :func:`iit_4_small_phi`, but without the absolute value: returns
the signed pointwise information density
:math:`p_i \log_2 (p_i / q_i)` (see :func:`information_density`)
evaluated at ``state``. The result may be negative.
Parameters
----------
p : np.ndarray
The first probability distribution.
q : np.ndarray
The second probability distribution.
state : State
Index into the squeezed pointwise array selecting the value returned.
Returns
-------
DistanceResult
The signed information density at ``state``, with
``method="IIT_4.0_SMALL_PHI_NO_ABSOLUTE_VALUE"``.
"""
return DistanceResult(
information_density(p, q).squeeze()[state],
method="IIT_4.0_SMALL_PHI_NO_ABSOLUTE_VALUE",
asymmetric=True,
state=state,
)
[docs]
@composite_measures.register("GENERALIZED_INTRINSIC_DIFFERENCE", asymmetric=True)
@composite_measures.register("INTRINSIC_SPECIFICATION", asymmetric=True)
def generalized_intrinsic_difference(
forward_repertoire: ArrayLike,
partitioned_forward_repertoire: ArrayLike,
selectivity_repertoire: ArrayLike,
state: State | None = None,
) -> Repertoire | float:
r"""Generalized intrinsic difference (GID) between two repertoires.
Pointwise, this is the selectivity repertoire weighting the
pointwise mutual information between the unpartitioned (forward) and
partitioned repertoires:
.. math::
\text{gid}_i = s_i \, \log_2 \left( \frac{f_i}{f^{\text{part}}_i} \right)
where :math:`f` is ``forward_repertoire``, :math:`f^{\text{part}}` is
``partitioned_forward_repertoire``, and :math:`s` is
``selectivity_repertoire``. Terms where the ratio is undefined are set to
zero (see :func:`pointwise_mutual_information_vector`). This is the
default mechanism- and system-level φ measure, and is also registered
under the name ``INTRINSIC_SPECIFICATION``.
Parameters
----------
forward_repertoire : np.ndarray
The unpartitioned repertoire.
partitioned_forward_repertoire : np.ndarray
The repertoire under the partition.
selectivity_repertoire : np.ndarray
The per-state weighting applied to the pointwise mutual information.
state : State or None
If ``None``, the full pointwise GID array is returned. Otherwise the
single value at ``state`` is returned as a
:class:`DistanceResult`.
Returns
-------
np.ndarray or DistanceResult
The full pointwise GID array when ``state`` is ``None``; otherwise
the value at ``state`` with ``method="GENERALIZED_INTRINSIC_DIFFERENCE"``.
"""
selectivity_repertoire = np.asarray(selectivity_repertoire)
informativeness = pointwise_mutual_information_vector(
forward_repertoire, partitioned_forward_repertoire
)
gid = selectivity_repertoire * informativeness
if state is None:
return gid
return DistanceResult(
gid[state],
method="GENERALIZED_INTRINSIC_DIFFERENCE",
asymmetric=True,
state=state,
)
intrinsic_specification = generalized_intrinsic_difference # alias
[docs]
def pointwise_intrinsic_differentiation(p):
r"""Pointwise surprisal :math:`-\log_2 p_i` of a distribution.
Entries where :math:`p_i = 0` are set to ``0`` rather than ``+inf``.
Parameters
----------
p : np.ndarray
A probability distribution.
Returns
-------
np.ndarray
The elementwise surprisal, with zeros where ``p`` is zero.
"""
p = np.asarray(p, dtype=float)
out = np.zeros_like(p)
return -np.log2(p, out=out, where=(p > 0))
[docs]
@state_aware_measures.register("INTRINSIC_DIFFERENTIATION")
def intrinsic_differentiation(p, state):
r"""Smallest positive surprisal in a repertoire slice.
Selects ``p.squeeze()[state]``, takes its pointwise intrinsic
differentiation (the surprisal :math:`-\log_2 p_i`; see
:func:`pointwise_intrinsic_differentiation`), and returns the minimum
over the entries that are positive up to ``config.numerics.precision``.
If no entry is positive (every selected probability equals ``1`` up to
floating-point noise), the result is ``0.0``.
Parameters
----------
p : np.ndarray
A probability distribution.
state : State
Index into the squeezed distribution selecting the slice to score.
Returns
-------
DistanceResult
The minimum positive surprisal of the selected slice (of the whole
repertoire when ``state`` is ``None``), with
``method="INTRINSIC_DIFFERENTIATION"``.
"""
p = p.squeeze()
if state is not None:
p = p[state]
surprisal = pointwise_intrinsic_differentiation(p)
positive_entries = surprisal[numerics.positive_mask(surprisal)]
return DistanceResult(
np.min(positive_entries) if positive_entries.size > 0 else 0.0,
method="INTRINSIC_DIFFERENTIATION",
asymmetric=False,
state=state,
)
# ---------------------------------------------------------------------------
# Resolver helpers: name -> typed measure callable
# ---------------------------------------------------------------------------
[docs]
def resolve_mechanism_measure(
name: str,
alphabet_sizes: tuple[int, ...] | None = None,
) -> (
DistributionMeasure
| StateAwareMeasure
| StatefulDistributionMeasure
| CompositeMeasure
):
"""Look up a measure usable at the mechanism level.
Mechanism-level integration accepts distribution measures (IIT 3.0
EMD/L1/KLD/...), state-aware pointwise measures, stateful-distribution
measures (IIT 4.0 small-phi variants), or composite measures (GID at
the partition layer).
When ``alphabet_sizes`` is provided, the resolved measure's
``supports_alphabet`` predicate is evaluated and
:class:`NotImplementedError` is raised for any measure that declares
itself incompatible with those cardinalities.
"""
from typing import cast
if name in distribution_measures:
measure = cast(DistributionMeasure, distribution_measures[name])
elif name in state_aware_measures:
measure = cast(StateAwareMeasure, state_aware_measures[name])
elif name in stateful_distribution_measures:
measure = cast(StatefulDistributionMeasure, stateful_distribution_measures[name])
elif name in composite_measures:
measure = cast(CompositeMeasure, composite_measures[name])
else:
available = sorted(
set(distribution_measures)
| set(state_aware_measures)
| set(stateful_distribution_measures)
| set(composite_measures)
)
raise ValueError(f"Unknown mechanism measure {name!r}. Available: {available}")
if alphabet_sizes is not None:
supports = getattr(measure, "supports_alphabet", None)
if supports is not None and not supports(alphabet_sizes):
raise NotImplementedError(
f"Measure {name!r} does not support alphabet sizes "
f"{alphabet_sizes}. "
f"For multi-valued substrates, use an alphabet-generic "
f"measure (AID, GID, INTRINSIC_INFORMATION, "
f"GENERALIZED_INTRINSIC_DIFFERENCE). "
f"See Gomez et al. 2021 §2.3 "
f"(https://doi.org/10.3390/e23010006) for the theoretical "
f"rationale."
)
return measure
[docs]
def resolve_system_measure(name: str) -> CompositeMeasure:
"""Look up a measure usable at the system level.
Only composite measures are valid system-level measures; the return
type is :class:`CompositeMeasure`.
"""
from typing import cast
if name in composite_measures:
return cast(CompositeMeasure, composite_measures[name])
raise ValueError(
f"Unknown system measure {name!r}. Available: {sorted(composite_measures)}"
)
[docs]
def resolve_distribution_measure(name: str) -> DistributionMeasure:
"""Look up a distribution measure for EMD ground distance and IIT-side dispatch.
Only distribution measures (two-distribution distances) are valid
here; the return type is :class:`DistributionMeasure`.
"""
from typing import cast
if name in distribution_measures:
return cast(DistributionMeasure, distribution_measures[name])
raise ValueError(
f"Unknown distribution measure {name!r}. "
f"Available: {sorted(distribution_measures)}"
)
[docs]
def resolve_actual_causation_measure(name: str) -> DistributionMeasure:
"""Look up a measure registered in :data:`actual_causation_measures`.
The actual-causation alpha computation uses ``(p, q) -> float``
distribution-shape callables from :data:`actual_causation_measures`
(e.g., ``PMI``, ``WPMI``). The return type is :class:`DistributionMeasure`.
"""
from typing import cast
if name in actual_causation_measures:
return cast(DistributionMeasure, actual_causation_measures[name])
raise ValueError(
f"Unknown actual-causation measure {name!r}. "
f"Available: {sorted(actual_causation_measures)}"
)
[docs]
def repertoire_distance(
r1: ArrayLike,
r2: ArrayLike,
direction: Direction | None = None,
repertoire_distance: (
DistributionMeasure
| StateAwareMeasure
| StatefulDistributionMeasure
| CompositeMeasure
| None
) = None,
**kwargs,
) -> float:
"""Compute the distance between two repertoires for the given direction.
Dispatches on the measure's signature: a stateful measure ``(p, q, state)``
is given ``state`` from ``kwargs``; a distribution measure ``(p, q)`` is
given ``direction`` and any other keyword arguments its signature accepts.
Parameters
----------
r1 : np.ndarray
The first repertoire.
r2 : np.ndarray
The second repertoire.
direction : Direction
``CAUSE`` or ``EFFECT``.
repertoire_distance : DistributionMeasure or StateAwareMeasure or \
StatefulDistributionMeasure or CompositeMeasure
A Protocol-typed measure callable. Required for callers below the
formalism-class boundary; public-API callers (``System.cause_info``,
etc.) resolve from config at their method boundary and pass the object
through.
Returns
-------
float
The distance between ``r1`` and ``r2``, rounded to
``config.numerics.precision``.
Raises
------
ValueError
If ``repertoire_distance`` is ``None``.
TypeError
If a stateful measure is passed without a ``state`` keyword argument,
or the measure's signature matches no dispatchable Protocol.
"""
if repertoire_distance is None:
raise ValueError(
"repertoire_distance must be provided explicitly; callers below "
"the formalism boundary thread the measure object as a kwarg."
)
func = repertoire_distance
if satisfies_stateful_distribution_measure(func):
# (p, q, state) — caller threads ``state`` via kwargs. ``direction``
# is not part of the stateful measure signature.
try:
state = kwargs.pop("state")
except KeyError as exc:
raise TypeError(
f"StatefulDistributionMeasure "
f"{getattr(func, 'name', repr(func))!r} requires a 'state' "
f"keyword argument, but none was provided."
) from exc
distance = func(r1, r2, state) # type: ignore[call-arg]
elif satisfies_distribution_measure(func):
# (p, q) — may accept an optional ``direction`` parameter (e.g. EMD)
# and/or additional keyword arguments. Pass only those that appear
# in the measure's signature so unrecognized kwargs do not produce
# signature-mismatch TypeErrors masquerading as bugs.
sig_params = inspect.signature(func).parameters
call_kwargs: dict[str, Any] = {
name: value for name, value in kwargs.items() if name in sig_params
}
if "direction" in sig_params:
call_kwargs["direction"] = direction
distance = func(r1, r2, **call_kwargs) # type: ignore[call-arg]
else:
observed_params = list(inspect.signature(func).parameters)
raise TypeError(
f"Cannot dispatch repertoire_distance with measure "
f"{getattr(func, 'name', repr(func))!r}: signature parameters "
f"{observed_params} do not match any registered measure Protocol "
f"that takes two repertoires (DistributionMeasure or "
f"StatefulDistributionMeasure)."
)
return round(distance, config.numerics.precision) # type: ignore[arg-type]
