diff --git a/docs/examples.rst b/docs/examples.rst index 05d2da2..7b32cd8 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -40,11 +40,90 @@ Other literature processes Tent map (Misiurewicz point) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Read through the two-letter kneading partition, split at the critical point +``c = 1/2`` (James et al. :cite:`James2013`, supplement Figs.~6--8): + * :func:`tent_map_misiurewicz_hmm` — non-unifilar HMM * :func:`tent_map_misiurewicz_forward`, :func:`tent_map_misiurewicz_reverse` * :func:`tent_map_misiurewicz_bidirectional` — information anatomy reference * :func:`tent_map_misiurewicz_information_expected` — expected measure dict +The four generating partitions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The critical point can be joined by either or both of its order-1 preimages, +``L = 1/(2a)`` and ``R = 1 - 1/(2a)``, giving four generating partitions of the +*same* dynamics. All four are handled together by one family of constructors, +keyed by the cuts they make: + +* :func:`tent_map_misiurewicz_partition_forward` — the ε-machine +* :func:`tent_map_misiurewicz_partition_symbol_matrices` — the ``T^(x)`` matrices +* :func:`tent_map_misiurewicz_partition_information_expected` — expected measure dict +* :func:`tent_map_misiurewicz_partition_cuts` — the cut points +* :data:`TENT_MAP_MISIUREWICZ_PARTITIONS` — the four keys, in refinement order + +Because every partition is generating they share the entropy rate +``h_mu = log2(a) ≈ 0.8232``, and each is strictly sofic with infinite Markov and +cryptic order. Only the anatomy split moves: + +.. list-table:: + :header-rows: 1 + :widths: 12 16 8 10 32 12 + + * - Partition + - Cuts + - States + - Alphabet + - ``r_mu`` + - ``r_mu`` ≈ + * - ``"c"`` + - ``c`` + - 4 + - 2 + - ``(59 + 7a - 11a**2)/57`` + - 0.6483 + * - ``"Lc"`` + - ``L, c`` + - 5 + - 3 + - ``(56 + 25a - 23a**2)/57`` + - 0.4953 + * - ``"cR"`` + - ``c, R`` + - 5 + - 3 + - ``(1 - 6a + 4a**2)/19`` + - 0.1529 + * - ``"LcR"`` + - ``L, c, R`` + - 5 + - 4 + - ``0`` + - 0.0000 + +Reducing modulo the parameter's minimal polynomial ``a**3 = 2a + 2`` makes every +transition probability and every ephemeral rate a quadratic in ``a`` with +rational coefficients. Two exact identities fall out. The rate is *modular* over +the two cuts, + +.. math:: + + r_\mu(\{c\}) - r_\mu(\{L, c\}) - r_\mu(\{c, R\}) + r_\mu(\{L, c, R\}) = 0, + +so each cut is worth a fixed number of bits whether or not the other has been +made; and the ``L`` cut's share, ``(1 - 6a + 4a**2)/19``, is exactly the +invariant measure of the two cells of the interval Markov chain that it +separates. Adding both cuts drives ``r_mu`` to zero, leaving the whole entropy +rate as bound information. + +``"c"`` is the kneading partition above, so +:func:`tent_map_misiurewicz_partition_forward` with ``"c"`` reproduces +:func:`tent_map_misiurewicz_forward` up to state names: the family names states +by decreasing stationary probability, making them comparable across partitions, +whereas :func:`tent_map_misiurewicz_forward` keeps the published figure's labels. +The three refinements are derived from the exact interval Markov chain, since the +2013 supplement's figures cover only the kneading partition. + Sofic-Dyck shifts ~~~~~~~~~~~~~~~~~ @@ -144,6 +223,11 @@ API .. autofunction:: tent_map_misiurewicz_bidirectional .. autofunction:: tent_map_misiurewicz_a .. autofunction:: tent_map_misiurewicz_information_expected +.. autofunction:: tent_map_misiurewicz_partition_cuts +.. autofunction:: tent_map_misiurewicz_partition_forward +.. autofunction:: tent_map_misiurewicz_partition_symbol_matrices +.. autofunction:: tent_map_misiurewicz_partition_information_expected +.. autodata:: TENT_MAP_MISIUREWICZ_PARTITIONS .. autofunction:: dyck_shift_order .. autofunction:: motzkin_shift .. autofunction:: sofic_dyck_fig1_shift diff --git a/sofic/examples/__init__.py b/sofic/examples/__init__.py index c74872b..b7b2156 100644 --- a/sofic/examples/__init__.py +++ b/sofic/examples/__init__.py @@ -3,6 +3,7 @@ from importlib import import_module as _import_module from sofic.examples.epsilon_machines import ( + TENT_MAP_MISIUREWICZ_PARTITIONS, alternating_biased_coins, bernoulli, butterfly_process, @@ -25,6 +26,10 @@ tent_map_misiurewicz_forward, tent_map_misiurewicz_hmm, tent_map_misiurewicz_information_expected, + tent_map_misiurewicz_partition_cuts, + tent_map_misiurewicz_partition_forward, + tent_map_misiurewicz_partition_information_expected, + tent_map_misiurewicz_partition_symbol_matrices, tent_map_misiurewicz_reverse, ) from sofic.examples.processes import * @@ -63,11 +68,16 @@ "sofic_dyck_fig1_shift", "sofic_dyck_nondeterminizable_shift", "sofic_dyck_zeta_example_shift", + "TENT_MAP_MISIUREWICZ_PARTITIONS", "tent_map_misiurewicz_a", "tent_map_misiurewicz_bidirectional", "tent_map_misiurewicz_forward", "tent_map_misiurewicz_hmm", "tent_map_misiurewicz_information_expected", + "tent_map_misiurewicz_partition_cuts", + "tent_map_misiurewicz_partition_forward", + "tent_map_misiurewicz_partition_information_expected", + "tent_map_misiurewicz_partition_symbol_matrices", "tent_map_misiurewicz_reverse", ] __all__ += _process_all diff --git a/sofic/examples/epsilon_machines.py b/sofic/examples/epsilon_machines.py index 73b8bf7..c999185 100644 --- a/sofic/examples/epsilon_machines.py +++ b/sofic/examples/epsilon_machines.py @@ -12,7 +12,11 @@ - Golden-mean shift (forbid ``11``, Parry max-entropy): standard symbolic dynamics; see e.g. Ellison et al., arXiv:1107.2168 Fig.~2. - Tent map (Misiurewicz point): James, Burke & Crutchfield (2013), supplement - to *Chaos Forgets and Remembers*; Figs.~6--8. + to *Chaos Forgets and Remembers*; Figs.~6--8. The ``partition`` family reads + the same dynamics through all four generating partitions built from the + critical point and its two order-1 preimages; only the kneading partition + appears in that paper's figures, so the three refinements are derived from the + interval Markov chain instead. """ from __future__ import annotations @@ -438,6 +442,11 @@ def tent_map_misiurewicz_fig7_symbol_matrices( Fig.~7. State ``A`` emits only ``1``; ``D`` has a nontrivial ``0`` branch to ``C`` and a ``1`` self-loop. When ``a`` is a sympy expression the matrices use object dtype with exact entries. + + This is the ``"c"`` member of the four-partition family of + :func:`tent_map_misiurewicz_partition_symbol_matrices`, keeping the figure's + state names. That function instead names states by decreasing stationary + probability, so its ``A, B, C, D`` are this function's ``D, C, B, A``. """ from sofic.generators.prob import is_symbolic, zeros @@ -465,6 +474,252 @@ def tent_map_misiurewicz_forward(a: Any | None = None) -> EpsilonMachine: return from_symbol_matrices(states, symbols, matrices) +#: Keys of the four generating partitions of the tent map at the Misiurewicz +#: point, in refinement order. Each names the cuts added to the critical point +#: ``c = 1/2``: nothing, the left preimage ``L = 1/(2a)``, the right preimage +#: ``R = 1 - 1/(2a)``, or both. See +#: :func:`tent_map_misiurewicz_partition_symbol_matrices`. +TENT_MAP_MISIUREWICZ_PARTITIONS: tuple[str, ...] = ("c", "Lc", "cR", "LcR") + +#: ``partition -> (states, alphabet, edges)`` where each edge is +#: ``(source, symbol, target, (k0, k1, k2), d)`` standing for the transition +#: probability ``(k0 + k1 * a + k2 * a**2) / d``. Every entry is reduced modulo +#: the parameter's minimal polynomial ``a**3 = 2a + 2``, which is why no +#: probability carries an ``a``-dependent denominator. States are named by +#: decreasing stationary probability, uniformly across the four partitions. +_TENT_MAP_PARTITION_EDGES: dict[ + str, + tuple[tuple[str, ...], tuple[int, ...], tuple[tuple[str, int, str, tuple[int, int, int], int], ...]], +] = { + "c": ( + ("A", "B", "C", "D"), + (0, 1), + ( + ("A", 0, "B", (4, 0, -1), 2), + ("A", 1, "A", (-2, 0, 1), 2), + ("B", 0, "D", (2, -2, 1), 6), + ("B", 1, "A", (4, 2, -1), 6), + ("C", 0, "B", (0, -1, 1), 2), + ("C", 1, "D", (2, 1, -1), 2), + ("D", 1, "C", (1, 0, 0), 1), + ), + ), + "Lc": ( + ("A", "B", "C", "D", "E"), + (0, 1, 2), + ( + ("A", 0, "E", (2, -1, 0), 2), + ("A", 1, "B", (2, 1, -1), 2), + ("A", 2, "A", (-2, 0, 1), 2), + ("B", 2, "A", (1, 0, 0), 1), + ("C", 0, "E", (-1, -1, 1), 2), + ("C", 1, "B", (1, 0, 0), 2), + ("C", 2, "D", (2, 1, -1), 2), + ("D", 2, "C", (1, 0, 0), 1), + ("E", 1, "D", (1, 0, 0), 1), + ), + ), + "cR": ( + ("A", "B", "C", "D", "E"), + (0, 1, 2), + ( + ("A", 0, "B", (1, 0, 0), 1), + ("B", 0, "D", (2, -2, 1), 6), + ("B", 1, "C", (-2, -1, 2), 6), + ("B", 2, "A", (2, 1, -1), 2), + ("C", 1, "C", (-2, 0, 1), 2), + ("C", 2, "A", (4, 0, -1), 2), + ("D", 1, "E", (2, 1, -1), 2), + ("D", 2, "A", (0, -1, 1), 2), + ("E", 1, "D", (1, 0, 0), 1), + ), + ), + "LcR": ( + ("A", "B", "C", "D", "E"), + (0, 1, 2, 3), + ( + ("A", 2, "A", (-2, 0, 1), 2), + ("A", 3, "B", (4, 0, -1), 2), + ("B", 0, "D", (2, -2, 1), 6), + ("B", 1, "A", (4, 2, -1), 6), + ("C", 2, "E", (2, 1, -1), 2), + ("C", 3, "B", (0, -1, 1), 2), + ("D", 1, "C", (1, 0, 0), 1), + ("E", 2, "C", (1, 0, 0), 1), + ), + ), +} + +#: ``partition -> ((k0, k1, k2), d)`` for the ephemeral information rate +#: ``r_mu = (k0 + k1 * a + k2 * a**2) / d``, again reduced modulo +#: ``a**3 = 2a + 2``. For ``"c"`` this is the published rate of James et al. +#: (2013) in reduced form; see +#: :func:`tent_map_misiurewicz_partition_information_expected`. +_TENT_MAP_PARTITION_EPHEMERAL: dict[str, tuple[tuple[int, int, int], int]] = { + "c": ((59, 7, -11), 57), + "Lc": ((56, 25, -23), 57), + "cR": ((1, -6, 4), 19), + "LcR": ((0, 0, 0), 1), +} + + +def _tent_map_partition_check(partition: str) -> None: + if partition not in _TENT_MAP_PARTITION_EDGES: + raise ValueError(f"unknown tent-map partition {partition!r}; expected one of {TENT_MAP_MISIUREWICZ_PARTITIONS}") + + +def _tent_map_quadratic(coeffs: tuple[int, int, int], denom: int, a: Any, symbolic: bool) -> Any: + """Evaluate ``(k0 + k1 * a + k2 * a**2) / denom`` exactly or in floats.""" + k0, k1, k2 = coeffs + if symbolic: + import sympy as sp + + return sp.Rational(k0, denom) + sp.Rational(k1, denom) * a + sp.Rational(k2, denom) * a**2 + return (k0 + k1 * a + k2 * a**2) / denom + + +def tent_map_misiurewicz_partition_cuts(partition: str, a: Any | None = None) -> tuple[Any, ...]: + """Return the ascending cut points of one of the four tent-map partitions. + + ``"c"`` cuts only at the critical point; ``"Lc"`` and ``"cR"`` add one + order-1 preimage of it; ``"LcR"`` adds both. With a symbolic ``a`` the cuts + are exact sympy expressions. + """ + _tent_map_partition_check(partition) + from sofic.generators.prob import is_symbolic + + if a is None: + a = tent_map_misiurewicz_a() + if is_symbolic(a): + import sympy as sp + + c = sp.Rational(1, 2) + left = 1 / (2 * a) + else: + c = 0.5 + left = 1.0 / (2.0 * a) + right = 1 - left + return {"c": (c,), "Lc": (left, c), "cR": (c, right), "LcR": (left, c, right)}[partition] + + +def tent_map_misiurewicz_partition_symbol_matrices( + partition: str, + a: Any | None = None, +) -> tuple[tuple[str, ...], tuple[int, ...], dict[int, np.ndarray]]: + """Symbol matrices ``T^(x)`` for one of the tent map's four generating partitions. + + At the Misiurewicz parameter the interval ``[0, 1]`` can be cut at the + critical point ``c = 1/2`` and, optionally, at either or both of its + order-1 preimages ``L = 1/(2a)`` and ``R = 1 - 1/(2a)``. All four choices + are generating, so all four read out the *same* dynamics at the same entropy + rate ``h_mu = log2(a)``; they differ in how many letters they spend and in + how much of that rate survives as bound information: + + ========== ========================= ====== ======== ================ + Partition Cells States Alphabet ``r_mu`` + ========== ========================= ====== ======== ================ + ``"c"`` ``c`` 4 2 ``(59 + 7a - 11a**2)/57`` + ``"Lc"`` ``L, c`` 5 3 ``(56 + 25a - 23a**2)/57`` + ``"cR"`` ``c, R`` 5 3 ``(1 - 6a + 4a**2)/19`` + ``"LcR"`` ``L, c, R`` 5 4 ``0`` + ========== ========================= ====== ======== ================ + + Symbols number the cells left to right, so ``"LcR"`` emits ``0`` on + ``[0, L)``, ``1`` on ``[L, c)``, ``2`` on ``[c, R)`` and ``3`` on + ``[R, 1]``. Reducing by the parameter's minimal polynomial + ``a**3 = 2a + 2`` makes every transition probability a quadratic in ``a`` + with rational coefficients, so none of them carries an ``a``-dependent + denominator. States are named by decreasing stationary probability in every + partition, which makes them comparable across the family; for ``"c"`` that + relabels the published figure, whose ``A, B, C, D`` are this function's + ``D, C, B, A`` (see :func:`tent_map_misiurewicz_fig7_symbol_matrices`). + + Derived from the exact interval Markov chain on the forward-orbit closure of + ``{c, L, R}``. The tent map, the Misiurewicz parameter and the ``"c"`` + presentation are from James, Burke & Crutchfield, *Chaos Forgets and + Remembers* (2013) :cite:`James2013`; that paper's figures cover only the + kneading partition, so the three refinements have no published figure to + cite. + """ + _tent_map_partition_check(partition) + from sofic.generators.prob import is_symbolic, zeros + + if a is None: + a = tent_map_misiurewicz_a() + states, alphabet, edges = _TENT_MAP_PARTITION_EDGES[partition] + symbolic = is_symbolic(a) + index = {state: i for i, state in enumerate(states)} + size = len(states) + matrices = {symbol: zeros((size, size), symbolic=symbolic) for symbol in alphabet} + for source, symbol, target, coeffs, denom in edges: + matrices[symbol][index[source], index[target]] = _tent_map_quadratic(coeffs, denom, a, symbolic) + return states, alphabet, matrices + + +def tent_map_misiurewicz_partition_forward(partition: str, a: Any | None = None) -> EpsilonMachine: + """Forward ε-machine of the tent map under one of its four generating partitions. + + See :func:`tent_map_misiurewicz_partition_symbol_matrices` for the partitions + and their presentations. Every one of the four is unifilar and strictly + sofic -- Markov and cryptic orders are infinite throughout -- so refining the + partition never buys finite memory. What it buys is bound information: + ``r_mu`` falls from ``0.6483`` bits/symbol at ``"c"`` to exactly zero at + ``"LcR"``. + + ``"c"`` is :func:`tent_map_misiurewicz_forward` up to the state relabeling + noted in :func:`tent_map_misiurewicz_partition_symbol_matrices`. + """ + states, symbols, matrices = tent_map_misiurewicz_partition_symbol_matrices(partition, a) + return from_symbol_matrices(states, symbols, matrices) + + +def tent_map_misiurewicz_partition_information_expected( + partition: str, + a: Any | None = None, +) -> dict[str, Any]: + """Closed-form anatomy of one of the tent map's four generating partitions. + + All four are generating, so ``entropy_rate = log2(a)`` throughout and only + the split into ``bound_mu`` and ``ephemeral_mu`` changes. Each ephemeral + rate is a quadratic in ``a`` with rational coefficients, tabulated in + :func:`tent_map_misiurewicz_partition_symbol_matrices`. + + Two exact facts about the family are worth noting. First, ``"LcR"`` has + ``r_mu = 0``: its machine is unifilar, no two edges share both a source and + a target, and every branch leads to a state with a distinguishable future, + so the past fixes the causal state, the future fixes the successor, and the + two together name the emitted symbol -- leaving nothing for + ``r_mu = H[X_0 | past, future]`` to measure. Second, the ephemeral rate is + *modular* over the two available cuts, + + ``r_mu("c") - r_mu("Lc") - r_mu("cR") + r_mu("LcR") = 0`` + + identically in ``a``, so each cut is worth a fixed number of bits regardless + of whether the other has been made; the ``L`` cut is worth + ``r_mu("cR") = (1 - 6a + 4a**2)/19``, which is also exactly the invariant + measure of the two fine cells it separates. + + For ``"c"`` this reproduces :func:`tent_map_misiurewicz_information_expected`, + which states the same rate in the unreduced form published by James, Burke & + Crutchfield (2013) :cite:`James2013`. + """ + _tent_map_partition_check(partition) + from sofic.generators.prob import is_symbolic + + if a is None: + a = tent_map_misiurewicz_a() + coeffs, denom = _TENT_MAP_PARTITION_EPHEMERAL[partition] + if is_symbolic(a): + import sympy as sp + + h_mu = sp.log(a, 2) + r_mu = sp.simplify(_tent_map_quadratic(coeffs, denom, a, True)) + return {"bound_mu": sp.simplify(h_mu - r_mu), "ephemeral_mu": r_mu, "entropy_rate": h_mu} + h_mu = math.log2(a) + r_mu = _tent_map_quadratic(coeffs, denom, a, False) + return {"bound_mu": h_mu - r_mu, "ephemeral_mu": r_mu, "entropy_rate": h_mu} + + def tent_map_misiurewicz_hmm(a: Any | None = None) -> MealyHMM: """Non-unifilar HMM from supplement Fig.~6 (right). @@ -777,6 +1032,10 @@ def tent_map_misiurewicz_information_expected(a: Any | None = None) -> dict[str, Returns floats when ``a`` is numeric, or sympy expressions when ``a`` is symbolic. The ephemeral rate is ``r_μ = (1/4)*(3 - 2/(a+1) - 4/(a+2) + 9/(2a+3))``. + + This is the kneading partition, i.e. the ``"c"`` member of the family of + :func:`tent_map_misiurewicz_partition_information_expected`, which states the + same rate reduced modulo ``a**3 = 2a + 2`` to ``(59 + 7a - 11a**2)/57``. """ from sofic.generators.prob import is_symbolic diff --git a/tests/test_information_anatomy.py b/tests/test_information_anatomy.py index 21700e9..3d4a350 100644 --- a/tests/test_information_anatomy.py +++ b/tests/test_information_anatomy.py @@ -2,10 +2,13 @@ from __future__ import annotations +import math + import pytest from sofic.examples import ( NRPS, + TENT_MAP_MISIUREWICZ_PARTITIONS, bernoulli, butterfly_process, even_process, @@ -13,9 +16,13 @@ golden_mean_forward, golden_mean_reverse, nemo_process, + tent_map_misiurewicz_a, tent_map_misiurewicz_bidirectional, tent_map_misiurewicz_forward, tent_map_misiurewicz_information_expected, + tent_map_misiurewicz_partition_cuts, + tent_map_misiurewicz_partition_forward, + tent_map_misiurewicz_partition_information_expected, ) from sofic.generators.bidirectional_epsilon_machine import BidirectionalEpsilonMachine @@ -481,6 +488,279 @@ def test_tent_map_misiurewicz_bidirectional_regression(): assert bidir.crypticity() == pytest.approx(bidir.statistical_complexity() - bidir.excess_entropy(), abs=1e-9) +# The four generating partitions of the tent map at the Misiurewicz point, each +# built from the critical point plus zero, one or both of its order-1 preimages. +TENT_MAP_PARTITION_SHAPE = { + "c": (["A", "B", "C", "D"], [0, 1]), + "Lc": (["A", "B", "C", "D", "E"], [0, 1, 2]), + "cR": (["A", "B", "C", "D", "E"], [0, 1, 2]), + "LcR": (["A", "B", "C", "D", "E"], [0, 1, 2, 3]), +} + +TENT_MAP_PARTITION_WEIGHTS = { + "c": {"A": 0.4870384416, "B": 0.2882345739, "C": 0.1123634923, "D": 0.1123634923}, + "Lc": { + "A": 0.4870384416, + "B": 0.2117654261, + "C": 0.1123634923, + "D": 0.1123634923, + "E": 0.0764691477, + }, + "cR": { + "A": 0.2882345739, + "B": 0.2882345739, + "C": 0.2752730155, + "D": 0.1123634923, + "E": 0.0358943445, + }, + "LcR": { + "A": 0.4870384416, + "B": 0.2882345739, + "C": 0.1123634923, + "D": 0.0764691477, + "E": 0.0358943445, + }, +} + +# Two-letter blocks the partition forbids; the coarsest partition forbids none. +TENT_MAP_PARTITION_FORBIDDEN_PAIRS = { + "c": set(), + "Lc": {(0, 0), (0, 2), (1, 0), (1, 1)}, + "cR": {(1, 0), (2, 1), (2, 2)}, + "LcR": {(0, 0), (0, 2), (0, 3), (1, 0), (1, 1), (2, 0), (2, 1), (3, 2), (3, 3)}, +} + +TENT_MAP_PARTITION_EPHEMERAL = { + "c": 0.6482578367935150, + "Lc": 0.4953195413389188, + "cR": 0.1529382954545961, + "LcR": 0.0, +} + +# (statistical complexity, excess entropy). C_μ - E is the crypticity. +TENT_MAP_PARTITION_COMPLEXITY = { + "c": (1.7315174883104, 0.3175426940962), + "Lc": (1.9720901450910, 0.5581153508764), + "cR": (2.0735446304436, 0.9001424930096), + "LcR": (1.8330694405681, 1.1407151497908), +} + + +def test_tent_map_partition_keys_are_ordered_by_refinement(): + assert TENT_MAP_MISIUREWICZ_PARTITIONS == ("c", "Lc", "cR", "LcR") + + a = tent_map_misiurewicz_a() + left, right = 1 / (2 * a), 1 - 1 / (2 * a) + assert tent_map_misiurewicz_partition_cuts("c") == pytest.approx([0.5]) + assert tent_map_misiurewicz_partition_cuts("Lc") == pytest.approx([left, 0.5]) + assert tent_map_misiurewicz_partition_cuts("cR") == pytest.approx([0.5, right]) + assert tent_map_misiurewicz_partition_cuts("LcR") == pytest.approx([left, 0.5, right]) + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_tent_map_partition_presentation(partition): + """Each partition gives a unifilar, strictly sofic machine of the stated shape.""" + forward = tent_map_misiurewicz_partition_forward(partition) + forward.validate_stochastic() + + states, alphabet = TENT_MAP_PARTITION_SHAPE[partition] + assert sorted(forward.states()) == states + assert sorted(forward.observation_alphabet) == alphabet + assert forward.is_unifilar() + # Refining never buys finite memory: r_μ falls to zero but the process stays + # strictly sofic, with infinite Markov and cryptic order throughout. + assert forward.is_strictly_sofic() + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_tent_map_partition_stationary_distribution(partition): + forward = tent_map_misiurewicz_partition_forward(partition) + index = forward.reindex() + pi = forward.stationary_distribution() + weights = {index.state(i): float(pi[i]) for i in range(len(index.states))} + + assert weights == pytest.approx(TENT_MAP_PARTITION_WEIGHTS[partition], abs=1e-9) + + +def test_tent_map_partition_edge_probabilities_are_quadratics_in_a(): + """Reduced by ``a**3 = 2a + 2`` every branching probability is a quadratic in ``a``.""" + a = tent_map_misiurewicz_a() + expected = { + "c": { + ("A", 0, "B"): (4 - a**2) / 2, + ("A", 1, "A"): (a**2 - 2) / 2, + ("B", 0, "D"): (2 - 2 * a + a**2) / 6, + ("B", 1, "A"): (4 + 2 * a - a**2) / 6, + ("C", 0, "B"): (a**2 - a) / 2, + ("C", 1, "D"): (2 + a - a**2) / 2, + ("D", 1, "C"): 1.0, + }, + "Lc": { + ("A", 0, "E"): (2 - a) / 2, + ("A", 1, "B"): (2 + a - a**2) / 2, + ("A", 2, "A"): (a**2 - 2) / 2, + ("B", 2, "A"): 1.0, + ("C", 0, "E"): (a**2 - a - 1) / 2, + ("C", 1, "B"): 0.5, + ("C", 2, "D"): (2 + a - a**2) / 2, + ("D", 2, "C"): 1.0, + ("E", 1, "D"): 1.0, + }, + "cR": { + ("A", 0, "B"): 1.0, + ("B", 0, "D"): (2 - 2 * a + a**2) / 6, + ("B", 1, "C"): (2 * a**2 - a - 2) / 6, + ("B", 2, "A"): (2 + a - a**2) / 2, + ("C", 1, "C"): (a**2 - 2) / 2, + ("C", 2, "A"): (4 - a**2) / 2, + ("D", 1, "E"): (2 + a - a**2) / 2, + ("D", 2, "A"): (a**2 - a) / 2, + ("E", 1, "D"): 1.0, + }, + "LcR": { + ("A", 2, "A"): (a**2 - 2) / 2, + ("A", 3, "B"): (4 - a**2) / 2, + ("B", 0, "D"): (2 - 2 * a + a**2) / 6, + ("B", 1, "A"): (4 + 2 * a - a**2) / 6, + ("C", 2, "E"): (2 + a - a**2) / 2, + ("C", 3, "B"): (a**2 - a) / 2, + ("D", 1, "C"): 1.0, + ("E", 2, "C"): 1.0, + }, + } + + for partition, edges in expected.items(): + forward = tent_map_misiurewicz_partition_forward(partition) + actual = {(t.source, t.data["emission"], t.target): float(t.data["prob"]) for t in forward.graph.transitions()} + assert actual.keys() == edges.keys(), partition + for edge, probability in edges.items(): + assert actual[edge] == pytest.approx(probability, abs=1e-12), (partition, edge) + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_tent_map_partition_forbidden_blocks(partition): + forward = tent_map_misiurewicz_partition_forward(partition) + _states, alphabet = TENT_MAP_PARTITION_SHAPE[partition] + words = forward.word_probabilities(2) + forbidden = TENT_MAP_PARTITION_FORBIDDEN_PAIRS[partition] + + for word in forbidden: + assert float(words.get(word, 0.0)) == pytest.approx(0.0, abs=1e-12) + allowed = {word for word, p in words.items() if float(p) > 1e-12} + assert allowed.isdisjoint(forbidden) + assert len(allowed) == len(alphabet) ** 2 - len(forbidden) + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_tent_map_partitions_are_all_generating(partition): + """Every partition is generating, so all four share the entropy rate ``log2(a)``.""" + pytest.importorskip("dit") + a = tent_map_misiurewicz_a() + forward = tent_map_misiurewicz_partition_forward(partition) + + assert forward.entropy_rate() == pytest.approx(math.log2(a), abs=1e-9) + assert forward.entropy_rate() == pytest.approx(tent_map_misiurewicz_forward().entropy_rate(), abs=1e-9) + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_tent_map_partition_anatomy_matches_closed_form(partition): + """The machine's measured anatomy matches the tabulated closed form.""" + pytest.importorskip("dit") + forward = tent_map_misiurewicz_partition_forward(partition) + expected = tent_map_misiurewicz_partition_information_expected(partition) + r_mu = TENT_MAP_PARTITION_EPHEMERAL[partition] + + assert expected["ephemeral_mu"] == pytest.approx(r_mu, abs=1e-12) + assert forward.ephemeral_information() == pytest.approx(r_mu, abs=1e-9) + assert forward.bound_information() == pytest.approx(expected["bound_mu"], abs=1e-9) + assert forward.entropy_rate() == pytest.approx(expected["entropy_rate"], abs=1e-9) + assert expected["bound_mu"] + expected["ephemeral_mu"] == pytest.approx(expected["entropy_rate"], abs=1e-12) + + +def test_tent_map_partition_ephemeral_rate_is_modular_over_the_two_cuts(): + """Each preimage cut is worth a fixed number of bits, independent of the other.""" + rates = { + p: tent_map_misiurewicz_partition_information_expected(p)["ephemeral_mu"] for p in ("c", "Lc", "cR", "LcR") + } + + # The ``L`` cut removes the same amount whether or not ``R`` has been cut. + assert rates["c"] - rates["Lc"] == pytest.approx(rates["cR"] - rates["LcR"], abs=1e-12) + # Equivalently, the interaction term vanishes and the two shares sum to the whole. + assert rates["c"] - rates["Lc"] - rates["cR"] + rates["LcR"] == pytest.approx(0.0, abs=1e-12) + assert rates["Lc"] + rates["cR"] == pytest.approx(rates["c"], abs=1e-12) + + # The ``L`` cut's share is the measure of the two fine cells it separates. + a = tent_map_misiurewicz_a() + assert rates["cR"] == pytest.approx(2 * (4 * a**2 - 6 * a + 1) / 38, abs=1e-12) + + +def test_tent_map_partition_c_is_the_published_figure_relabeled(): + """``"c"`` is the Fig.~7 machine with states renamed by decreasing weight.""" + pytest.importorskip("dit") + family = tent_map_misiurewicz_partition_forward("c") + published = tent_map_misiurewicz_forward() + relabel = {"A": "D", "B": "C", "C": "B", "D": "A"} + + family_edges = { + (relabel[t.source], t.data["emission"], relabel[t.target]): float(t.data["prob"]) + for t in family.graph.transitions() + } + published_edges = { + (t.source, t.data["emission"], t.target): float(t.data["prob"]) for t in published.graph.transitions() + } + assert family_edges.keys() == published_edges.keys() + for edge, probability in published_edges.items(): + assert family_edges[edge] == pytest.approx(probability, abs=1e-12) + + # And the two closed forms of the kneading rate agree. + reduced = tent_map_misiurewicz_partition_information_expected("c") + as_published = tent_map_misiurewicz_information_expected() + for key in ("entropy_rate", "ephemeral_mu", "bound_mu"): + assert reduced[key] == pytest.approx(as_published[key], abs=1e-12) + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_tent_map_partition_statistical_complexity_exceeds_excess_entropy(partition): + """Regression on C_μ and E; the gap is the machine's crypticity.""" + pytest.importorskip("dit") + forward = tent_map_misiurewicz_partition_forward(partition) + c_mu, excess = TENT_MAP_PARTITION_COMPLEXITY[partition] + + assert forward.statistical_complexity() == pytest.approx(c_mu, abs=1e-9) + assert forward.excess_entropy() == pytest.approx(excess, abs=1e-9) + assert forward.statistical_complexity() > forward.excess_entropy() + + +def test_tent_map_partition_excess_entropy_is_modular_but_complexity_is_not(): + """``E`` inherits the cuts' additivity; stored history does not.""" + pytest.importorskip("dit") + excess = {} + complexity = {} + for partition in TENT_MAP_MISIUREWICZ_PARTITIONS: + forward = tent_map_misiurewicz_partition_forward(partition) + excess[partition] = forward.excess_entropy() + complexity[partition] = forward.statistical_complexity() + + interaction = excess["c"] - excess["Lc"] - excess["cR"] + excess["LcR"] + assert interaction == pytest.approx(0.0, abs=1e-9) + + interaction = complexity["c"] - complexity["Lc"] - complexity["cR"] + complexity["LcR"] + assert abs(interaction) > 0.4 + + +@pytest.mark.parametrize( + "call", + [ + tent_map_misiurewicz_partition_cuts, + tent_map_misiurewicz_partition_forward, + tent_map_misiurewicz_partition_information_expected, + ], +) +def test_tent_map_partition_rejects_unknown_key(call): + with pytest.raises(ValueError, match="unknown tent-map partition"): + call("Rc") + + def test_tent_forward_matches_generator_path(): pytest.importorskip("dit") from sofic.examples.epsilon_machines import tent_map_misiurewicz_hmm diff --git a/tests/test_symbolic_hmm.py b/tests/test_symbolic_hmm.py index 142ec95..a9a766c 100644 --- a/tests/test_symbolic_hmm.py +++ b/tests/test_symbolic_hmm.py @@ -12,11 +12,15 @@ import sympy as sp from sofic.examples.epsilon_machines import ( + TENT_MAP_MISIUREWICZ_PARTITIONS, tent_map_misiurewicz_a, tent_map_misiurewicz_bidirectional, tent_map_misiurewicz_forward, tent_map_misiurewicz_hmm, tent_map_misiurewicz_information_expected, + tent_map_misiurewicz_partition_forward, + tent_map_misiurewicz_partition_information_expected, + tent_map_misiurewicz_partition_symbol_matrices, ) from sofic.generators.epsilon_machine import EpsilonMachine from sofic.generators.prob import ( @@ -28,6 +32,61 @@ from sofic.generators.words import hmm_word_probability +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_symbolic_partition_rows_sum_to_one_identically(partition): + """Each partition's rows normalize for a free ``a``, not just at the root.""" + a = sp.symbols("a", positive=True) + states, symbols, matrices = tent_map_misiurewicz_partition_symbol_matrices(partition, a) + for i in range(len(states)): + row = sum(matrices[x][i, j] for x in symbols for j in range(len(states))) + assert sp.simplify(row) == 1 + + +@pytest.mark.parametrize("partition", TENT_MAP_MISIUREWICZ_PARTITIONS) +def test_symbolic_partition_substitution_matches_numeric(partition): + """Substituting the Misiurewicz root into the symbolic machine recovers the floats.""" + a = sp.symbols("a", positive=True) + a_num = tent_map_misiurewicz_a() + symbolic = tent_map_misiurewicz_partition_forward(partition, a) + numeric = tent_map_misiurewicz_partition_forward(partition) + + def edges(machine, substitute): + # Certain transitions are stored as an exact ``1``, so sympify before substituting. + return { + (t.source, t.data["emission"], t.target): float( + sp.sympify(t.data["prob"]).subs(a, a_num) if substitute else t.data["prob"] + ) + for t in machine.graph.transitions() + } + + assert edges(symbolic, True) == pytest.approx(edges(numeric, False), abs=1e-12) + assert float(symbolic.entropy_rate().subs(a, a_num)) == pytest.approx(math.log2(a_num), abs=1e-9) + + +def test_symbolic_partition_ephemeral_rate_is_modular_identically_in_a(): + """The cuts' additivity is an identity in ``a``, not a numerical coincidence.""" + a = sp.symbols("a", positive=True) + rates = { + partition: tent_map_misiurewicz_partition_information_expected(partition, a)["ephemeral_mu"] + for partition in TENT_MAP_MISIUREWICZ_PARTITIONS + } + + interaction = rates["c"] - rates["Lc"] - rates["cR"] + rates["LcR"] + assert sp.simplify(interaction) == 0 + + +def test_symbolic_kneading_rate_reduces_to_the_published_form(): + """``(59 + 7a - 11a**2)/57`` is the published rate modulo ``a**3 = 2a + 2``.""" + a = sp.symbols("a", positive=True) + published = tent_map_misiurewicz_information_expected(a)["ephemeral_mu"] + reduced = tent_map_misiurewicz_partition_information_expected("c", a)["ephemeral_mu"] + + difference = sp.together(sp.simplify(published - reduced)) + numerator, _denominator = sp.fraction(sp.cancel(difference)) + remainder = sp.rem(sp.Poly(sp.expand(numerator), a), sp.Poly(a**3 - 2 * a - 2, a)) + assert remainder.as_expr() == 0 + + def test_symbolic_edge_probabilities_preserved(): a = sp.symbols("a", positive=True) hmm = tent_map_misiurewicz_forward(a)