#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
M14 Dichte-Sweep: Wo ueberleben lokal konsistente Teilbelegungen am laengsten?

These: Haerte entsteht dort, wo lokal konsistente Zustaende am laengsten
ueberleben; das sollte knapp oberhalb der Schwelle sein und zu hoher Dichte
hin abfallen (Ueberbestimmung -> Widerspruch ueberall und redundant).

Zwei Messgroessen je Instanz:
  (1) DPLL-Knoten (Unit-Propagation, keine Lernklauseln) -- Suchaufwand,
      fuer loesbare und unloesbare Instanzen (klassisches Easy-Hard-Easy).
  (2) Kiesel-Zahl k*: kleinstes k, fuer das starke k-Konsistenz die leere
      Belegung entfernt (Spoiler gewinnt das existenzielle k-Kiesel-Spiel).
      Nach Atserias-Dalmau entspricht das bis auf +-1 der Resolutionsweite.
      Nur fuer unloesbare Instanzen definiert.

Dichte relativ zur eigenen Schwelle: faktor * SCHWELLE[typ].
Seed: seed*1000003 + n*7919 + i*104729 + Typ-Offset (wie M13) + 1000*faktor-Index

Abhaengigkeit: m13_durchlass.py (Generator, CNF) im selben Ordner.

Aufruf:
  nohup python3 m14_dichte.py --jobs 16 > m14.log 2>&1 &
"""
import argparse, csv, itertools, signal, statistics, sys, time
from multiprocessing import Pool

import m13_durchlass as m13

FAKTOREN = "0.5,0.7,0.85,1.0,1.15,1.3,1.6,2,3,4,6,8"


# ---------------------------------------------------------------- DPLL
def dpll(klauseln, n, knotenmax):
    """Liefert (sat, knoten). knoten = Zahl der Verzweigungen. -1 bei Kappe."""
    kl = [tuple(c) for c in klauseln]
    zaehl = [0]

    def prop(kl, bel):
        while True:
            neu = []
            einheit = None
            for c in kl:
                offen = []
                erf = False
                for l in c:
                    v = bel.get(abs(l))
                    if v is None:
                        offen.append(l)
                    elif v == (l > 0):
                        erf = True; break
                if erf: continue
                if not offen: return None
                if len(offen) == 1 and einheit is None: einheit = offen[0]
                neu.append(offen)
            if einheit is None: return neu
            bel[abs(einheit)] = einheit > 0
            kl = neu

    def rek(kl, bel):
        kl = prop(kl, bel)
        if kl is None: return False
        if not kl: return True
        zaehl[0] += 1
        if zaehl[0] > knotenmax: raise OverflowError
        # Verzweigung: haeufigste Variable in kuerzesten Klauseln
        kmin = min(len(c) for c in kl)
        h = {}
        for c in kl:
            if len(c) == kmin:
                for l in c: h[abs(l)] = h.get(abs(l), 0) + 1
        v = max(h, key=lambda x: (h[x], -x))
        for w in (True, False):
            b = dict(bel); b[v] = w
            if rek(kl, b): return True
        return False

    try:
        s = rek(kl, {})
        return int(s), zaehl[0]
    except OverflowError:
        return -1, -1


# ---------------------------------------------------------------- k-Konsistenz
def k_konsistenz(klauseln, n, k):
    """True, wenn starke k-Konsistenz die leere Belegung entfernt (Spoiler gewinnt)."""
    vs = list(range(n))
    erlaubt = {}
    # Klauseln als (vars-Tupel sortiert, verbotene Belegung als Dict)
    kls = []
    for c in klauseln:
        vv = tuple(sorted({abs(l) - 1 for l in c}))
        if len(vv) > k: continue  # breitere Klauseln wirken erst ueber Teilmengen nicht
        verb = {abs(l) - 1: (l < 0) for l in c}  # Belegung, die die Klausel verletzt
        kls.append((frozenset(vv), verb))
    for i in range(k + 1):
        for S in itertools.combinations(vs, i):
            Sset = set(S)
            maske = 0
            rel = [verb for vv, verb in kls if vv <= Sset]
            for a in range(1 << i):
                ok = True
                for verb in rel:
                    if all(((a >> S.index(v)) & 1) == int(w) for v, w in verb.items()):
                        ok = False; break
                if ok: maske |= 1 << a
            erlaubt[S] = maske
    # Projektionstabellen: proj[i][j][a] = a ohne Bit j
    proj = {}
    for i in range(1, k + 1):
        proj[i] = [[((a & ((1 << j) - 1)) | ((a >> (j + 1)) << j)) for a in range(1 << i)]
                   for j in range(i)]
    geaendert = True
    while geaendert:
        geaendert = False
        for i in range(k, 0, -1):
            for T in itertools.combinations(vs, i):
                mt = erlaubt[T]
                for j in range(i):
                    S = T[:j] + T[j + 1:]
                    ms = erlaubt[S]
                    pj = proj[i][j]
                    # runter: a in T nur, wenn Projektion in S erlaubt
                    neu_t = 0; stuetze = 0
                    a_bits = mt
                    while a_bits:
                        lsb = a_bits & -a_bits; a = lsb.bit_length() - 1; a_bits ^= lsb
                        b = pj[a]
                        if (ms >> b) & 1:
                            neu_t |= lsb; stuetze |= 1 << b
                    if neu_t != mt:
                        mt = neu_t; geaendert = True
                    # vor: b in S nur, wenn eine Erweiterung in T erlaubt ist
                    neu_s = ms & stuetze
                    if neu_s != ms:
                        erlaubt[S] = neu_s; geaendert = True
                erlaubt[T] = mt
        if erlaubt[()] == 0:
            return True
    return erlaubt[()] == 0


def kiesel_zahl(klauseln, n, kmin, kmax):
    for k in range(kmin, kmax + 1):
        if k_konsistenz(klauseln, n, k): return k
    return kmax + 1  # bedeutet "> kmax"


# ---------------------------------------------------------------- Aufgaben
class Zeit(Exception): pass
def _alarm(*_): raise Zeit()


def aufgabe(a):
    art, typ, n, fi, faktor, seed, i, par, tmax = a
    signal.signal(signal.SIGALRM, _alarm); signal.alarm(tmax)
    t0 = time.time()
    z = dict(art=art, typ=typ, n=n, faktor=faktor, i=i, sat="", knoten="", kiesel="", status="ok")
    try:
        cons = m13.erzeuge(typ, n, faktor, seed + 1000 * fi, i)
        kl = [frozenset(c) for c in m13.cnf(typ, cons)]
        s, kn = dpll(kl, n, par["knotenmax"])
        if s == -1:
            z["status"] = "kappe"
        else:
            z.update(sat=s, knoten=kn)
            if art == "kiesel" and s == 0:
                z["kiesel"] = kiesel_zahl(kl, n, 1, par["kmax"])
    except Zeit:
        z["status"] = "zeit"
    finally:
        signal.alarm(0)
    z["sek"] = round(time.time() - t0, 3)
    return z


# ---------------------------------------------------------------- Selbsttest
def selbsttest():
    fehler = 0
    for typ in ("3sat", "xor3", "rel5"):
        for n in (6, 8, 10):
            for i in range(12):
                f = [0.7, 1.5, 3.0][i % 3]
                cons = m13.erzeuge(typ, n, f, 77, i)
                kl = m13.cnf(typ, cons)
                def ok(x):
                    return all(any(((x >> (abs(l) - 1)) & 1) == (l > 0) for l in c) for c in kl)
                wahr = any(ok(x) for x in range(1 << n))
                s, _ = dpll(kl, n, 10**7)
                if s != int(wahr): fehler += 1
                # k = n + 1 muss Unloesbarkeit genau erkennen (volle Konsistenz bei kleinem n)
                if n <= 8:
                    if k_konsistenz(kl, n, n) != (not wahr): fehler += 1
                # Monotonie: wenn k widerlegt, dann auch k+1
                if not wahr and n <= 8:
                    ks = kiesel_zahl(kl, n, 1, n)
                    if ks <= n - 1 and not k_konsistenz(kl, n, ks + 1): fehler += 1
    # Kontrolle: 2-SAT-Widerspruch x->y->-x, -x->z->x ist mit k=3 widerlegbar, nicht mit k=1
    kl = [frozenset(c) for c in ([-1, 2], [-2, -1], [1, 3], [-3, 1])]
    if k_konsistenz(kl, 3, 1) or not k_konsistenz(kl, 3, 3): fehler += 1
    return fehler


# ---------------------------------------------------------------- Bericht
def med(xs):
    return statistics.median(xs) if xs else float("nan")


def bericht(zs, args):
    out = []; p = out.append
    p("M14 Dichte-Sweep — DPLL-Knoten (alle) und Kiesel-Zahl k* (nur unloesbar)")
    p(f"seed={args.seed} inst={args.inst} kmax={args.kmax} knotenmax={args.knotenmax}")
    st = {}
    for z in zs: st[(z['art'], z['status'])] = st.get((z['art'], z['status']), 0) + 1
    p(f"Status: {st}")
    for typ in args.typen:
        p("")
        p(f"=== {typ}  (Schwelle {m13.SCHWELLE[typ]})")
        p("DPLL: je Faktor  SAT-Anteil | Knoten-Median   (Spalten: n)")
        p("  faktor  alpha  | " + " | ".join(f"n={n:>3}" for n in args.ns_dpll))
        for f in args.faktoren:
            teile = []
            for n in args.ns_dpll:
                r = [z for z in zs if z["art"] == "dpll" and z["typ"] == typ and z["n"] == n
                     and z["faktor"] == f and z["status"] == "ok"]
                teile.append(f"{sum(z['sat'] for z in r)/len(r):.2f} {med([z['knoten'] for z in r]):>6.0f}"
                             if r else "   -       ")
            p(f"  {f:>5}  {f*m13.SCHWELLE[typ]:>6.2f} | " + " | ".join(teile))
        p("Kiesel-Zahl k* auf unloesbaren Instanzen:  Median (Anteil >kmax) [N]")
        p("  faktor  alpha  | " + " | ".join(f"n={n:>3}" for n in args.ns_kiesel))
        for f in args.faktoren:
            teile = []
            for n in args.ns_kiesel:
                r = [z for z in zs if z["art"] == "kiesel" and z["typ"] == typ and z["n"] == n
                     and z["faktor"] == f and z["status"] == "ok" and z["sat"] == 0]
                if r:
                    ks = [z["kiesel"] for z in r]
                    gross = sum(1 for k in ks if k > args.kmax) / len(ks)
                    teile.append(f"{med(ks):>4.1f} ({gross:.2f}) [{len(r):>2}]")
                else:
                    teile.append("      -        ")
            p(f"  {f:>5}  {f*m13.SCHWELLE[typ]:>6.2f} | " + " | ".join(teile))
    p("")
    p("Lesart: Die These sagt ein Maximum von k* (und der DPLL-Knoten) knapp oberhalb")
    p("faktor 1 voraus und einen Abfall zu hoher Dichte. k*=kmax+1 heisst nur '> kmax'.")
    p("Bei kleinem n ist 'hohe Dichte' nur hohe Dichte in alpha, nicht das Regime n^1.4.")
    return "\n".join(out)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--typen", default="3sat,xor3,rel5")
    ap.add_argument("--faktoren", default=FAKTOREN)
    ap.add_argument("--ns_dpll", default="20,30,40,50")
    ap.add_argument("--ns_kiesel", default="10,12,14,16")
    ap.add_argument("--inst", type=int, default=60)
    ap.add_argument("--seed", type=int, default=14)
    ap.add_argument("--kmax", type=int, default=5)
    ap.add_argument("--knotenmax", type=int, default=2_000_000)
    ap.add_argument("--jobs", type=int, default=16)
    ap.add_argument("--tmax", type=int, default=900)
    ap.add_argument("--csv", default="m14_dichte.csv")
    ap.add_argument("--report", default="m14_report.txt")
    ap.add_argument("--nur_selbsttest", action="store_true")
    args = ap.parse_args()
    args.typen = args.typen.split(",")
    args.faktoren = [float(x) for x in args.faktoren.split(",")]
    args.ns_dpll = [int(x) for x in args.ns_dpll.split(",")] if args.ns_dpll else []
    args.ns_kiesel = [int(x) for x in args.ns_kiesel.split(",")] if args.ns_kiesel else []

    t = time.time(); f = selbsttest()
    print(f"Selbsttest (DPLL und k-Konsistenz gegen Vollaufzaehlung): Fehler {f}  [{time.time()-t:.1f}s]", flush=True)
    if f:
        print("ABBRUCH: Selbsttest fehlgeschlagen."); sys.exit(1)
    if args.nur_selbsttest: return

    par = dict(kmax=args.kmax, knotenmax=args.knotenmax)
    aufg = []
    for art, ns in (("dpll", args.ns_dpll), ("kiesel", args.ns_kiesel)):
        for n in ns:
            for typ in args.typen:
                for fi, fa in enumerate(args.faktoren):
                    for i in range(args.inst):
                        aufg.append((art, typ, n, fi, fa, args.seed, i, par, args.tmax))
    aufg.sort(key=lambda a: (a[0] != "kiesel", -a[2]))  # teure zuerst
    zs = []; t0 = time.time()
    with Pool(args.jobs) as pool, open(args.csv, "w", newline="") as fh:
        w = None
        for k, z in enumerate(pool.imap_unordered(aufgabe, aufg), 1):
            if w is None:
                w = csv.DictWriter(fh, fieldnames=list(z.keys())); w.writeheader()
            w.writerow(z); fh.flush(); zs.append(z)
            if k % 200 == 0 or k == len(aufg):
                print(f"[{k}/{len(aufg)}] {(time.time()-t0)/60:.1f} Min", flush=True)
    txt = bericht(zs, args)
    with open(args.report, "w") as fh: fh.write(txt + "\n")
    print(txt)


if __name__ == "__main__":
    main()