#!/usr/bin/env python3
"""M7: Bewegung des Loesungsraums pro Formelaenderung  vs.  Expansion des Hypergraphen.

Gepaarte Instanzen: zu jeder 3-SAT-Instanz (bitgleich zu ebene1_dp.py, Typ sat3) wird eine 3-XOR-Instanz
auf EXAKT DENSELBEN Variablen-Tripeln gebaut (xor3p); nur die Relation unterscheidet sich (Paritaet statt ODER).
-> Expansion ist bei sat3 und xor3p per Konstruktion identisch. Jeder Unterschied in der Bewegung kommt aus der Relation.
Optional --xor-unabh: zusaetzlich unabhaengige XOR-Instanzen (xor3), bitgleich zu ebene1_dp.py.

Bewegung (exakt, Brute Force ueber alle 2^n Belegungen):
  d_j = Zahl der Loesungen, die beim Entfernen von Einheit j hinzukommen
      = Zahl der Belegungen, die NUR Einheit j verletzen  (ein einziger Durchlauf).
  Einheit = Klausel (sat3) bzw. ganze XOR-Gleichung (xor3p/xor3), nie CNF-Viertel.
  cv    = std(d)/mean(d)                       Ungleichfoermigkeit der Bewegung
  H     = Entropie der Werteverteilung von d    Beschreibungslaenge der Bewegung in Bit (XOR: wenige Werte -> klein)
  dist  = Anzahl verschiedener d-Werte;  frac0 = Anteil Einheiten ohne Bewegung
Expansion des Hypergraphen (Randvariablen = Variablen, die in genau einer Kante der Teilmenge S vorkommen):
  e2, e3   exaktes Minimum |Rand(S)|/|S| ueber zusammenhaengende S (Minimum wird immer von zusammenhaengenden S erreicht)
  eXX      lokale Suche fuer |S| = XX % von n  -> OBERE Schranke des Minimums (mehr --restarts = schaerfer)
  kern     kleinste gefundene Kantenmenge mit leerem Rand (gierig, obere Schranke)
  ACHTUNG: bei n <= ~30 enthaelt jede dichte Instanz ~3*alpha^2/n identische Tripel -> e2 = 0, kern = 2.
           Expansion ist im Brute-Force-Bereich degeneriert. Deshalb Option --dp:
--dp (Ebene-1-Kosten, braucht ebene1_dp.py im selben Ordner): DP-Variablenelimination (Reihenfolge minres)
  l2cnf = log2 max Klauselzahl, weite = breiteste Klausel  (Ben-Sasson/Wigderson-Kette: Expansion -> Weite -> Laenge)
Nutzung:  python3 m7_bewegung.py --ns 12,14,16,18,20 --inst 60 --alphas 0.918,4.267 --jobs 15
"""
import argparse, csv, math, sys, time
import numpy as np
from multiprocessing import Pool

# ---------------- Instanzen ----------------
def gen_sat3(n, alpha, seed):
    """bitgleich zu ebene1_dp.gen(n,'sat3',default_rng(seed))"""
    rng = np.random.default_rng(seed); m = int(round(alpha * n)); tri, cls = [], []
    for _ in range(m):
        vs = [int(v) for v in rng.choice(n, 3, replace=False)]
        ss = rng.choice([-1, 1], 3)
        tri.append(tuple(vs)); cls.append(tuple(int(s) * (v + 1) for v, s in zip(vs, ss)))
    return tri, cls

def gen_xor_gepaart(tri, seed):
    rng = np.random.default_rng([seed, 2])
    return [(t, int(rng.integers(2))) for t in tri]

def gen_xor_unabh(n, alpha, seed):
    """bitgleich zu ebene1_dp.gen(n,'xor3',...) mit dessen Seed (Offset 55511), als Gleichungen statt CNF"""
    rng = np.random.default_rng(seed + 55511); m = int(round(alpha * n)); eq = []
    for _ in range(m):
        vs = tuple(int(v) for v in rng.choice(n, 3, replace=False))
        eq.append((vs, int(rng.integers(2))))
    return eq

# ---------------- Bewegung ----------------
_BITS = {}
def bits(n):
    if n not in _BITS:
        _BITS.clear()
        idx = np.arange(1 << n, dtype=np.uint32)
        _BITS[n] = [((idx >> v) & 1).astype(bool) for v in range(n)]
    return _BITS[n]

def bewegung(n, verletzt):
    """verletzt: Liste boolescher Arrays (je Einheit: welche Belegungen sie verletzt)"""
    m = len(verletzt); N = 1 << n
    cnt = np.zeros(N, np.uint8 if m < 255 else np.uint16)
    last = np.full(N, -1, np.int32)
    for j, v in enumerate(verletzt):
        cnt += v; last[v] = j
    base = int((cnt == 0).sum())
    d = np.bincount(last[cnt == 1], minlength=m)
    return base, d

def verletzt_sat3(B, cls):
    out = []
    for c in cls:
        v = None
        for l in c:
            x = ~B[l - 1] if l > 0 else B[-l - 1]
            v = x if v is None else (v & x)
        out.append(v)
    return out

def verletzt_xor(B, eq):
    return [(B[a] ^ B[b] ^ B[c]) ^ bool(p) for (a, b, c), p in eq]

def kennzahlen(base, d):
    d = d.astype(float); mean = d.mean(); sd = d.std()
    vals, cts = np.unique(d, return_counts=True); p = cts / cts.sum()
    return dict(sol=base, l2sol=math.log2(base + 1), mean=mean,
                cv=sd / mean if mean > 0 else float("nan"),
                H=float(-(p * np.log2(p)).sum()), dist=len(vals), frac0=float((d == 0).mean()))

# ---------------- Expansion ----------------
def _nachbarn(tri, inc):
    return [sorted({k for v in e for k in inc[v]} - {j}) for j, e in enumerate(tri)]

def _rand(edges, tri):
    c = {}
    for j in edges:
        for v in tri[j]: c[v] = c.get(v, 0) + 1
    return sum(1 for x in c.values() if x == 1)

def exakt_klein(tri, nb):
    m = len(tri); b2 = b3 = float("inf")
    for j in range(m):
        for k in nb[j]:
            if k > j: b2 = min(b2, _rand((j, k), tri) / 2)
            for l in set(nb[j]) | set(nb[k]):
                if l != j and l != k and l > j:
                    b3 = min(b3, _rand((j, k, l), tri) / 3)
    return b2, b3

def lokal(tri, inc, n, s, restarts, rng, swaps_faktor=3):
    m = len(tri); best = float("inf")
    if s > m: return float("nan")
    for _ in range(restarts):
        cnt = [0] * n; inS = [False] * m; S = []; bd = 0
        def delta_add(j): return sum(1 if cnt[v] == 0 else (-1 if cnt[v] == 1 else 0) for v in tri[j])
        def delta_rem(j): return sum(-1 if cnt[v] == 1 else (1 if cnt[v] == 2 else 0) for v in tri[j])
        def kand():
            K = {k for j in S for v in tri[j] if cnt[v] == 1 for k in inc[v] if not inS[k]}
            return list(K) if K else [k for k in range(m) if not inS[k]]
        j0 = int(rng.integers(m)); bd += delta_add(j0)
        for v in tri[j0]: cnt[v] += 1
        inS[j0] = True; S.append(j0)
        while len(S) < s:
            K = kand(); ds = [delta_add(k) for k in K]; mn = min(ds)
            k = K[int(rng.choice([i for i, x in enumerate(ds) if x == mn]))]
            bd += mn
            for v in tri[k]: cnt[v] += 1
            inS[k] = True; S.append(k)
        for _ in range(swaps_faktor * s):          # Tausch: raus/rein, akzeptiere nicht-schlechter (Plateaus erlaubt)
            i = int(rng.integers(len(S))); out = S[i]
            dr = delta_rem(out)
            for v in tri[out]: cnt[v] -= 1
            inS[out] = False; S.pop(i)
            K = kand(); ds = [delta_add(k) for k in K]; mn = min(ds)
            if dr + mn <= 0:
                k = K[int(rng.choice([q for q, x in enumerate(ds) if x == mn]))]; bd += dr + mn
            else:
                k = out
            for v in tri[k]: cnt[v] += 1
            inS[k] = True; S.append(k)
        best = min(best, bd / s)
    return best

def kern(tri, inc, n, restarts, rng):
    m = len(tri); best = float("inf")
    for _ in range(restarts):
        cnt = [0] * n; inS = [False] * m; S = []; bd = 0
        j0 = int(rng.integers(m))
        while True:
            if not S: k = j0
            else:
                K = list({k for j in S for v in tri[j] if cnt[v] == 1 for k in inc[v] if not inS[k]})
                if not K: break
                ds = [sum(1 if cnt[v] == 0 else (-1 if cnt[v] == 1 else 0) for v in tri[q]) for q in K]
                mn = min(ds); k = K[int(rng.choice([i for i, x in enumerate(ds) if x == mn]))]
            for v in tri[k]:
                bd += 1 if cnt[v] == 0 else (-1 if cnt[v] == 1 else 0); cnt[v] += 1
            inS[k] = True; S.append(k)
            if bd == 0 or len(S) >= best: break
        if bd == 0: best = min(best, len(S))
    return best if best < float("inf") else float("nan")

def expansion(n, tri, anteile, restarts, seed):
    rng = np.random.default_rng([seed, 7])
    inc = [[] for _ in range(n)]
    for j, e in enumerate(tri):
        for v in e: inc[v].append(j)
    nb = _nachbarn(tri, inc)
    e2, e3 = exakt_klein(tri, nb)
    r = dict(e2=e2, e3=e3)
    for f in anteile:
        r[f"e{int(round(f * 100))}"] = lokal(tri, inc, n, max(2, int(round(f * n))), restarts, rng)
    r["kern"] = kern(tri, inc, n, restarts, rng)
    return r

# ---------------- Job ----------------
G = {}
G_DP = [False]
def _init(a): G.update(a)

def ebene1(n, cnf):
    import ebene1_dp
    mx, capped, _, w = ebene1_dp.dp(n, [frozenset(c) for c in cnf], "minres", G["cap"], G["tmax"])
    return dict(l2cnf=math.log2(max(mx, 1)), weite=w, gekappt=int(capped))

def xor_cnf(eq):
    out = []
    for vs, b in eq:
        for bits_ in range(8):
            xs = [(bits_ >> j) & 1 for j in range(3)]
            if (xs[0] ^ xs[1] ^ xs[2]) != b:
                out.append(tuple((-(v + 1) if x else (v + 1)) for v, x in zip(vs, xs)))
    return out

LEER = dict(l2cnf=float("nan"), weite=float("nan"), gekappt=0)

def job(key):
    n, alpha, i, seed = key
    B = bits(n); tri, cls = gen_sat3(n, alpha, seed)
    ex = expansion(n, tri, G["anteile"], G["restarts"], seed)
    rows = []
    base, d = bewegung(n, verletzt_sat3(B, cls))
    rows.append(dict(n=n, alpha=alpha, i=i, typ="sat3", **kennzahlen(base, d), **ex,
                     **(ebene1(n, cls) if G["dp"] else LEER)))
    eq = gen_xor_gepaart(tri, seed)
    base, d = bewegung(n, verletzt_xor(B, eq))
    rows.append(dict(n=n, alpha=alpha, i=i, typ="xor3p", **kennzahlen(base, d), **ex,
                     **(ebene1(n, xor_cnf(eq)) if G["dp"] else LEER)))
    if G["xor_unabh"]:
        eq = gen_xor_unabh(n, alpha, seed)
        exu = expansion(n, [t for t, _ in eq], G["anteile"], G["restarts"], seed + 55511)
        base, d = bewegung(n, verletzt_xor(B, eq))
        rows.append(dict(n=n, alpha=alpha, i=i, typ="xor3", **kennzahlen(base, d), **exu,
                         **(ebene1(n, xor_cnf(eq)) if G["dp"] else LEER)))
    return rows

# ---------------- Statistik ----------------
def _raenge(x):
    x = np.asarray(x, float); o = np.argsort(x, kind="mergesort"); r = np.empty(len(x)); r[o] = np.arange(len(x))
    for v in np.unique(x):                      # Bindungen mitteln
        w = x == v
        if w.sum() > 1: r[w] = r[w].mean()
    return r

def _pearson(x, y):
    x, y = np.asarray(x, float), np.asarray(y, float)
    if len(x) < 4 or x.std() == 0 or y.std() == 0: return float("nan")
    return float(np.corrcoef(x, y)[0, 1])

def spearman(x, y): return _pearson(_raenge(x), _raenge(y))

def partiell(x, y, z):
    rxy, rxz, ryz = spearman(x, y), spearman(x, z), spearman(y, z)
    den = math.sqrt(max(1e-12, (1 - rxz ** 2) * (1 - ryz ** 2)))
    return (rxy - rxz * ryz) / den

def med(v):
    v = [x for x in v if not (isinstance(x, float) and math.isnan(x))]
    return float(np.median(v)) if v else float("nan")

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--ns", default="12,14,16,18,20")
    ap.add_argument("--inst", type=int, default=60)
    ap.add_argument("--alphas", default="0.918,4.267")
    ap.add_argument("--jobs", type=int, default=15)
    ap.add_argument("--seed", type=int, default=11)
    ap.add_argument("--anteile", default="0.25,0.5,1.0", help="|S| fuer lokale Expansionssuche als Anteil von n")
    ap.add_argument("--restarts", type=int, default=60)
    ap.add_argument("--xor-unabh", action="store_true")
    ap.add_argument("--dp", action="store_true", help="Ebene-1-Kosten (Weite, log2 cnf) mitmessen")
    ap.add_argument("--cap", type=int, default=400000)
    ap.add_argument("--tmax", type=float, default=30.0)
    ap.add_argument("--csv", default="m7_bewegung.csv")
    ap.add_argument("--report", default="report_m7.txt")
    a = ap.parse_args()
    ns = [int(x) for x in a.ns.split(",")]; als = [float(x) for x in a.alphas.split(",")]
    anteile = [float(x) for x in a.anteile.split(",")]
    keys = [(n, al, i, a.seed * 1000003 + n * 7919 + i * 104729) for al in als for n in ns for i in range(a.inst)]
    keys.sort(key=lambda k: -k[0])                 # grosse n zuerst -> bessere Auslastung
    R, t0 = [], time.time()
    with Pool(a.jobs, initializer=_init, initargs=(dict(anteile=anteile, restarts=a.restarts, xor_unabh=a.xor_unabh, dp=a.dp, cap=a.cap, tmax=a.tmax),)) as p:
        for k, rows in enumerate(p.imap_unordered(job, keys), 1):
            R.extend(rows)
            if k % max(1, len(keys) // 40) == 0 or k == len(keys):
                el = time.time() - t0
                print(f"[{k}/{len(keys)}] {el/60:6.1f} min, Rest ~{el/k*(len(keys)-k)/60:6.1f} min", file=sys.stderr, flush=True)
    with open(a.csv, "w", newline="") as f:
        w = csv.DictWriter(f, fieldnames=list(R[0].keys())); w.writeheader(); w.writerows(R)

    G_DP[0] = a.dp
    ecols = ["e2", "e3"] + [f"e{int(round(f*100))}" for f in anteile] + ["kern"]
    emain = f"e{int(round(sorted(anteile)[len(anteile)//2]*100))}"
    typen = ["sat3", "xor3p"] + (["xor3"] if a.xor_unabh else [])
    L = [f"M7  seed={a.seed}  inst={a.inst}  restarts={a.restarts}  Haupt-Expansion={emain}", ""]
    for al in als:
        L.append(f"===== alpha = {al} =====")
        L.append("typ    n | SAT  | l2sol | cv(med) H(med) dist(med/max) frac0 | l2cnf weite gek | " + " ".join(f"{c:>6s}" for c in ecols))
        for t in typen:
            for n in ns:
                rr = [r for r in R if r["alpha"] == al and r["n"] == n and r["typ"] == t]
                m = int(round(al * n))
                L.append(f"{t:5s} {n:3d} | {np.mean([r['sol'] > 0 for r in rr]):.2f} | {med([r['l2sol'] for r in rr]):5.2f} | "
                         f"{med([r['cv'] for r in rr]):7.3f} {med([r['H'] for r in rr]):6.3f} {med([r['dist'] for r in rr]):5.1f}/{max(r['dist'] for r in rr):3d} "
                         f"{med([r['frac0'] for r in rr]):5.2f} | {med([r['l2cnf'] for r in rr]):5.2f} {med([r['weite'] for r in rr]):5.1f} "
                         f"{sum(r['gekappt'] for r in rr):3d} | " + " ".join(f"{med([r[c] for r in rr]):6.3f}" for c in ecols))
            L.append("")
        L.append("Gepaart (gleicher Hypergraph): Anteil Seeds mit H(sat3) > H(xor3p) | Median H(sat3)-H(xor3p)")
        for n in ns:
            S = {r["i"]: r for r in R if r["alpha"] == al and r["n"] == n and r["typ"] == "sat3"}
            X = {r["i"]: r for r in R if r["alpha"] == al and r["n"] == n and r["typ"] == "xor3p"}
            dH = [S[i]["H"] - X[i]["H"] for i in S if i in X]
            L.append(f"  n={n:3d}:  {np.mean([x > 0 for x in dH]):.2f}  |  {med(dH):+.3f}")
        L.append("")
        L.append("Gleiche Loesungszahl (l2sol gerundet): dist sat3 vs xor3p  (Median [Anzahl])")
        for b in sorted({round(r["l2sol"]) for r in R if r["alpha"] == al}):
            q = {t: [r["dist"] for r in R if r["alpha"] == al and r["typ"] == t and round(r["l2sol"]) == b] for t in ("sat3", "xor3p")}
            if q["sat3"] and q["xor3p"]:
                L.append(f"  l2sol~{b:2d}:  sat3 {med(q['sat3']):5.1f} [{len(q['sat3'])}]   xor3p {med(q['xor3p']):5.1f} [{len(q['xor3p'])}]")
        L.append("")
        if G_DP[0]:
            L.append("Bewegung vs Ebene-1-Kosten innerhalb sat3 (Spearman je n; nicht gekappte Instanzen):")
            L.append("   n | rho(dist,l2cnf) rho(H,l2cnf) rho(dist,weite) | rho(l2cnf,l2sol) | partiell rho(dist,l2cnf | l2sol)  [Anzahl]")
            for n in ns:
                rr = [r for r in R if r["alpha"] == al and r["n"] == n and r["typ"] == "sat3" and not r["gekappt"]]
                if len(rr) < 6:
                    L.append(f" {n:3d} | zu wenige ({len(rr)})"); continue
                D = [r["dist"] for r in rr]; H = [r["H"] for r in rr]; C = [r["l2cnf"] for r in rr]
                W = [r["weite"] for r in rr]; S = [r["l2sol"] for r in rr]
                L.append(f" {n:3d} | {spearman(D, C):+.3f}  {spearman(H, C):+.3f}  {spearman(D, W):+.3f} | {spearman(C, S):+.3f} | "
                         f"{partiell(D, C, S):+.3f}  [{len(rr)}]")
            L.append("")
        L.append(f"Korrelation innerhalb sat3 (Spearman, je n; nur Instanzen mit definiertem cv):")
        L.append(f"   n | rho(H,{emain}) rho(cv,{emain}) | rho(H,l2sol) rho(cv,l2sol) | partiell rho(H,{emain} | l2sol) | rho(H,kern)  [Anzahl]")
        for n in ns:
            rr = [r for r in R if r["alpha"] == al and r["n"] == n and r["typ"] == "sat3"
                  and not math.isnan(r["cv"]) and not math.isnan(r[emain])]
            if len(rr) < 6:
                L.append(f" {n:3d} | zu wenige Instanzen ({len(rr)})"); continue
            H = [r["H"] for r in rr]; cv = [r["cv"] for r in rr]; e = [r[emain] for r in rr]; ls = [r["l2sol"] for r in rr]
            kr = [r["kern"] if not math.isnan(r["kern"]) else 1e9 for r in rr]
            L.append(f" {n:3d} | {spearman(H, e):+.3f}  {spearman(cv, e):+.3f} | {spearman(H, ls):+.3f}  {spearman(cv, ls):+.3f} | "
                     f"{partiell(H, e, ls):+.3f} | {spearman(H, kr):+.3f}  [{len(rr)}]")
        L.append("")
    L += ["Lesart:",
          " - sat3 und xor3p haben identische Expansion (gleiche Tripel). Unterschied in H/cv/dist = allein die Relation.",
          "   dist = Zahl verschiedener Bewegungswerte ist das Hauptmass (nicht cv): xor3p theoretisch hoechstens ~2",
          "   (Gleichung streichen verdoppelt den affinen Raum oder aendert nichts), sat3 viele und wachsend.",
          "   Fairer Vergleich nur bei gleicher Loesungszahl (Abschnitt 'Gleiche Loesungszahl').",
          " - Expansion kann sat3 und xor3p per Konstruktion NICHT trennen. Passt dazu, dass Resolution auf beiden",
          "   exponentiell ist (Ben-Sasson/Wigderson gilt fuer beide): Expansion ist eine Ebene-1-Groesse.",
          " - Expansion (eXX, kern) ist bei kleinem n durch doppelte Tripel degeneriert (~3*alpha^2/n Stueck je Instanz).",
          "   Die Bruecke zur Literatur laeuft deshalb ueber --dp: korreliert dist mit Weite/Laenge, bei herausgerechneter",
          "   Loesungszahl (partiell)? Positiv = die Bewegung misst dieselbe Haerte, die Resolution bezahlt.",
          " - Korrelation innerhalb sat3: H und cv haengen vermutlich an der Loesungszahl. Deshalb die partielle Korrelation",
          "   (H gegen Expansion bei herausgerechneter l2sol). Nur sie sagt, ob Bewegung und Expansion DASSELBE messen.",
          " - eXX und kern sind obere Schranken aus lokaler Suche; Rauschen dort schwaecht Korrelationen ab (konservativ)."]
    txt = "\n".join(L); print(txt); open(a.report, "w").write(txt + "\n")

if __name__ == "__main__":
    main()