SAT im Phasenraum

"""
SAT als kontinuierliches dynamisches System (Ercsey-Ravasz/Toroczkai 2011).

DAS SYSTEM. Variablen werden entspannt: s_i in [-1,1] statt {-1,+1}.
Die Verletzung einer Klausel m ist

    K_m(s) = 2^-k * prod_{i in m} (1 - c_mi s_i)      in [0,1]

K_m = 0 genau dann, wenn ein Literal voll erfuellt ist. Jede Klausel bekommt
ein Gewicht a_m, das exponentiell waechst, solange sie verletzt ist:

    ds_i/dt = 2 * sum_m a_m K_m K_mi c_mi        (K_mi = K_m ohne den Faktor i)
    da_m/dt = a_m K_m

Der s-Teil ist Gradientenabstieg auf V = sum_m a_m K_m^2; der a-Teil hebt
genau die Klauseln an, die haengen bleiben, und macht damit jedes lokale
Minimum instabil. Folge: **es gibt keine Attraktoren ausser den Loesungen.**
Der Wuerfel [-1,1]^n ist invariant.

Der Preis: der Weg dorthin ist ein CHAOTISCHER TRANSIENT. Die Zeit bis zur
Loesung ist unbeschraenkt, die Einzugsgebiete der Loesungen haben FRAKTALE
Raender, und die Rate, mit der Trajektorien das chaotische Gebiet verlassen
(die Fluchtrate kappa), geht an der SAT/UNSAT-Schwelle gegen null. Haerte
wird damit zu einer dynamischen Groesse: 1/kappa.

Alles hier ist ueber das Gitter vektorisiert -- ein RK4-Schritt bewegt
hunderttausend Trajektorien gleichzeitig.
"""
import struct
import zlib

import numpy as np


# ------------------------------------------------------------------ Aufbau

def bau(klauseln, n):
    """lit[m,k] Variablenindex, sgn[m,k] Vorzeichen, streu[3m,n] Streumatrix."""
    k = len(klauseln[0])
    m = len(klauseln)
    lit = np.zeros((m, k), np.int32)
    sgn = np.zeros((m, k), np.float32)
    for j, c in enumerate(klauseln):
        for l, x in enumerate(c):
            lit[j, l] = abs(x) - 1
            sgn[j, l] = 1.0 if x > 0 else -1.0
    streu = np.zeros((m * k, n), np.float32)
    streu[np.arange(m * k), lit.ravel()] = 1.0
    return lit, sgn, streu


def ableitung(s, a, lit, sgn, streu):
    """ds/dt und dK. s: (P,n), a: (P,M). Gibt (ds, K) zurueck."""
    P = s.shape[0]
    M, k = lit.shape
    f = 1.0 - sgn[None, :, :] * s[:, lit]          # (P,M,k), jeder Faktor in [0,2]
    K = f.prod(axis=2) * (0.5 ** k)                # (P,M)
    # K_mi = K ohne den Faktor i  ->  Produkt der uebrigen
    Kmi = np.empty_like(f)
    for l in range(k):
        andere = [x for x in range(k) if x != l]
        p = f[:, :, andere[0]]
        for x in andere[1:]:
            p = p * f[:, :, x]
        Kmi[:, :, l] = p * (0.5 ** (k - 1))
    beitrag = (2.0 * a * K)[:, :, None] * Kmi * sgn[None, :, :]
    ds = beitrag.reshape(P, M * k) @ streu
    return ds, K


def erfuellt(s, lit, sgn):
    """Ist sign(s) eine Loesung? (P,) bool"""
    z = np.sign(s)
    z[z == 0] = 1.0
    ok = (sgn[None, :, :] * z[:, lit]) > 0        # (P,M,k)
    return ok.any(axis=2).all(axis=1)


def kode(s):
    """Vorzeichenvektor als Ganzzahl -- Kennung der erreichten Loesung."""
    b = (s > 0).astype(np.int64)
    g = np.zeros(s.shape[0], np.int64)
    for j in range(s.shape[1]):
        g |= b[:, j] << j
    return g


# ------------------------------------------------------------------ Integration

def laufe(s0, klauseln, n, tmax=200.0, eta=0.05, dtmax=0.5, pruef=1,
          maxschritt=3000, adeckel=1e8):
    """Integriert alle Startpunkte gleichzeitig (RK4, punktweise Schrittweite).

    ZWEI Abbruchbedingungen, und beide werden gebraucht: die analoge Zeit t
    UND die Schrittzahl. Denn die Schrittweite ist ~ 1/|ds| ~ 1/a, und a
    waechst exponentiell -- im chaotischen Transienten friert t praktisch ein,
    waehrend die Trajektorie in s weiterwandert. Ohne die Schrittgrenze laeuft
    so ein Punkt endlos. Die Schrittzahl misst die BOGENLAENGE der Bahn und ist
    damit das ehrlichere Mass fuer den Rechenaufwand.

    Gibt (zeit, loesung, fertig) zurueck.
    """
    lit, sgn, streu = bau(klauseln, n)
    M = len(klauseln)
    P = s0.shape[0]
    s = s0.astype(np.float32).copy()
    a = np.ones((P, M), np.float32)
    t = np.zeros(P, np.float32)
    bogen = np.zeros(P, np.int32)
    zeit = np.full(P, np.nan, np.float32)
    schritte = np.full(P, -1, np.int32)
    loes = np.full(P, -1, np.int64)
    offen = np.ones(P, bool)
    runde = 0
    while offen.any() and runde < maxschritt:
        runde += 1
        idx = np.flatnonzero(offen)
        ss, aa = s[idx], a[idx]
        d1, K1 = ableitung(ss, aa, lit, sgn, streu)
        dt = np.clip(eta / (1e-6 + np.abs(d1).max(axis=1)), 1e-5, dtmax)[:, None]
        d2, K2 = ableitung(ss + 0.5 * dt * d1, aa * np.exp(0.5 * dt * K1), lit, sgn, streu)
        d3, K3 = ableitung(ss + 0.5 * dt * d2, aa * np.exp(0.5 * dt * K2), lit, sgn, streu)
        d4, K4 = ableitung(ss + dt * d3, aa * np.exp(dt * K3), lit, sgn, streu)
        ss = np.clip(ss + (dt / 6.0) * (d1 + 2 * d2 + 2 * d3 + d4), -1.0, 1.0)
        aa = np.minimum(aa * np.exp((dt / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)), adeckel)
        s[idx], a[idx] = ss, aa
        t[idx] += dt[:, 0]
        bogen[idx] += 1
        if runde % pruef == 0:
            fertig = erfuellt(ss, lit, sgn)
            if fertig.any():
                g = idx[fertig]
                zeit[g] = t[g]
                schritte[g] = bogen[g]
                loes[g] = kode(s[g])
                offen[g] = False
        raus = idx[t[idx] > tmax]
        offen[raus] = False
    return zeit, loes, np.isfinite(zeit), schritte


# ------------------------------------------------------------------ PNG

def png(pfad, bild):
    H, W, _ = bild.shape
    roh = b"".join(b"\x00" + bild[y].tobytes() for y in range(H))

    def brocken(typ, daten):
        return (struct.pack(">I", len(daten)) + typ + daten
                + struct.pack(">I", zlib.crc32(typ + daten) & 0xFFFFFFFF))

    d = b"\x89PNG\r\n\x1a\n"
    d += brocken(b"IHDR", struct.pack(">IIBBBBB", W, H, 8, 2, 0, 0, 0))
    d += brocken(b"IDAT", zlib.compress(roh, 6))
    d += brocken(b"IEND", b"")
    with open(pfad, "wb") as f:
        f.write(d)
    return pfad

"""
Woher der Fraktal kommt: die Streckrate lambda, und die Probe kappa/lambda.

DIE KETTE. Der Rand zwischen zwei Einzugsgebieten ist unter dem Fluss
INVARIANT -- wer genau auf ihm startet, konvergiert nie, denn er muesste den
Rand verlassen, und Raender gehen unter einem Fluss wieder in Raender ueber.
Auf dem Rand lebt also eine invariante Menge ohne Konvergenz: der chaotische
SATTEL. Der Rand selbst ist dessen STABILE MANNIGFALTIGKEIT.

Ein Sattel hat streckende und stauchende Richtungen. Weil der Wuerfel
[-1,1]^n beschraenkt und invariant ist, kann die Streckung nicht davonlaufen
-- die Bahn muss zurueckgefaltet werden. Strecken und Falten, wiederholt, ist
das Hufeisen: quer zur stabilen Richtung entsteht eine Cantormenge. Genau die
sieht man im 500-fachen Zoom als Lamination.

DIE PROBE. Fuer einen chaotischen Sattel verknuepft die Beziehung von
Kantz und Grassberger (1985) die drei Groessen:

    alpha  =  kappa / lambda          und   D_Rand = D_Raum - alpha

  kappa   Fluchtrate -- wie schnell Bahnen den Sattel verlassen
  lambda  Lyapunov-Exponent -- wie schnell Nachbarn auseinanderlaufen

Anschaulich: die Streckung erzeugt in jeder Zeiteinheit neue Unsicherheit,
die Flucht raeumt sie ab. Das Verhaeltnis ist die Fraktalitaet. Ist die
Streckung viel schneller als die Flucht, ist alpha klein und der Rand fast
flaechenfuellend.

Gemessen wird lambda nach Benettin: Paare mit Anfangsabstand delta0, gleicher
Schrittweite, regelmaessig auf delta0 zurueckskaliert; lambda ist das Mittel
von ln(Zuwachs) je Zeiteinheit -- gezaehlt nur, solange beide Bahnen noch im
Transienten sind.

VORHERSAGE, VOR DEM LAUF NOTIERT. Fuer diese Instanz wurden alpha = 0,273 und
(fuer alpha_Dichte = 4,6) kappa ungefaehr 0,67 gemessen. Dann muss lambda bei
etwa 0,67/0,273 = 2,5 liegen. Trifft das, ist die Herkunft des Fraktals nicht
nur erzaehlt, sondern geschlossen.
"""
import argparse

import numpy as np

import chaos_sat as C
import gruppen as G


def paarlauf(s0, klauseln, n, delta0=1e-6, eta=0.05, schritte=4000,
             normiere=20, rng=None):
    """Benettin: P Paare, gemeinsame Schrittweite je Paar."""
    lit, sgn, streu = C.bau(klauseln, n)
    M = len(klauseln)
    P = s0.shape[0]
    rng = rng or np.random.default_rng(0)

    d0 = rng.normal(size=(P, n)).astype(np.float32)
    d0 /= np.linalg.norm(d0, axis=1, keepdims=True)
    s = np.repeat(s0.astype(np.float32), 2, axis=0)
    s[1::2] = np.clip(s[1::2] + delta0 * d0, -1, 1)
    a = np.ones((2 * P, M), np.float32)

    summe = np.zeros(P)          # aufsummiertes ln(Zuwachs)
    zeit = np.zeros(P)           # zugehoerige analoge Zeit
    lebt = np.ones(P, bool)

    for k in range(schritte):
        if not lebt.any():
            break
        d1, K1 = C.ableitung(s, a, lit, sgn, streu)
        # Schrittweite aus der Referenzbahn, fuer beide Partner dieselbe
        dtr = np.clip(eta / (1e-6 + np.abs(d1[0::2]).max(axis=1)), 1e-5, 0.5)
        dt = np.repeat(dtr, 2)[:, None]
        d2, K2 = C.ableitung(s + 0.5 * dt * d1, a * np.exp(0.5 * dt * K1), lit, sgn, streu)
        d3, K3 = C.ableitung(s + 0.5 * dt * d2, a * np.exp(0.5 * dt * K2), lit, sgn, streu)
        d4, K4 = C.ableitung(s + dt * d3, a * np.exp(dt * K3), lit, sgn, streu)
        s = np.clip(s + (dt / 6.0) * (d1 + 2 * d2 + 2 * d3 + d4), -1.0, 1.0)
        a = np.minimum(a * np.exp((dt / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)), 1e8)

        fertig = C.erfuellt(s[0::2], lit, sgn) | C.erfuellt(s[1::2], lit, sgn)
        lebt &= ~fertig
        zeit += np.where(lebt, dtr, 0.0)

        if (k + 1) % normiere == 0:
            diff = s[1::2] - s[0::2]
            d = np.linalg.norm(diff, axis=1)
            gut = lebt & (d > 0)
            summe[gut] += np.log(d[gut] / delta0)
            skal = np.where(d > 0, delta0 / np.maximum(d, 1e-30), 1.0)[:, None]
            s[1::2] = np.clip(s[0::2] + diff * skal, -1, 1)

    brauchbar = zeit > 1.0
    lam = summe[brauchbar] / zeit[brauchbar]
    return lam, zeit[brauchbar]


def lauf(args):
    n, alpha, seed = args.n, args.alpha, args.seed
    kl = G.zufalls_cnf(n, int(round(alpha * n)), 3, np.random.default_rng(1000 * n + seed))
    rng = np.random.default_rng(args.rng)
    print(f"\n  Instanz der Bilder: n = {n}, alpha = {alpha}, seed = {seed}\n")

    # kappa fuer GENAU diese Instanz
    s0 = rng.uniform(-1, 1, (args.punkte, n))
    z, l, f, sch = C.laufe(s0, kl, n, tmax=1e9, eta=args.eta, pruef=1, maxschritt=6000)
    zs = np.sort(z[f])
    p = 1.0 - np.arange(1, len(zs) + 1) / len(f)
    m = (p < 0.5) & (p > 0.02) & (zs > 0)
    k, b = np.polyfit(zs[m], np.log(p[m]), 1)
    kappa = -k
    r = np.log(p[m]) - (k * zs[m] + b)
    r2 = 1 - r.var() / np.log(p[m]).var()
    print(f"  Fluchtrate      kappa  = {kappa:.4f}   (R^2 = {r2:.3f}, "
          f"{len(zs)} von {len(f)} konvergiert)")

    # lambda auf Startpunkten im Randgebiet (lange Transienten)
    schwelle = np.percentile(sch[f], args.perzentil)
    kand = s0[f & (sch >= schwelle)]
    if len(kand) > args.paare:
        kand = kand[rng.choice(len(kand), args.paare, replace=False)]
    lam, zt = paarlauf(kand.astype(np.float32), kl, n, eta=args.eta,
                       schritte=args.schritte, rng=rng)
    L = float(np.median(lam))
    print(f"  Streckrate      lambda = {L:.4f}   (Median ueber {len(lam)} Paare, "
          f"Quartile {np.percentile(lam,25):.3f} / {np.percentile(lam,75):.3f})")
    print(f"\n  Vorhersage      alpha  = kappa/lambda = {kappa/L:.4f}")
    print(f"  Gemessen (Bild) alpha  = 0,273           ->  D_Rand = {2 - kappa/L:.3f} "
          f"gegen gemessene 1,727")


if __name__ == "__main__":
    p = argparse.ArgumentParser()
    p.add_argument("--n", type=int, default=16)
    p.add_argument("--alpha", type=float, default=4.6)
    p.add_argument("--seed", type=int, default=3)
    p.add_argument("--punkte", type=int, default=20000)
    p.add_argument("--paare", type=int, default=500)
    p.add_argument("--schritte", type=int, default=3000)
    p.add_argument("--perzentil", type=float, default=90.0)
    p.add_argument("--eta", type=float, default=0.05)
    p.add_argument("--rng", type=int, default=99)
    lauf(p.parse_args())

  Instanz der Bilder: n = 16, alpha = 4.6, seed = 3

  Fluchtrate      kappa  = 0.1308   (R^2 = 0.937, 20000 von 20000 konvergiert)
  Streckrate      lambda = 0.6126   (Median ueber 500 Paare, Quartile 0.455 / 0.822)

  Vorhersage      alpha  = kappa/lambda = 0.2135
  Gemessen (Bild) alpha  = 0,273           ->  D_Rand = 1.786 gegen gemessene 1,727