#!/usr/bin/env python3
"""Dependency-free audit for the matrix quotient calibration.

This verifies the finite partition census, the dominance/rank direction, the
small hostile collisions, the Jordan scaling formula, and the sharp finite
observability horizon.  The all-dimension claims still rest on the displayed
proofs, not this finite regression suite.
"""

from fractions import Fraction


def partitions(n, ceiling=None):
    if n == 0:
        yield ()
        return
    ceiling = n if ceiling is None else min(n, ceiling)
    for first in range(ceiling, 0, -1):
        for tail in partitions(n - first, first):
            yield (first,) + tail


def conjugate(partition):
    return tuple(sum(part >= k for part in partition)
                 for k in range(1, max(partition, default=0) + 1))


def nullities(partition, n):
    return tuple(sum(min(k, block) for block in partition) for k in range(1, n + 1))


def recover_conjugate(nullity):
    previous = 0
    out = []
    for value in nullity:
        out.append(value - previous)
        previous = value
    while out and out[-1] == 0:
        out.pop()
    return tuple(out)


def centralizer_dimension(partition):
    return sum(column * column for column in conjugate(partition))


def dominates(left, right):
    width = max(len(left), len(right))
    return all(sum(left[:k]) >= sum(right[:k]) for k in range(1, width + 1))


def nilpotent_ranks(partition, n):
    return tuple(n - value for value in nullities(partition, n))


def matrix_multiply(left, right):
    return [[sum(left[i][k] * right[k][j] for k in range(len(right)))
             for j in range(len(right[0]))] for i in range(len(left))]


def row_rank(rows):
    rows = [[Fraction(x) for x in row] for row in rows]
    rank = 0
    for column in range(len(rows[0]) if rows else 0):
        pivot = next((i for i in range(rank, len(rows)) if rows[i][column]), None)
        if pivot is None:
            continue
        rows[rank], rows[pivot] = rows[pivot], rows[rank]
        scale = rows[rank][column]
        rows[rank] = [x / scale for x in rows[rank]]
        for i in range(len(rows)):
            if i != rank and rows[i][column]:
                scale = rows[i][column]
                rows[i] = [x - scale * y for x, y in zip(rows[i], rows[rank])]
        rank += 1
    return rank


def jordan_nilpotent(n):
    return [[1 if j == i + 1 else 0 for j in range(n)] for i in range(n)]


def observability_ranks(n):
    a = jordan_nilpotent(n)
    power = [[int(i == j) for j in range(n)] for i in range(n)]
    rows = []
    ranks = []
    for _ in range(n):
        rows.append(power[0][:])  # C=e_1^T
        ranks.append(row_rank(rows))
        power = matrix_multiply(power, a)
    return tuple(ranks)


def main():
    all_partitions = [(n, p) for n in range(1, 8) for p in partitions(n)]
    assert len(all_partitions) == 44
    for n, partition in all_partitions:
        nullity = nullities(partition, n)
        assert recover_conjugate(nullity) == conjugate(partition)
        assert centralizer_dimension(partition) == sum(x * x for x in recover_conjugate(nullity))

    comparisons = 0
    for n in range(1, 8):
        parts = list(partitions(n))
        for larger in parts:
            for smaller in parts:
                if dominates(larger, smaller):
                    comparisons += 1
                    left = nilpotent_ranks(larger, n)
                    right = nilpotent_ranks(smaller, n)
                    assert all(a >= b for a, b in zip(left, right))
    assert comparisons == 233

    # 2I_2 and J_2(2): equal characteristic polynomial and power traces,
    # different shifted nullities and centralizer dimensions.
    assert nullities((1, 1), 2) == (2, 2)
    assert nullities((2,), 2) == (1, 2)
    assert centralizer_dimension((1, 1)) == 4
    assert centralizer_dimension((2,)) == 2
    assert tuple(2 * (2 ** k) for k in range(1, 7)) == tuple((2 ** k) + (2 ** k) for k in range(1, 7))

    # Exact diagonal-scaling witnesses D_t J_m(lambda) D_t^{-1}.
    scaling_rows = []
    for m, t in ((2, Fraction(1, 2)), (2, Fraction(1, 4)),
                 (3, Fraction(1, 2)), (3, Fraction(1, 4)),
                 (4, Fraction(1, 2)), (4, Fraction(1, 4))):
        excess = (m - 1) * t * t
        condition = t ** (-(m - 1))
        assert excess > 0 and condition >= 1
        scaling_rows.append((m, str(t), str(excess), str(condition)))

    # Cayley--Hamilton horizon sharpness: successive observation rows for
    # (J_n(0), C=e_1^T) increase rank at every horizon through n-1.
    for n in range(1, 8):
        assert observability_ranks(n) == tuple(range(1, n + 1))

    print("partitions_through_dimension_7=44")
    print(f"dominance_comparisons={comparisons}")
    print("nullity_receipt_recoveries=44")
    print("centralizer_formula_checks=44")
    print(f"jordan_scaling_rows={len(scaling_rows)}")
    print("sharp_observability_chains=7")
    print("two_by_two_collision=PASS")
    print("status=PASS")


if __name__ == "__main__":
    main()