check_relative_degree_exhaustive.py

#!/usr/bin/env python3
"""Exhaustive exact check for the relative-degree criterion in dimensions 1 and 2.

The domain is every integer linear representation over {a,b} whose entries lie
in {-1,0,1}.  The direct oracle tests the definition using all suffix words of
length below d; this is exact because those columns span the reachable space.
"""

from fractions import Fraction
from itertools import product
import hashlib
import json
import sys


VALUES = (-1, 0, 1)


def dot(row, col):
    return sum(x * y for x, y in zip(row, col))


def row_times(row, matrix):
    d = len(row)
    return tuple(sum(row[i] * matrix[i][j] for i in range(d)) for j in range(d))


def matrix_times(matrix, col):
    return tuple(dot(row, col) for row in matrix)


def add_basis(basis, vector):
    """Append vector iff it increases row rank; dimensions here are at most two."""
    rows = [[Fraction(x) for x in row] for row in (*basis, vector)]
    rank = 0
    columns = len(vector)
    for column in range(columns):
        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]:
                factor = rows[i][column]
                rows[i] = [x - factor * y for x, y in zip(rows[i], rows[rank])]
        rank += 1
    if rank > len(basis):
        basis.append(vector)
        return True
    return False


def reachable_basis(gamma, matrices):
    basis = []
    add_basis(basis, gamma)
    cursor = 0
    while cursor < len(basis):
        vector = basis[cursor]
        cursor += 1
        for matrix in matrices:
            add_basis(basis, matrix_times(matrix, vector))
    return basis


def safe_row(row, gamma, b_matrix, basis):
    row_b = row_times(row, b_matrix)
    return dot(row, gamma) == 0 and all(dot(row_b, col) == 0 for col in basis)


def krylov_algorithm(lam, gamma, a_matrix, b_matrix):
    d = len(lam)
    basis = reachable_basis(gamma, (a_matrix, b_matrix))
    row = lam
    earlier_safe = True
    for degree in range(d):
        if earlier_safe and dot(row_times(row, b_matrix), gamma) != 0:
            return degree
        earlier_safe = earlier_safe and safe_row(row, gamma, b_matrix, basis)
        row = row_times(row, a_matrix)
    return None


def suffix_columns(gamma, matrices, max_length):
    levels = [(gamma,)]
    for _ in range(max_length):
        levels.append(tuple(matrix_times(matrix, col) for col in levels[-1] for matrix in matrices))
    return tuple(col for level in levels for col in level)


def direct_candidates(lam, gamma, a_matrix, b_matrix):
    """Test the definition independently for every candidate 0 <= r < d."""
    d = len(lam)
    suffixes = suffix_columns(gamma, (a_matrix, b_matrix), d - 1)
    candidates = []
    row = lam
    earlier_ok = True
    for degree in range(d):
        witness = dot(row_times(row, b_matrix), gamma) != 0
        if earlier_ok and witness:
            candidates.append(degree)
        row_b = row_times(row, b_matrix)
        stop_zero = dot(row, gamma) == 0
        deviation_zero = all(dot(row_b, suffix) == 0 for suffix in suffixes)
        earlier_ok = earlier_ok and stop_zero and deviation_zero
        row = row_times(row, a_matrix)
    return tuple(candidates)


def matrices(d):
    for entries in product(VALUES, repeat=d * d):
        yield tuple(tuple(entries[i * d : (i + 1) * d]) for i in range(d))


def representations(d):
    vectors = tuple(product(VALUES, repeat=d))
    mats = tuple(matrices(d))
    for lam in vectors:
        for gamma in vectors:
            for a_matrix in mats:
                for b_matrix in mats:
                    yield lam, gamma, a_matrix, b_matrix


def run():
    by_dimension = {}
    total = 0
    degrees = {}
    transcript = hashlib.sha256()
    for d in (1, 2):
        count = 0
        with_degree = 0
        local_degrees = {str(r): 0 for r in range(d)}
        for index, representation in enumerate(representations(d)):
            algorithm = krylov_algorithm(*representation)
            direct = direct_candidates(*representation)
            if len(direct) > 1:
                raise AssertionError(("nonunique", d, index, direct, representation))
            expected = direct[0] if direct else None
            if algorithm != expected:
                raise AssertionError(("mismatch", d, index, algorithm, direct, representation))
            transcript.update(f"{d}:{index}:{algorithm}\n".encode())
            count += 1
            if algorithm is not None:
                with_degree += 1
                local_degrees[str(algorithm)] += 1
        by_dimension[str(d)] = {
            "representations": count,
            "with_relative_degree": with_degree,
            "degree_counts": local_degrees,
        }
        total += count
        degrees.update({f"d{d}_r{r}": n for r, n in local_degrees.items()})

    return {
        "protocol": "fprd-relative-degree-exhaustive.v1",
        "domain": {
            "dimensions": [1, 2],
            "entry_set": list(VALUES),
            "alphabet": ["a", "b"],
            "enumeration": "all lambda, gamma, M_a, M_b over the entry set",
        },
        "oracle": "direct definition test for every r<d using every suffix word of length <d",
        "result": "pass",
        "total_representations": total,
        "by_dimension": by_dimension,
        "degree_counts": degrees,
        "mismatches": 0,
        "nonunique_candidate_sets": 0,
        "transcript_sha256": transcript.hexdigest(),
        "python": sys.version.split()[0],
    }


if __name__ == "__main__":
    print(json.dumps(run(), indent=2, sort_keys=True))