#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ M13 Durchlass: Wie viel vom Formelzufall kommt im Loesungsraum an? Fuer jede loesbare Instanz (3-SAT, 3-XOR, 2-SAT) wird das exakte reduzierte BDD der GANZEN Loesungsmenge gebaut, minimal ueber mehrere Variablenordnungen. Vergleichsbasis: eine echt zufaellige Menge GLEICHER Groesse (gleiches Volumen), unter denselben Ordnungen. Das ist die maximal inkomprimierbare Menge. quote = BDD(Loesungsmenge) / BDD(Zufallsmenge gleicher Groesse) quote ~ 1 -> Loesungsraum so ungeordnet wie reiner Zufall ("Zufall im Quadrat") quote << 1 -> Loesungsraum hat eine kurze Regel (z.B. affin bei XOR) Dichten relativ zur eigenen Schwelle jedes Typs: 3-SAT 4.267, 3-XOR 0.918, 2-SAT 1.0 (faktor * schwelle) Seed: seed*1000003 + n*7919 + i*104729 + Typ-Offset (3sat 0, xor3 55511, sat2 33331) Aufruf (gross): nohup python3 m13_durchlass.py --jobs 16 > m13.log 2>&1 & (Standard: n=12..28, inst 60, Faktoren 0.3/0.5/0.7/0.9, Typen 3sat,rel6,rel5,rel4,xor3,sat2, DP-Block n=12..20 bei Faktor 1.5) Hinweis: nahe der Schwelle sind die Loesungsmengen winzig (wenige Bits), dort sagt die Quote wenig. Aussagekraeftig ist der Vergleich gleicher Volumenklasse (log2S) im Bericht. """ import argparse, csv, math, os, random, signal, statistics, sys, time from bisect import bisect_left from multiprocessing import Pool # relK: Constraint auf 3 Variablen mit K erlaubten von 8 Mustern (8-K zufaellig # verbotene). rel7 waere genau 3-SAT. Schwellen fuer relK sind SCHAETZUNGEN: # Erstes-Moment-Schranke ln2/ln(8/K), skaliert mit 4.267/5.191 (Faktor bei 3-SAT). def _schaetz(k): return round(4.267 / 5.191 * math.log(2) / math.log(8 / k), 3) SCHWELLE = {"3sat": 4.267, "xor3": 0.918, "sat2": 1.0, "rel6": _schaetz(6), "rel5": _schaetz(5), "rel4": _schaetz(4), "horn3": 4.267} # horn3: gleiche Klauselzahl wie 3sat (immer loesbar) OFFSET = {"3sat": 0, "xor3": 55511, "sat2": 33331, "rel6": 66661, "rel5": 44441, "rel4": 22221, "horn3": 88881} TYPEN = ["3sat", "rel6", "rel5", "rel4", "xor3", "sat2", "horn3"] # ---------------------------------------------------------------- BDD class BDD: """Reduziertes geordnetes BDD. Knoten 0/1 = Terminale, Ebene n.""" def __init__(self, n): self.n = n self.lev = [n, n]; self.lo = [0, 1]; self.hi = [0, 1] self.uniq = {} def mk(self, l, lo, hi): if lo == hi: return lo k = (l, lo, hi) u = self.uniq.get(k) if u is None: u = len(self.lev) self.lev.append(l); self.lo.append(lo); self.hi.append(hi) self.uniq[k] = u return u def und(self, a, b, memo): if a == 0 or b == 0: return 0 if a == 1: return b if b == 1 or a == b: return a if a > b: a, b = b, a k = (a, b) r = memo.get(k) if r is not None: return r la, lb = self.lev[a], self.lev[b] l = la if la < lb else lb a0, a1 = (self.lo[a], self.hi[a]) if la == l else (a, a) b0, b1 = (self.lo[b], self.hi[b]) if lb == l else (b, b) r = self.mk(l, self.und(a0, b0, memo), self.und(a1, b1, memo)) memo[k] = r return r def klausel(self, lits): """lits: Liste (ebene, wert_der_erfuellt). ODER der Literale.""" lits = sorted(lits, reverse=True) u = 0 for l, w in lits: u = self.mk(l, u, 1) if w else self.mk(l, 1, u) return u def xor(self, ebenen, par): """Paritaet der Variablen (Ebenen) == par.""" ebenen = sorted(ebenen, reverse=True) # z[p] = Knoten fuer "Restparitaet p noch zu erreichen" z = {0: 1, 1: 0} for l in ebenen: z = {p: self.mk(l, z[p], z[1 - p]) for p in (0, 1)} return z[par] def groesse(self, u): gesehen = set(); st = [u] while st: v = st.pop() if v < 2 or v in gesehen: continue gesehen.add(v); st.append(self.lo[v]); st.append(self.hi[v]) return len(gesehen) def zaehle(self, u): memo = {0: 0, 1: 1} def c(v): if v in memo: return memo[v] l = self.lev[v] lo, hi = self.lo[v], self.hi[v] r = c(lo) * (1 << (self.lev[lo] - l - 1)) + c(hi) * (1 << (self.lev[hi] - l - 1)) memo[v] = r return r return c(u) * (1 << self.lev[u]) if u > 1 else (u << self.n) def ausmenge(self, items, n): """ROBDD aus sortierter Liste von Ganzzahlen (Bit n-1-l = Ebene l).""" def b(l, lo_i, hi_i, basis): if lo_i == hi_i: return 0 if l == n: return 1 mitte = basis | (1 << (n - 1 - l)) m = bisect_left(items, mitte, lo_i, hi_i) return self.mk(l, b(l + 1, lo_i, m, basis), b(l + 1, m, hi_i, mitte)) return b(0, 0, len(items), 0) # ---------------------------------------------------------------- Instanzen def erzeuge(typ, n, faktor, seed, i): r = random.Random(seed * 1000003 + n * 7919 + i * 104729 + OFFSET[typ]) m = max(1, round(faktor * SCHWELLE[typ] * n)) k = 2 if typ == "sat2" else 3 cons = [] for _ in range(m): vs = r.sample(range(n), k) if typ == "xor3": cons.append((tuple(vs), r.randint(0, 1))) elif typ.startswith("rel"): kk = int(typ[3:]) verb = r.sample(range(8), 8 - kk) # verbotene Muster (Bit j = Wert vs[j]) cons.append((tuple(vs), tuple(sorted(verb)))) elif typ == "horn3": # Horn: hoechstens ein positives Literal; 1 von 8 Mustern verboten wie 3-SAT pos = r.randint(-1, 2) # -1 = alle negativ, sonst Index des positiven Literals cons.append(tuple((v, 1 if j == pos else 0) for j, v in enumerate(vs))) else: cons.append(tuple((v, r.randint(0, 1)) for v in vs)) return cons def cvars(typ, c): return c[0] if (typ == "xor3" or typ.startswith("rel")) else [v for v, _ in c] def cnf(typ, cons): """Klauseln als Liste von frozensets von Literalen (+v+1 / -(v+1)).""" out = [] for c in cons: if typ == "xor3" or typ.startswith("rel"): vs = c[0] if typ == "xor3": verb = [m for m in range(8) if bin(m).count("1") % 2 != c[1]] else: verb = c[1] for m in verb: # Klausel, die genau Muster m verbietet out.append(frozenset((vs[j] + 1) if not (m >> j) & 1 else -(vs[j] + 1) for j in range(3))) else: out.append(frozenset((v + 1) if w else -(v + 1) for v, w in c)) return out def ordnungen(typ, cons, n, anzahl, r): """Natuerlich, Cuthill-McKee-artig aus dem Interaktionsgraphen, Rest zufaellig.""" ords = [list(range(n))] nb = [set() for _ in range(n)] for c in cons: vs = cvars(typ, c) for a in vs: for b in vs: if a != b: nb[a].add(b) rest = set(range(n)); o = [] while rest: s = min(rest, key=lambda v: (len(nb[v]), v)) q = [s]; rest.discard(s) while q: v = q.pop(0); o.append(v) for w in sorted(nb[v] & rest, key=lambda x: len(nb[x])): rest.discard(w); q.append(w) ords.append(o) while len(ords) < anzahl: p = list(range(n)); r.shuffle(p); ords.append(p) return ords def bdd_loesung(typ, cons, n, ordn): lvl = [0] * n for l, v in enumerate(ordn): lvl[v] = l B = BDD(n) if typ == "xor3": teile = [B.xor([lvl[v] for v in vs], p) for vs, p in cons] schl = [min(lvl[v] for v in vs) for vs, _ in cons] elif typ.startswith("rel"): teile = [] for vs, verb in cons: u = 1 for m in verb: k = B.klausel([(lvl[vs[j]], 1 - ((m >> j) & 1)) for j in range(3)]) u = B.und(u, k, {}) teile.append(u) schl = [min(lvl[v] for v in vs) for vs, _ in cons] else: teile = [B.klausel([(lvl[v], w) for v, w in c]) for c in cons] schl = [min(lvl[v] for v, _ in c) for c in cons] # tief liegende Klauseln zuerst verunden (kleinere Zwischenstufen) teile = [t for _, t in sorted(zip(schl, teile), reverse=True)] u = 1 for t in teile: u = B.und(u, t, {}) if u == 0: break return B, u def basis(n, menge_int, ordn): """BDD-Groesse der Zufallsmenge unter Ordnung ordn.""" items = [] for x in menge_int: y = 0 for l, v in enumerate(ordn): if (x >> v) & 1: y |= 1 << (n - 1 - l) items.append(y) items.sort() B = BDD(n) return B.groesse(B.ausmenge(items, n)) def dp_weite(klauseln, n, kappe): """Davis-Putnam-Elimination (min pos*neg). Liefert (unsat, max Weite, log2 max Klauselzahl).""" S = set(klauseln); wmax = max(len(c) for c in S); mmax = len(S) rest = set(range(1, n + 1)) while rest: if frozenset() in S: return 1, wmax, round(math.log2(mmax), 3) best = None for v in rest: p = sum(1 for c in S if v in c); q = sum(1 for c in S if -v in c) if best is None or p * q < best[0]: best = (p * q, v) v = best[1]; rest.discard(v) P = [c for c in S if v in c]; N = [c for c in S if -v in c] S = {c for c in S if v not in c and -v not in c} for a in P: for b in N: r = (a - {v}) | (b - {-v}) if any(-x in r for x in r): continue S.add(frozenset(r)) if len(r) > wmax: wmax = len(r) if len(S) > mmax: mmax = len(S) if mmax > kappe: return -1, wmax, round(math.log2(mmax), 3) return (1 if frozenset() in S else 0), wmax, round(math.log2(max(mmax, 1)), 3) def aufgabe_dp(a): typ, n, faktor, seed, i, kappe, tmax = a signal.signal(signal.SIGALRM, _alarm); signal.alarm(tmax) t0 = time.time() z = dict(typ=typ, n=n, faktor=faktor, i=i, unsat="", weite="", log2max="", status="ok") try: cons = erzeuge(typ, n, faktor, seed + 500, i) u, w, lm = dp_weite(cnf(typ, cons), n, kappe) if u == -1: z["status"] = "kappe" z.update(unsat=max(u, 0), weite=w, log2max=lm) except Zeit: z["status"] = "zeit" finally: signal.alarm(0) z["sek"] = round(time.time() - t0, 3) return z class Zeit(Exception): pass def _alarm(*_): raise Zeit() def aufgabe(a): typ, n, faktor, seed, i, nord, basis_max, tmax = a signal.signal(signal.SIGALRM, _alarm); signal.alarm(tmax) t0 = time.time() zeile = dict(typ=typ, n=n, faktor=faktor, i=i, sat=0, log2S=0.0, bdd_min=0, bdd_nat=0, basis_min="", quote="", sek=0.0, status="ok") try: cons = erzeuge(typ, n, faktor, seed, i) r = random.Random(seed * 7 + n * 13 + i * 17 + OFFSET[typ]) ords = ordnungen(typ, cons, n, nord, r) groessen = []; S = None for o in ords: B, u = bdd_loesung(typ, cons, n, o) if u == 0: S = 0; break if S is None: S = B.zaehle(u) groessen.append(B.groesse(u)) if S: zeile.update(sat=1, log2S=round(math.log2(S), 4), bdd_min=min(groessen), bdd_nat=groessen[0]) if S <= basis_max: rb = random.Random(seed * 31 + n * 37 + i * 41 + OFFSET[typ] + 1) menge = rb.sample(range(1 << n), S) bm = min(basis(n, menge, o) for o in ords) zeile.update(basis_min=bm, quote=round(min(groessen) / bm, 4)) except Zeit: zeile["status"] = "zeit" except MemoryError: zeile["status"] = "speicher" finally: signal.alarm(0) zeile["sek"] = round(time.time() - t0, 3) return zeile # ---------------------------------------------------------------- Pruefungen def selbsttest(): """Exakte Kontrolle gegen Vollaufzaehlung auf kleinen n (BDD, CNF, DP).""" fehler = 0 # rel7 muss genau 3-SAT-Klauseln entsprechen: ein verbotenes Muster assert all(len(c) == 3 for c in cnf("3sat", erzeuge("3sat", 8, 1.0, 1, 0))) for typ in TYPEN: for n in (6, 8, 10): for i in range(15): cons = erzeuge(typ, n, 0.9 if i < 8 else 1.6, 99, i) kl = cnf(typ, cons) def ok(x): for c in kl: # gegen CNF pruefen: testet zugleich die CNF-Uebersetzung if not any((x >> (abs(l) - 1)) & 1 == (l > 0) for l in c): return False for c in cons: if typ.startswith("rel"): m = sum(((x >> c[0][j]) & 1) << j for j in range(3)) if m in c[1]: return False elif typ == "xor3": vs, p = c if sum((x >> v) & 1 for v in vs) % 2 != p: return False elif not any(((x >> v) & 1) == w for v, w in c): return False return True loes = [x for x in range(1 << n) if ok(x)] r = random.Random(i) for o in ordnungen(typ, cons, n, 3, r): B, u = bdd_loesung(typ, cons, n, o) S = B.zaehle(u) if u else 0 if S != len(loes): fehler += 1 if loes: # BDD der Loesungsmenge muss = BDD aus Menge sein if basis(n, loes, o) != B.groesse(u): fehler += 1 # DP: unsat genau dann, wenn keine Loesung u, _, _ = dp_weite(kl, n, 10**7) if u != (0 if loes else 1): fehler += 1 # XOR-Kontrolle: affiner Raum -> Groesse unabhaengig von Ordnung beschraenkt return fehler # ---------------------------------------------------------------- Bericht def med(xs): return statistics.median(xs) if xs else float("nan") def bericht(zeilen, dpz, args): out = [] p = out.append typen = args.typen p("M13 Durchlass — BDD der Loesungsmenge gegen Zufallsmenge gleicher Groesse") p(f"seed={args.seed} inst={args.inst} ordnungen={args.ordnungen} basis_max={args.basis_max}") p("Kopplungsachse: 3sat = 7 von 8 Mustern erlaubt, rel6/rel5/rel4 = 6/5/4 erlaubt (zufaellig),") p(" xor3 = 4 erlaubt, aber affin geordnet; sat2 = 2-SAT;") p(" horn3 = 7 erlaubt wie 3sat, gleiche Klauselzahl, aber geordnet (Horn).") p("Schwellen: " + ", ".join(f"{t} {SCHWELLE[t]}" for t in typen) + " (rel = Schaetzung)") st = {} for z in zeilen: st[z["status"]] = st.get(z["status"], 0) + 1 p(f"Status BDD: {st}") p("") p("Spalten: SAT-Anteil | Volumen log2S/n | BDD-min Median | Basis Median | Quote Median (S>=16) | IQR | N") for typ in typen: for f in args.faktoren: p(f"--- {typ} faktor {f} (Dichte {f*SCHWELLE[typ]:.3f})") for n in args.ns: zs = [z for z in zeilen if z["typ"] == typ and z["n"] == n and z["faktor"] == f and z["status"] == "ok"] if not zs: continue sat = [z for z in zs if z["sat"]] q = sorted(z["quote"] for z in sat if z["quote"] != "" and z["log2S"] >= 4) iqr = f"{q[len(q)//4]:.3f}-{q[(3*len(q))//4]:.3f}" if len(q) >= 4 else "-" bas = [z["basis_min"] for z in sat if z["basis_min"] != ""] p(f"n={n:>2} sat {len(sat)/len(zs):.2f} | vol {med([z['log2S']/n for z in sat]):.3f} | " f"bdd {med([z['bdd_min'] for z in sat]):.0f} | basis {med(bas):.0f} | " f"quote {med(q):.3f} | iqr {iqr} | N {len(q)}") p("") p("HAUPTBLOCK: Quote nach Volumenklasse log2S (alle n, alle Faktoren) — Median/N") klassen = [(4, 6), (6, 8), (8, 10), (10, 12), (12, 14), (14, 17), (17, 29)] for typ in typen: teile = [] for a, b in klassen: q = [z["quote"] for z in zeilen if z["typ"] == typ and z["status"] == "ok" and z["sat"] and z["quote"] != "" and a <= z["log2S"] < b] teile.append(f"[{a},{b}) {med(q):.3f}/{len(q)}" if q else f"[{a},{b}) -") p(f" {typ:>5}: " + " ".join(teile)) p("") p("Quote nach n innerhalb Volumenklasse [8,14) — waechst der Abstand mit n?") for typ in typen: teile = [] for n in args.ns: q = [z["quote"] for z in zeilen if z["typ"] == typ and z["n"] == n and z["status"] == "ok" and z["sat"] and z["quote"] != "" and 8 <= z["log2S"] < 14] if q: teile.append(f"n{n} {med(q):.3f}/{len(q)}") p(f" {typ:>5}: " + " ".join(teile)) p("") for lo, hi in ((0.3, 0.5), (0.5, 0.7)): p(f"HAUPTBLOCK 2: Quote bei festem Raumanteil log2S/n in [{lo},{hi}) nach n — Median/N") for typ in typen: teile = [] for n in args.ns: q = [z["quote"] for z in zeilen if z["typ"] == typ and z["n"] == n and z["status"] == "ok" and z["sat"] and z["quote"] != "" and lo <= z["log2S"] / n < hi] teile.append(f"n{n} {med(q):.3f}/{len(q)}" if len(q) >= 10 else f"n{n} -") p(f" {typ:>5}: " + " ".join(teile)) p("") p("Basis fehlt (S > basis_max): " + str(sum(1 for z in zeilen if z['sat'] and z['quote'] == ""))) if zeilen: p("Laufzeit BDD Median/Max s: " + f"{med([z['sek'] for z in zeilen]):.2f}/{max(z['sek'] for z in zeilen):.1f}") if dpz: p("") p(f"DP-BLOCK (Ebene 1, Beweisseite): Davis-Putnam auf unloesbaren Instanzen, Faktor {args.dp_faktor}") p("Weite = max. Klauselbreite (obere Schranke fuer Resolutionsweite); log2max = log2 groesste Klauselmenge") st = {} for z in dpz: st[z["status"]] = st.get(z["status"], 0) + 1 p(f"Status DP: {st}") for typ in typen: teile = [] for n in args.dp_ns: zs = [z for z in dpz if z["typ"] == typ and z["n"] == n] uz = [z for z in zs if z["status"] == "ok" and z["unsat"] == 1] kap = sum(1 for z in zs if z["status"] in ("kappe", "zeit")) if zs: teile.append(f"n{n} U{len(uz)}/{len(zs)} w{med([z['weite'] for z in uz]):.0f} " f"m{med([z['log2max'] for z in uz]):.1f} k{kap}") p(f" {typ:>5}: " + " ".join(teile)) p(" (U = unloesbar/gesamt, w = Weite Median, m = log2max Median, k = Kappe/Zeit)") return "\n".join(out) def main(): ap = argparse.ArgumentParser() ap.add_argument("--ns", default="12,14,16,18,20,22,24,26,28") ap.add_argument("--inst", type=int, default=60) ap.add_argument("--faktoren", default="0.3,0.5,0.7,0.9") ap.add_argument("--typen", default=",".join(TYPEN)) ap.add_argument("--seed", type=int, default=13) ap.add_argument("--ordnungen", type=int, default=6) ap.add_argument("--basis_max", type=int, default=400000) ap.add_argument("--dp_ns", default="12,14,16,18,20") ap.add_argument("--dp_faktor", type=float, default=1.5) ap.add_argument("--dp_kappe", type=int, default=300000) ap.add_argument("--jobs", type=int, default=16) ap.add_argument("--tmax", type=int, default=900) ap.add_argument("--csv", default="m13_durchlass.csv") ap.add_argument("--dpcsv", default="m13_dp.csv") ap.add_argument("--report", default="m13_report.txt") ap.add_argument("--nur_selbsttest", action="store_true") args = ap.parse_args() args.ns = [int(x) for x in args.ns.split(",")] args.dp_ns = [int(x) for x in args.dp_ns.split(",")] if args.dp_ns else [] args.faktoren = [float(x) for x in args.faktoren.split(",")] args.typen = args.typen.split(",") t = time.time(); f = selbsttest() print(f"Selbsttest (BDD, CNF, DP gegen Vollaufzaehlung, {len(TYPEN)} Typen x n=6/8/10): Fehler {f} [{time.time()-t:.1f}s]", flush=True) if f: print("ABBRUCH: Selbsttest fehlgeschlagen."); sys.exit(1) if args.nur_selbsttest: return aufg = [(typ, n, fa, args.seed, i, args.ordnungen, args.basis_max, args.tmax) for n in args.ns for typ in args.typen for fa in args.faktoren for i in range(args.inst)] aufg.sort(key=lambda a: -a[1]) dpa = [(typ, n, args.dp_faktor, args.seed, i, args.dp_kappe, args.tmax) for n in args.dp_ns for typ in args.typen for i in range(args.inst)] dpa.sort(key=lambda a: -a[1]) zeilen, dpz = [], []; t0 = time.time() with Pool(args.jobs) as pool: for name, fn, liste, ziel, pfad in (("BDD", aufgabe, aufg, zeilen, args.csv), ("DP", aufgabe_dp, dpa, dpz, args.dpcsv)): if not liste: continue with open(pfad, "w", newline="") as fh: w = None for k, z in enumerate(pool.imap_unordered(fn, liste), 1): if w is None: w = csv.DictWriter(fh, fieldnames=list(z.keys())); w.writeheader() w.writerow(z); fh.flush(); ziel.append(z) if k % 100 == 0 or k == len(liste): print(f"{name} [{k}/{len(liste)}] {(time.time()-t0)/60:.1f} Min", flush=True) txt = bericht(zeilen, dpz, args) with open(args.report, "w") as fh: fh.write(txt + "\n") print(txt) if __name__ == "__main__": main()