Bahnen — fünf Tafeln

"""
Die BAHN statt der Momentaufnahme.

MATHIAS' KORREKTUR. Uebersetzt man Klausel fuer Klausel in Geometrie, so
kann die Loesungsmenge erst leer sein, wenn die letzte relevante Klausel
die Erfuellbarkeit kippt. Bis dahin ist sie per Definition nichtleer. Der
Gegenstand ist also die VERLAUFSKURVE, nicht der Endpunkt — und
restgeometrie.py hat nur eine Momentaufnahme gemessen.

Das macht den dortigen Befund erst scharf: einen Schritt vor dem Kippen
sind nur noch 5,3 Variablen frei. Der Zusammenbruch ist also nicht das Werk
der letzten Klausel. WO passiert er?

VERFAHREN. Klauseln in der Erzeugungsreihenfolge hinzufuegen, den Kipppunkt
k* per Bisektion finden, und dann die letzten Schritte davor vermessen:
   spann     Zahl der Variablen, die ueber 20 Loesungen variieren
   hamming   mittlerer paarweiser Abstand
   loesungen wieviele der 20 verlangten ueberhaupt existieren

Loesungen exakt ueber Blockierklauseln, kein Absuchen.

VORHERSAGE, VOR DEM LAUF NOTIERT. Die Kurve wird von der Dichte bestimmt
und daher zwischen Instanzen aehnlich sein; instanzspezifische Information
steckt allenfalls in der ABWEICHUNG von der mittleren Kurve. Ob die die
Haerte vorhersagt, ist die eigentliche Frage.
"""
import os, sys, pickle, resource
import numpy as np
from restgeometrie import loesungen_sammeln, rang


def kipppunkt(cl, n, budget=8000000):
    """Kleinstes k, sodass die ersten k Klauseln unerfuellbar sind."""
    from pysat.solvers import Cadical153
    lo, hi = 1, len(cl)
    while lo < hi:
        mid = (lo + hi) // 2
        so = Cadical153(bootstrap_with=cl[:mid]); so.conf_budget(budget)
        r = so.solve_limited(); so.delete()
        if r is False:
            hi = mid
        else:
            lo = mid + 1
    return lo


def geometrie(cl, n, k=20):
    L = loesungen_sammeln(cl, n, k)
    if len(L) < 2:
        return {"spann": 0, "hamming": 0.0, "loesungen": len(L)}
    A = np.array(L)
    d = [np.sum(A[a] != A[b]) for a in range(len(A)) for b in range(a + 1, len(A))]
    return {"spann": int((A.max(0) != A.min(0)).sum()),
            "hamming": float(np.mean(d)), "loesungen": len(L)}


if __name__ == "__main__":
    resource.setrlimit(resource.RLIMIT_AS, (3500 * 1024 * 1024,) * 2)
    resource.setrlimit(resource.RLIMIT_CPU, (3300, 3300))
    from pysat.solvers import Cadical153
    from math import sqrt, erfc
    SCR = os.environ["SCR"]
    with open(os.path.join(SCR, "instanzen.pkl"), "rb") as f:
        CLS = pickle.load(f)
    N, ZIEL = 100, 30
    RUECK = [40, 30, 20, 15, 10, 6, 3, 2, 1]
    print("Kipppunkt suchen und die Bahn davor vermessen, %d Instanzen, n = %d" % (ZIEL, N))
    K, BAHN, KST = [], [], []
    for i, (cl, R) in enumerate(CLS[:ZIEL]):
        so = Cadical153(bootstrap_with=cl); so.conf_budget(20000000)
        so.solve_limited(); K.append(so.accum_stats().get("conflicts", 0)); so.delete()
        ks = kipppunkt(cl, N)
        KST.append(ks)
        z = {}
        for r in RUECK:
            if ks - r < 3:
                continue
            z[r] = geometrie(cl[:ks - r], N)
        BAHN.append(z)
        if (i + 1) % 10 == 0:
            sys.stdout.write(" %d" % (i + 1)); sys.stdout.flush()
    print()
    print()
    print("=" * 88)
    print("DIE BAHN: Geometrie r Klauseln VOR dem Kipppunkt")
    print("k* liegt im Median bei %d von %d Klauseln (%.1f %%)"
          % (np.median(KST), len(CLS[0][0]), 100 * np.median(KST) / len(CLS[0][0])))
    print("=" * 88)
    print("%-10s %10s %12s %12s %14s" % ("r", "Instanzen", "Spann", "Hamming", "Loesungen"))
    for r in RUECK:
        v = [b[r] for b in BAHN if r in b]
        if not v:
            continue
        print("%-10d %10d %12.1f %12.2f %14.1f"
              % (r, len(v), np.mean([x["spann"] for x in v]),
                 np.mean([x["hamming"] for x in v]),
                 np.mean([x["loesungen"] for x in v])))
    print("=" * 88)
    print()
    y = np.log2(np.maximum(np.array(K, float), 1))
    print("Sagt die BAHN die Widerlegungskosten vorher?  (Massstab d(50 %) = +0,56)")
    print("%-22s %11s %11s %10s" % ("Groesse", "Pearson", "Spearman", "p"))
    kand = {}
    for r in (20, 10, 3):
        for nm in ("spann", "hamming"):
            w = [(b[r][nm] if r in b else np.nan) for b in BAHN]
            kand["%s bei r=%d" % (nm, r)] = w
    # Steilheit des Zusammenbruchs: Spann bei r=20 minus Spann bei r=3
    kand["Abfall Spann 20->3"] = [(b[20]["spann"] - b[3]["spann"])
                                  if (20 in b and 3 in b) else np.nan for b in BAHN]
    kand["k* selbst"] = [float(x) for x in KST]
    for nm, w in kand.items():
        w = np.array(w, float)
        gut = ~np.isnan(w)
        if gut.sum() < 10 or w[gut].std() < 1e-9:
            print("%-22s   (zu wenige/konstant)" % nm); continue
        pr = np.corrcoef(w[gut], y[gut])[0, 1]
        sp_ = np.corrcoef(rang(w[gut]), rang(y[gut]))[0, 1]
        p_ = erfc(abs(pr) * sqrt(gut.sum() - 3) / sqrt(2))
        print("%-22s %+10.4f %+11.4f %10.4f%s"
              % (nm, pr, sp_, p_, "  <-" if p_ < 0.05 else ""))
        sys.stdout.flush()

Vollständige Kennzahlen aller 10 Instanzen (2 je n ∈ {30,50,80,100,140}) — die einzige Ebene der 692 MB, die klein genug zum Einbetten ist.

n,m,alpha,cdcl_konflikte,d50,R_anzahl,R_anteil,bC_mittel,bC_min,bC_max,sp_konv,sp_runden,sp_eta_mittel,sp_joker,sp_bias,sortierungen_zufall,sortierungen_gezielt
100,427,4.27,843,11.5180,15,0.03513,0.9033,0.8500,0.9400,0,2000,0.20968,0.10769,0.52097,3000,6
100,427,4.27,670,11.2904,6,0.01405,0.8300,0.6400,0.9400,1,170,0.20873,0.12726,0.39213,3000,6
140,598,4.27,2998,15.2872,0,0.00000,,,,0,2000,0.23622,0.09529,0.37168,500,6
140,598,4.27,2740,15.4087,26,0.04348,0.9333,0.9214,0.9429,1,122,0.19351,0.14551,0.40951,500,6
30,128,4.27,31,4.6589,9,0.07031,0.9111,0.8667,0.9667,1,309,0.24064,0.09424,0.35464,20000,6
30,128,4.27,29,4.4776,4,0.03125,0.9778,0.9667,1.0000,1,395,0.26895,0.06843,0.26829,20000,6
50,213,4.27,63,6.6248,3,0.01408,0.9467,0.9200,0.9600,0,2000,0.23420,0.09666,0.49654,15000,6
50,213,4.27,66,6.5477,1,0.00469,0.7800,0.7800,0.7800,1,200,0.22146,0.11572,0.40301,15000,6
80,342,4.27,482,9.7937,19,0.05556,0.9042,0.8500,0.9625,1,138,0.22085,0.10957,0.42749,8000,6
80,342,4.27,316,9.5109,3,0.00877,0.9042,0.9000,0.9125,0,2000,0.23126,0.08497,0.50138,8000,6

sortierung_id,sortierung_art,k,klausel_index,klausel_in_R,backbone_abs,backbone_anteil,dbackbone,modelle,unsat
0,erzeugung,1,0,0,0,0.0000,0,2,0
0,erzeugung,2,1,0,0,0.0000,0,29,0
0,erzeugung,3,2,0,0,0.0000,0,16,0
0,erzeugung,4,3,0,0,0.0000,0,16,0
0,erzeugung,5,4,0,0,0.0000,0,16,0
0,erzeugung,6,5,0,0,0.0000,0,17,0
0,erzeugung,7,6,1,0,0.0000,0,16,0
0,erzeugung,8,7,0,0,0.0000,0,16,0
0,erzeugung,9,8,0,0,0.0000,0,16,0
0,erzeugung,10,9,0,0,0.0000,0,16,0
0,erzeugung,11,10,0,0,0.0000,0,14,0
0,erzeugung,12,11,0,0,0.0000,0,4,0
0,erzeugung,13,12,0,0,0.0000,0,5,0
0,erzeugung,14,13,0,0,0.0000,0,7,0
0,erzeugung,15,14,0,0,0.0000,0,8,0