diagonal_embedding_search.py

#!/usr/bin/env python3
"""Exhaustive search for arrival-compatible diagonal embeddings of 2D posets.

Elements are labeled by their arrival order 0,1,...,n-1.  A pair of
permutations (lx, ly) is a 2-realizer of a naturally labeled poset precisely
when their inversion sets (relative to arrival order) are disjoint.  Its
comparabilities are the pairs inverted by neither permutation.

For a fixed realizer we seek positive gaps along lx and ly such that x_i+y_i
is nondecreasing in arrival order.  Homogeneity lets us normalize every gap
to be at least one, leaving an ordinary linear feasibility problem.
"""

from __future__ import annotations

import argparse
from dataclasses import dataclass
from fractions import Fraction
from itertools import combinations, permutations

import numpy as np
from scipy.optimize import linprog


Permutation = tuple[int, ...]


@dataclass(frozen=True)
class PermData:
    order: Permutation
    rank: tuple[int, ...]
    inversions: int


def pair_index(n: int) -> dict[tuple[int, int], int]:
    return {pair: k for k, pair in enumerate(combinations(range(n), 2))}


def permutation_data(n: int) -> list[PermData]:
    index = pair_index(n)
    result = []
    for order in permutations(range(n)):
        rank = [0] * n
        for i, element in enumerate(order):
            rank[element] = i
        inversions = sum(
            1 << index[i, j]
            for i, j in combinations(range(n), 2)
            if rank[i] > rank[j]
        )
        result.append(PermData(order, tuple(rank), inversions))
    return result


def difference_row(rank: tuple[int, ...], before: int, after: int) -> list[int]:
    """Coefficients of coord(after)-coord(before) in consecutive-order gaps."""

    lo, hi = sorted((rank[before], rank[after]))
    sign = 1 if rank[before] < rank[after] else -1
    return [sign if lo <= k < hi else 0 for k in range(len(rank) - 1)]


def diagonal_matrix(px: PermData, py: PermData) -> np.ndarray:
    n = len(px.order)
    rows = []
    for before in range(n - 1):
        rows.append(
            difference_row(px.rank, before, before + 1)
            + difference_row(py.rank, before, before + 1)
        )
    return np.asarray(rows, dtype=float).reshape((n - 1, 2 * (n - 1)))


def feasible_realizer(
    px: PermData,
    py: PermData,
    *,
    strict_diagonal: bool = False,
    return_solution: bool = False,
):
    """Test existence of positive coordinate gaps by homogeneous LP.

    Non-strict diagonal order uses Dg >= 0 with g >= 1.
    Strict diagonal order is equivalently Dg >= 1 with g >= 1, because any
    strictly feasible homogeneous solution can be scaled to that normalization.
    """

    dmat = diagonal_matrix(px, py)
    if not dmat.shape[1]:
        return (True, tuple()) if return_solution else True
    rhs = -np.ones(dmat.shape[0]) if strict_diagonal else np.zeros(dmat.shape[0])
    result = linprog(
        np.zeros(dmat.shape[1]),
        A_ub=-dmat,
        b_ub=rhs,
        bounds=[(1, None)] * dmat.shape[1],
        method="highs",
    )
    if return_solution:
        if not result.success:
            return False, None
        # HiGHS is used only to propose a rational point.  The returned point
        # is reconstructed and checked with exact Fraction arithmetic before
        # it is accepted as a positive certificate.
        gaps = tuple(Fraction(float(value)).limit_denominator(1_000_000) for value in result.x)
        exact_dmat = diagonal_matrix(px, py).astype(int).tolist()
        target = 1 if strict_diagonal else 0
        assert all(gap >= 1 for gap in gaps)
        assert all(
            sum(Fraction(coeff) * gap for coeff, gap in zip(row, gaps)) >= target
            for row in exact_dmat
        )
        return True, gaps
    return result.success


def coordinates_from_gaps(
    order: Permutation, gaps: tuple[Fraction, ...]
) -> tuple[Fraction, ...]:
    coordinates = [Fraction(0)] * len(order)
    for k, gap in enumerate(gaps):
        coordinates[order[k + 1]] = coordinates[order[k]] + gap
    return tuple(coordinates)


def check_exact_strictification(
    px: PermData, py: PermData, gaps: tuple[Fraction, ...]
) -> None:
    """Turn an exact weak-diagonal witness into an exact strict one.

    All gaps have been normalized to at least one.  Adding i/(2n) to both
    coordinates of arrival element i changes any old coordinate difference
    by less than 1/2, hence preserves both realizer orders, while adding i/n
    to the diagonal sum.  Thus every weak arrival inequality becomes strict.
    """

    n = len(px.order)
    if n <= 1:
        return
    x = coordinates_from_gaps(px.order, gaps[: n - 1])
    y = coordinates_from_gaps(py.order, gaps[n - 1 :])
    epsilon = Fraction(1, 2 * n)
    xx = tuple(value + epsilon * i for i, value in enumerate(x))
    yy = tuple(value + epsilon * i for i, value in enumerate(y))
    assert all(xx[px.order[k]] < xx[px.order[k + 1]] for k in range(n - 1))
    assert all(yy[py.order[k]] < yy[py.order[k + 1]] for k in range(n - 1))
    sums = tuple(xx[i] + yy[i] for i in range(n))
    assert all(sums[i] < sums[i + 1] for i in range(n - 1))


def all_natural_presentations(n: int, *, strict_diagonal: bool = False):
    """Return all naturally labeled 2D posets and those with an embedding."""

    pdata = permutation_data(n)
    full = (1 << (n * (n - 1) // 2)) - 1
    all_posets: set[int] = set()
    good_posets: set[int] = set()
    witnesses: dict[int, tuple[Permutation, Permutation, tuple[Fraction, ...]]] = {}
    pair_count = 0
    lp_count = 0

    for px in pdata:
        for py in pdata:
            if px.inversions & py.inversions:
                continue
            pair_count += 1
            poset = full ^ (px.inversions | py.inversions)
            all_posets.add(poset)
            if poset in good_posets:
                continue
            lp_count += 1
            feasible, gaps = feasible_realizer(
                px, py, strict_diagonal=strict_diagonal, return_solution=True
            )
            if feasible:
                if not strict_diagonal:
                    check_exact_strictification(px, py, gaps)
                good_posets.add(poset)
                witnesses[poset] = (px.order, py.order, tuple(gaps))

    return all_posets, good_posets, witnesses, pair_count, lp_count


def audit_six_point_example() -> None:
    original_lx = (0, 1, 2, 3, 4, 5)
    original_rank_y = (2, 5, 4, 1, 0, 3)
    original_ly = tuple(sorted(range(6), key=original_rank_y.__getitem__))
    audit_lx = (0, 1, 2, 4, 3, 5)
    audit_ly = (3, 4, 0, 5, 2, 1)
    x = (-9, -8, -1, 1, 0, 2)
    y = (-1, 5, 4, -12, -2, 0)
    arrival = (3, 0, 1, 4, 5, 2)

    def relation(lx: Permutation, ly: Permutation):
        rx = {e: i for i, e in enumerate(lx)}
        ry = {e: i for i, e in enumerate(ly)}
        return {
            (a, b)
            for a in range(6)
            for b in range(6)
            if a != b and rx[a] < rx[b] and ry[a] < ry[b]
        }

    assert original_ly == (4, 3, 0, 5, 2, 1)
    assert relation(original_lx, original_ly) == relation(audit_lx, audit_ly)
    assert all(x[audit_lx[i]] < x[audit_lx[i + 1]] for i in range(5))
    assert all(y[audit_ly[i]] < y[audit_ly[i + 1]] for i in range(5))
    sums = tuple(x[e] + y[e] for e in arrival)
    assert all(sums[i] < sums[i + 1] for i in range(5))
    assert sums == (-11, -10, -3, -2, 2, 3)
    print("audit six-point alternate realizer: verified")
    print(f"  original Ly order: {original_ly}")
    print(f"  arrival diagonal sums: {sums}")


def relations_of_mask(n: int, mask: int) -> tuple[tuple[int, int], ...]:
    return tuple(
        pair for bit, pair in enumerate(combinations(range(n), 2)) if mask >> bit & 1
    )


def width_of_mask(n: int, mask: int) -> int:
    relations = set(relations_of_mask(n, mask))
    return max(
        len(subset)
        for size in range(n + 1)
        for subset in combinations(range(n), size)
        if all(pair not in relations for pair in combinations(subset, 2))
    )


def realizers_of_mask(n: int, mask: int) -> list[tuple[PermData, PermData]]:
    """Enumerate all ordered 2-realizers of one naturally labeled poset."""

    pdata = permutation_data(n)
    full = (1 << (n * (n - 1) // 2)) - 1
    return [
        (px, py)
        for px in pdata
        for py in pdata
        if not px.inversions & py.inversions
        and full ^ (px.inversions | py.inversions) == mask
    ]


def check_seven_point_farkas(mask: int) -> None:
    """Verify the exact certificate for the first seven-point obstruction."""

    assert relations_of_mask(7, mask) == (
        (0, 1),
        (0, 5),
        (0, 6),
        (1, 5),
        (2, 3),
        (2, 5),
        (2, 6),
        (4, 5),
    )
    assert width_of_mask(7, mask) == 4
    realizers = realizers_of_mask(7, mask)
    assert len(realizers) == 2
    expected = {
        (
            (2, 3, 0, 6, 1, 4, 5),
            (4, 0, 1, 2, 5, 6, 3),
        ),
        (
            (4, 0, 1, 2, 5, 6, 3),
            (2, 3, 0, 6, 1, 4, 5),
        ),
    }
    assert {(px.order, py.order) for px, py in realizers} == expected

    # For the displayed order of coordinates, rows 1,3,5 of D correspond to
    # arrival inequalities 1->2, 3->4, and 5->6.  Their exact sum is strictly
    # negative on every positive gap vector.
    px, py = next(
        (px, py)
        for px, py in realizers
        if px.order == (2, 3, 0, 6, 1, 4, 5)
    )
    dmat = diagonal_matrix(px, py).astype(int)
    certificate = dmat[1] + dmat[3] + dmat[5]
    assert tuple(certificate) == (
        -1,
        0,
        0,
        -1,
        0,
        -1,
        -1,
        -1,
        0,
        -1,
        0,
        -1,
    )
    assert not feasible_realizer(px, py, strict_diagonal=False)
    assert not feasible_realizer(py, px, strict_diagonal=False)
    print("seven-point obstruction: exact certificate verified")
    print(f"  relations: {relations_of_mask(7, mask)}")
    print(f"  unique unordered realizer: {px.order} / {py.order}")
    print("  Farkas rows: arrival steps 1->2, 3->4, 5->6")
    print("  sum: -(A0+A3+A5+B0+B1+B3+B5) >= 0, impossible")


def check_dual_seven_point_farkas(mask: int) -> None:
    """Verify the arrival-reversed/order-dual obstruction."""

    assert relations_of_mask(7, mask) == (
        (0, 4),
        (0, 6),
        (1, 2),
        (1, 4),
        (1, 5),
        (1, 6),
        (3, 4),
        (5, 6),
    )
    assert width_of_mask(7, mask) == 4
    realizers = realizers_of_mask(7, mask)
    assert len(realizers) == 2
    expected = {
        (
            (1, 2, 5, 0, 6, 3, 4),
            (3, 0, 1, 4, 5, 6, 2),
        ),
        (
            (3, 0, 1, 4, 5, 6, 2),
            (1, 2, 5, 0, 6, 3, 4),
        ),
    }
    assert {(px.order, py.order) for px, py in realizers} == expected
    px, py = next(
        (px, py)
        for px, py in realizers
        if px.order == (1, 2, 5, 0, 6, 3, 4)
    )
    dmat = diagonal_matrix(px, py).astype(int)
    certificate = dmat[0] + dmat[2] + dmat[4]
    assert tuple(certificate) == (
        -1,
        0,
        -1,
        0,
        0,
        -1,
        -1,
        0,
        -1,
        0,
        -1,
        -1,
    )
    assert not feasible_realizer(px, py, strict_diagonal=False)
    assert not feasible_realizer(py, px, strict_diagonal=False)
    print("dual seven-point obstruction: exact certificate verified")
    print(f"  relations: {relations_of_mask(7, mask)}")
    print(f"  unique unordered realizer: {px.order} / {py.order}")
    print("  Farkas rows: arrival steps 0->1, 2->3, 4->5")
    print("  sum: -(A0+A2+A5+B0+B2+B4+B5) >= 0, impossible")


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--max-n", type=int, default=6)
    parser.add_argument("--strict", action="store_true")
    args = parser.parse_args()
    audit_six_point_example()
    for n in range(1, args.max_n + 1):
        all_posets, good_posets, _, pair_count, lp_count = all_natural_presentations(
            n, strict_diagonal=args.strict
        )
        print(
            f"n={n}: realizer pairs={pair_count}, presentations={len(all_posets)}, "
            f"LPs={lp_count}, diagonal={len(good_posets)}, "
            f"obstructions={len(all_posets-good_posets)}"
        )
        if all_posets != good_posets:
            bad = sorted(all_posets - good_posets)
            print("  obstruction masks:", bad)
            if n == 7 and 289329 in bad:
                check_seven_point_farkas(289329)
            if n == 7 and 1083240 in bad:
                check_dual_seven_point_farkas(1083240)
            break


if __name__ == "__main__":
    main()