"""
Coordinate systems on the Gray-Wyner region.
The Gray-Wyner rate region, the *region of tension* of Prabhakaran &
Prabhakaran :cite:`prabhakaran2014assisted`, the *extension profile* of
Matveev & Romashchenko :cite:`matveev2026beyond`, and the *mutual information
region* of Li & El Gamal :cite:`li2018extended` are all images of one another
under invertible affine maps. Every one of them is generated by the same
auxiliary variable ``W``; they differ only in which affine functions of the
entropies of ``(X_0, ..., X_{n-1}, W)`` are recorded.
Writing ``R_0 = I[X_{0:n} : W]`` and ``R_i = H[X_i | W]`` for the Gray-Wyner
rates, the maps are
.. math::
\\tau_i &= R_0 + R_i - H[X_i] &&= I[X_{-i} : W | X_i] \\\\
\\tau_{res} &= R_0 + \\sum_i R_i - H[X_{0:n}] &&= T[X_{0:n} | W] \\\\
v_i &= H[X_i] - R_i &&= I[X_i : W] \\\\
v_{joint} &= R_0 &&= I[X_{0:n} : W]
The tension identities are the chain rule ``I[X_{0:n}:W] = I[X_i:W] +
I[X_{-i}:W|X_i]`` and the definition of the conditional total correlation
``T[X_{0:n}|W] = sum_i H[X_i|W] - H[X_{0:n}|W]``. The matrix carrying
``(R_0, ..., R_n)`` to ``(tau_0, ..., tau_{n-1}, tau_res)`` is unimodular, so
nothing is lost in either direction.
At ``n = 2`` the residual coordinate is ``I[X:Y|W]`` and these reduce to the
map ``g_{X,Y}`` of :cite:`prabhakaran2014assisted` (their Theorem 5.2).
Conventions
-----------
``dit`` follows the coordinate ordering of :cite:`prabhakaran2014assisted` and
:cite:`li2018extended`: the tension coordinate indexed by ``i`` is
``I[X_{-i} : W | X_i]``, so for a pair it is ``(I[Y:W|X], I[X:W|Y],
I[X:Y|W])``. Csirmaz :cite:`csirmaz2023short` swaps the first two coordinates
and omits the increasing hull (and correspondingly uses a cardinality bound of
``|X||Y| + 3`` rather than ``+ 2``). Formulas transported between these papers
must be checked against the convention in force.
The tension region proper is the *increasing hull* of the attainable tension
points -- every point dominating an attainable one is in the region. The
routines here report attainable points on the lower boundary; the hull is
implicit.
Besides the coordinate maps, this module holds the three quantities that
:cite:`li2018extended` identify as extreme points of the mutual information
region: the maximal interaction information and the two private interaction
informations. They are the reason `GrayWynerOptimizer` accepts signed
weights.
"""
from collections import namedtuple
import numpy as np
from ...exceptions import ditException
from ...utils import unitful
from .optimizer import GrayWynerPoint
__all__ = (
"ExtensionPoint",
"MutualInformationPoint",
"TensionPoint",
"asymmetric_private_interaction_information",
"extension_to_rates",
"maximal_interaction_information",
"mutual_information_to_rates",
"rates_to_extension",
"rates_to_mutual_information",
"rates_to_tension",
"symmetric_private_interaction_information",
"tension_to_rates",
)
TensionPoint = namedtuple("TensionPoint", ["tensions", "residual"])
TensionPoint.__doc__ = """
A point of the region of tension.
Attributes
----------
tensions : tuple of float
The per-source tensions ``tau_i = I[X_{-i} : W | X_i]``.
residual : float
The residual dependence ``tau_res = T[X_{0:n} | W]``, which for a pair is
the conditional mutual information ``I[X:Y|W]``.
"""
ExtensionPoint = namedtuple("ExtensionPoint", ["conditional_entropies", "residual"])
ExtensionPoint.__doc__ = """
A point of the extension profile.
Attributes
----------
conditional_entropies : tuple of float
The conditional entropies ``H[X_i | W]``.
residual : float
The residual dependence ``T[X_{0:n} | W]``.
"""
MutualInformationPoint = namedtuple("MutualInformationPoint", ["marginals", "joint"])
MutualInformationPoint.__doc__ = """
A point of the mutual information region of :cite:`li2018extended`.
Attributes
----------
marginals : tuple of float
The per-source informations ``I[X_i : W]``.
joint : float
The joint information ``I[X_{0:n} : W]``.
"""
def _check_lengths(values, marginal_entropies, what):
"""
Validate that a coordinate tuple matches the number of sources.
Parameters
----------
values : iterable of float
The per-source coordinates.
marginal_entropies : iterable of float
The marginal entropies ``H[X_i]``.
what : str
A name for the coordinate tuple, used in the error message.
Returns
-------
values, marginal_entropies : np.ndarray, np.ndarray
The validated inputs as arrays.
"""
values = np.asarray(values, dtype=float)
marginal_entropies = np.asarray(marginal_entropies, dtype=float)
if values.size != marginal_entropies.size:
msg = f"`{what}` has length {values.size} but there are {marginal_entropies.size} sources."
raise ditException(msg)
return values, marginal_entropies
[docs]
def rates_to_tension(point, marginal_entropies, joint_entropy):
"""
Convert a Gray-Wyner rate point to tension coordinates.
Parameters
----------
point : GrayWynerPoint
The rate point ``(R_0, (R_1, ..., R_n))``.
marginal_entropies : iterable of float
The marginal entropies ``H[X_i]``.
joint_entropy : float
The joint entropy ``H[X_{0:n}]``.
Returns
-------
tension : TensionPoint
The corresponding point of the region of tension.
"""
private, marginals = _check_lengths(point.private, marginal_entropies, "point.private")
tensions = point.common + private - marginals
residual = point.common + private.sum() - joint_entropy
return TensionPoint(tensions=tuple(tensions), residual=float(residual))
[docs]
def tension_to_rates(tension, marginal_entropies, joint_entropy):
"""
Convert tension coordinates back to a Gray-Wyner rate point.
Inverts :func:`rates_to_tension`. Summing the per-source tensions and
subtracting the residual eliminates the private rates and leaves
``(n - 1) R_0``, which is why at least two sources are required.
Parameters
----------
tension : TensionPoint
The point of the region of tension.
marginal_entropies : iterable of float
The marginal entropies ``H[X_i]``.
joint_entropy : float
The joint entropy ``H[X_{0:n}]``.
Returns
-------
point : GrayWynerPoint
The corresponding Gray-Wyner rate point.
"""
tensions, marginals = _check_lengths(tension.tensions, marginal_entropies, "tension.tensions")
n = tensions.size
if n < 2:
msg = "The region of tension requires at least two sources."
raise ditException(msg)
common = (tensions.sum() - tension.residual + marginals.sum() - joint_entropy) / (n - 1)
private = tensions + marginals - common
return GrayWynerPoint(common=float(common), private=tuple(private))
[docs]
def rates_to_extension(point, joint_entropy):
"""
Convert a Gray-Wyner rate point to extension-profile coordinates.
The private rates *are* the conditional entropies ``H[X_i | W]``, so this
map only computes the residual coordinate.
Parameters
----------
point : GrayWynerPoint
The rate point ``(R_0, (R_1, ..., R_n))``.
joint_entropy : float
The joint entropy ``H[X_{0:n}]``.
Returns
-------
extension : ExtensionPoint
The corresponding point of the extension profile.
"""
private = np.asarray(point.private, dtype=float)
residual = point.common + private.sum() - joint_entropy
return ExtensionPoint(conditional_entropies=tuple(private), residual=float(residual))
[docs]
def extension_to_rates(extension, joint_entropy):
"""
Convert extension-profile coordinates back to a Gray-Wyner rate point.
Inverts :func:`rates_to_extension`.
Parameters
----------
extension : ExtensionPoint
The point of the extension profile.
joint_entropy : float
The joint entropy ``H[X_{0:n}]``.
Returns
-------
point : GrayWynerPoint
The corresponding Gray-Wyner rate point.
"""
private = np.asarray(extension.conditional_entropies, dtype=float)
common = extension.residual + joint_entropy - private.sum()
return GrayWynerPoint(common=float(common), private=tuple(private))
def _pair_network(dist, rvs, crvs, bound):
"""
Build the Gray-Wyner network backing the interaction informations.
Parameters
----------
dist : Distribution
The distribution of interest.
rvs : list of lists, None
The two sources.
crvs : list, None
Variables to condition on.
bound : int, None
Optional cap on the cardinality of the probe.
Returns
-------
network : GrayWynerNetwork
The network, guaranteed to have exactly two sources.
"""
from .network import GrayWynerNetwork
network = GrayWynerNetwork(dist, rvs=rvs, crvs=crvs, bound=bound)
if network.n != 2:
msg = f"The interaction informations are defined for 2 sources, not {network.n}."
raise ditException(msg)
return network
def _no_path_of_length_three(network):
"""
Whether every outcome has ``p(x|y) = 1`` or ``p(y|x) = 1``.
This is exactly when the interaction informations vanish
:cite:`li2018extended`. Checked only without conditioning variables.
"""
if network.crvs:
return False
d = network.dist.coalesce(network.rvs)
d = d.copy(base="linear")
d.make_dense()
p = d.pmf.reshape([len(a) for a in d.alphabet])
with np.errstate(divide="ignore", invalid="ignore"):
x_given_y = p / p.sum(axis=0, keepdims=True)
y_given_x = p / p.sum(axis=1, keepdims=True)
determined = np.isclose(x_given_y, 1.0) | np.isclose(y_given_x, 1.0)
return bool(np.all(determined[p > 0]))