r(d) Curve Shape

"""The shape of the r(d) curve carries more information than the single
crossing point d(50%). The evidence behind ordnung_kurvenform.html.

Builds on d50_messung.py (same instance generation, same
propagation measurement r(d)), but samples the whole curve instead
of only searching for the 50% point, and derives four features from it:

  d(50%)    inflection point (where r(d) = 0.5), as before
  Area      integral of r(d) over the measured range
  Early     r(d) at the smallest measured d (instability)
  Slope     numerical derivative of r(d) at the inflection point

Call: python3 kurvenform_messung.py [n-list=40,60,80,100] [instances=200] [seed=7]
Writes daten_kurvenform.json.
"""
import json
import sys

import numpy as np

from d50_messung import instanzen, r_von_d


def kurve(klauseln, n, dgrid, proben, rng):
    return np.array([r_von_d(klauseln, n, d, proben, rng) for d in dgrid])


def merkmale(dgrid, rd):
    # Area: trapezoidal rule
    flaeche = float(np.trapezoid(rd, dgrid)) if hasattr(np, "trapezoid") else float(np.trapz(rd, dgrid))
    fruehe = float(rd[0])
    # interpolate d(50%) through the measured grid
    if rd[-1] < 0.5:
        d50 = float(dgrid[-1])
    elif rd[0] > 0.5:
        d50 = float(dgrid[0])
    else:
        i = int(np.searchsorted(rd, 0.5))
        d_lo, d_hi = dgrid[i - 1], dgrid[i]
        r_lo, r_hi = rd[i - 1], rd[i]
        frac = (0.5 - r_lo) / (r_hi - r_lo) if r_hi != r_lo else 0.0
        d50 = float(d_lo + frac * (d_hi - d_lo))
    # slope: central difference at the point closest to d(50%)
    i = int(np.argmin(np.abs(dgrid - d50)))
    i = min(max(i, 1), len(dgrid) - 2)
    steigung = float((rd[i + 1] - rd[i - 1]) / (dgrid[i + 1] - dgrid[i - 1]))
    return d50, flaeche, fruehe, steigung


def partielle_korrelation(x, y, kontrolle):
    """Residual of x and y after linear regression on `kontrolle`, then Pearson."""
    A = np.column_stack([np.ones(len(kontrolle)), kontrolle])
    bx, *_ = np.linalg.lstsq(A, x, rcond=None)
    by, *_ = np.linalg.lstsq(A, y, rcond=None)
    rx = x - A @ bx
    ry = y - A @ by
    if rx.std() < 1e-12 or ry.std() < 1e-12:
        return 0.0
    return float(np.corrcoef(rx, ry)[0, 1])


def lauf(ns=(40, 60, 80, 100), anzahl=200, seed=7, proben=150):
    ergebnis = {}
    for n in ns:
        R = instanzen(n, anzahl, seed + n)
        unsat = [r for r in R if r["sat"] == 0]
        dgrid = np.unique(np.linspace(2, min(n, 20), 10).round().astype(int))
        rng = np.random.default_rng(3000 + n)

        d50s, flaechen, fruehe, steigungen, konflikte = [], [], [], [], []
        for r in unsat:
            rd = kurve(r["klauseln"], n, dgrid, proben, rng)
            d50, fl, fr, st = merkmale(dgrid, rd)
            d50s.append(d50); flaechen.append(fl); fruehe.append(fr); steigungen.append(st)
            konflikte.append(np.log2(max(r["konflikte"], 1)))
        d50s, flaechen, fruehe, steigungen, y = map(np.array, (d50s, flaechen, fruehe, steigungen, konflikte))

        def r(x):
            return float(np.corrcoef(x, y)[0, 1]) if x.std() > 1e-12 else 0.0

        r_d50, r_fl, r_fr, r_st = r(d50s), r(flaechen), r(fruehe), r(steigungen)
        r_st_partial = partielle_korrelation(steigungen, y, d50s)

        # best single d (grid point) as predictor
        rng2 = np.random.default_rng(4000 + n)
        beste_d, beste_r = None, 0.0
        for j, d in enumerate(dgrid):
            rds = np.array([r_von_d(rr["klauseln"], n, int(d), 80, rng2) for rr in unsat])
            if rds.std() < 1e-9:
                continue
            rr_ = float(np.corrcoef(rds, y)[0, 1])
            if abs(rr_) > abs(beste_r):
                beste_r, beste_d = rr_, int(d)

        ergebnis[n] = {
            "n": n, "unsat_instanzen": len(unsat), "dgrid": [int(x) for x in dgrid],
            "r_d50": r_d50, "r_flaeche": r_fl, "r_fruehe": r_fr,
            "r_steigung": r_st, "r_steigung_partial": r_st_partial,
            "bestes_d": beste_d, "bestes_d_r": beste_r,
        }
        print(f"n={n:3d}  UNSAT={len(unsat):3d}  d(50%) r={r_d50:+.2f}  Area r={r_fl:+.2f}  "
              f"Early r={r_fr:+.2f}  Slope r={r_st:+.2f} (partial {r_st_partial:+.2f})  "
              f"best d={beste_d} (r={beste_r:+.2f})")

    with open("daten_kurvenform.json", "w") as f:
        json.dump(ergebnis, f, indent=1)
    print("\nsaved -> daten_kurvenform.json")


if __name__ == "__main__":
    a = sys.argv[1:]
    ns = tuple(int(x) for x in a[0].split(",")) if len(a) > 0 else (40, 60, 80, 100)
    lauf(ns=ns,
        anzahl=int(a[1]) if len(a) > 1 else 200,
        seed=int(a[2]) if len(a) > 2 else 7)

Fresh instances (α=4.267), r(d) over a 10-point grid d∈[2,20] with 150 samples each. Same signs and order of magnitude as in the text (area negatively correlated, d(50%) positively correlated), individual values differ — own sample, coarser grid, different seed. File is genuinely runnable, not reproduced identically.

{
 "40": {
  "n": 40,
  "unsat_instanzen": 73,
  "dgrid": [
   2,
   4,
   6,
   8,
   10,
   12,
   14,
   16,
   18,
   20
  ],
  "r_d50": 0.3479606073029222,
  "r_flaeche": -0.36251842134884277,
  "r_fruehe": -0.17370088667120148,
  "r_steigung": -0.05819098798997463,
  "r_steigung_partial": -0.11979069357446034,
  "bestes_d": 6,
  "bestes_d_r": -0.3020972629964603
 },
 "60": {
  "n": 60,
  "unsat_instanzen": 94,
  "dgrid": [
   2,
   4,
   6,
   8,
   10,
   12,
   14,
   16,
   18,
   20
  ],
  "r_d50": 0.2243167757534666,
  "r_flaeche": -0.15259385658657534,
  "r_fruehe": 0.15194037920638617,
  "r_steigung": 0.0402972934352017,
  "r_steigung_partial": -0.007313354762488981,
  "bestes_d": 16,
  "bestes_d_r": -0.14619166732875394
 },
 "80": {
  "n": 80,
  "unsat_instanzen": 85,
  "dgrid": [
   2,
   4,
   6,
   8,
   10,
   12,
   14,
   16,
   18,
   20
  ],
  "r_d50": 0.050865332976337885,
  "r_flaeche": -0.23906609449956126,
  "r_fruehe": -0.14568779293536763,
  "r_steigung": 0.015305609677908677,
  "r_steigung_partial": 0.008388184736164003,
  "bestes_d": 10,
  "bestes_d_r": -0.25859175329911815
 },
 "100": {
  "n": 100,
  "unsat_instanzen": 98,
  "dgrid": [
   2,
   4,
   6,
   8,
   10,
   12,
   14,
   16,
   18,
   20
  ],
  "r_d50": 0.18863854659675688,
  "r_flaeche": -0.2818229821799546,
  "r_fruehe": 0.0,
  "r_steigung": -0.005817917565504143,
  "r_steigung_partial": -0.03533360218296536,
  "bestes_d": 12,
  "bestes_d_r": -0.2792504496684058
 }
}