verify-flat-koszul-blocks.py

#!/usr/bin/env python3
"""Dependency-free exact calibration for flatwise GMA Koszul blocks.

For the graded Moebius algebra A(M), the basis vectors of the Koszul
complex are pairs (F, J), with F a flat and J an exterior subset.  The
script verifies that

    X(F, J) = closure(F union J)

is preserved by every nonzero differential term.  It then computes H_2
over Q both globally and as the sum of its X-labelled blocks.

The examples are deliberately structural controls:

* U(2,6), PG(2,2), and AG(2,3) are Sylvester examples used in the source;
* M(K4) is quadratic but not Sylvester, and should have beta_(2,4) = 3.

All matrix ranks use sparse Gaussian elimination with Fraction, so the
reported values are exact over Q and require only the Python standard
library.
"""

from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction
from itertools import combinations
from typing import Callable, Iterable


Subset = frozenset[int]
BasisVector = tuple[Subset, tuple[int, ...]]


def gf_rank(rows: list[list[int]], prime: int) -> int:
    """Rank of a small matrix over the prime field GF(prime)."""
    if not rows:
        return 0
    matrix = [[entry % prime for entry in row] for row in rows]
    row_count, column_count = len(matrix), len(matrix[0])
    rank = 0
    for column in range(column_count):
        pivot = next(
            (row for row in range(rank, row_count) if matrix[row][column]),
            None,
        )
        if pivot is None:
            continue
        matrix[rank], matrix[pivot] = matrix[pivot], matrix[rank]
        inverse = pow(matrix[rank][column], -1, prime)
        matrix[rank] = [(inverse * entry) % prime for entry in matrix[rank]]
        for row in range(row_count):
            if row == rank or not matrix[row][column]:
                continue
            factor = matrix[row][column]
            matrix[row] = [
                (entry - factor * pivot_entry) % prime
                for entry, pivot_entry in zip(matrix[row], matrix[rank])
            ]
        rank += 1
        if rank == row_count:
            break
    return rank


def rational_sparse_rank(rows: Iterable[dict[int, int]]) -> int:
    """Exact row rank over Q for a sparse integer matrix."""
    pivots: dict[int, dict[int, Fraction]] = {}
    for integer_row in rows:
        row = {
            column: Fraction(value)
            for column, value in integer_row.items()
            if value
        }
        while row:
            column = min(row)
            if column not in pivots:
                leading = row[column]
                row = {key: value / leading for key, value in row.items()}
                pivots[column] = row
                break
            factor = row[column]
            for key, value in pivots[column].items():
                updated = row.get(key, Fraction(0)) - factor * value
                if updated:
                    row[key] = updated
                else:
                    row.pop(key, None)
    return len(pivots)


@dataclass(frozen=True)
class Matroid:
    name: str
    size: int
    rank_function: Callable[[Subset], int]

    def rank(self, subset: Subset) -> int:
        return self.rank_function(subset)

    def closure(self, subset: Subset) -> Subset:
        rank = self.rank(subset)
        return frozenset(
            element
            for element in range(self.size)
            if self.rank(subset | {element}) == rank
        )

    def flats(self) -> list[Subset]:
        all_flats = {
            self.closure(frozenset(subset))
            for length in range(self.size + 1)
            for subset in combinations(range(self.size), length)
        }
        return sorted(all_flats, key=lambda flat: (self.rank(flat), tuple(flat)))


def uniform_rank_two(size: int) -> Matroid:
    return Matroid(f"U(2,{size})", size, lambda subset: min(2, len(subset)))


def vector_matroid(
    name: str, columns: list[tuple[int, ...]], prime: int
) -> Matroid:
    def rank_function(subset: Subset) -> int:
        rows = [
            [columns[column][coordinate] for column in subset]
            for coordinate in range(len(columns[0]))
        ]
        return gf_rank(rows, prime)

    return Matroid(name, len(columns), rank_function)


def fano_plane() -> Matroid:
    columns = [tuple((value >> bit) & 1 for bit in range(3)) for value in range(1, 8)]
    return vector_matroid("PG(2,2)", columns, 2)


def affine_plane_order_three() -> Matroid:
    columns = [(x, y, 1) for x in range(3) for y in range(3)]
    return vector_matroid("AG(2,3)", columns, 3)


def complete_graph_four() -> Matroid:
    edges = list(combinations(range(4), 2))
    columns = []
    for left, right in edges:
        column = tuple(
            int(vertex == left) + int(vertex == right) for vertex in range(3)
        )
        columns.append(column)
    return vector_matroid("M(K4)", columns, 2)


def koszul_basis(
    matroid: Matroid,
    flats_by_rank: dict[int, list[Subset]],
    homological_degree: int,
    internal_degree: int,
) -> list[BasisVector]:
    coefficient_degree = internal_degree - homological_degree
    if coefficient_degree < 0:
        return []
    return [
        (flat, wedge)
        for flat in flats_by_rank.get(coefficient_degree, [])
        for wedge in combinations(range(matroid.size), homological_degree)
    ]


def total_flat(matroid: Matroid, vector: BasisVector) -> Subset:
    flat, wedge = vector
    return matroid.closure(flat | frozenset(wedge))


def differential_rows(
    matroid: Matroid,
    source: list[BasisVector],
    target: list[BasisVector],
) -> tuple[list[dict[int, int]], list[dict[int, int]]]:
    """Return a differential as sparse rows and sparse columns."""
    target_index = {vector: index for index, vector in enumerate(target)}
    rows: list[dict[int, int]] = [dict() for _ in target]
    columns: list[dict[int, int]] = [dict() for _ in source]
    for source_index, (flat, wedge) in enumerate(source):
        flat_rank = matroid.rank(flat)
        source_label = total_flat(matroid, (flat, wedge))
        for position, element in enumerate(wedge):
            joined = matroid.closure(flat | {element})
            if matroid.rank(joined) != flat_rank + 1:
                continue
            output = (joined, wedge[:position] + wedge[position + 1 :])
            assert total_flat(matroid, output) == source_label
            target_row = target_index[output]
            value = -1 if position % 2 else 1
            rows[target_row][source_index] = value
            columns[source_index][target_row] = value
    return rows, columns


def assert_chain_zero(
    first_columns: list[dict[int, int]],
    second_columns: list[dict[int, int]],
) -> None:
    """Assert that the product of two consecutive differentials is zero."""
    for source_column in second_columns:
        composite: dict[int, int] = {}
        for middle_index, second_value in source_column.items():
            for target_index, first_value in first_columns[middle_index].items():
                composite[target_index] = (
                    composite.get(target_index, 0)
                    + second_value * first_value
                )
        assert not {index: value for index, value in composite.items() if value}


def restricted_rank(
    rows: list[dict[int, int]],
    target_labels: list[Subset],
    source_labels: list[Subset],
    label: Subset,
) -> int:
    selected_columns = [
        index for index, source_label in enumerate(source_labels) if source_label == label
    ]
    column_map = {old: new for new, old in enumerate(selected_columns)}
    selected_rows = []
    for row, target_label in zip(rows, target_labels):
        if target_label != label:
            continue
        selected_rows.append(
            {column_map[column]: value for column, value in row.items()}
        )
    return rational_sparse_rank(selected_rows)


def three_circuit_counts(matroid: Matroid) -> list[int]:
    counts = []
    for left, right in combinations(range(matroid.size), 2):
        counts.append(
            sum(
                1
                for third in range(matroid.size)
                if third not in {left, right}
                and matroid.rank(frozenset({left, right, third})) == 2
            )
        )
    return counts


def audit(matroid: Matroid) -> list[str]:
    flats = matroid.flats()
    flats_by_rank: dict[int, list[Subset]] = {}
    for flat in flats:
        flats_by_rank.setdefault(matroid.rank(flat), []).append(flat)
    matroid_rank = matroid.rank(frozenset(range(matroid.size)))
    circuit_counts = three_circuit_counts(matroid)
    sylvester = min(circuit_counts) > 0
    output = [
        (
            f"{matroid.name}: n={matroid.size}, rank={matroid_rank}, "
            f"flats={len(flats)}, Sylvester={sylvester}, "
            f"third-points-per-pair=[{min(circuit_counts)},{max(circuit_counts)}]"
        )
    ]

    for internal_degree in range(3, matroid_rank + 3):
        k1 = koszul_basis(matroid, flats_by_rank, 1, internal_degree)
        k2 = koszul_basis(matroid, flats_by_rank, 2, internal_degree)
        k3 = koszul_basis(matroid, flats_by_rank, 3, internal_degree)
        d2_rows, d2_columns = differential_rows(matroid, k2, k1)
        d3_rows, d3_columns = differential_rows(matroid, k3, k2)
        assert_chain_zero(d2_columns, d3_columns)

        rank_d2 = rational_sparse_rank(d2_rows)
        rank_d3 = rational_sparse_rank(d3_rows)
        h2 = len(k2) - rank_d2 - rank_d3

        labels1 = [total_flat(matroid, vector) for vector in k1]
        labels2 = [total_flat(matroid, vector) for vector in k2]
        labels3 = [total_flat(matroid, vector) for vector in k3]
        all_labels = sorted(
            set(labels1 + labels2 + labels3),
            key=lambda flat: (matroid.rank(flat), tuple(flat)),
        )
        block_h2 = 0
        nonzero_blocks = 0
        for label in all_labels:
            block_k2 = sum(block_label == label for block_label in labels2)
            block_rank_d2 = restricted_rank(
                d2_rows, labels1, labels2, label
            )
            block_rank_d3 = restricted_rank(
                d3_rows, labels2, labels3, label
            )
            block_dimension = block_k2 - block_rank_d2 - block_rank_d3
            assert block_dimension >= 0
            block_h2 += block_dimension
            nonzero_blocks += int(block_dimension > 0)
        assert h2 == block_h2
        output.append(
            f"  degree {internal_degree}: dim K2={len(k2)}, "
            f"rank(d2)={rank_d2}, rank(d3)={rank_d3}, "
            f"H2={h2}, block-sum={block_h2}, "
            f"nonzero-blocks={nonzero_blocks}"
        )
    return output


if __name__ == "__main__":
    examples = [
        uniform_rank_two(6),
        fano_plane(),
        affine_plane_order_three(),
        complete_graph_four(),
    ]
    for position, example in enumerate(examples):
        if position:
            print()
        print("\n".join(audit(example)))