# tropmat.py - tropical (min-plus / max-plus) matrix algebra
# Copyright (c) 2026 Carey Witkov. MIT License (see LICENSE).
#
# Provided "as is", without warranty of any kind. No guarantee of
# correctness or fitness for any particular purpose; use at your own risk.
# Verify results against your own model before relying on them.
#
# Target: NumWorks (MicroPython 1.17, Python 3.4-compatible, single-precision float)
# Pure Python, no imports. Matrices are lists of lists of floats.
# zero = additive identity of the semiring (INF for min-plus, NEG for max-plus).
# one  = multiplicative identity = 0.0 in both.
#
# Notes for NumWorks single precision:
#   float('inf') and float('-inf') are supported on epsilon.
#   EPS is loosened to 1e-4 because 32-bit floats carry ~7 significant digits,
#   so a 1e-6 tolerance can give false negatives in the periodicity / critical
#   tests. Use exact integer weights where you can.

INF = float('inf')
NEG = float('-inf')
EPS = 1e-4            # tolerance for critical-cycle / periodicity tests (single precision)

class Semiring:
    def __init__(self, oplus, otimes, zero, one, name):
        self.oplus = oplus
        self.otimes = otimes
        self.zero = zero
        self.one = one
        self.name = name

def _min(a, b):
    return a if a < b else b
def _max(a, b):
    return a if a > b else b
def _add(a, b):
    return a + b

MINPLUS = Semiring(_min, _add, INF, 0.0, "min-plus")   # shortest paths
MAXPLUS = Semiring(_max, _add, NEG, 0.0, "max-plus")    # longest / critical paths

# ---------- constructors ----------

def dims(M):
    return (len(M), len(M[0]) if M else 0)

def zeros(n, m, S):
    return [[S.zero for j in range(m)] for i in range(n)]

def identity(n, S):
    I = zeros(n, n, S)
    for i in range(n):
        I[i][i] = S.one
    return I

def copy(M):
    return [row[:] for row in M]

def from_edges(n, edges, S):
    M = zeros(n, n, S)
    for e in edges:
        i = e[0]; j = e[1]; w = e[2]
        M[i][j] = S.oplus(M[i][j], w)   # keep best if duplicates given
    return M

def _num(tok):
    tok = tok.strip().lower()
    if tok == "inf" or tok == "+inf":
        return INF
    if tok == "-inf":
        return NEG
    return float(tok)

def from_native(s):
    t = s.replace("], [", "];[").replace("],[", "];[")
    t = t.replace("] [", "];[").replace("][", "];[")
    t = t.replace("[", " ").replace("]", " ")
    t = t.replace("\n", ";")
    rows = []
    for part in t.split(";"):
        part = part.strip()
        if part == "":
            continue
        part = part.replace(",", " ")
        row = []
        for tok in part.split():
            row.append(_num(tok))
        if row:
            rows.append(row)
    return rows

# ---------- core operations ----------

def madd(A, B, S):
    n = len(A); m = len(A[0])
    return [[S.oplus(A[i][j], B[i][j]) for j in range(m)] for i in range(n)]

def mmul(A, B, S):
    n = len(A); p = len(B); m = len(B[0])
    C = zeros(n, m, S)
    for i in range(n):
        Ai = A[i]; Ci = C[i]
        for k in range(p):
            a = Ai[k]
            if a == S.zero:          # absorbing: skip
                continue
            Bk = B[k]
            for j in range(m):
                Ci[j] = S.oplus(Ci[j], S.otimes(a, Bk[j]))
    return C

def mvmul(A, x, S):
    # matrix (X) vector product: y = A (X) x. Use to iterate x(k) = A (X) x(k-1).
    n = len(A); p = len(x)
    y = [S.zero] * n
    for i in range(n):
        Ai = A[i]; acc = S.zero
        for k in range(p):
            a = Ai[k]
            if a == S.zero:
                continue
            acc = S.oplus(acc, S.otimes(a, x[k]))
        y[i] = acc
    return y

def mpow(A, e, S):
    n = len(A)
    R = identity(n, S)
    P = copy(A)
    while e > 0:
        if e & 1:
            R = mmul(R, P, S)
        e = e >> 1
        if e > 0:
            P = mmul(P, P, S)
    return R

def mtrace(A, S):
    # semiring trace: (+)-sum of the diagonal (max of diag in max-plus).
    t = S.zero
    for i in range(len(A)):
        t = S.oplus(t, A[i][i])
    return t

def closure(A, S):
    # Kleene star A* = I (+) A (+) A^2 (+) ...  (finite when no positive cycle for
    # max-plus / no negative cycle for min-plus). Floyd-Warshall style.
    n = len(A)
    D = copy(A)
    for i in range(n):
        D[i][i] = S.oplus(D[i][i], S.one)   # allow the empty (zero-cost) path
    for k in range(n):
        Dk = D[k]
        for i in range(n):
            Di = D[i]
            aik = Di[k]
            if aik == S.zero:
                continue
            for j in range(n):
                Di[j] = S.oplus(Di[j], S.otimes(aik, Dk[j]))
    return D

# ---------- convenience wrappers ----------

def shortest_paths(A):
    return closure(A, MINPLUS)
def critical_paths(A):
    return closure(A, MAXPLUS)
def minmul(A, B):
    return mmul(A, B, MINPLUS)
def maxmul(A, B):
    return mmul(A, B, MAXPLUS)

# ---------- display ----------

def _f(x):
    if x == INF:
        return "inf"
    if x == NEG:
        return "-inf"
    if x == int(x):          # safe: x is finite here
        return str(int(x))
    return str(x)

def show(M, label=None):
    if label is not None:
        print(label)
    for row in M:
        print(" ".join([_f(v) for v in row]))

def showv(v, label=None):
    if label is not None:
        print(label)
    print(" ".join([_f(x) for x in v]))

def lit(M):
    def cell(x):
        if x == INF:
            return "float('inf')"
        if x == NEG:
            return "float('-inf')"
        if x == int(x):
            return str(int(x))
        return str(x)
    rows = []
    for row in M:
        rows.append("[" + ",".join([cell(v) for v in row]) + "]")
    return "[" + ",".join(rows) + "]"

# ---------- spectral (max cycle mean = max-plus eigenvalue) ----------

def max_cycle_mean(A):
    # Karp's algorithm. lambda = max over cycles of (cycle weight / cycle length).
    n = len(A)
    if n == 0:
        return NEG
    D = [[NEG] * n for k in range(n + 1)]
    for v in range(n):
        D[0][v] = 0.0
    for k in range(n):
        Dk = D[k]; Dk1 = D[k + 1]
        for u in range(n):
            du = Dk[u]
            if du == NEG:
                continue
            Au = A[u]
            for v in range(n):
                w = Au[v]
                if w == NEG:
                    continue
                s = du + w
                if s > Dk1[v]:
                    Dk1[v] = s
    best = NEG
    for v in range(n):
        dnv = D[n][v]
        if dnv == NEG:
            continue
        worst = INF                      # min over k; finite (k=0 term exists)
        for k in range(n):
            dkv = D[k][v]
            if dkv == NEG:
                continue
            val = (dnv - dkv) / (n - k)
            if val < worst:
                worst = val
        if worst > best:
            best = worst
    return best

eigenvalue = max_cycle_mean

def _neg_for_min(M):
    n = len(M)
    R = [[NEG] * n for i in range(n)]
    for i in range(n):
        for j in range(n):
            x = M[i][j]
            R[i][j] = NEG if x == INF else -x
    return R

def min_cycle_mean(A):
    m = max_cycle_mean(_neg_for_min(A))
    if m == NEG:
        return INF
    return -m

def _normalize(A, lam):
    n = len(A)
    return [[(A[i][j] - lam) if A[i][j] != NEG else NEG for j in range(n)]
            for i in range(n)]

def _spectral(A):
    lam = max_cycle_mean(A)
    if lam == NEG:
        return (NEG, None, [])
    Al = _normalize(A, lam)
    star = closure(Al, MAXPLUS)          # A_lambda *
    plus = mmul(Al, star, MAXPLUS)       # A_lambda +  (max cycle weight through i)
    n = len(A)
    crit = []
    for i in range(n):
        if abs(plus[i][i]) <= EPS:       # on a critical (0-mean) cycle
            crit.append(i)
    return (lam, star, crit)

def critical_nodes(A):
    lam, star, crit = _spectral(A)
    return crit

def eigenvector(A, node=None):
    # Returns (lam, v) with A (X) v = lam (X) v. Column of A_lambda* over a
    # critical node. node selects which critical class; default = first.
    lam, star, crit = _spectral(A)
    if lam == NEG or not crit:
        return None
    if node is None:
        node = crit[0]
    elif node not in crit:
        return None
    n = len(A)
    v = [star[j][node] for j in range(n)]
    return (lam, v)

def check_eig(A, tol=None):
    # Numerically verify A (X) v == lam (X) v. Returns (ok, lam, v, residual).
    if tol is None:
        tol = EPS
    res = eigenvector(A)
    if res is None:
        return (False, None, None, None)
    lam, v = res
    Av = mvmul(A, v, MAXPLUS)
    worst = 0.0
    for i in range(len(v)):
        lhs = Av[i]
        rhs = (lam + v[i]) if v[i] != NEG else NEG
        if lhs == NEG and rhs == NEG:
            continue
        d = abs(lhs - rhs)
        if d > worst:
            worst = d
    return (worst <= tol, lam, v, worst)

# ---------- structure: reachability / irreducibility ----------

def reach(A, S):
    # Boolean reachability (paths of length >= 1). R[i][j] True iff j reachable
    # from i. Underlies the Perron-Frobenius (irreducibility) hypothesis.
    n = len(A)
    R = [[(A[i][j] != S.zero) for j in range(n)] for i in range(n)]
    for k in range(n):
        for i in range(n):
            if R[i][k]:
                Rk = R[k]; Ri = R[i]
                for j in range(n):
                    if Rk[j]:
                        Ri[j] = True
    return R

def is_irreducible(A, S=MAXPLUS):
    # True iff the graph is strongly connected (every event eventually depends on
    # every other). This is the hinge for max-plus Perron-Frobenius.
    n = len(A)
    if n <= 1:
        return True
    R = reach(A, S)
    for i in range(n):
        for j in range(n):
            if i != j and not R[i][j]:
                return False
    return True

# ---------- residuation: greatest subsolution of A (X) x <= b (max-plus) ----------

def ldiv(A, b):
    # "A \ b": greatest x with A (X) x <= b in max-plus.
    #   x_j = min over i (with A[i][j] finite) of (b_i - A[i][j]); +inf if none.
    n = len(A); m = len(A[0])
    x = [INF] * m
    for j in range(m):
        best = INF
        for i in range(n):
            a = A[i][j]
            if a == NEG:
                continue
            bi = b[i]
            if bi == NEG:          # constraint (-inf) forces x_j -> -inf
                best = NEG
                break
            val = bi - a
            if val < best:
                best = val
        x[j] = best
    return x

# ---------- dynamics ----------

def orbit(A, x0, steps, S=MAXPLUS):
    # Iterate x(k) = A (X) x(k-1). Returns [x0, x1, ..., x_steps].
    traj = [x0[:]]
    x = x0[:]
    for _ in range(steps):
        x = mvmul(A, x, S)
        traj.append(x)
    return traj

def transient(A, maxiter=100):
    # Cyclicity: find (lam, c, T) with A^(k+c) == lam^c (X) A^k for all k >= T.
    lam = max_cycle_mean(A)
    if lam == NEG:
        return None
    seq = []                 # seq[t] = A^(t+1)
    P = copy(A)
    for k in range(1, maxiter + 1):
        for t in range(len(seq)):
            ok, const = _periodic_match(P, seq[t])
            if ok:
                T = t + 1
                c = k - T
                return (lam, c, T)
        seq.append(copy(P))
        P = mmul(A, P, MAXPLUS)
    return None

def _periodic_match(P, Q):
    n = len(P)
    const = None
    for i in range(n):
        Pi = P[i]; Qi = Q[i]
        for j in range(n):
            a = Pi[j]; b = Qi[j]
            ainf = (a == INF or a == NEG)
            binf = (b == INF or b == NEG)
            if ainf != binf:
                return (False, None)
            if ainf:
                if a != b:
                    return (False, None)
                continue
            d = a - b
            if const is None:
                const = d
            elif abs(d - const) > EPS:
                return (False, None)
    return (True, const if const is not None else 0.0)

# =====================================================================
# Kinesin coordination model (GRC Cytoskeletal Motors 2026)
#   tf = swing (detachment) time  tau_f
#   tg = stance (bound) time      tau_g
#   tD = must-wait gating overlap tau_Delta
# Composites:
#   tgam = tf + tg          intrinsic head cycle
#   tdel = tf + tD          gating composite; coordination cycle = 2*tdel
# Recurrent 2x2 block (eliminating the swing phase), poster appendix form:
#   B = [[tgam,        tdel                 ],
#        [tgam + tdel, max(tgam, 2*tdel)     ]]
# Its max-plus eigenvalue (plain max cycle mean) is lambda = max(tgam, 2*tdel).
# =====================================================================

def kinesin_B(tf, tg, tD):
    tgam = tf + tg
    tdel = tf + tD
    return [[tgam, tdel], [tgam + tdel, max(tgam, 2.0 * tdel)]]

def kinesin_companion(tf, tg, tD):
    # Per-step companion form of t[n] = max(t[n-1]+tdel, t[n-2]+tgam):
    #   [t_n; t_{n-1}] = C (X) [t_{n-1}; t_{n-2}],  C = [[tdel, tgam],[0, -inf]].
    # Per-step eigenvalue is max(tdel, tgam/2); the full-cycle period is 2x that,
    # i.e. max(2*tdel, tgam) = lambda(B). Provided so both derivations can be
    # cross-checked on the calculator.
    tgam = tf + tg
    tdel = tf + tD
    return [[tdel, tgam], [0.0, NEG]]

def kinesin_gait(tf, tg, tD, steps=12):
    # Simulate the alternating step sequence and return the settled cadence.
    # Returns (times, step_interval, cycle_period).
    tgam = tf + tg
    tdel = tf + tD
    t = [0.0, tdel]                      # first handoff at tdel
    for n in range(2, steps):
        t.append(_max(t[n - 1] + tdel, t[n - 2] + tgam))
    step_interval = t[-1] - t[-2]        # settles to max(tdel, tgam/2)
    cycle_period = t[-1] - t[-3]         # two steps: -> max(2*tdel, tgam) = lambda
    return (t, step_interval, cycle_period)

def crossover_tg(tf, tD):
    # Value of tg where the regime flips (the knee): tgam == 2*tdel.
    #   tf + tg = 2(tf + tD)  =>  tg* = tf + 2*tD.
    return tf + 2.0 * tD

def kinesin_report(tf, tg, tD):
    tgam = tf + tg
    tdel = tf + tD
    B = kinesin_B(tf, tg, tD)
    lam = max_cycle_mean(B)
    gap = abs(tgam - 2.0 * tdel)
    if 2.0 * tdel > tgam + EPS:
        regime = "gating-limited (coordination cycle critical)"
    elif tgam > 2.0 * tdel + EPS:
        regime = "intrinsic-limited (head cycle critical)"
    else:
        regime = "knee (tgam == 2*tdel)"
    ev = eigenvector(B)
    offset = None
    if ev is not None:
        _, v = ev
        offset = max(v) - min(v)
    print("tf=%s tg=%s tD=%s" % (_f(tf), _f(tg), _f(tD)))
    print("tgam=tf+tg=%s   tdel=tf+tD=%s   2*tdel=%s" % (_f(tgam), _f(tdel), _f(2.0 * tdel)))
    show(B, "B =")
    print("lambda (period) = %s  = max(tgam,2*tdel) = %s" % (_f(lam), _f(max(tgam, 2.0 * tdel))))
    print("regime: " + regime)
    print("Gamma (robustness gap) = %s" % _f(gap))
    if offset is not None:
        print("head offset (from eigenvector) = %s   [= tdel in gating-limited regime]" % _f(offset))
    return (lam, regime, gap, offset)

# ---------- self-test: verifies each poster claim ----------

def poster_selftest():
    # Type  import tropmat ; tropmat.poster_selftest()  on the calculator.
    tally = [0, 0]           # [passed, failed]; list avoids nonlocal
    def check(name, cond):
        if cond:
            tally[0] += 1; print("PASS " + name)
        else:
            tally[1] += 1; print("FAIL " + name)

    # 1) lambda = max(tgam, 2*tdel) across regimes (D4).
    cases = [(1.0, 5.0, 0.5), (3.0, 1.0, 2.0), (2.0, 2.0, 1.0)]
    ok = True
    for (tf, tg, tD) in cases:
        B = kinesin_B(tf, tg, tD)
        lam = max_cycle_mean(B)
        want = max(tf + tg, 2.0 * (tf + tD))
        if abs(lam - want) > EPS:
            ok = False
    check("D4  lambda = max(tgam, 2*tdel)", ok)

    # 2) gating-limited head offset == tdel (D3).
    tf, tg, tD = 3.0, 1.0, 2.0
    B = kinesin_B(tf, tg, tD)
    _, v = eigenvector(B)
    off = max(v) - min(v)
    check("D3  gating-limited offset = tdel", abs(off - (tf + tD)) <= EPS)

    # 3) event graph strongly connected -> irreducible (D2).
    C = kinesin_companion(2.0, 2.0, 1.0)
    check("D2  companion irreducible", is_irreducible(C))

    # 4) eigen-relation holds numerically for B.
    okk, lam, v, resid = check_eig(B)
    check("     B (X) v = lambda (X) v", okk)

    # 5) simulated cadence matches lambda (D2/D4: start-independent period).
    _, _, per = kinesin_gait(3.0, 1.0, 2.0, 16)
    check("     simulated period -> lambda", abs(per - 10.0) <= EPS)

    # 6) crossover tg* = tf + 2*tD flips the regime (D6 knee).
    tf, tD = 2.0, 1.0
    tgs = crossover_tg(tf, tD)
    Bk = kinesin_B(tf, tgs, tD)
    check("D6  knee at tg = tf + 2*tD", abs(max_cycle_mean(Bk) - (tf + tgs)) <= EPS
          and abs((tf + tgs) - 2.0 * (tf + tD)) <= EPS)

    # 7) general tropical sanity: 2-cycle mean 5/2.
    A = [[NEG, 2.0], [3.0, NEG]]
    check("     tropical cycle mean 5/2", abs(max_cycle_mean(A) - 2.5) <= EPS)

    # 8) regime structure -- guards the intrinsic corner where the reduction
    #    matters (a wrong B can match lambda everywhere yet mispredict here):
    #      gating-limited    -> ONE critical circuit  -> unique two-beat (D2)
    #      intrinsic-limited -> TWO critical circuits  -> gait not pinned (D7)
    Bg = kinesin_B(3.0, 1.0, 2.0)                 # gating: 2*tdel=10 > tgam=4
    check("D2  gating: single critical circuit (unique gait)",
          len(critical_nodes(Bg)) == 1)
    Bi = kinesin_B(1.0, 6.0, 1.0)                 # intrinsic: tgam=7 > 2*tdel=4
    li = max_cycle_mean(Bi)
    check("D7  intrinsic: two critical circuits (gait degenerate)",
          len(critical_nodes(Bi)) == 2)
    check("     intrinsic: both self-loops carry lambda",
          abs(Bi[0][0] - li) <= EPS and abs(Bi[1][1] - li) <= EPS)

    print("%d passed, %d failed" % (tally[0], tally[1]))
    return tally[1] == 0