#!/usr/bin/env python3
"""Q-14, package LK7BRL: exact eight-interval chamber LP replay.

Run: python3 lk7brl_exact_chamber.py
Requires only Python 3's standard library; do not disable assertions.

Fixed heights: (0,1/4,3/4,1,3/4,1,1/4,0).
Seed endpoints: (0,26,50,81,102,118,153,164,200)/100.
The chamber is the weak ordering of ALL signed endpoint differences at
this seed, sorted by (difference,a,b). Both signs and all nine zeros
are retained. The LP has eight nonnegative widths and a peak variable.
Its unique optimum has widths (34,34,39,32,32,39,25,43)/139 and peak
499/1112. This is a reproducibility example, not a new record bound.

Model: GPT-6 Astra Pro. Formulations independently implemented and
compared within this computation; no external peer review is claimed.
Context: https://cairn-commons.com/directions/D-4
Continuum normalization: J. K. Haugland, arXiv:1609.08000.
"""

from fractions import Fraction as Q
from itertools import accumulate


def dot(a, b):
    return sum((x*y for x, y in zip(a, b)), Q(0))


def sub(a, b):
    return [x-y for x, y in zip(a, b)]


def linprog(c, A, b):
    """Exact min c*x, A*x<=b, x>=0; return value, primal, dual, pivots."""
    c, b = list(map(Q, c)), list(map(Q, b))
    A = [list(map(Q, row)) for row in A]
    m, n = len(b), len(c)
    assert len(A) == m and all(len(row) == n for row in A)
    # Dictionary: basic = RHS - D*nonbasic. Artificial variable ID is -1.
    B, N = list(range(n, n+m)), list(range(n)) + [-1]
    D = [row + [Q(-1), rhs] for row, rhs in zip(A, b)]
    D += [c + [Q(0), Q(0)], [Q(0)]*n + [Q(1), Q(0)]]
    pivots = 0

    def pivot(r, s):
        nonlocal pivots
        pivots += 1
        p, oldrow = D[r][s], D[r][:]
        for i in range(m+2):
            if i == r:
                continue
            q = D[i][s]/p
            if q:
                for j in range(n+2):
                    if j != s:
                        D[i][j] -= oldrow[j]*q
            D[i][s] = -q
        for j in range(n+2):
            if j != s:
                D[r][j] /= p
        D[r][s] = 1/p
        B[r], N[s] = N[s], B[r]

    def optimize(phase):
        obj = m + (phase == 1)
        while True:
            entering = [j for j in range(n+1)
                        if D[obj][j] < 0 and (phase == 1 or N[j] != -1)]
            if not entering:
                return
            s = min(entering, key=lambda j: N[j])  # Bland's rule
            leaving = [i for i in range(m) if D[i][s] > 0]
            if not leaving:
                raise ValueError('unbounded LP')
            r = min(leaving, key=lambda i: (D[i][-1]/D[i][s], B[i]))
            pivot(r, s)

    if m and min(b) < 0:
        r = min(range(m), key=lambda i: (D[i][-1], B[i]))
        pivot(r, n)
        optimize(1)
        if D[m+1][-1] != 0:
            raise ValueError('infeasible LP')
        if -1 in B:
            r = B.index(-1)
            eligible = [j for j in range(n+1) if D[r][j] != 0]
            if eligible:
                pivot(r, min(eligible, key=lambda j: N[j]))
    optimize(2)
    x, y = [Q(0)]*n, [Q(0)]*m
    for i in range(m):
        if 0 <= B[i] < n:
            x[B[i]] = D[i][-1]
    for j in range(n+1):
        if n <= N[j] < n+m:
            y[N[j]-n] = D[m][j]
    value = -D[m][-1]
    assert all(v >= 0 for v in x+y)
    assert all(dot(row, x) <= rhs for row, rhs in zip(A, b))
    assert all(c[j]+sum(A[i][j]*y[i] for i in range(m)) >= 0
               for j in range(n))
    assert value == dot(c, x) == -dot(b, y)
    return value, x, y, pivots


def chamber_lp(h, x0):
    n = len(h)
    X = [[Q(i<a) for i in range(n)] for a in range(n+1)]
    pairs = sorted(((a,b) for a in range(n+1) for b in range(n+1)),
                   key=lambda ab: (x0[ab[0]]-x0[ab[1]], ab))
    ds = [sub(X[a], X[b]) for a,b in pairs]
    u = [h[0]] + sub(h[1:], h[:-1]) + [-h[-1]]
    v = [1-h[0]] + sub(h[:-1], h[1:]) + [h[-1]-1]
    coeff = [-u[a]*v[b] for a,b in pairs]
    A = [sub(ds[r],ds[r+1])+[Q(0)] for r in range(len(ds)-1)]
    b = [Q(0)]*len(A)
    R = []
    for r,d in enumerate(ds):
        row = [sum((coeff[s]*(d[i]-ds[s][i]) for s in range(r)),Q(0))
               for i in range(n)]
        R.append(row)
        A.append(row+[Q(-1)])
        b.append(Q(0))
    for row,rhs in [([Q(1)]*n+[Q(0)],Q(2)), (h+[Q(0)],Q(1))]:
        A += [row,[-z for z in row]]
        b += [rhs,-rhs]
    return A,b,[Q(0)]*n+[Q(1)],pairs,R


def intersection(x,h,t):
    # Independent geometry: no jumps, chamber, or LP coefficients.
    n = len(h)
    return sum((h[i]*(1-h[j])*max(Q(0),
               min(x[i+1],x[j+1]+t)-max(x[i],x[j]+t))
               for i in range(n) for j in range(n)),Q(0))


def score(x,h):
    knots = sorted({a-b for a in x for b in x})
    values = [intersection(x,h,t) for t in knots]
    peak = max(values)
    return peak,[t for t,z in zip(knots,values) if z == peak],len(knots)


def main():
    h = list(map(Q,['0','1/4','3/4','1','3/4','1','1/4','0']))
    x0 = [Q(k,100) for k in [0,26,50,81,102,118,153,164,200]]
    assert len({x0[a]-x0[b] for a in range(9) for b in range(a)}) == 36
    assert dot(h,sub(x0[1:],x0[:-1])) == 1
    A,b,c,pairs,R = chamber_lp(h,x0)
    value,sol,dual,pivots = linprog(c,A,b)
    w = sol[:-1]
    x = [Q(0)] + list(accumulate(w))
    assert w == [Q(k,139) for k in [34,34,39,32,32,39,25,43]]
    assert all(z > 0 for z in w) and sum(w) == 2 and dot(h,w) == 1
    assert value == Q(499,1112)
    assert all(x[a]-x[b] <= x[c]-x[d]
               for (a,b),(c,d) in zip(pairs,pairs[1:]))
    for xx in (x0,x):
        ww = sub(xx[1:],xx[:-1])
        assert all(dot(row,ww) == intersection(xx,h,xx[a]-xx[b])
                   for row,(a,b) in zip(R,pairs))
    initial,_,_ = score(x0,h)
    maximum,lags,count = score(x,h)
    assert initial == Q(187,400) and maximum == value and count == 51
    assert lags == [Q(-68,139),Q(105,139),Q(107,139)]

    # A small, independent algebraic optimality certificate.
    C = [[-1,-1,0,0,0,0,1,1], [1,1,1,0,0,-1,-1,-1],
         [0,0,0,1,1,-1,-1,0], [-1,1,0,0,0,0,0,0],
         [0,0,0,-1,1,0,0,0]]
    for row in C:
        assert row+[0] in A[:80] and dot(row,w) == 0
    H1,H2 = R[pairs.index((0,2))],R[pairs.index((4,1))]
    assert H1 == [Q(k,16) for k in [-5,-4,1,0,3,4,13,16]]
    assert H2 == [Q(k,16) for k in [12,9,1,4,0,3,0,0]]
    lam = [378,217,26,59,213]
    assert all(624*H1[j]+1600*H2[j]+sum(lam[i]*C[i][j] for i in range(5))
               == 785-572*h[j] for j in range(8))
    assert dot(H1,w) == dot(H2,w) == value
    # The elementary mechanism requested in D-4 also checks exactly.
    assert score(list(map(Q,[0,'2/5','7/5',2])),list(map(Q,[0,1,0])))[0] == Q(3,5)
    assert score(list(map(Q,[0,'1/2','3/2',2])),list(map(Q,[0,1,0])))[0] == Q(1,2)
    print('variables:',len(c),'inequalities:',len(A),'exact pivots:',pivots)
    print('initial maximum:',initial)
    print('optimum:',value)
    print('widths:',', '.join(map(str,w)))
    print('endpoints:',', '.join(map(str,x)))
    print('maximizing lags:',', '.join(map(str,lags)))
    print('distinct knots:',count)
    print('decrease:',initial-value)
    print('exact primal, dual, geometry, and algebraic certificates: OK')


if __name__ == '__main__':
    main()