"""The second-order test of the Kuhn lattice: Theorem 4.3, the vertex symmetry fails at second order.

Paper: "Quantum Gravity from the Recognition Ledger: A Finite One-Loop Path Sum, Einstein Gravity at Long Wavelengths, and the
Massless-Graviton Problem", J. Washburn (2026), version 1.1, https://doi.org/10.5281/zenodo.23813297. Section and theorem numbers
below are that version's.

The question
  Expand a solution of the Regge equations in its amplitude, l = l0 + t h + (t^2/2) d + ..., about the flat lattice l0. The first
  order says that h solves the linearized equations, D2S(h, .) = 0. The second order asks for an edge field d with
      D2S(d, .) + D3S(h, h, .) = 0.
  The flat Hessian D2S annihilates every vertex move g, because flat configurations stay flat when a vertex is moved. Pairing the
  second-order equation with g therefore kills its first term, so a solution d can exist only if D3S(h, h, g) = 0 for every vertex
  move g: the second-order form of the vertex symmetry. Theorem 4.3 shows that on the Kuhn lattice at kappa = 1 this fails, with
  D3S(h, h, g) = 756 sqrt3 for one explicit h and g. This program recomputes that number from the Regge action.

What it computes
  * The lattice. The Kuhn triangulation of Z^4 at kappa = 1: each unit cell is cut into 24 four-simplices, one for each order in
    which the four counters advance; an edge of class d in {0,1}^4, d != 0 (fifteen classes) has squared length |d|(1 - |d|).
  * The action's derivatives. The Regge action is S = sum_h A_h delta_h over the hinges (triangles) h, with A_h the area and
    delta_h the deficit angle. Schlaefli's identity, sum_h A_h d(delta_h) = 0, removes every term carrying the constant 2 pi at the
    flat lattice (where every deficit is zero), and leaves forms that are sums over four-simplices s of
        D2S = - sum_s sum_{h in s} theta'_{h,s} (x) A'_h,
        D3S(u, v, w) = - sum_s sum_{h in s} [theta''_{h,s}(u, v) A'_h(w) + theta'_{h,s}(v) A''_h(u, w) + theta'_{h,s}(u) A''_h(v, w)],
    symmetrized, where theta_{h,s} is the dihedral angle of s at h and the primes are derivatives with respect to the ten squared
    edge lengths of s. The program evaluates them by third-order Taylor arithmetic on the exact flat simplex, in double precision.
    Every entry of 16 sqrt3 times the quadratic form and of 96 sqrt3 times the cubic form is an integer (checked); the program
    rounds them and from then on works in exact rational arithmetic. (Section 4 says the proof uses "the cubic cell form of the
    Regge action at kappa = 1, whose numerators are integers over 1152 sqrt3". Integers over 96 sqrt3 are also integers over
    1152 sqrt3, since 1152 = 12 x 96, so these forms, and their sums over a cell, are consistent with it.)
  * The field h of Theorem 4.3: the sampled weak-field metric of the static potential Phi = xi1^3 - 3 xi1 xi2^2, with
    xi1 = (x0 + x1 - x2 - x3)/2 and xi2 = (x0 - x1 + x2 - x3)/2, read at the edge midpoints, plus its slaved correction. The paper
    does not print the weak-field metric; the program reads it as the linearized field of a static potential with -h_tt/2 = Phi,
    the normalization of Section 4: h_ab = -2 Phi (eta_ab + 2 u_a u_b), with u = (1,1,1,1)/sqrt12 the unit timelike vector
    (eta(u, u) = -1), so that an edge of class d gets h = -2 Phi (|d| + |d|^2/2).
  * The slaved correction. Section 4 describes it as lying in the five non-metric combinations of edge classes and as linear in
    the tidal field. The metric part of the class equations is the symmetric matrix sum_d d d^T Phi_d (Phi_d the equation of class
    d), ten of the fifteen combinations (Section 4); the program takes the non-metric combinations to be the class patterns c with
    sum_d c_d d d^T = 0. It constructs the correction from that description: for each harmonic quadratic potential, whose tidal
    field is constant (xi3 = (x0 - x1 - x2 + x3)/2 completes the spatial coordinates), it finds the constant field in the
    non-metric combinations that makes the sampled field an exact solution of the linearized equations on the thirty edges at a
    record. That solution exists and is unique; it is zero for the two diagonal tidal components and nonzero for the three
    off-diagonal ones. The program takes the correction of the cubic potential to be this linear map applied to its tidal field at
    each edge midpoint. The construction is the program's, made from the paper's description; the paper does not print the
    correction.
  * The record y = (2, 1, -1, 0) and the vertex move g of y along the first counter, in the convention of the paper's Appendix A
    (a displacement eta changes the squared length of an edge of class d by 2 d^T eta at its ends): an edge p -> q of class d
    changes by 2 sum_a d_a (zeta_a(q) - zeta_a(p)), with zeta = e_0 at y and zero elsewhere. Theorem 4.3 states the theorem on the
    torus of side 8 with every coordinate read in {-4, ..., 3}; the statement is local at y, and the edges and simplices it uses
    carry the same values on Z^4, where the program works (Section 4: it holds on the infinite lattice too).

Checks (printed, each with PASS or FAIL; all exact except the first two)
  * The flat lattice: the dihedral angles at every hinge add to 2 pi, so every flat deficit vanishes; and the 24 simplex types
    have the same forms in their own edge order.
  * The forms are integers over 16 sqrt3 (quadratic) and 96 sqrt3 (cubic).
  * The slaved correction: the five non-metric combinations give rank five on the thirty equations; the correction found
    satisfies all thirty exactly for each of the five harmonic quadratics.
  * Theorem 4.3, with the program's correction: h solves the linearized equations on the thirty edges at y exactly (without
    the correction it does not); the flat Hessian annihilates g; D3S(h, h, g) = 756 sqrt3 exactly. Hence no edge field d solves
    the second-order equations on the thirty edges at y.

Output
  With --json FILE, writes the thirty linearized equations at y for the sampled field without and with its correction, the
  correction's class pattern, and D3S(h, h, g) for the moves of y along each counter.

Requires Python 3.9+ and numpy. Run: python3 vertex_symmetry.py   [--json vertex.json]
License: MIT. Copyright (c) 2026 Recognition Physics Institute.
"""
import argparse
import itertools
import json
import sys
from fractions import Fraction as Fr

import numpy as np

PERMS = list(itertools.permutations(range(4)))
PAIRS = [(0, 1), (0, 2), (0, 3), (0, 4), (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]
TRI = list(itertools.combinations(range(5), 3))
CLASSES = [frozenset(c) for r in range(1, 5) for c in itertools.combinations(range(4), r)]
RESULTS = []


def check(name, ok, detail):
    RESULTS.append(ok)
    print(f"[{'PASS' if ok else 'FAIL'}] {name}: {detail}")


def eta(v, w):
    """the flat metric I - 11^T at kappa = 1"""
    return sum(a * b for a, b in zip(v, w)) - sum(v) * sum(w)


def verts(cell, s):
    v = [tuple(cell)]
    cur = list(cell)
    for a in PERMS[s]:
        cur = cur[:]
        cur[a] += 1
        v.append(tuple(cur))
    return v


# ---------------------------------------------------------------- third-order Taylor arithmetic in the ten squared lengths
def sym3(t):
    return (t + t.transpose(0, 2, 1) + t.transpose(1, 0, 2) + t.transpose(1, 2, 0) + t.transpose(2, 0, 1)
            + t.transpose(2, 1, 0)) / 6


class Jet:
    """f(l0 + e) = a0 + a1.e + e.a2.e + a3[e,e,e] to third order (a2, a3 symmetric)."""
    __slots__ = ('a0', 'a1', 'a2', 'a3')

    def __init__(self, a0, a1=None, a2=None, a3=None):
        self.a0 = float(a0)
        self.a1 = np.zeros(10) if a1 is None else a1
        self.a2 = np.zeros((10, 10)) if a2 is None else a2
        self.a3 = np.zeros((10, 10, 10)) if a3 is None else a3

    def __add__(self, o):
        if not isinstance(o, Jet):
            return Jet(self.a0 + o, self.a1, self.a2, self.a3)
        return Jet(self.a0 + o.a0, self.a1 + o.a1, self.a2 + o.a2, self.a3 + o.a3)

    def __neg__(self):
        return Jet(-self.a0, -self.a1, -self.a2, -self.a3)

    def __sub__(self, o):
        return self + (-o)

    def __mul__(self, o):
        if not isinstance(o, Jet):
            return Jet(self.a0 * o, self.a1 * o, self.a2 * o, self.a3 * o)
        a2 = self.a0 * o.a2 + o.a0 * self.a2 + 0.5 * (np.outer(self.a1, o.a1) + np.outer(o.a1, self.a1))
        a3 = self.a0 * o.a3 + o.a0 * self.a3 + sym3(np.einsum('i,jk->ijk', self.a1, o.a2) + np.einsum('i,jk->ijk', o.a1, self.a2))
        return Jet(self.a0 * o.a0, self.a0 * o.a1 + o.a0 * self.a1, a2, a3)

    def apply(self, f0, f1, f2, f3):
        d = Jet(0.0, self.a1, self.a2, self.a3)
        d2 = d * d
        return Jet(f0) + d * f1 + d2 * (f2 / 2) + (d2 * d) * (f3 / 6)

    def inv(self):
        x = self.a0
        return self.apply(1 / x, -1 / x ** 2, 2 / x ** 3, -6 / x ** 4)

    def sqrt(self):
        x = self.a0
        r = np.sqrt(x)
        return self.apply(r, 1 / (2 * r), -1 / (4 * r * x), 3 / (8 * r * x * x))

    def acos(self):
        x = self.a0
        u = 1 - x * x
        return self.apply(np.arccos(x), -1 / np.sqrt(u), -x / u ** 1.5, -(1 + 2 * x * x) / u ** 2.5)


def simplex_forms(s):
    """For simplex type s at the flat lattice: its dihedral angles and the Schlaefli forms D2S_s (10 x 10), D3S_s (10 x 10 x 10)."""
    V = verts((0, 0, 0, 0), s)
    L = {}
    for k, (i, j) in enumerate(PAIRS):
        e = [V[j][a] - V[i][a] for a in range(4)]
        a1 = np.zeros(10)
        a1[k] = 1.0
        L[(i, j)] = Jet(eta(e, e), a1)

    def ell(i, j):
        return Jet(0.0) if i == j else L[(min(i, j), max(i, j))]

    q2 = np.zeros((10, 10))
    q3 = np.zeros((10, 10, 10))
    angles = {}
    for tri in TRI:
        o, i1, i2 = tri
        a3, a4 = [k for k in range(5) if k not in tri]

        def gram(i, j):
            return (ell(o, i) + ell(o, j) - ell(i, j)) * 0.5
        g = [[gram(u, v) for v in (i1, i2)] for u in (i1, i2)]
        det = g[0][0] * g[1][1] - g[0][1] * g[1][0]
        if det.a0 >= 0:
            raise SystemExit('a flat hinge is not timelike')
        area = (-det).sqrt() * 0.5
        idet = det.inv()
        ginv = [[g[1][1] * idet, -g[0][1] * idet], [-g[1][0] * idet, g[0][0] * idet]]
        base = (i1, i2)

        def proj(a, b):
            r = gram(a, b)
            for u in range(2):
                for v in range(2):
                    r = r - gram(a, base[u]) * ginv[u][v] * gram(base[v], b)
            return r
        theta = (proj(a3, a4) * (proj(a3, a3) * proj(a4, a4)).sqrt().inv()).acos()
        angles[frozenset(V[k] for k in tri)] = theta.a0
        dA, d2A, dT, d2T = area.a1, 2 * area.a2, theta.a1, 2 * theta.a2
        q2 -= np.outer(dT, dA)
        q3 -= np.einsum('uv,w->uvw', d2T, dA) + np.einsum('v,uw->uvw', dT, d2A) + np.einsum('u,vw->uvw', dT, d2A)
    return angles, (q2 + q2.T) / 2, sym3(q3)


# ---------------------------------------------------------------- exact linear algebra over the rationals
def nullspace(rows, n):
    """basis of {x in Q^n : r.x = 0 for every r in rows}"""
    m = [[Fr(x) for x in r] for r in rows]
    piv, r = [], 0
    for c in range(n):
        p = next((i for i in range(r, len(m)) if m[i][c] != 0), None)
        if p is None:
            continue
        m[r], m[p] = m[p], m[r]
        m[r] = [x / m[r][c] for x in m[r]]
        for i in range(len(m)):
            if i != r and m[i][c] != 0:
                f = m[i][c]
                m[i] = [a - f * b for a, b in zip(m[i], m[r])]
        piv.append(c)
        r += 1
    free = [c for c in range(n) if c not in piv]
    basis = []
    for fc in free:
        x = [Fr(0)] * n
        x[fc] = Fr(1)
        for i, pc in enumerate(piv):
            x[pc] = -m[i][fc]
        basis.append(x)
    return basis


def solve_exact(A, b):
    """the unique solution of A x = b (A: m x n with rank n, m >= n), or None; also returns the rank"""
    m, n = len(A), len(A[0])
    M = [[Fr(x) for x in A[i]] + [Fr(b[i])] for i in range(m)]
    r, piv = 0, []
    for c in range(n):
        p = next((i for i in range(r, m) if M[i][c] != 0), None)
        if p is None:
            continue
        M[r], M[p] = M[p], M[r]
        M[r] = [x / M[r][c] for x in M[r]]
        for i in range(m):
            if i != r and M[i][c] != 0:
                f = M[i][c]
                M[i] = [a - f * bb for a, bb in zip(M[i], M[r])]
        piv.append(c)
        r += 1
    if r < n or any(M[i][n] != 0 for i in range(r, m)):
        return None, r
    return [M[i][n] for i in range(n)], r


# ---------------------------------------------------------------- fields on edges
def cls(p, q):
    return frozenset(a for a in range(4) if q[a] != p[a])


def midpoint(p, q):
    return [Fr(p[a] + q[a], 2) for a in range(4)]


def xi(x):
    return (Fr(x[0] + x[1] - x[2] - x[3], 2), Fr(x[0] - x[1] + x[2] - x[3], 2), Fr(x[0] - x[1] - x[2] + x[3], 2))


def memo(fun):
    """an edge field with its exact values cached"""
    cache = {}

    def f(p, q):
        key = (p, q)
        if key not in cache:
            cache[key] = Fr(fun(p, q))
        return cache[key]
    return f


def sampled(phi):
    def h(p, q):
        n = len(cls(p, q))
        return -2 * phi(xi(midpoint(p, q))) * (n + Fr(n * n, 2))
    return memo(h)


def edges_at(z):
    out = []
    for S in CLASSES:
        v = [1 if a in S else 0 for a in range(4)]
        out.append((tuple(z), tuple(z[a] + v[a] for a in range(4))))
        out.append((tuple(z[a] - v[a] for a in range(4)), tuple(z)))
    return out


PREFIX = [set(verts((0, 0, 0, 0), s)) for s in range(24)]
_STAR = {}


def simplices_containing(points):
    """the four-simplices (as vertex lists) that contain every given lattice point"""
    key = tuple(sorted(tuple(p) for p in points))
    if key not in _STAR:
        lo = [max(p[a] for p in key) - 1 for a in range(4)]
        hi = [min(p[a] for p in key) for a in range(4)]
        out = []
        for cell in itertools.product(*[range(lo[a], hi[a] + 1) for a in range(4)]):
            offs = [tuple(p[a] - cell[a] for a in range(4)) for p in key]
            for s in range(24):
                if all(o in PREFIX[s] for o in offs):
                    out.append(verts(cell, s))
        _STAR[key] = out
    return _STAR[key]


class Forms:
    """The exact Schlaefli forms: D2S = N2 / (16 sqrt3) and D3S = N3 / (96 sqrt3) per simplex, in its own edge order."""

    def __init__(self, n2, n3):
        self.n2, self.n3 = n2, n3
        self.nz3 = [(e, f, k, int(n3[e, f, k])) for e in range(10) for f in range(10) for k in range(10) if n3[e, f, k] != 0]

    def row(self, field, edge):
        """16 sqrt3 times D2S(field, delta_edge), exactly"""
        tot = Fr(0)
        for V in simplices_containing(edge):
            local = [(V[i], V[j]) for (i, j) in PAIRS]
            k = local.index(edge)
            tot += sum(int(self.n2[k, f]) * field(*local[f]) for f in range(10) if self.n2[k, f] != 0)
        return tot

    def rows(self, field, z):
        return [self.row(field, e) for e in edges_at(z)]

    def cubic(self, field, move, z):
        """D3S(field, field, move) / sqrt3, exactly, for a move supported on the edges at z"""
        tot = 0
        for V in simplices_containing([z]):
            local = [(V[i], V[j]) for (i, j) in PAIRS]
            h = [field(*e) * 64 for e in local]          # every value of these fields is a multiple of 1/64
            if any(x.denominator != 1 for x in h):
                raise SystemExit('a field value is not a multiple of 1/64')
            h = [int(x) for x in h]
            g = [int(move(*e)) for e in local]
            tot += sum(c * h[e] * h[f] * g[k] for (e, f, k, c) in self.nz3 if g[k] != 0)
        return Fr(tot, 288 * 64 * 64)          # (1 / (96 sqrt3)) / sqrt3 = 1/288


def vertex_move(y, axis):
    """the move of the record y along one counter, in the convention of the paper's Appendix A: 2 d^T (zeta(q) - zeta(p))"""
    def g(p, q):
        d = [q[a] - p[a] for a in range(4)]
        zq = 1 if tuple(q) == tuple(y) else 0
        zp = 1 if tuple(p) == tuple(y) else 0
        return 2 * d[axis] * (zq - zp)
    return memo(g)


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('--json', help='write the data of the page figure to this file')
    args = ap.parse_args()
    r3 = np.sqrt(3.0)

    # 1. the flat simplex forms
    per_type = [simplex_forms(s) for s in range(24)]
    angle_of = {}
    for s in range(24):
        for key, th in per_type[s][0].items():
            angle_of[(s, tuple(sorted(key)))] = th
    hinge_sum = {}
    origin = (0, 0, 0, 0)
    for V in simplices_containing([origin]):
        rel_v = [tuple(v[a] - V[0][a] for a in range(4)) for v in V]
        s = next(t for t in range(24) if rel_v == verts(origin, t))
        for tri in TRI:
            pts = [V[k] for k in tri]
            if origin in pts:      # every simplex containing such a hinge contains the origin, so its sum is complete
                rel = tuple(sorted(tuple(p[a] - V[0][a] for a in range(4)) for p in pts))
                key = frozenset(pts)
                hinge_sum[key] = hinge_sum.get(key, 0.0) + angle_of[(s, rel)]
    worst = max(abs(v - 2 * np.pi) for v in hinge_sum.values())
    check('Flat lattice', worst < 1e-12,
          f'around each of the {len(hinge_sum)} hinges at a record the dihedral angles add to 2 pi (largest deviation {worst:.1e}), '
          'so every flat deficit vanishes')
    q2s = np.array([t[1] for t in per_type])
    q3s = np.array([t[2] for t in per_type])
    same = max(np.max(np.abs(q2s - q2s[0])), np.max(np.abs(q3s - q3s[0])))
    check('Symmetry', same < 1e-12, f'the 24 simplex types have the same forms in their own edge order (largest difference {same:.1e})')

    # 2. exact forms
    a2, a3 = 16 * r3 * q2s[0], 96 * r3 * q3s[0]
    n2, n3 = np.round(a2), np.round(a3)
    e2, e3 = np.max(np.abs(a2 - n2)), np.max(np.abs(a3 - n3))
    check('Exact forms', e2 < 1e-9 and e3 < 1e-9,
          f'16 sqrt3 D2S and 96 sqrt3 D3S per simplex are integer arrays (distances {e2:.1e}, {e3:.1e}; largest entries '
          f'{int(np.max(np.abs(n2)))} and {int(np.max(np.abs(n3)))}); exact arithmetic from here on')
    forms = Forms(n2.astype(int), n3.astype(int))

    # 3. the slaved correction, derived
    metric = [[(1 if a == b else 2) if (a in S and b in S) else 0 for S in CLASSES] for a in range(4) for b in range(a, 4)]
    nonmetric = nullspace(metric, 15)
    basis_rows = [forms.rows(memo(lambda p, q, v=v: v[CLASSES.index(cls(p, q))]), origin) for v in nonmetric]
    A = [[basis_rows[k][i] for k in range(5)] for i in range(30)]
    quadratics = [('xi1 xi2', lambda x: x[0] * x[1], (0, 1)), ('xi1 xi3', lambda x: x[0] * x[2], (0, 2)),
                  ('xi2 xi3', lambda x: x[1] * x[2], (1, 2)), ('xi1^2 - xi2^2', lambda x: x[0] ** 2 - x[1] ** 2, None),
                  ('xi1^2 + xi2^2 - 2 xi3^2', lambda x: x[0] ** 2 + x[1] ** 2 - 2 * x[2] ** 2, None)]
    patterns, ranks, report = {}, [], []
    for name, phi, pair in quadratics:
        r = forms.rows(sampled(phi), origin)
        gamma, rank = solve_exact(A, [-x for x in r])
        ranks.append(rank)
        if gamma is None:
            check('Slaved correction', False, f'no constant non-metric correction solves the equations for {name}')
            continue
        pattern = [sum(gamma[k] * nonmetric[k][i] for k in range(5)) for i in range(15)]
        if pair:
            patterns[pair] = pattern
        nz = sorted(set(abs(x) for x in pattern if x != 0))
        report.append(f"{name}: {'none needed' if not nz else 'values ' + ', '.join('+-' + str(x) for x in nz)}")
    check('Slaved correction', len(patterns) == 3 and all(rk == 5 for rk in ranks),
          'the five non-metric combinations give rank 5 on the thirty equations at a record, so the constant correction of a '
          'constant tidal field is unique; it solves all thirty exactly: ' + '; '.join(report) +
          ' (per unit of the tidal component d_i d_j Phi)')

    # 4. Theorem 4.3
    def phi(x):
        return x[0] ** 3 - 3 * x[0] * x[1] ** 2

    def tidal(x):          # the second derivatives of Phi at x: (11, 22, 33, 12, 13, 23)
        return {(0, 1): -6 * x[1], (0, 2): Fr(0), (1, 2): Fr(0)}

    def correction(p, q):
        t = tidal(xi(midpoint(p, q)))
        i = CLASSES.index(cls(p, q))
        return sum(t[pair] * patterns[pair][i] for pair in patterns)

    bare = sampled(phi)

    field = memo(lambda p, q: bare(p, q) + correction(p, q))
    y = (2, 1, -1, 0)
    rows_bare = forms.rows(bare, y)
    rows_full = forms.rows(field, y)
    scale = 16 * r3
    check('Theorem 4.3 (i), with the program\'s correction', all(x == 0 for x in rows_full) and any(x != 0 for x in rows_bare),
          f'h solves the linearized Regge equations on the thirty edges at y = (2, 1, -1, 0) exactly; the sampled field alone '
          f'does not (largest equation {max(abs(float(x)) for x in rows_bare) / scale:.4f})')
    g = vertex_move(y, 0)
    star_edges = sorted({(V[i], V[j]) for V in simplices_containing([y]) for (i, j) in PAIRS})
    hg = all(forms.row(g, e) == 0 for e in star_edges)
    check('Theorem 4.3 (ii)', hg, f'the flat Hessian annihilates the vertex move g of y along the first counter: its equation '
          f'vanishes exactly on each of the {len(star_edges)} edges of the simplices at y, the only edges it reaches')
    values = [forms.cubic(field, vertex_move(y, a), y) for a in range(4)]
    check('Theorem 4.3 (iii)', values[0] == 756,
          f'D3S(h, h, g) = {values[0]} sqrt3 exactly, so no edge field d solves the second-order equations on the thirty edges '
          f'at y; the moves along the other counters give {values[1]}, {values[2]} and {values[3]} times sqrt3')
    print()
    print(f"{sum(RESULTS)} of {len(RESULTS)} checks passed. The forms come from a double-precision evaluation rounded to "
          'integers, which the program checks; everything after that is exact. The slaved correction is constructed by the '
          "program from the paper's description of it.")
    if args.json:
        out = {'program': 'vertex_symmetry.py', 'paper': 'https://doi.org/10.5281/zenodo.23813297', 'record': list(y),
               'edges': [[list(p), list(q)] for (p, q) in edges_at(y)],
               'rows_sampled': [float(x) / scale for x in rows_bare], 'rows_corrected': [float(x) / scale for x in rows_full],
               'correction_per_unit_tidal': {f'{i + 1}{j + 1}': [str(x) for x in pat] for (i, j), pat in patterns.items()},
               'classes': [''.join(str(a) for a in sorted(S)) for S in CLASSES],
               'd3s_over_sqrt3': [str(v) for v in values]}
        with open(args.json, 'w') as fh:
            json.dump(out, fh, indent=1)
        print(f'wrote {args.json}')
    return 0 if all(RESULTS) else 1


if __name__ == '__main__':
    sys.exit(main())
