#!/usr/bin/env python3
"""Ebene 1: reines Ausschlussverfahren (Davis-Putnam-Variablenelimination) gegen gewonnene Information.

Typen:  sat3 (3-OR, alpha=4.267)   xor3 (3-XOR als je 4 Klauseln, alpha=0.918)   -> gleicher Schwellenabstand
Reihenfolgen:  idx (fest nach Index, rein idiomatisch)   minres (wenigste Resolventen zuerst)
DP: nur Resolution + Streichen von Tautologien und Duplikaten. Keine Subsumption, keine Heuristik sonst.

Je Instanz (exakt, Brute Force):
  need  = n - log2(#Loesungen)      Bits, die die Formel tatsaechlich ausschliesst (UNSAT: n)
  indep = Summe der Einzelklausel-Bits (jede k-Klausel: log2(2^k/(2^k-1)))
  syn   = need - indep              Information, die nur aus dem Zusammenspiel entsteht
  cnf   = max Klauselzahl waehrend DP  (Ebene-1-Aufwand; gekappt bei --cap)
Fuer SAT-Instanzen zusaetzlich, je Eliminationsschritt fuer die Projektion  EXISTS(eliminierte) F :
  bdd   = max ROBDD-Groesse der Projektion (generische komprimierte Darstellung)
  aff   = Anteil der Schritte, in denen die Projektion ein affiner Raum ist (Gauss-Sprache)
Hinweis: bdd/aff messen, ob das Zwischenwissen PRINZIPIELL komprimierbar ist (semantisch),
nicht, ob ein Verfahren es billig erreicht. Sie werden aus der Brute-Force-Loesungsmenge gelesen.
Nutzung:  python3 ebene1_dp.py --ns 10,12,14,16,18 --inst 40 --jobs 20
"""
import argparse, math, numpy as np
from multiprocessing import Pool

ALPHA = {"sat3": 4.267, "xor3": 0.918}

def _set_alpha(a3, ax):
    # Worker-Initialisierung: noetig, weil ab Python 3.14 'forkserver' Standard ist und
    # Worker die Modul-Globals sonst mit den Standardwerten neu laden.
    ALPHA["sat3"], ALPHA["xor3"] = a3, ax

def gen(n, typ, rng):
    m = int(round(ALPHA[typ] * n)); cls = []
    for _ in range(m):
        vs = [int(v) for v in rng.choice(n, 3, replace=False)]
        if typ == "sat3":
            ss = rng.choice([-1, 1], 3)
            cls.append(frozenset(int(s) * (v + 1) for v, s in zip(vs, ss)))
        else:
            b = int(rng.integers(2))                     # x1^x2^x3 = b
            for bits in range(8):
                xs = [(bits >> j) & 1 for j in range(3)]
                if (xs[0] ^ xs[1] ^ xs[2]) != b:         # verbotene Belegung -> Klausel, die sie ausschliesst
                    cls.append(frozenset((-(v + 1) if x else (v + 1)) for v, x in zip(vs, xs)))
    return cls

def solutions(n, cls):
    idx = np.arange(1 << n, dtype=np.uint32)
    ok = np.ones(1 << n, bool)
    for c in cls:
        sat = np.zeros(1 << n, bool)
        for l in c:
            bit = ((idx >> (abs(l) - 1)) & 1).astype(bool)
            sat |= bit if l > 0 else ~bit
        ok &= sat
    return ok

def dp(n, cls, order, cap, tmax):
    import time
    S = set(cls); mx = len(S); elim = []; remaining = set(range(1, n + 1))
    capped = False; t0 = time.time()
    wmax = max((len(c) for c in S), default=0)   # Resolutionsweite: breiteste Klausel im Lauf
    while remaining:
        if tmax and time.time() - t0 > tmax: capped = True; break
        if order == "idx":
            v = min(remaining)
        else:
            def cost(u):
                p = sum(1 for c in S if u in c); q = sum(1 for c in S if -u in c)
                return p * q - p - q
            v = min(remaining, key=cost)
        P = [c for c in S if v in c]; N = [c for c in S if -v in c]
        R = set(); cnt = 0
        for a in P:
            for b in N:
                cnt += 1
                if tmax and (cnt & 1023) == 0 and time.time() - t0 > tmax:
                    capped = True; break                 # Zeit-Deckel auch mitten in einer Elimination
                r = (a - {v}) | (b - {-v})
                if any(-l in r for l in r): continue     # Tautologie
                R.add(frozenset(r))
                if len(r) > wmax: wmax = len(r)
                if len(R) > cap: capped = True; break
            if capped: break
            if tmax and time.time() - t0 > tmax: capped = True; break
        S = (S - set(P) - set(N)) | R
        mx = max(mx, len(S)); remaining.discard(v); elim.append(v - 1)
        if tmax and time.time() - t0 > tmax: capped = True; break   # Zeit-Deckel pro Instanz
        if capped or frozenset() in S: break
    return mx, capped, elim, wmax

def robdd_size(tt):
    m = int(round(math.log2(len(tt)))); tot = 0
    for i in range(m):                                   # oberste Variable = hoechstes Bit
        ch = tt.reshape(1 << i, 1 << (m - i))
        h = 1 << (m - i - 1)
        diff = (ch[:, :h] != ch[:, h:]).any(axis=1)
        if diff.any():
            tot += len(np.unique(np.packbits(ch[diff], axis=1), axis=0))
    return tot

def is_affine(tt, rng, trials=1500):
    pts = np.flatnonzero(tt)
    if len(pts) <= 2: return True
    if len(pts) & (len(pts) - 1): return False           # Groesse keine Zweierpotenz
    a, b, c = (pts[rng.integers(len(pts), size=trials)] for _ in range(3))
    return bool(tt[a ^ b ^ c].all())

def job(args):
    n, typ, order, i, seed, cap, tmax = args
    rng = np.random.default_rng(seed)
    cls = gen(n, typ, rng)
    ok = solutions(n, cls); C = int(ok.sum())
    need = n if C == 0 else n - math.log2(C)
    indep = sum(math.log2((1 << len(c)) / ((1 << len(c)) - 1)) for c in cls)
    mx, capped, elim, wmax = dp(n, cls, order, cap, tmax)
    bdd, aff = -1, float("nan")
    if C > 0:
        cube = ok.reshape([2] * n)                       # Achse 0 = Variable n-1 (hoechstes Bit)
        sizes, affs = [], []
        for k in range(1, n):
            axes = tuple(n - 1 - v for v in elim[:k])
            proj = cube.any(axis=axes).reshape(-1)
            sizes.append(robdd_size(proj)); affs.append(is_affine(proj, rng))
        bdd = max(sizes) if sizes else 0; aff = float(np.mean(affs)) if affs else float("nan")
    return dict(n=n, typ=typ, order=order, i=i, sat=C > 0, need=need, indep=indep,
                cnf=mx, capped=capped, width=wmax, bdd=bdd, aff=aff)

def slope(x, y):
    x, y = np.asarray(x, float), np.asarray(y, float)
    return np.polyfit(x, y, 1)[0] if len(set(x)) > 1 else float("nan")

def _fortschritt(p, job, keys, cs=1):
    import sys, time
    R, t0, tot = [], time.time(), len(keys)
    for k, r in enumerate(p.imap_unordered(job, keys, chunksize=cs), 1):
        R.append(r)
        if k % max(1, tot // 40) == 0 or k == tot:
            el = time.time() - t0
            print(f"[{k}/{tot}] {el/60:6.1f} min, Rest ~{el/k*(tot-k)/60:6.1f} min", file=sys.stderr, flush=True)
    return R

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--ns", default="8,10,12,14,16")
    ap.add_argument("--inst", type=int, default=40)
    ap.add_argument("--jobs", type=int, default=20)
    ap.add_argument("--cap", type=int, default=400000)
    ap.add_argument("--tmax", type=float, default=30.0,
                    help="Zeit-Deckel pro Instanz in Sekunden (0 = aus); gekappte zaehlen voll bei need/indep/syn/bdd/aff, cnf wird zur unteren Schranke")
    ap.add_argument("--seed", type=int, default=11)
    ap.add_argument("--alpha-sat3", type=float, default=4.267)
    ap.add_argument("--alpha-xor3", type=float, default=0.918,
                    help="gleiches Geruest: --alpha-sat3 0.918 oder --alpha-xor3 4.267")
    ap.add_argument("--report", default="report_ebene1.txt")
    a = ap.parse_args()
    ALPHA["sat3"], ALPHA["xor3"] = a.alpha_sat3, a.alpha_xor3
    ns = [int(x) for x in a.ns.split(",")]
    keys = [(n, t, o, i, a.seed * 1000003 + n * 7919 + i * 104729 + (0 if t == "sat3" else 55511), a.cap)
            for n in ns for t in ALPHA for o in ("idx", "minres") for i in range(a.inst)]
    keys = [k + (a.tmax,) for k in keys]
    print(f"alpha: sat3={ALPHA['sat3']}  xor3={ALPHA['xor3']}", flush=True)
    with Pool(a.jobs, initializer=_set_alpha, initargs=(a.alpha_sat3, a.alpha_xor3)) as p:
        R = _fortschritt(p, job, keys, 2)
    L = [f"alpha: sat3={a.alpha_sat3}  xor3={a.alpha_xor3}", "typ  order   n | SAT | need  indep  syn  (Bits, Median) | log2 cnf [ungekappt] (Anzahl gekappt) | Weite [ungek.] | log2 bdd SAT | affin"]
    for t in ALPHA:
        for o in ("idx", "minres"):
            xs, ycnf, ybdd, yneed, ywd = [], [], [], [], []
            for n in ns:
                rr = [r for r in R if r["n"] == n and r["typ"] == t and r["order"] == o]
                sat = [r for r in rr if r["sat"]]
                need = np.median([r["need"] for r in rr]); ind = np.median([r["indep"] for r in rr])
                syn = np.median([r["need"] - r["indep"] for r in rr])
                lc = np.median([math.log2(r["cnf"]) for r in rr]); cp = sum(r["capped"] for r in rr)
                lc_un = np.median([math.log2(r["cnf"]) for r in rr if not r["capped"]]) if cp < len(rr) else float("nan")
                wd = np.median([r["width"] for r in rr])
                wd_un = np.median([r["width"] for r in rr if not r["capped"]]) if cp < len(rr) else float("nan")
                lb = np.median([math.log2(max(r["bdd"], 1)) for r in sat]) if sat else float("nan")
                af = np.nanmean([r["aff"] for r in sat]) if sat else float("nan")
                L.append(f"{t} {o:6s} {n:3d} | {len(sat)/len(rr):.2f} | {need:5.2f} {ind:5.2f} {syn:+5.2f} | "
                         f"{lc:6.2f} [ungek. {lc_un:5.2f}] ({cp}/{len(rr)} gekappt) | {wd:5.1f} [{wd_un:5.1f}] | {lb:6.2f} | {af:.2f}")
                xs.append(n); ywd.append(wd_un if not math.isnan(wd_un) else wd); ycnf.append(lc_un if not math.isnan(lc_un) else lc); ybdd.append(lb); yneed.append(need)
            L.append(f"   Steigung je Variable:  need {slope(xs, yneed):.3f}   log2 cnf {slope(xs, ycnf):.3f}   "
                     f"Weite {slope(xs, ywd):.3f}   log2 bdd {slope(xs, ybdd):.3f}")
            L.append("")
    L.append("Hinweis: gekappte Instanzen (Zeit-Deckel) zaehlen voll bei need/indep/syn/bdd/affin; ihr cnf ist nur untere Schranke.")
    L.append("  cnf-Median [ungek.] und die cnf-Steigung nutzen nur NICHT-gekappte Instanzen; viele Kappungen = Verfahren explodiert hier (selbst ein Befund).")
    L.append("Weite = breiteste Klausel, die DP in dieser Reihenfolge erzeugt (obere Schranke fuer die minimale Resolutionsweite;")
    L.append("  Ben-Sasson/Wigderson: grosse Mindestweite erzwingt exponentielle Laenge). Gekappt: nur bis zum Abbruch gesehen.")
    L.append("Lesart: cnf-Steigung >> need-Steigung -> Ausschlussverfahren zahlt exponentiell fuer linear viel Information.")
    L.append("xor3: Projektionen affin (aff=1) -> das explodierende Wissen ist in Gauss-Sprache kompakt (Hyperspace).")
    L.append("bdd: waechst es bei sat3 UND xor3, ist BDD nicht der Hyperspace, der Gauss erklaert.")
    txt = "\n".join(L); print(txt); open(a.report, "w").write(txt + "\n")

if __name__ == "__main__":
    main()