Jump to content

The Principia Machina/TPC-PY/tpc eom.py: Difference between revisions

From The Invariance Itself
TPC-PY source tpc_eom.py (0.1.0), SHA-256 215e78b3ddcdb0c4; conformed TPS 0.3.0; PROPOSED, not declared (CO-8)
 
m 1 revision imported
(No difference)

Revision as of 04:37, 28 September 2026

# -*- coding: utf-8 -*-
"""
TPC-PY --- the structure check of a kernel's equations of motion
(tpc_eom 0.1.0, PROPOSED; D176.1 and D177.1, taken 2026-09-26; standing
rulings S22-S31 of tpm_eom_round_selections 0.2.0 and 0.3.0)
=====================================================================
Standard library only; Python 3.8+; Windows/Linux. Spec and report:
tpm_eom_check 0.1.0 (agent R2, EOM round pass 4), section 3. Vectors:
vectors/eom_corpus.json (gen_eom_corpus.py). TPC-PHP mirrors it operation
for operation (tpc_php/tpc_eom.php).

THE INPUT is a TFK in the schema 0.5.0 shape (BRIEF_COMMON_P4 section 2),
read as plain objects; tpc_zone does not validate it here. READ, in this
order, the first problem refusing (["REFUSED", reason], never raised):
  the TFK; the TPK; DOF (each NAME, an identifier; ROLE absent or "INPUT");
  CONSTANT (each NAME; VALUE a string when present); LAW (each NAME); the
  NAMEs unique across DOF, CONSTANT and LAW (the least duplicate named);
  CONSTRAINT (each TARGET and ACTION, strings; ACTION FIX or REMOVE; TARGET
  a DOF; a FIX's VALUE present and a decimal); LAW in list order (ROLE
  required: EOM(<NAME>), CONSTRAINT, CONSTITUTIVE, INVARIANT or DEFINITION,
  spelled exactly so, EOM's NAME a DOF; EQUATION, when present, a string
  that parses under tpc_form and names only DOFs and CONSTANTs, the least
  other name refused; without EQUATION, REFS a list of strings). A kind
  present is a list whose entries are all objects, checked before any
  entry's fields; an absent kind holds no entry. Not read, not checked: ASSUMPTION, BOUNDARY, the TPK
  PARAMETERS (D_VARIABLES), DIMENSION, UNIT, UNCERTAINTY, STATEMENT, a
  formed law's REFS, a CONSTRAINT's NAME and a REMOVE's VALUE, the TrM;
  the TMS only in eom_class.

KNOWNS. Every CONSTANT. A DOF is FIXED (a FIX targets it), else REMOVED (a
REMOVE targets it), else INPUT (ROLE INPUT), else INDEPENDENT (a D
variable), in that precedence. Every other DOF is an UNKNOWN. The D
VARIABLES are the second arguments of every D in every EQUATION read
(every ROLE; in eom_register_test also the entering records' forms).

ZERO PROPAGATION, exact, on each law's form before matching: D(x, s, n) of
a FIXED x is 0; a FIXED DOF every FIX of which has the VALUE 0 is 0; a
literal 0 is 0; then a*0, 0*a, 0/a (a not 0), -0, (0), 0+0, 0-0, 0^(p/q)
with p/q > 0, abs(0), sqrt(0), sin(0), tan(0) are 0, and a+0, 0+a, a-0,
0-a keep a's names; exp(0), cos(0), ln(0), a/0 and 0^(p/q <= 0) are not 0
and name nothing of the zero. A zero subtree names nothing. A law's names
after propagation are the identifiers (D's operand and variable included)
outside zero subtrees; an unknown's ORDER in a law is the highest n of a
D(unknown, s, n) there, 0 when it appears underived.

THE ANALYSIS. ANALYSIS is "NO FORM" when no law read carries EQUATION,
"PDE" when the forms use two or more distinct D variables (not analyzed),
else "RUN". Under RUN, the ROWS are the laws with EQUATION whose ROLE
enters the matching (EOM, CONSTRAINT, CONSTITUTIVE, DEFINITION) and that
name an unknown after propagation; one that names none is CLOSED: ZERO
(both sides 0), BOUND (every name a CONSTANT with VALUE or a FIXED DOF:
T3's ground fragment) or UNBOUND. The COLUMNS are every unknown. A maximum
matching (augmenting paths) gives the Dulmage-Mendelsohn coarse parts:
UNDER (the unknowns alternating paths reach from an unmatched unknown, and
the rows beside them), OVER (likewise from an unmatched row), SQUARE (the
rest). In each part a matching saturating its smaller side (both in
SQUARE) is chosen, lexicographically: the greatest sum of ORDER, then the
most ROLE agreements (a law EOM(q) matched to q), then the least sequence
of partner NAMEs reading the saturated side (the laws in UNDER and SQUARE,
the unknowns in OVER) in NAME order; integer weights ORDER x K + agreement,
K = rows + columns of the part + 1, a Hungarian assignment (O(n^3)) on the
padded square cost matrix, then the least sequence among the assignments
tight under its potentials. On SQUARE, Pryce's signature method: that
assignment is a highest-value transversal; the offsets c (laws), d
(unknowns) are the least with d_j - c_i >= ORDER_ij and equality on it, by
the fixed point from c = 0 (at most rows + 2 passes); INDEX = max c + (1 if
some d_j = 0), INITIAL_VALUES = sum d - sum c (information; 0 and 0 when
SQUARE is empty). Every tie by NAME in code point order.

THE FINDINGS, in this order (subject, then kind): UNDER (per unknown the
chosen UNDER matching leaves unmatched, with its connected block), OVER-EOM
(per q: the laws with EQUATION and ROLE EOM(q), when two or more, and the
laws with a side of the equation-of-motion shape for an unknown q and one
order, when two or more), OVER (per connected block of the over part that
remains when every OVER-EOM set but its first law is set aside), CLOSED
(UNBOUND laws, with their unbound names), ROLE MISMATCH (an EOM(q) law
whose form has no D(q, s, n): NO DERIVATIVE; an EOM(q) law, prose or not,
whose q is known: KNOWN), PROSE (a law without EQUATION whose ROLE enters
the matching and whose REFS name an unknown), NO FORM, PDE. Under PDE and
NO FORM only ROLE MISMATCH and PROSE run beside the line. INVARIANT laws
leave the matching and are listed as information. The EQUATION-OF-MOTION
SHAPE of a side: after parentheses and unary minus, a product and quotient
whose numerator factors hold exactly one D(q, s, n), q an unknown, every
other factor (the denominators included) naming only CONSTANTs.

eom_class reads STRUCTURE_CHECK (TMS PARAMETERS; absent = ON): OFF gives
"Structure not Validated at Inception" with CHECK OFF alone and reads no
TPK; ON gives "structurally regular" when there is no finding at all, else
the label with its reasons FINDINGS (any of UNDER, OVER-EOM, OVER, CLOSED,
ROLE MISMATCH), PROSE, PDE, NO FORM, in that order (S30). A SQUARE part is
structurally regular by construction whenever the analysis runs.

eom_register_test adds the standing records (STATUS a status word other
than STRUCK, HELD, REFUTED, VOID) that declare a ROLE entering the matching
and an EQUATION, as laws named by their IDs ("<prefix>-<n>"), and returns
the map and findings of the whole, the IDs used (by prefix, then n
numerically), SETTLED (the TFK alone's UNDER unknowns now outside the
combined UNDER part; none when the whole is not RUN) and OPEN (the rest).
Every record is read whole, standing or not.

PARSER: this module reads a form through tpc_form's public form_tree (S32,
tpc_form 0.2.0; C-R2-1 closed), never its private parser; a refusal is
form_tree's reason, which is parse_form's. It relies on the tree shapes
form_tree documents: ["eq", l, r], ["num", s], ["id", n], ["neg", a],
["paren", a], ["bin", op, a, b], ["pow", a, [p, q]], ["fn", f, a],
["D", x, s, n]. Every walk runs on an explicit stack, so no form depth
reaches Python's recursion limit.
"""

import re
from collections import deque

import tpc_zone as Z
import tpc_form as F

TPC_EOM_VERSION = "0.1.0"      # proposed; unreleased working label (D176.1, D177.1)

ENTER = ("EOM", "CONSTRAINT", "CONSTITUTIVE", "DEFINITION")
NONSTANDING = ("STRUCK", "HELD", "REFUTED", "VOID")
REGULAR = "structurally regular"
NOT_VALIDATED = "Structure not Validated at Inception"
CLASS_ORDER = ["UNDER", "OVER-EOM", "OVER", "CLOSED", "ROLE MISMATCH", "PROSE", "NO FORM", "PDE"]
FINDING_CLASSES = ["UNDER", "OVER-EOM", "OVER", "CLOSED", "ROLE MISMATCH"]
MISMATCH_KINDS = ["NO DERIVATIVE", "KNOWN"]
_ROLE = re.compile(r"^(?:EOM\(([A-Za-z_][A-Za-z0-9_]*)\)|CONSTRAINT|CONSTITUTIVE|INVARIANT|DEFINITION)\Z")
_RID = re.compile(r"^[A-Za-z_][A-Za-z0-9_]*-[1-9][0-9]*\Z")
_STATUS = re.compile(r"^[A-Z][A-Z_]*\Z")


class _Refused(Exception):
    pass


def _refuse(msg):
    return ["REFUSED", Z._show(msg)]


# ---- reading ----------------------------------------------------------------------------
def _kind(tpk, kind):
    if kind not in tpk:
        return []
    v = tpk[kind]
    if not isinstance(v, list):
        raise _Refused("TPK kind " + kind + " is not a list")
    for e in v:
        if not isinstance(e, dict):
            raise _Refused("TPK kind " + kind + ": entry is not an object")
    return v


def _field_str(tag, e, f):
    if f not in e:
        raise _Refused(tag + ": missing field " + f)
    if not isinstance(e[f], str):
        raise _Refused(tag + ": field " + f + " has the wrong type")
    return e[f]


def _name_of(tag, e):
    n = _field_str(tag, e, "NAME")
    if not Z._is_ident(n):
        raise _Refused(tag + ": NAME " + n + " is not an identifier")
    return n


def _role(tag, r, dofs):
    """(kind, q) of a ROLE string; refused when it is not a role."""
    m = _ROLE.match(r)
    if m is None:
        raise _Refused(tag + ": ROLE " + r + " is not a role")
    q = m.group(1)
    if q is None:
        return r, None
    if q not in dofs:
        raise _Refused(tag + ": ROLE " + r + ": " + q + " is not a DOF")
    return "EOM", q


def _kids(node):
    k = node[0]
    if k == "bin":
        return [node[2], node[3]]
    if k in ("neg", "paren", "pow"):
        return [node[1]]
    if k == "fn":
        return [node[2]]
    if k == "eq":
        return [node[1], node[2]]
    return []


def _walk(tree):
    """Every node, parents first; explicit stack."""
    out, stack = [], [tree]
    while stack:
        n = stack.pop()
        out.append(n)
        ks = _kids(n)
        for i in range(len(ks) - 1, -1, -1):
            stack.append(ks[i])
    return out


def _idents(tree):
    acc = set()
    for n in _walk(tree):
        if n[0] == "id":
            acc.add(n[1])
        elif n[0] == "D":
            acc.add(n[1])
            acc.add(n[2])
    return acc


def _form(tag, text, declared):
    if not isinstance(text, str):
        raise _Refused(tag + ": field EQUATION has the wrong type")
    r = F.form_tree(text, "equation")
    if r[0] != "OK":
        raise _Refused(tag + ": EQUATION does not parse: " + r[1])
    tree = r[1]
    for n in sorted(_idents(tree)):
        if n not in declared:
            raise _Refused(tag + ": EQUATION names " + n + ", which is not a DOF or CONSTANT")
    return tree


def _read(tfk):
    if not isinstance(tfk, dict):
        raise _Refused("TFK is not an object")
    tpk = tfk.get("TPK")
    if not isinstance(tpk, dict):
        raise _Refused("TPK: zone is not an object")
    dofs, inputs = [], set()
    for e in _kind(tpk, "DOF"):
        tag = "DOF:" + Z._name(e)
        n = _name_of(tag, e)
        if "ROLE" in e:
            r = e["ROLE"]
            if not isinstance(r, str):
                raise _Refused(tag + ": field ROLE has the wrong type")
            if r != "INPUT":
                raise _Refused(tag + ": ROLE " + r + " is not INPUT")
            inputs.add(n)
        dofs.append(n)
    cnames, valued = [], set()
    for e in _kind(tpk, "CONSTANT"):
        tag = "CONSTANT:" + Z._name(e)
        n = _name_of(tag, e)
        if "VALUE" in e:
            if not isinstance(e["VALUE"], str):
                raise _Refused(tag + ": field VALUE has the wrong type")
            valued.add(n)
        cnames.append(n)
    raw = []
    for e in _kind(tpk, "LAW"):
        tag = "LAW:" + Z._name(e)
        raw.append((_name_of(tag, e), tag, e))
    seen, dup = set(), set()
    for n in dofs + cnames + [r[0] for r in raw]:
        if n in seen:
            dup.add(n)
        seen.add(n)
    if dup:
        raise _Refused("TPK: duplicate name " + sorted(dup)[0])
    dofset = set(dofs)
    fixes, removed = {}, set()
    for e in _kind(tpk, "CONSTRAINT"):
        tag = "CONSTRAINT:" + Z._name(e)
        t = _field_str(tag, e, "TARGET")
        a = _field_str(tag, e, "ACTION")
        if a not in ("FIX", "REMOVE"):
            raise _Refused(tag + ": ACTION " + a + " is not REMOVE or FIX")
        if t not in dofset:
            raise _Refused(tag + ": TARGET " + t + " is not a DOF")
        if a == "FIX":
            if "VALUE" not in e:
                raise _Refused(tag + ": FIX requires VALUE")
            v = e["VALUE"]
            if not isinstance(v, str):
                raise _Refused(tag + ": field VALUE has the wrong type")
            cv = Z.canon_decimal(v)
            if cv is None:
                raise _Refused(tag + ": VALUE " + v + " is not a decimal")
            fixes.setdefault(t, []).append(cv == "0")
        else:
            removed.add(t)
    declared = dofset | set(cnames)
    laws = []
    for n, tag, e in raw:
        r = _field_str(tag, e, "ROLE")
        kind, q = _role(tag, r, dofset)
        tree, refs = None, None
        if "EQUATION" in e:
            tree = _form(tag, e["EQUATION"], declared)
        else:
            if "REFS" not in e:
                raise _Refused(tag + ": missing field REFS")
            refs = e["REFS"]
            if not isinstance(refs, list) or not all(isinstance(x, str) for x in refs):
                raise _Refused(tag + ": field REFS has the wrong type")
        laws.append({"name": n, "kind": kind, "q": q, "tree": tree, "refs": refs})
    return {"dofs": dofs, "inputs": inputs, "consts": set(cnames), "valued": valued,
            "fixed": set(fixes), "zerofix": set(t for t, zs in fixes.items() if all(zs)),
            "removed": removed, "declared": declared, "dofset": dofset, "laws": laws}


def _read_records(records, K):
    if not isinstance(records, list):
        raise _Refused("records is not a list")
    seen, entering = set(), []
    for r in records:
        if not isinstance(r, dict):
            raise _Refused("records: entry is not an object")
        tag = "RECORD:" + (r["ID"] if isinstance(r.get("ID"), str) else "?")
        rid = _field_str(tag, r, "ID")
        if _RID.match(rid) is None:
            raise _Refused(tag + ": ID " + rid + " is not a record identifier")
        if rid in seen:
            raise _Refused("records: duplicate ID " + rid)
        seen.add(rid)
        st = _field_str(tag, r, "STATUS")
        if _STATUS.match(st) is None:
            raise _Refused(tag + ": STATUS " + st + " is not a status word")
        kind, q, tree = None, None, None
        if "ROLE" in r:
            if not isinstance(r["ROLE"], str):
                raise _Refused(tag + ": field ROLE has the wrong type")
            kind, q = _role(tag, r["ROLE"], K["dofset"])
        if "EQUATION" in r:
            tree = _form(tag, r["EQUATION"], K["declared"])
        if st not in NONSTANDING and kind in ENTER and tree is not None:
            entering.append({"name": rid, "kind": kind, "q": q, "tree": tree, "refs": None})
    return entering


# ---- the per-law readings -------------------------------------------------------------
def _dvars(tree, acc):
    for n in _walk(tree):
        if n[0] == "D":
            acc.add(n[2])


def _has_d_of(tree, q):
    for n in _walk(tree):
        if n[0] == "D" and n[1] == q:
            return True
    return False


def _merge(a, b):
    """Union of two {name: order} maps, the higher order kept; the larger map
    is updated in place (the children's maps are used once)."""
    if len(a) < len(b):
        a, b = b, a
    for k, v in b.items():
        if k not in a or a[k] < v:
            a[k] = v
    return a


def _prop(tree, K):
    """(zero, {name: order}) of a subtree after zero propagation; post-order on
    an explicit stack."""
    stack, vals = [(tree, False)], []
    while stack:
        node, done = stack.pop()
        k = node[0]
        if not done:
            ks = _kids(node)
            if ks:
                stack.append((node, True))
                for i in range(len(ks) - 1, -1, -1):
                    stack.append((ks[i], False))
                continue
            if k == "num":
                vals.append((node[1] == "0", {}))
            elif k == "id":
                if node[1] in K["zerofix"]:
                    vals.append((True, {}))
                else:
                    vals.append((False, {node[1]: 0}))
            else:                                          # D
                x, s, n = node[1], node[2], node[3]
                if x in K["fixed"]:
                    vals.append((True, {}))
                else:
                    occ = {s: 0}
                    occ[x] = n if x not in occ or occ[x] < n else occ[x]
                    vals.append((False, occ))
            continue
        if k in ("neg", "paren"):
            continue                                       # its one value stands
        if k == "bin":
            zb, ob = vals.pop()
            za, oa = vals.pop()
            op = node[1]
            if op in ("+", "-"):
                z = za and zb
            elif op == "*":
                z = za or zb
            else:
                z = za and not zb
            vals.append((True, {}) if z else (False, _merge(oa, ob)))
        elif k == "pow":
            za, oa = vals.pop()
            if za:
                vals.append((node[2][0] > 0, {}))
            else:
                vals.append((False, oa))
        elif k == "fn":
            za, oa = vals.pop()
            if za:
                vals.append((node[1] in ("abs", "sqrt", "sin", "tan"), {}))
            else:
                vals.append((False, oa))
        else:                                              # eq: kept as the two sides
            pass
    return vals


def _factors(side):
    num, den, stack = [], [], [side]
    while stack:
        n = stack.pop()
        k = n[0]
        if k in ("paren", "neg"):
            stack.append(n[1])
        elif k == "bin" and n[1] == "*":
            stack.append(n[3])
            stack.append(n[2])
        elif k == "bin" and n[1] == "/":
            den.append(n[3])
            stack.append(n[2])
        else:
            num.append(n)
    return num, den


def _shapes(tree, K, unknowns):
    out = set()
    for side in (tree[1], tree[2]):
        num, den = _factors(side)
        ds = [f for f in num if f[0] == "D"]
        if len(ds) != 1 or ds[0][1] not in unknowns:
            continue
        rest = [f for f in num if f[0] != "D"] + den
        if all(all(n in K["consts"] for n in _idents(f)) for f in rest):
            out.add((ds[0][1], ds[0][3]))
    return out


# ---- matching ---------------------------------------------------------------------------
def _max_matching(rows, adj):
    """A greedy pass (each row, in order, takes its first free column), then
    augmenting paths (Kuhn) from each row still unmatched, in order, each
    row's columns in order; an explicit stack. Any maximum matching gives the
    same Dulmage-Mendelsohn parts. -> (row -> col, col -> row)."""
    mr, mc = {}, {}
    for r in rows:
        for cc in adj[r]:
            if cc not in mc:
                mr[r] = cc
                mc[cc] = r
                break
    for r in rows:
        if r in mr:
            continue
        visited, parent = set(), {}
        stack, found = [[r, 0]], None
        while stack:
            top = stack[-1]
            x, lst = top[0], adj[top[0]]
            advanced = False
            while top[1] < len(lst):
                cc = lst[top[1]]
                top[1] += 1
                if cc in visited:
                    continue
                visited.add(cc)
                parent[cc] = x
                if cc not in mc:
                    found = cc
                else:
                    stack.append([mc[cc], 0])
                    advanced = True
                break
            if found is not None:
                break
            if not advanced:
                stack.pop()
        if found is not None:
            cc = found
            while True:
                x = parent[cc]
                prev = mr.get(x)
                mr[x] = cc
                mc[cc] = x
                if x == r:
                    break
                cc = prev
    return mr, mc


def _dm(rows, cols, adj, cadj):
    mr, mc = _max_matching(rows, adj)
    cu = set(c for c in cols if c not in mc)
    ru, queue = set(), deque(sorted(cu))
    while queue:
        cc = queue.popleft()
        for r in cadj[cc]:
            if r not in ru:
                ru.add(r)
                m = mr.get(r)
                if m is not None and m not in cu:
                    cu.add(m)
                    queue.append(m)
    ro = set(r for r in rows if r not in mr)
    co, queue = set(), deque(sorted(ro))
    while queue:
        r = queue.popleft()
        for cc in adj[r]:
            if cc not in co:
                co.add(cc)
                m = mc.get(cc)
                if m is not None and m not in ro:
                    ro.add(m)
                    queue.append(m)
    return cu, ru, ro, co


def _hungarian(n, cost):
    """Min-cost perfect assignment of an n x n matrix whose entries are
    integers or None (no edge), a perfect assignment over the edges existing.
    -> (row -> col, u, v) with u[i] + v[j] <= cost[i][j] on every edge and
    equality on the assignment (1-based potentials, e-maxx form). Of the
    columns of least reduced cost a free one is taken first, then the least
    index; the answer's use (the tight graph of an optimal dual) does not
    depend on which optimum is found."""
    u = [0] * (n + 1)
    v = [0] * (n + 1)
    p = [0] * (n + 1)
    way = [0] * (n + 1)
    for i in range(1, n + 1):
        p[0] = i
        j0 = 0
        minv = [None] * (n + 1)
        used = [False] * (n + 1)
        while True:
            used[j0] = True
            i0 = p[j0]
            delta, j1 = None, -1
            row = cost[i0 - 1]
            for j in range(1, n + 1):
                if used[j]:
                    continue
                cij = row[j - 1]
                if cij is not None:
                    cur = cij - u[i0] - v[j]
                    if minv[j] is None or cur < minv[j]:
                        minv[j] = cur
                        way[j] = j0
                if minv[j] is not None and (delta is None or minv[j] < delta
                                            or (minv[j] == delta and p[j] == 0 and p[j1] != 0)):
                    delta = minv[j]
                    j1 = j
            if j1 == -1:
                raise _Refused("internal: no assignment in a part that has one")
            for j in range(n + 1):
                if used[j]:
                    u[p[j]] += delta
                    v[j] -= delta
                elif minv[j] is not None:
                    minv[j] -= delta
            j0 = j1
            if p[j0] == 0:
                break
        while True:
            j1 = way[j0]
            p[j0] = p[j1]
            j0 = j1
            if j0 == 0:
                break
    assign = [0] * n
    for j in range(1, n + 1):
        assign[p[j] - 1] = j - 1
    return assign, u, v


def _best(A, B, w):
    """A (sorted) saturated into B (sorted), |A| <= |B|; w[a][b] the integer
    weight of each edge. The greatest total weight, then the least sequence
    of partner NAMEs reading A in order. -> {a: b}."""
    na, n = len(A), len(B)
    if na == 0:
        return {}
    bidx = dict((b, j) for j, b in enumerate(B))
    cost = []
    for a in A:
        row = [None] * n
        for b, x in w[a].items():
            row[bidx[b]] = -x
        cost.append(row)
    for _ in range(n - na):
        cost.append([0] * n)
    assign, u, v = _hungarian(n, cost)
    tight = []
    for i in range(n):
        tight.append([j for j in range(n) if cost[i][j] is not None and u[i + 1] + v[j + 1] == cost[i][j]])
    owner = [0] * n
    for i in range(n):
        owner[assign[i]] = i
    fixed = [False] * n
    for i in range(na):
        c0 = assign[i]
        for j in tight[i]:
            if j == c0:
                break
            r = owner[j]
            if fixed[r]:
                continue
            via, queue, found = {r: (i, j)}, deque([r]), None
            while queue and found is None:
                x = queue.popleft()
                for y in tight[x]:
                    if y == assign[x]:
                        continue
                    if y == c0:
                        found = x
                        break
                    o = owner[y]
                    if fixed[o] or o == i or o in via:
                        continue
                    via[o] = (x, y)
                    queue.append(o)
            if found is None:
                continue
            cur, take = found, c0
            while True:
                assign[cur] = take
                owner[take] = cur
                pr, col = via[cur]
                if pr == i:
                    assign[i] = col
                    owner[col] = i
                    break
                cur, take = pr, col
            break
        fixed[i] = True
    return dict((A[i], B[assign[i]]) for i in range(na))


def _components(rows, cols, adj):
    """Connected blocks of the subgraph on rows and cols: [(rows, cols)] in
    order of their least member (columns first, then rows)."""
    radj, cadj = {}, {}
    for r in rows:
        for cc in adj[r]:
            if cc in cols:
                radj.setdefault(r, []).append(cc)
                cadj.setdefault(cc, []).append(r)
    seen, out = set(), []
    for start in [("c", x) for x in sorted(cols)] + [("r", x) for x in sorted(rows)]:
        if start in seen:
            continue
        br, bc, queue = set(), set(), deque([start])
        seen.add(start)
        while queue:
            t, x = queue.popleft()
            if t == "r":
                br.add(x)
                nxt = [("c", y) for y in radj.get(x, [])]
            else:
                bc.add(x)
                nxt = [("r", y) for y in cadj.get(x, [])]
            for y in nxt:
                if y not in seen:
                    seen.add(y)
                    queue.append(y)
        out.append((br, bc))
    return out


def _pryce(rows, cols, sig, T):
    """Offsets (c, d) by the fixed point from c = 0; refused if it does not
    settle (it settles for a highest-value transversal)."""
    c = dict((r, 0) for r in rows)
    crow = {}
    for r in rows:
        for j in sig[r]:
            crow.setdefault(j, []).append(r)
    for _ in range(len(rows) + 2):
        d = dict((j, max(sig[r][j] + c[r] for r in crow[j])) for j in cols)
        c2 = dict((r, d[T[r]] - sig[r][T[r]]) for r in rows)
        if c2 == c:
            return c, d
        c = c2
    raise _Refused("internal: the square part's offsets do not settle")


# ---- phrases (the questions for H-DECL) ----------------------------------------------
def _is_rec(n):
    return "-" in n


def _subj(n):
    return ("RECORD " if _is_rec(n) else "LAW ") + n


def _phrase(names):
    laws = sorted(n for n in names if not _is_rec(n))
    recs = sorted(n for n in names if _is_rec(n))
    out = []
    if laws:
        out.append(("LAW " if len(laws) == 1 else "LAWs ") + ", ".join(laws))
    if recs:
        out.append(("RECORD " if len(recs) == 1 else "RECORDs ") + ", ".join(recs))
    return " and ".join(out)


def _unk(names):
    names = sorted(names)
    return ("the unknown " if len(names) == 1 else "the unknowns ") + ", ".join(names)


# ---- the analysis ----------------------------------------------------------------------
def _analyse(K, laws):
    """-> (R, findings, combined UNDER unknowns or None, UNDER list)."""
    dofs = sorted(K["dofs"])
    dv = set()
    for l in laws:
        if l["tree"] is not None:
            _dvars(l["tree"], dv)
    known = {}
    for d in dofs:
        if d in K["fixed"]:
            known[d] = "FIXED"
        elif d in K["removed"]:
            known[d] = "REMOVED"
        elif d in K["inputs"]:
            known[d] = "INPUT"
        elif d in dv:
            known[d] = "INDEPENDENT"
    unknowns = set(d for d in dofs if d not in known)
    formed = [l for l in laws if l["tree"] is not None]
    invariant = sorted(l["name"] for l in laws if l["kind"] == "INVARIANT")
    prose = []
    for l in laws:
        if l["tree"] is None and l["kind"] in ENTER:
            u = sorted(set(l["refs"]) & unknowns)
            if u:
                prose.append([l["name"], u])
    prose.sort()
    if not formed:
        analysis = "NO FORM"
    elif len(dv) >= 2:
        analysis = "PDE"
    else:
        analysis = "RUN"
    fnd = []
    for l in sorted(laws, key=lambda l: l["name"]):
        if l["kind"] != "EOM":
            continue
        q = l["q"]
        if l["tree"] is not None and not _has_d_of(l["tree"], q):
            fnd.append(["ROLE MISMATCH", l["name"], q, "NO DERIVATIVE",
                        "%s declares ROLE EOM(%s), but its EQUATION has no derivative of %s: which ROLE does it "
                        "have?" % (_subj(l["name"]), q, q)])
        if q in known:
            fnd.append(["ROLE MISMATCH", l["name"], q, "KNOWN",
                        "%s declares ROLE EOM(%s), but %s is known (%s): which ROLE does it have?"
                        % (_subj(l["name"]), q, q, known[q])])
    kind_of = dict((l["name"], l) for l in laws)
    for name, u in prose:
        l = kind_of[name]
        rl = "EOM(%s)" % l["q"] if l["kind"] == "EOM" else l["kind"]
        fnd.append(["PROSE", name, u, "%s (ROLE %s) has no EQUATION and names %s: its equation is outside the "
                    "analysis; can it be declared in form?" % (_subj(name), rl, _unk(u))])
    mp, closed_info = [], []
    parts = [["UNDER", [], []], ["SQUARE", [], []], ["OVER", [], []]]
    index, iv, offs = "NONE", "NONE", [[], []]
    under_cols, under_list = None, []
    if analysis == "NO FORM":
        for d in dofs:
            mp.append([d, known.get(d, "NO LAW"), ""])
        u = sorted(unknowns)
        if u:
            q = "no LAW carries an EQUATION: %s %s?" % (
                "which law gives the equation of" if len(u) == 1 else "which laws give the equations of", _unk(u))
        else:
            q = "no LAW carries an EQUATION and the kernel has no unknown: is that intended?"
        fnd.append(["NO FORM", u, q])
    elif analysis == "PDE":
        for d in dofs:
            mp.append([d, known.get(d, "NOT ANALYZED"), ""])
        fnd.append(["PDE", sorted(dv), "PDE: not analyzed"])
    else:
        bound = K["valued"] | K["fixed"]
        rows, sig, role, shape_of = [], {}, {}, {}
        for l in formed:
            if l["kind"] not in ENTER:
                continue
            (zl, ol), (zr, orr) = _prop(l["tree"], K)
            occ = _merge(ol, orr)
            us = dict((n, o) for n, o in occ.items() if n in unknowns)
            shape_of[l["name"]] = _shapes(l["tree"], K, unknowns)
            if not us:
                if zl and zr:
                    closed_info.append([l["name"], "ZERO"])
                elif all(n in bound for n in occ):
                    closed_info.append([l["name"], "BOUND"])
                else:
                    ub = sorted(n for n in occ if n not in bound)
                    closed_info.append([l["name"], "UNBOUND"])
                    fnd.append(["CLOSED", l["name"], ub,
                                "%s names no unknown once the fixed DOFs are applied, and depends on %s, which the "
                                "kernel does not bind: does it hold for every value of %s, or is a force or an "
                                "unknown missing?" % (_subj(l["name"]), ", ".join(ub), ", ".join(ub))])
                continue
            rows.append(l["name"])
            sig[l["name"]] = us
            role[l["name"]] = l["q"] if l["kind"] == "EOM" else None
        closed_info.sort()
        rows.sort()
        cols = sorted(unknowns)
        adj = dict((r, sorted(sig[r])) for r in rows)
        cadj = dict((cc, []) for cc in cols)
        for r in rows:
            for cc in adj[r]:
                cadj[cc].append(r)
        cu, ru, ro, co = _dm(rows, cols, adj, cadj)
        rs = set(rows) - ru - ro
        cs = set(cols) - cu - co
        parts = [["UNDER", sorted(cu), sorted(ru)], ["SQUARE", sorted(cs), sorted(rs)],
                 ["OVER", sorted(co), sorted(ro)]]

        def weights(A, B, by_row):
            k = len(A) + len(B) + 1
            w = {}
            for a in A:
                w[a] = {}
                for b in B:
                    r, cc = (a, b) if by_row else (b, a)
                    if cc in sig[r]:
                        w[a][b] = sig[r][cc] * k + (1 if role[r] == cc else 0)
            return w

        RU, CU, RS, CS, RO, CO = sorted(ru), sorted(cu), sorted(rs), sorted(cs), sorted(ro), sorted(co)
        mU = _best(RU, CU, weights(RU, CU, True))
        mS = _best(RS, CS, weights(RS, CS, True))
        mO = _best(CO, RO, weights(CO, RO, False))
        law_of = {}
        for r, cc in mU.items():
            law_of[cc] = ("UNDER", r)
        for r, cc in mS.items():
            law_of[cc] = ("DETERMINED", r)
        for cc, r in mO.items():
            law_of[cc] = ("OVER", r)
        for d in dofs:
            if d in known:
                mp.append([d, known[d], ""])
            elif d in law_of:
                mp.append([d, law_of[d][0], law_of[d][1]])
            else:
                mp.append([d, "NO LAW", ""])
        if RS:
            sS = dict((r, dict((cc, o) for cc, o in sig[r].items() if cc in cs)) for r in RS)
            cc_, dd_ = _pryce(RS, CS, sS, mS)
            index = max(cc_.values()) + (1 if any(x == 0 for x in dd_.values()) else 0)
            iv = sum(dd_.values()) - sum(cc_.values())
            offs = [[[r, cc_[r]] for r in RS], [[j, dd_[j]] for j in CS]]
        else:
            index, iv = 0, 0
        under_cols = cu
        matched_u = set(mU.values())
        comps = _components(ru, cu, adj)
        for u in CU:
            if u in matched_u:
                continue
            under_list.append(u)
            br, bc = [b for b in comps if u in b[1]][0]
            if br:
                q = ("unknown %s is determined by no law with a declared form (%s %s %s under-determined): which "
                     "equation determines %s?" % (u, _phrase(br), "leaves" if len(br) == 1 else "leave",
                                                  ", ".join(sorted(bc)), u))
            else:
                q = "unknown %s is named by no law with a declared form: which equation determines %s?" % (u, u)
            fnd.append(["UNDER", u, sorted(bc), sorted(br), q])
        sets = {}
        for l in formed:
            if l["kind"] == "EOM":
                sets.setdefault(l["q"], []).append(l["name"])
        oe = {}
        for q, g in sets.items():
            if len(g) >= 2:
                oe.setdefault(q, set()).update(g)
        byshape = {}
        for name in sorted(shape_of):
            for sh in shape_of[name]:
                byshape.setdefault(sh, []).append(name)
        for sh, g in byshape.items():
            if len(g) >= 2:
                oe.setdefault(sh[0], set()).update(g)
        aside, kept = set(), set()
        for q in sorted(oe):
            ls = sorted(oe[q])
            kept.add(ls[0])
            aside.update(ls[1:])
            fnd.append(["OVER-EOM", q, ls, "%s each state the motion of %s: are they one equation, or is one to be "
                        "removed?" % (_phrase(ls), q)])
        aside -= kept
        rrows = [r for r in rows if r not in aside]
        rcadj = dict((cc, [r for r in cadj[cc] if r not in aside]) for cc in cols)
        _, _, ro2, co2 = _dm(rrows, cols, adj, rcadj)
        for br, bc in _components(ro2, co2, adj):
            ls, us = sorted(br), sorted(bc)
            fnd.append(["OVER", ls, us, "%s give more equations than %s they determine: which law is implied by "
                        "the others, or is one an independent claim?" % (_phrase(ls), _unk(us))])

    def fkey(f):
        k = CLASS_ORDER.index(f[0])
        if f[0] == "OVER":
            return (k, f[1][0], 0)
        if f[0] == "ROLE MISMATCH":
            return (k, f[1], MISMATCH_KINDS.index(f[3]))
        if f[0] in ("NO FORM", "PDE"):
            return (k, "", 0)
        return (k, f[1], 0)
    fnd.sort(key=fkey)
    R = [["ANALYSIS", analysis], ["D_VARIABLES", sorted(dv)], ["MAP", mp],
         ["PARTS", parts[0], parts[1], parts[2]], ["INDEX", index], ["INITIAL_VALUES", iv],
         ["OFFSETS", offs[0], offs[1]], ["CLOSED", closed_info], ["INVARIANT", invariant],
         ["PROSE", prose]]
    return R, fnd, (under_cols if analysis == "RUN" else None), under_list


# ---- the four public functions -------------------------------------------------------
def eom_map(tfk):
    """["OK", [ANALYSIS, D_VARIABLES, MAP, PARTS, INDEX, INITIAL_VALUES,
    OFFSETS, CLOSED, INVARIANT, PROSE]] or ["REFUSED", reason]. MAP: one
    [DOF, label, law] per DOF by NAME; label DETERMINED, UNDER or OVER (with
    the law the part's matching gives), NO LAW, FIXED, REMOVED, INPUT,
    INDEPENDENT or NOT ANALYZED; law "" when none."""
    try:
        K = _read(tfk)
        return ["OK", _analyse(K, K["laws"])[0]]
    except _Refused as e:
        return _refuse(str(e))


def eom_findings(tfk):
    """["OK", [finding, ...]] in canonical order, each ending with its question
    for H-DECL; or ["REFUSED", reason]. Reads no switch (S26)."""
    try:
        K = _read(tfk)
        return ["OK", _analyse(K, K["laws"])[1]]
    except _Refused as e:
        return _refuse(str(e))


def eom_class(tfk):
    """S30. ["STRUCTURE", "structurally regular"] or ["STRUCTURE",
    "Structure not Validated at Inception", [reasons]], or ["REFUSED",
    reason]. Reads the TFK, the TMS, its PARAMETERS and STRUCTURE_CHECK
    (absent = ON) first; under OFF nothing else."""
    try:
        if not isinstance(tfk, dict):
            raise _Refused("TFK is not an object")
        tms = tfk.get("TMS")
        if not isinstance(tms, dict):
            raise _Refused("TMS: zone is not an object")
        params = tms.get("PARAMETERS")
        if not isinstance(params, dict):
            raise _Refused("TMS: part PARAMETERS is not an object")
        sc = params.get("STRUCTURE_CHECK", "ON") if "STRUCTURE_CHECK" in params else "ON"
        if not isinstance(sc, str):
            raise _Refused("TMS: field STRUCTURE_CHECK has the wrong type")
        if sc not in ("ON", "OFF"):
            raise _Refused("TMS: STRUCTURE_CHECK " + sc + " is not ON or OFF")
        if sc == "OFF":
            return ["STRUCTURE", NOT_VALIDATED, ["CHECK OFF"]]
        K = _read(tfk)
        fnd = _analyse(K, K["laws"])[1]
    except _Refused as e:
        return _refuse(str(e))
    present = set(f[0] for f in fnd)
    reasons = []
    if any(c in present for c in FINDING_CLASSES):
        reasons.append("FINDINGS")
    for r in ("PROSE", "PDE", "NO FORM"):
        if r in present:
            reasons.append(r)
    if not reasons:
        return ["STRUCTURE", REGULAR]
    return ["STRUCTURE", NOT_VALIDATED, reasons]


def _id_key(i):
    p, n = i.rsplit("-", 1)
    return (p, len(n), n)


def eom_register_test(tfk, records):
    """["OK", R, findings, ["USED", ids], ["SETTLED", unknowns], ["OPEN",
    unknowns]] or ["REFUSED", reason]. records: [{"ID", "STATUS", "ROLE"?,
    "EQUATION"?}], the provisional shape C-16 will carry (fixtures only);
    the TFK is read first, then every record whole."""
    try:
        K = _read(tfk)
        entering = _read_records(records, K)
        _, _, _, under0 = _analyse(K, K["laws"])
        R, fnd, cu, _ = _analyse(K, K["laws"] + entering)
    except _Refused as e:
        return _refuse(str(e))
    settled = [] if cu is None else sorted(u for u in under0 if u not in cu)
    still = sorted(u for u in under0 if u not in settled)
    used = sorted((l["name"] for l in entering), key=_id_key)
    return ["OK", R, fnd, ["USED", used], ["SETTLED", settled], ["OPEN", still]]