form_geometry.py

"""Regional form observations and exact affine-law admission.

The supplied ordered triples fix the frame
``z=(x0-x1)/sqrt(2) + i*(x0+x1-2*x2)/sqrt(6)``. Its orientation is
derived from scalar EPI, not the primitive nodal phase. Exact rational fields
refer to admitted binary64 inputs and the separately materialized canonical
nodal product. Numerical polar estimates have explicit availability statuses.
The affine admission helper checks a separately supplied fine law against
the same mean/Gram observation. Neither selects a pressure law or evolution.
"""

from __future__ import annotations

import math
from collections.abc import Mapping, Set
from dataclasses import dataclass
from decimal import ROUND_HALF_EVEN, Context, Decimal, localcontext
from fractions import Fraction
from numbers import Real
from typing import Any

from .._exact_time import finite_represented_real
from ..alias import get_attr
from ..constants.aliases import ALIAS_DNFR, ALIAS_EPI, ALIAS_VF
from ..dynamics.canonical import compute_canonical_nodal_derivative
from ..mathematics.epi import BEPIElement
from ..types import require_finite_real_scalar_epi
from ._cycle_algebra import ordered_vector

__all__ = [
    "RegionalFormRegion",
    "RegionalFormObservation",
    "RegionalAffineClosure",
    "observe_regional_form",
    "derive_regional_affine_closure",
]

_MISSING = object()


@dataclass(frozen=True)
class RegionalFormRegion:
    """One fixed real frame; imaginary Gram numerators elsewhere use sqrt(12)."""

    nodes: tuple[Any, Any, Any]
    mean: Fraction
    contrast_a: Fraction
    contrast_b: Fraction
    mean_rate: Fraction
    contrast_a_rate: Fraction
    contrast_b_rate: Fraction
    intensity: Fraction
    intensity_rate: Fraction
    amplitude_estimate: float | None
    phase_estimate: float | None
    amplitude_rate_estimate: float | None
    phase_rate_estimate: float | None
    estimate_status: tuple[tuple[str, str], ...]


@dataclass(frozen=True)
class RegionalFormObservation:
    """Immutable represented-state observation, not a closed reduced model.

    ``gram_real + i*gram_imag_numerator/sqrt(12)`` is ``z*z.conjugate().T``.
    Its rate uses the same convention. Rational products retain tiny nonzero
    information even when a polar estimate is unavailable. No result authenticates
    pressure refresh, a measured derivative, a trajectory or physical emergence.
    """

    nodes: tuple
    epi: tuple[Fraction, ...]
    capacity: tuple[Fraction, ...]
    stored_pressure: tuple[Fraction, ...]
    nodal_rate: tuple[Fraction, ...]
    nodal_rate_rounding_defect: tuple[Fraction, ...]
    regions: tuple[RegionalFormRegion, ...]
    gram_real: tuple[tuple[Fraction, ...], ...]
    gram_imag_numerator: tuple[tuple[Fraction, ...], ...]
    gram_rate_real: tuple[tuple[Fraction, ...], ...]
    gram_rate_imag_numerator: tuple[tuple[Fraction, ...], ...]
    scope: tuple[str, ...] = (
        "supplied_ordered_triples_and_fixed_orthonormal_form_frames",
        "exact_rational_algebra_after_binary64_input_and_nodal_rate_materialization",
        "unforced_nu_times_stored_pressure_prediction_not_measured_derivative",
        "no_pressure_refresh_callbacks_cache_writes_or_primitive_phase_consumption",
        "gamma_operator_jumps_clipping_and_other_runtime_rates_excluded",
        "instantaneous_pushforward_not_closed_reduced_dynamics_or_future_certificate",
        "polar_estimates_not_rigorous_transcendental_enclosures",
    )


@dataclass(frozen=True)
class RegionalAffineClosure:
    """Exact mean/Gram projectability of a supplied fixed affine fine law.

    The law is ``x_dot = generator*x + source`` in ``nodes`` order. Source
    has form/time units; for a nodal pressure source F it is ``nu*F``.
    Each region fixes the same orthonormal frame as the form observer.
    When admitted, the complex contrast generator is
    ``contrast_generator_real + i*sqrt(3)*contrast_generator_imag_over_sqrt3``.
    Its mean law is ``mu_dot = mean_generator*mu + mean_source`` and its
    Gram law is ``H_dot = A*H + H*A.conjugate().T``.

    Inputs and residuals are retained even on failed admission; complex
    reduced generators are unavailable then. The property tests the full
    independent real mean/contrast domain, not a restricted preparation.
    Each block defect subtracts the circulant generated by that block's first
    row, in original node order. Source contrasts use unscaled (a,b) coordinates.
    This detached algebra does not authenticate a graph, pressure realization,
    frozen support/capacity, callback behavior or future runtime execution.
    """

    nodes: tuple
    regions: tuple[tuple[Any, Any, Any], ...]
    generator: tuple[tuple[Fraction, ...], ...]
    source: tuple[Fraction, ...]
    block_circulant_defect: tuple[tuple[Fraction, ...], ...]
    source_contrast_a: tuple[Fraction, ...]
    source_contrast_b: tuple[Fraction, ...]
    mean_generator: tuple[tuple[Fraction, ...], ...]
    mean_source: tuple[Fraction, ...]
    contrast_generator_real: tuple[tuple[Fraction, ...], ...] | None
    contrast_generator_imag_over_sqrt3: tuple[tuple[Fraction, ...], ...] | None
    scope: tuple[str, ...] = (
        "fixed_affine_fine_law_on_full_independent_real_mean_and_contrast_domain",
        "exact_rational_inputs_or_materialized_binary64_coefficients",
        "common_contrast_rotation_fibers_with_source_held_fixed",
        "zero_contrast_and_zero_capacity_need_no_polar_coordinate_or_division",
        "no_pressure_inference_graph_authentication_or_runtime_certificate",
        "no_operator_jumps_clipping_moving_partition_or_new_constitutive_law",
    )

    @property
    def all_state_closed(self) -> bool:
        """Return whether every exact generator/source defect vanishes."""
        return not (
            any(any(row) for row in self.block_circulant_defect)
            or any(self.source_contrast_a)
            or any(self.source_contrast_b)
        )


def _partition(nodes, regions):
    if not nodes:
        raise ValueError("regional form requires a nonempty graph")
    forbidden = (str, bytes, bytearray, Mapping, Set)
    if isinstance(regions, forbidden):
        raise TypeError("regions must be an ordered collection of ordered triples")
    try:
        raw = tuple(regions)
        if any(isinstance(region, forbidden) for region in raw):
            raise TypeError("each region must be an ordered triple")
        blocks = tuple(tuple(region) for region in raw)
    except TypeError as exc:
        raise TypeError("regions must contain ordered triples") from exc
    if not blocks or any(len(block) != 3 for block in blocks):
        raise ValueError("each region must contain exactly three nodes")
    flat = tuple(node for block in blocks for node in block)
    try:
        valid = len(set(flat)) == len(flat) and set(flat) == set(nodes)
    except TypeError as exc:
        raise TypeError("region labels must be graph node identities") from exc
    if not valid:
        raise ValueError("regions must be disjoint and cover exactly the graph nodes")
    return blocks


def derive_regional_affine_closure(
    nodes, regions, *, generator, source
) -> RegionalAffineClosure:
    """Derive exact mean/Gram admission for ``x_dot=G*x+b`` on fixed triples.

    ``nodes`` binds matrix/vector order; ``regions`` is the observer's disjoint
    ordered partition. G and b must be supplied independently of an evaluated
    response. For held nodal pressure ``p=e*G_W*x+F``, pass
    ``G=diag(nu)*e*G_W`` and ``b=diag(nu)*F``. This helper neither reconstructs
    pressure from a derivative nor proves that the caller holds these data.

    In the equally oriented triangle frames, homogeneous mean/Gram closure
    is equivalent to every three-by-three block of G being real circulant.
    An affine source additionally must be constant within each region; its
    mean need not vanish. Rejected residuals describe failed reduction, not
    an invalid fine law. Polar phase is never consumed, including at zero.
    Rationals retain exact values; other admitted reals use binary64 once.
    """
    forbidden = (str, bytes, bytearray, Mapping, Set)
    if isinstance(nodes, forbidden):
        raise TypeError("nodes must be an ordered collection of distinct labels")
    nodes = tuple(nodes)
    if len(set(nodes)) != len(nodes):
        raise ValueError("nodes must be distinct")
    blocks = _partition(nodes, regions)
    if isinstance(generator, forbidden):
        raise TypeError("generator must be an ordered square matrix")
    matrix = tuple(ordered_vector(row, "generator row") for row in generator)
    count = len(nodes)
    if len(matrix) != count or any(len(row) != count for row in matrix):
        raise ValueError("generator must be square in the supplied node order")
    source = ordered_vector(source, "affine rate source")
    if len(source) != count:
        raise ValueError("source must match the supplied node order")
    positions = {node: index for index, node in enumerate(nodes)}
    indices = tuple(tuple(positions[node] for node in block) for block in blocks)
    defects = [[Fraction(0) for _ in nodes] for _ in nodes]
    mean_matrix, real_matrix, imaginary_matrix = [], [], []
    for left in indices:
        mean_row, real_row, imaginary_row = [], [], []
        for right in indices:
            c0, c1, c2 = (matrix[left[0]][j] for j in right)
            for i, fine_i in enumerate(left):
                for j, fine_j in enumerate(right):
                    defects[fine_i][fine_j] = (
                        matrix[fine_i][fine_j] - matrix[left[0]][right[(j - i) % 3]]
                    )
            mean_row.append(sum((matrix[i][j] for i in left for j in right)) / 3)
            real_row.append(c0 - (c1 + c2) / 2)
            imaginary_row.append((c2 - c1) / 2)
        mean_matrix.append(tuple(mean_row))
        real_matrix.append(tuple(real_row))
        imaginary_matrix.append(tuple(imaginary_row))
    source_a = tuple(source[i] - source[j] for i, j, _ in indices)
    source_b = tuple(source[i] + source[j] - 2 * source[k] for i, j, k in indices)
    admitted = not (any(any(row) for row in defects) or any(source_a) or any(source_b))
    return RegionalAffineClosure(
        nodes,
        blocks,
        matrix,
        source,
        tuple(tuple(row) for row in defects),
        source_a,
        source_b,
        tuple(mean_matrix),
        tuple(sum(source[i] for i in block) / 3 for block in indices),
        tuple(real_matrix) if admitted else None,
        tuple(imaginary_matrix) if admitted else None,
    )


def _read(data, aliases, label):
    raw = get_attr(data, aliases, _MISSING, strict=True, conv=lambda value: value)
    if raw is _MISSING:
        raise ValueError(f"missing consumed {label}")
    if aliases == ALIAS_EPI:
        # Preserve raw scalar admission before the BEPI materialization path.
        if isinstance(raw, Real):
            raw = finite_represented_real(raw, label)[0]
        elif not isinstance(raw, BEPIElement):
            raise TypeError(
                f"{label} must be a real scalar or materialized BEPIElement"
            )
        raw = require_finite_real_scalar_epi(raw, label)
    return finite_represented_real(raw, label)


def _decimal(value):
    return Decimal(value.numerator) / Decimal(value.denominator)


def _estimate(value):
    result = float(value)
    if not math.isfinite(result):
        return None, "unrepresentable_overflow"
    if value and result == 0.0:
        return None, "unrepresentable_underflow"
    return result, "available"


def _polar(a, b, da, db, intensity, intensity_rate):
    names = (
        "amplitude_estimate",
        "phase_estimate",
        "amplitude_rate_estimate",
        "phase_rate_estimate",
    )
    if not intensity:
        return (0.0, None, None, None), tuple(
            (name, "available" if i == 0 else "zero_amplitude")
            for i, name in enumerate(names)
        )
    # The caller's Decimal precision, rounding, exponent bounds and traps must
    # not change this observer's estimates or turn them into exceptions.
    with localcontext(Context(prec=80, rounding=ROUND_HALF_EVEN)):
        radius = _decimal(intensity).sqrt()
        amplitude = _estimate(radius)
        radial = _estimate(_decimal(intensity_rate) / (2 * radius))
        angular = _estimate(
            _decimal(a * db - b * da) / (Decimal(12).sqrt() * _decimal(intensity))
        )
        u, v = _decimal(a) / Decimal(2).sqrt(), _decimal(b) / Decimal(6).sqrt()
        scale = max(abs(u), abs(v))
        angle = math.atan2(float(v / scale), float(u / scale))
        phase = (
            (None, "unrepresentable_underflow")
            if b and a > 0 and angle == 0.0
            else (angle, "available")
        )
    pairs = (amplitude, phase, radial, angular)
    return tuple(value for value, _ in pairs), tuple(
        (name, status) for name, (_, status) in zip(names, pairs)
    )


def observe_regional_form(graph, regions) -> RegionalFormObservation:
    """Read complete ordered triples and push stored unforced rates into form.

    All graph nodes must occur exactly once. Each triple fixes the displayed
    frame; no edge, phase or geometry-dependent partition is inferred. Missing
    consumed EPI/capacity/pressure rejects; unused primitive phase is not read.
    Zero capacity is admitted. An overflowing or underflowed nonzero nodal
    product rejects before any derived observation is returned.

    The pressure can be stale or operator-written. Rates are explicitly the
    canonical unforced product of that stored value, not full runtime motion.
    Exact scaled Gram coordinates remain available when a numerical polar
    estimate cannot be represented. Node labels retain their opaque identities.
    """
    nodes = tuple(graph.nodes())
    blocks = _partition(nodes, regions)
    positions = {node: index for index, node in enumerate(nodes)}
    epi, capacity, pressure, rates, defects = [], [], [], [], []
    for node in nodes:
        data = graph.nodes[node]
        _, x = _read(data, ALIAS_EPI, f"EPI[{node!r}]")
        nu, nu_exact = _read(data, ALIAS_VF, f"capacity[{node!r}]")
        p, p_exact = _read(data, ALIAS_DNFR, f"stored pressure[{node!r}]")
        rate = compute_canonical_nodal_derivative(nu, p).derivative
        if nu and p and rate == 0.0:
            raise ValueError("nonzero nodal product underflows to represented zero")
        rate_exact = finite_represented_real(rate, "nodal rate")[1]
        epi.append(x)
        capacity.append(nu_exact)
        pressure.append(p_exact)
        rates.append(rate_exact)
        defects.append(rate_exact - nu_exact * p_exact)

    rows = []
    for block in blocks:
        x0, x1, x2 = (epi[positions[node]] for node in block)
        v0, v1, v2 = (rates[positions[node]] for node in block)
        a, b = x0 - x1, x0 + x1 - 2 * x2
        da, db = v0 - v1, v0 + v1 - 2 * v2
        intensity = a * a / 2 + b * b / 6
        intensity_rate = a * da + b * db / 3
        estimates, statuses = _polar(a, b, da, db, intensity, intensity_rate)
        rows.append(
            RegionalFormRegion(
                block,
                (x0 + x1 + x2) / 3,
                a,
                b,
                (v0 + v1 + v2) / 3,
                da,
                db,
                intensity,
                intensity_rate,
                *estimates,
                statuses,
            )
        )

    def matrix(entry):
        return tuple(tuple(entry(left, right) for right in rows) for left in rows)

    return RegionalFormObservation(
        nodes,
        tuple(epi),
        tuple(capacity),
        tuple(pressure),
        tuple(rates),
        tuple(defects),
        tuple(rows),
        matrix(
            lambda left, right: left.contrast_a * right.contrast_a / 2
            + left.contrast_b * right.contrast_b / 6
        ),
        matrix(
            lambda left, right: left.contrast_b * right.contrast_a
            - left.contrast_a * right.contrast_b
        ),
        matrix(
            lambda left, right: (
                left.contrast_a_rate * right.contrast_a
                + left.contrast_a * right.contrast_a_rate
            )
            / 2
            + (
                left.contrast_b_rate * right.contrast_b
                + left.contrast_b * right.contrast_b_rate
            )
            / 6
        ),
        matrix(
            lambda left, right: left.contrast_b_rate * right.contrast_a
            + left.contrast_b * right.contrast_a_rate
            - left.contrast_a_rate * right.contrast_b
            - left.contrast_a * right.contrast_b_rate
        ),
    )