From 69061155ae2afac663011d9d4ef33e892f7df796 Mon Sep 17 00:00:00 2001 From: ParitKansal Date: Fri, 23 Oct 2026 08:16:28 +0530 Subject: [PATCH 1/3] Add symbolic (partition-based) entropy, MI and CMI estimators (#59) Implements the symbolic estimators of Cafaro, Lord, Sun and Bollt (Chaos, 2015): symbolize, symbolic_entropy, symbolic_mutual_information and symbolic_conditional_mutual_information, available as method='symbolic' in conditional_mutual_information and information='symbolic' in discover_network. Adds 24 tests, README and CHANGELOG entries. --- CHANGELOG.md | 1 + README.md | 2 +- causationentropy/core/discovery.py | 11 +- .../conditional_mutual_information.py | 96 ++++++++ causationentropy/core/information/entropy.py | 164 +++++++++++++ .../core/information/mutual_information.py | 69 +++++- .../tests/core/information/test_symbolic.py | 217 ++++++++++++++++++ 7 files changed, 557 insertions(+), 3 deletions(-) create mode 100644 causationentropy/tests/core/information/test_symbolic.py diff --git a/CHANGELOG.md b/CHANGELOG.md index fef357b..0942aca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,6 +23,7 @@ - `linear_gaussian_from_graph` for generating synthetic data from a chosen `networkx` graph (#39). - Block-bootstrap confidence intervals: `moving_block_bootstrap_indices`, `stationary_bootstrap_indices`, `bootstrap_confidence_interval` and `bootstrap_cmi_confidence_interval` (#40). - Optional futility early stopping in `shuffle_test` (`early_stop=True`). It is only available when calling `shuffle_test` directly (#41). +- Symbolic (partition-based) estimators following Cafaro, Lord, Sun and Bollt (Chaos, 2015): `symbolize`, `symbolic_entropy`, `symbolic_mutual_information` and `symbolic_conditional_mutual_information`, available as `method='symbolic'` in `conditional_mutual_information` and `information='symbolic'` in `discover_network` (#59). - Multiple-testing corrections `bonferroni_correction`, `benjamini_hochberg_correction`, `benjamini_yekutieli_correction`, `adaptive_bh_correction` and `estimate_null_proportion`, plus `apply_test_correction` to correct the p-values of a discovered network after discovery (adds a `P_Adjusted` column to `network_to_dataframe`). Discovery itself does not apply a correction (#37). - Tutorials and expanded theory documentation (#30, #49, #50). diff --git a/README.md b/README.md index 645caa9..5beec9b 100644 --- a/README.md +++ b/README.md @@ -95,7 +95,7 @@ from causationentropy import discover_network network = discover_network( data, method='standard', # 'standard', 'alternative', base lines: 'information_lasso', or 'lasso' - information='gaussian', # 'gaussian', 'knn', 'kde', 'geometric_knn', or 'poisson' + information='gaussian', # 'gaussian', 'knn', 'kde', 'geometric_knn', 'poisson', or 'symbolic' max_lag=5, # Maximum time lag to consider alpha_forward=0.05, # Forward selection significance alpha_backward=0.05, # Backward elimination significance diff --git a/causationentropy/core/discovery.py b/causationentropy/core/discovery.py index b658d7f..6a00ef4 100644 --- a/causationentropy/core/discovery.py +++ b/causationentropy/core/discovery.py @@ -84,6 +84,8 @@ def discover_network( - 'kde': Kernel density estimation - 'geometric_knn': Geometric mean k-NN estimator - 'poisson': Poisson distribution assumption + - 'symbolic': Plug-in estimate on symbolized data (Cafaro et al., 2015); + each variable is split at its median (two symbols) max_lag : int, default=5 Maximum time lag to consider in causal relationships. The algorithm examines lags from 1 to max_lag (inclusive). @@ -176,7 +178,14 @@ def discover_network( if method not in ["standard", "alternative", "information_lasso", "lasso"]: raise NotImplementedError(f"discover_network: method={method} not supported.") - supported_information_types = ["gaussian", "knn", "kde", "geometric_knn", "poisson"] + supported_information_types = [ + "gaussian", + "knn", + "kde", + "geometric_knn", + "poisson", + "symbolic", + ] if information not in supported_information_types: raise NotImplementedError( f"discover_network: information={information} not supported. " diff --git a/causationentropy/core/information/conditional_mutual_information.py b/causationentropy/core/information/conditional_mutual_information.py index b910271..5933ee4 100644 --- a/causationentropy/core/information/conditional_mutual_information.py +++ b/causationentropy/core/information/conditional_mutual_information.py @@ -1,3 +1,5 @@ +from typing import Optional + import numpy as np from scipy.spatial.distance import cdist from scipy.special import digamma @@ -7,12 +9,15 @@ kde_entropy, poisson_entropy, poisson_joint_entropy, + symbolic_entropy, + symbolize, ) from causationentropy.core.information.mutual_information import ( gaussian_mutual_information, geometric_knn_mutual_information, kde_mutual_information, knn_mutual_information, + symbolic_mutual_information, ) @@ -371,6 +376,84 @@ def poisson_conditional_mutual_information(X, Y, Z): return cmi +def symbolic_conditional_mutual_information( + X: np.ndarray, + Y: np.ndarray, + Z: Optional[np.ndarray] = None, + n_symbols: int = 2, + partition: str = "quantile", +) -> float: + r""" + Plug-in conditional mutual information of symbolized data, in nats. + + Every variable is symbolized with + :func:`~causationentropy.core.information.entropy.symbolize` (a partition of each + variable's range, as in Cafaro et al. [1]_), and the conditional mutual information + of Eq. (15) in [1]_ is estimated from empirical symbol frequencies: + + .. math:: + + \hat I(X; Y | Z) = \hat H(X, Z) + \hat H(Y, Z) - \hat H(Z) - \hat H(X, Y, Z). + + Parameters + ---------- + X : np.ndarray, shape (n_samples,) or (n_samples, n_features_x) + First variable. + Y : np.ndarray, shape (n_samples,) or (n_samples, n_features_y) + Second variable. + Z : np.ndarray, shape (n_samples,) or (n_samples, n_features_z), or None + Conditioning variable. If None, returns the mutual information I(X; Y). + n_symbols : int, default=2 + Number of symbols per variable. + partition : {'quantile', 'uniform', 'precomputed'}, default='quantile' + Partition type. + + Returns + ------- + float + Conditional mutual information estimate in nats. + + Notes + ----- + Cafaro et al. [1]_ show that the inferred causal structure can depend sensitively + on the partition, so ``n_symbols`` and ``partition`` should be chosen with care. + The plug-in estimate is biased upward when the number of occupied joint symbols is + not small compared with the sample size [2]_; in :func:`discover_network` the + permutation test uses the same estimator for its null distribution. + Time complexity is :math:`O(n \log n)` per call. + + Examples + -------- + >>> import numpy as np + >>> rng = np.random.default_rng(0) + >>> z = rng.integers(0, 2, 1000) + >>> x = z.copy(); y = z.copy() + >>> round(symbolic_conditional_mutual_information(x, y, z, partition="precomputed"), 6) + 0.0 + + References + ---------- + .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from + symbolic representations of dynamical systems," Chaos 25, 043106 (2015). + arXiv:1507.07262. + .. [2] L. Paninski, "Estimation of entropy and mutual information," Neural + Computation 15(6), 1191-1253 (2003). + """ + if Z is None: + return symbolic_mutual_information( + X, Y, n_symbols=n_symbols, partition=partition + ) + + SX = symbolize(X, n_symbols=n_symbols, partition=partition) + SY = symbolize(Y, n_symbols=n_symbols, partition=partition) + SZ = symbolize(Z, n_symbols=n_symbols, partition=partition) + + def H(*parts): + return symbolic_entropy(np.hstack(parts), partition="precomputed") + + return H(SX, SZ) + H(SY, SZ) - H(SZ) - H(SX, SY, SZ) + + def conditional_mutual_information( X, Y, @@ -380,6 +463,8 @@ def conditional_mutual_information( k=6, bandwidth="silverman", kernel="gaussian", + n_symbols: int = 2, + partition: str = "quantile", ): """ Compute conditional mutual information using specified estimation method. @@ -417,6 +502,7 @@ def conditional_mutual_information( - 'knn': k-nearest neighbor (KSG) estimator - 'geometric_knn': Geometric k-NN with manifold corrections - 'poisson': For discrete count data with Poisson assumptions + - 'symbolic': Plug-in estimate on symbolized data (Cafaro et al., 2015) metric : str or None, default=None Distance metric for k-NN based methods. @@ -426,6 +512,10 @@ def conditional_mutual_information( Bandwidth parameter for KDE methods. kernel : str, default='gaussian' Kernel function for KDE methods. + n_symbols : int, default=2 + Number of symbols per variable for the 'symbolic' method. + partition : {'quantile', 'uniform', 'precomputed'}, default='quantile' + Partition type for the 'symbolic' method. Returns ------- @@ -495,6 +585,11 @@ def conditional_mutual_information( elif method == "poisson": cmi = poisson_conditional_mutual_information(X, Y, Z) + elif method == "symbolic": + cmi = symbolic_conditional_mutual_information( + X, Y, Z, n_symbols=n_symbols, partition=partition + ) + else: supported_methods = [ "gaussian", @@ -503,6 +598,7 @@ def conditional_mutual_information( "knn", "geometric_knn", "poisson", + "symbolic", ] raise ValueError( f"Method '{method}' unavailable. Supported methods: {supported_methods}" diff --git a/causationentropy/core/information/entropy.py b/causationentropy/core/information/entropy.py index 7fc1790..14e7b39 100644 --- a/causationentropy/core/information/entropy.py +++ b/causationentropy/core/information/entropy.py @@ -344,3 +344,167 @@ def poisson_joint_entropy(Cov): Ent1 = np.sum(poisson_entropy(U)) Ent2 = np.sum(T) return Ent1 + Ent2 + + +def symbolize( + X: np.ndarray, + n_symbols: int = 2, + partition: str = "quantile", +) -> np.ndarray: + r""" + Convert continuous data into symbols by partitioning each variable's range. + + Following Cafaro et al. [1]_, a partition :math:`\mathcal{A} = \{A_0, \dots, A_m\}` + of the state space assigns to each observation :math:`x_t` the symbol + + .. math:: + + s_t = \sum_{i=0}^{m} \chi_{A_i}(x_t) \cdot i, + + where :math:`\chi_{A}` is the indicator function of the cell :math:`A`. Here the + partition is applied to each column (variable) separately, with ``n_symbols`` cells. + + Parameters + ---------- + X : np.ndarray, shape (n_samples,) or (n_samples, n_features) + Input data. A 1-D array is treated as a single column. + n_symbols : int, default=2 + Number of cells (symbols) per variable. Must be at least 2. Ignored when + ``partition='precomputed'``. + partition : {'quantile', 'uniform', 'precomputed'}, default='quantile' + How the cells are chosen for each column: + + - ``'quantile'``: approximately equal-frequency cells. Edges are placed at + midpoints between neighbouring distinct values, as close as possible to the + empirical quantiles, so tied values always share a cell (with + ``n_symbols=2`` this is a split near the median). Data with few distinct + values may therefore use fewer than ``n_symbols`` symbols. + - ``'uniform'``: equal-width cells between the column minimum and maximum. + - ``'precomputed'``: the data are already symbols (integer labels) and are + returned unchanged as integers. + + For ``'uniform'``, a value equal to a cell edge is assigned to the upper cell. + + Returns + ------- + np.ndarray, shape (n_samples, n_features) + Integer symbols in ``{0, ..., n_symbols - 1}`` (or the given labels for + ``'precomputed'``). + + Raises + ------ + ValueError + If ``n_symbols < 2``, the partition type is unknown, or precomputed data are + not integer-valued. + + Examples + -------- + >>> import numpy as np + >>> symbolize(np.array([0.1, 0.9, 0.4, 0.6]), n_symbols=2).ravel() + array([0, 1, 0, 1]) + + References + ---------- + .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from + symbolic representations of dynamical systems," Chaos 25, 043106 (2015). + arXiv:1507.07262. + """ + X = np.asarray(X) + if X.ndim == 1: + X = X.reshape(-1, 1) + + if partition == "precomputed": + if not np.all(np.equal(np.mod(X, 1), 0)): + raise ValueError("partition='precomputed' requires integer-valued symbols.") + return X.astype(np.int64) + + if n_symbols < 2: + raise ValueError(f"n_symbols must be at least 2, got {n_symbols}.") + + X = X.astype(float) + symbols = np.empty(X.shape, dtype=np.int64) + for j in range(X.shape[1]): + col = X[:, j] + if partition == "quantile": + inner_edges = _quantile_edges(col, n_symbols) + elif partition == "uniform": + inner_edges = np.linspace(col.min(), col.max(), n_symbols + 1)[1:-1] + else: + raise ValueError( + f"Unknown partition '{partition}'. " + "Use 'quantile', 'uniform' or 'precomputed'." + ) + symbols[:, j] = np.searchsorted(inner_edges, col, side="right") + return symbols + + +def _quantile_edges(col: np.ndarray, n_symbols: int) -> np.ndarray: + """Equal-frequency cell edges placed between distinct values of ``col``.""" + values, counts = np.unique(col, return_counts=True) + if values.size < 2: + return np.empty(0) + midpoints = (values[:-1] + values[1:]) / 2.0 + share_below = np.cumsum(counts)[:-1] / col.size + targets = np.arange(1, n_symbols) / n_symbols + nearest = np.abs(share_below[None, :] - targets[:, None]).argmin(axis=1) + return np.unique(midpoints[nearest]) + + +def _joint_labels(S: np.ndarray) -> np.ndarray: + """Map each row of a 2-D symbol array to one integer label for its joint symbol.""" + _, labels = np.unique(S, axis=0, return_inverse=True) + return labels.ravel() + + +def symbolic_entropy( + X: np.ndarray, + n_symbols: int = 2, + partition: str = "quantile", +) -> float: + r""" + Plug-in (maximum-likelihood) Shannon entropy of symbolized data, in nats. + + The data are symbolized with :func:`symbolize`, and the joint symbol of each row is + counted. With empirical frequencies :math:`\hat P(s)`, + + .. math:: + + \hat H(X) = -\sum_{s} \hat P(s) \log \hat P(s). + + Parameters + ---------- + X : np.ndarray, shape (n_samples,) or (n_samples, n_features) + Input data; several columns are treated as one joint variable. + n_symbols : int, default=2 + Number of symbols per variable (see :func:`symbolize`). + partition : {'quantile', 'uniform', 'precomputed'}, default='quantile' + Partition type (see :func:`symbolize`). + + Returns + ------- + float + Entropy estimate in nats. + + Notes + ----- + The plug-in estimator is biased downward for small samples relative to the number + of occupied joint symbols [2]_. Time complexity is :math:`O(n \log n)`. + + Examples + -------- + >>> import numpy as np + >>> x = np.array([0, 1, 2, 3]) + >>> round(symbolic_entropy(x, partition="precomputed"), 6) == round(np.log(4), 6) + True + + References + ---------- + .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from + symbolic representations of dynamical systems," Chaos 25, 043106 (2015). + .. [2] L. Paninski, "Estimation of entropy and mutual information," Neural + Computation 15(6), 1191-1253 (2003). + """ + labels = _joint_labels(symbolize(X, n_symbols=n_symbols, partition=partition)) + counts = np.bincount(labels) + p = counts[counts > 0] / labels.size + return float(-np.sum(p * np.log(p))) + 0.0 # + 0.0 turns -0.0 into 0.0 diff --git a/causationentropy/core/information/mutual_information.py b/causationentropy/core/information/mutual_information.py index 3131c85..fb0dc62 100644 --- a/causationentropy/core/information/mutual_information.py +++ b/causationentropy/core/information/mutual_information.py @@ -4,7 +4,12 @@ from scipy.spatial.distance import cdist from scipy.special import digamma -from causationentropy.core.information.entropy import geometric_knn_entropy, kde_entropy +from causationentropy.core.information.entropy import ( + geometric_knn_entropy, + kde_entropy, + symbolic_entropy, + symbolize, +) from causationentropy.core.linalg import correlation_log_determinant @@ -240,3 +245,65 @@ def geometric_knn_mutual_information(X, Y, metric="euclidean", k=1): # Ensure non-negativity return max(0.0, mi) + + +def symbolic_mutual_information( + X: np.ndarray, + Y: np.ndarray, + n_symbols: int = 2, + partition: str = "quantile", +) -> float: + r""" + Plug-in mutual information of symbolized data, in nats. + + Each variable is symbolized with + :func:`~causationentropy.core.information.entropy.symbolize`, and + + .. math:: + + \hat I(X; Y) = \hat H(X) + \hat H(Y) - \hat H(X, Y), + + which equals Eq. (14) of Cafaro et al. [1]_ with empirical probabilities. + + Parameters + ---------- + X : np.ndarray, shape (n_samples,) or (n_samples, n_features_x) + First variable. + Y : np.ndarray, shape (n_samples,) or (n_samples, n_features_y) + Second variable. + n_symbols : int, default=2 + Number of symbols per variable. + partition : {'quantile', 'uniform', 'precomputed'}, default='quantile' + Partition type. + + Returns + ------- + float + Mutual information estimate in nats (non-negative up to rounding). + + Notes + ----- + The plug-in estimate is biased upward when the number of occupied joint symbols + is not small compared with the sample size [2]_. + + Examples + -------- + >>> import numpy as np + >>> x = np.array([0, 1, 0, 1]) + >>> round(symbolic_mutual_information(x, x, partition="precomputed"), 6) + 0.693147 + + References + ---------- + .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from + symbolic representations of dynamical systems," Chaos 25, 043106 (2015). + .. [2] L. Paninski, "Estimation of entropy and mutual information," Neural + Computation 15(6), 1191-1253 (2003). + """ + SX = symbolize(X, n_symbols=n_symbols, partition=partition) + SY = symbolize(Y, n_symbols=n_symbols, partition=partition) + return ( + symbolic_entropy(SX, partition="precomputed") + + symbolic_entropy(SY, partition="precomputed") + - symbolic_entropy(np.hstack([SX, SY]), partition="precomputed") + ) diff --git a/causationentropy/tests/core/information/test_symbolic.py b/causationentropy/tests/core/information/test_symbolic.py new file mode 100644 index 0000000..05c4727 --- /dev/null +++ b/causationentropy/tests/core/information/test_symbolic.py @@ -0,0 +1,217 @@ +"""Tests for the symbolic (partition-based) information estimators (issue #59).""" + +import networkx as nx +import numpy as np +import pytest + +from causationentropy.core.discovery import discover_network +from causationentropy.core.information.conditional_mutual_information import ( + conditional_mutual_information, + symbolic_conditional_mutual_information, +) +from causationentropy.core.information.entropy import symbolic_entropy, symbolize +from causationentropy.core.information.mutual_information import ( + symbolic_mutual_information, +) + + +def _cmi_by_definition(x, y, z): + """I(X;Y|Z) from Eq. (15) of Cafaro et al. (2015), summed term by term (nats).""" + n = len(x) + total = 0.0 + for xv, yv, zv in set(zip(x, y, z)): + p_xyz = np.sum((x == xv) & (y == yv) & (z == zv)) / n + p_z = np.sum(z == zv) / n + p_xz = np.sum((x == xv) & (z == zv)) / n + p_yz = np.sum((y == yv) & (z == zv)) / n + total += p_xyz * np.log((p_xyz / p_z) / ((p_xz / p_z) * (p_yz / p_z))) + return total + + +class TestSymbolize: + def test_median_split(self): + x = np.array([0.1, 0.9, 0.4, 0.6]) + assert symbolize(x, n_symbols=2).ravel().tolist() == [0, 1, 0, 1] + + def test_quantile_cells_are_equal_frequency(self): + x = np.random.default_rng(0).normal(size=1000) + counts = np.bincount(symbolize(x, n_symbols=4).ravel()) + assert counts.tolist() == [250, 250, 250, 250] + + def test_uniform_cells(self): + x = np.array([0.0, 0.2, 0.5, 0.8, 1.0]) + assert symbolize(x, n_symbols=2, partition="uniform").ravel().tolist() == [ + 0, + 0, + 1, + 1, + 1, + ] + + def test_edge_value_goes_to_upper_cell(self): + x = np.array([0.0, 0.5, 1.0]) + assert symbolize(x, n_symbols=2, partition="uniform").ravel().tolist() == [ + 0, + 1, + 1, + ] + + def test_columns_are_symbolized_separately(self): + X = np.column_stack([[1.0, 2.0, 3.0, 4.0], [400.0, 300.0, 200.0, 100.0]]) + S = symbolize(X, n_symbols=2) + assert S.shape == (4, 2) + assert S[:, 0].tolist() == [0, 0, 1, 1] + assert S[:, 1].tolist() == [1, 1, 0, 0] + + def test_precomputed_returns_integers(self): + S = symbolize(np.array([2.0, 0.0, 1.0]), partition="precomputed") + assert S.dtype == np.int64 + assert S.ravel().tolist() == [2, 0, 1] + + def test_invalid_inputs(self): + with pytest.raises(ValueError): + symbolize(np.arange(5.0), n_symbols=1) + with pytest.raises(ValueError): + symbolize(np.arange(5.0), partition="unknown") + with pytest.raises(ValueError): + symbolize(np.array([0.5, 1.0]), partition="precomputed") + + +class TestSymbolicEntropy: + def test_uniform_symbols(self): + x = np.array([0, 1, 2, 3]) + assert np.isclose(symbolic_entropy(x, partition="precomputed"), np.log(4)) + + def test_constant_has_zero_entropy(self): + assert symbolic_entropy(np.zeros(10), partition="precomputed") == 0.0 + + def test_joint_entropy_of_two_columns(self): + X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]]) + assert np.isclose(symbolic_entropy(X, partition="precomputed"), np.log(4)) + + +class TestSymbolicMutualInformation: + def test_self_information_equals_entropy(self): + x = np.random.default_rng(1).normal(size=500) + assert np.isclose( + symbolic_mutual_information(x, x, n_symbols=3), + symbolic_entropy(x, n_symbols=3), + ) + + def test_independent_variables_near_zero(self): + rng = np.random.default_rng(2) + mi = symbolic_mutual_information(rng.normal(size=20000), rng.normal(size=20000)) + assert 0.0 <= mi < 1e-3 + + def test_invariant_to_monotone_transform(self): + rng = np.random.default_rng(3) + x = rng.normal(size=1000) + y = x + rng.normal(size=1000) + assert np.isclose( + symbolic_mutual_information(x, y, n_symbols=3), + symbolic_mutual_information(np.exp(x), y**3, n_symbols=3), + ) + + +class TestSymbolicConditionalMutualInformation: + def test_matches_definition(self): + rng = np.random.default_rng(4) + z = rng.integers(0, 3, 400) + x = (z + rng.integers(0, 2, 400)) % 3 + y = (x + rng.integers(0, 2, 400)) % 2 + got = symbolic_conditional_mutual_information(x, y, z, partition="precomputed") + assert np.isclose(got, _cmi_by_definition(x, y, z)) + + def test_chain_is_screened_off(self): + # X -> Z -> Y with Y a function of Z: I(X;Y) > 0 but I(X;Y|Z) = 0. + rng = np.random.default_rng(5) + x = rng.integers(0, 2, 2000) + z = np.where(rng.random(2000) < 0.9, x, 1 - x) + y = z.copy() + assert symbolic_mutual_information(x, y, partition="precomputed") > 0.2 + assert np.isclose( + symbolic_conditional_mutual_information(x, y, z, partition="precomputed"), + 0.0, + ) + + def test_none_z_returns_mutual_information(self): + rng = np.random.default_rng(6) + x = rng.normal(size=300) + y = x + rng.normal(size=300) + assert symbolic_conditional_mutual_information(x, y, None) == pytest.approx( + symbolic_mutual_information(x, y) + ) + + def test_logistic_map_generating_partition_is_memoryless(self): + # The logistic map at r = 4 is conjugate to the tent map; the partition at + # 1/2 is generating and its symbol sequence is a fair coin (Cafaro et al., + # 2015, Sec. IV), so one past symbol carries no information about the next. + x = np.empty(20001) + x[0] = 0.1234 + for t in range(20000): + x[t + 1] = 4.0 * x[t] * (1.0 - x[t]) + s = symbolize(x, n_symbols=2, partition="uniform").ravel() + assert symbolic_entropy(s, partition="precomputed") == pytest.approx( + np.log(2), abs=1e-3 + ) + assert ( + symbolic_mutual_information(s[:-1], s[1:], partition="precomputed") < 1e-3 + ) + + def test_dispatcher(self): + rng = np.random.default_rng(7) + x, y, z = rng.normal(size=(3, 500)) + direct = symbolic_conditional_mutual_information(x, y, z, n_symbols=3) + via = conditional_mutual_information(x, y, z, method="symbolic", n_symbols=3) + assert via == pytest.approx(max(0.0, direct)) + + def test_dispatcher_lists_symbolic_on_error(self): + with pytest.raises(ValueError, match="symbolic"): + conditional_mutual_information( + np.zeros(5), np.zeros(5), None, method="nonexistent" + ) + + +class TestSymbolicDiscovery: + def test_recovers_lagged_edge(self): + rng = np.random.default_rng(42) + T = 600 + x0 = rng.normal(size=T) + x1 = np.zeros(T) + x1[1:] = 0.9 * x0[:-1] + 0.3 * rng.normal(size=T - 1) + G = discover_network( + np.column_stack([x0, x1]), + information="symbolic", + max_lag=2, + n_shuffles=100, + n_jobs=1, + ) + assert isinstance(G, nx.MultiDiGraph) + edges = {(u, v, a.get("lag")) for u, v, a in G.edges(data=True)} + assert edges == {("X0", "X1", 1)} + + +class TestSymbolizeTies: + def test_heavily_tied_binary_data_keeps_both_symbols(self): + # 70% zeros: a median edge at 0 would put every point in one cell. + x = np.array([0.0] * 70 + [1.0] * 30) + assert symbolize(x, n_symbols=2).ravel().tolist() == [0] * 70 + [1] * 30 + assert symbolic_entropy(x) == pytest.approx( + -(0.7 * np.log(0.7) + 0.3 * np.log(0.3)) + ) + + def test_tied_values_share_a_cell(self): + x = np.array([0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3, 3], dtype=float) + s = symbolize(x, n_symbols=3).ravel() + for v in np.unique(x): + assert len(set(s[x == v])) == 1 + + def test_constant_column(self): + assert symbolize(np.ones(10), n_symbols=3).ravel().tolist() == [0] * 10 + assert symbolic_entropy(np.ones(10)) == 0.0 + + def test_count_data_dependence_is_detected(self): + rng = np.random.default_rng(8) + x = rng.poisson(0.5, 3000).astype(float) + y = x + rng.poisson(0.2, 3000) + assert symbolic_mutual_information(x, y, n_symbols=3) > 0.1 From 8032922f10fefd6f616f55866a160bdcd84e3b4e Mon Sep 17 00:00:00 2001 From: ParitKansal Date: Fri, 23 Oct 2026 08:37:03 +0530 Subject: [PATCH 2/3] Reject NaN/inf in symbolize, fix symbolic_entropy doctest on NumPy 2, clarify median wording --- causationentropy/core/discovery.py | 2 +- causationentropy/core/information/entropy.py | 10 ++++++---- .../tests/core/information/test_symbolic.py | 6 ++++++ 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/causationentropy/core/discovery.py b/causationentropy/core/discovery.py index 6a00ef4..eec27d9 100644 --- a/causationentropy/core/discovery.py +++ b/causationentropy/core/discovery.py @@ -85,7 +85,7 @@ def discover_network( - 'geometric_knn': Geometric mean k-NN estimator - 'poisson': Poisson distribution assumption - 'symbolic': Plug-in estimate on symbolized data (Cafaro et al., 2015); - each variable is split at its median (two symbols) + each variable is split near its median (two symbols) max_lag : int, default=5 Maximum time lag to consider in causal relationships. The algorithm examines lags from 1 to max_lag (inclusive). diff --git a/causationentropy/core/information/entropy.py b/causationentropy/core/information/entropy.py index 14e7b39..baf4fc7 100644 --- a/causationentropy/core/information/entropy.py +++ b/causationentropy/core/information/entropy.py @@ -394,8 +394,8 @@ def symbolize( Raises ------ ValueError - If ``n_symbols < 2``, the partition type is unknown, or precomputed data are - not integer-valued. + If ``n_symbols < 2``, the partition type is unknown, the data contain NaN or + infinite values, or precomputed data are not integer-valued. Examples -------- @@ -422,6 +422,8 @@ def symbolize( raise ValueError(f"n_symbols must be at least 2, got {n_symbols}.") X = X.astype(float) + if not np.all(np.isfinite(X)): + raise ValueError("symbolize requires finite values (no NaN or inf).") symbols = np.empty(X.shape, dtype=np.int64) for j in range(X.shape[1]): col = X[:, j] @@ -494,8 +496,8 @@ def symbolic_entropy( -------- >>> import numpy as np >>> x = np.array([0, 1, 2, 3]) - >>> round(symbolic_entropy(x, partition="precomputed"), 6) == round(np.log(4), 6) - True + >>> round(symbolic_entropy(x, partition="precomputed"), 6) # log(4) + 1.386294 References ---------- diff --git a/causationentropy/tests/core/information/test_symbolic.py b/causationentropy/tests/core/information/test_symbolic.py index 05c4727..4fc4050 100644 --- a/causationentropy/tests/core/information/test_symbolic.py +++ b/causationentropy/tests/core/information/test_symbolic.py @@ -76,6 +76,12 @@ def test_invalid_inputs(self): with pytest.raises(ValueError): symbolize(np.array([0.5, 1.0]), partition="precomputed") + @pytest.mark.parametrize("partition", ["quantile", "uniform", "precomputed"]) + @pytest.mark.parametrize("bad", [np.nan, np.inf]) + def test_non_finite_values_are_rejected(self, partition, bad): + with pytest.raises(ValueError): + symbolize(np.array([1.0, bad, 2.0, 3.0]), partition=partition) + class TestSymbolicEntropy: def test_uniform_symbols(self): From 9d8bf546d81a38d00a43ec68788b4ea258ef6413 Mon Sep 17 00:00:00 2001 From: ParitKansal Date: Fri, 23 Oct 2026 09:10:08 +0530 Subject: [PATCH 3/3] Check for NaN/inf before the precomputed branch in symbolize With partition='precomputed', NaN or inf input raised the misleading 'requires integer-valued symbols' error. The finite check now runs first, so every partition type gives the same error. The test now checks the message. --- causationentropy/core/information/entropy.py | 6 +++--- causationentropy/tests/core/information/test_symbolic.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/causationentropy/core/information/entropy.py b/causationentropy/core/information/entropy.py index baf4fc7..eefae73 100644 --- a/causationentropy/core/information/entropy.py +++ b/causationentropy/core/information/entropy.py @@ -412,6 +412,9 @@ def symbolize( X = np.asarray(X) if X.ndim == 1: X = X.reshape(-1, 1) + X = X.astype(float) + if not np.all(np.isfinite(X)): + raise ValueError("symbolize requires finite values (no NaN or inf).") if partition == "precomputed": if not np.all(np.equal(np.mod(X, 1), 0)): @@ -421,9 +424,6 @@ def symbolize( if n_symbols < 2: raise ValueError(f"n_symbols must be at least 2, got {n_symbols}.") - X = X.astype(float) - if not np.all(np.isfinite(X)): - raise ValueError("symbolize requires finite values (no NaN or inf).") symbols = np.empty(X.shape, dtype=np.int64) for j in range(X.shape[1]): col = X[:, j] diff --git a/causationentropy/tests/core/information/test_symbolic.py b/causationentropy/tests/core/information/test_symbolic.py index 4fc4050..65e36df 100644 --- a/causationentropy/tests/core/information/test_symbolic.py +++ b/causationentropy/tests/core/information/test_symbolic.py @@ -79,7 +79,7 @@ def test_invalid_inputs(self): @pytest.mark.parametrize("partition", ["quantile", "uniform", "precomputed"]) @pytest.mark.parametrize("bad", [np.nan, np.inf]) def test_non_finite_values_are_rejected(self, partition, bad): - with pytest.raises(ValueError): + with pytest.raises(ValueError, match="finite"): symbolize(np.array([1.0, bad, 2.0, 3.0]), partition=partition)