1import json
  2import math
  3import random
  4from pathlib import Path
  5
  6import numpy as np
  7
  8
  9def phi(x):
 10    x = np.asarray(x)
 11    return np.stack([
 12        0.55 * np.sin(2 * np.pi * x) + 0.25 * np.cos(4 * np.pi * x),
 13        0.55 * np.cos(2 * np.pi * x) - 0.25 * np.sin(6 * np.pi * x),
 14    ], axis=-1)
 15
 16
 17def generate(r, ell=2, n=7000, burn=45, seed=0):
 18    # Sample the stationary attractor without long forward floating-point
 19    # doubling orbits (which eventually lose all mantissa bits).  For a
 20    # current uniform x, choose inverse branches x'=(x+k)/ell and accumulate
 21    # h = sum_{j>=0} r^j phi(x_{-1-j}).
 22    rng = np.random.default_rng(seed)
 23    x = rng.random(n)
 24    h = np.zeros((n, 2))
 25    weight = 1.0
 26    for _ in range(burn):
 27        branch = rng.integers(0, ell, size=n)
 28        x = (x + branch) / ell
 29        h += weight * phi(x)
 30        weight *= r
 31    return np.column_stack([x, h])
 32
 33
 34def box_dimension(z, scales=(6, 8, 12, 16, 24, 32)):
 35    # Normalize coordinates so grid resolution has comparable meaning.
 36    lo, hi = z.min(0), z.max(0)
 37    q = (z - lo) / (hi - lo + 1e-12)
 38    counts = []
 39    for m in scales:
 40        cells = np.floor(np.minimum(q, 1 - 1e-12) * m).astype(np.int64)
 41        counts.append(np.unique(cells, axis=0).shape[0])
 42    x = np.log(np.asarray(scales, float))
 43    y = np.log(np.asarray(counts, float))
 44    # Fine-scale endpoint bins are most informative, while retaining robustness.
 45    slope = float(np.polyfit(x[-5:], y[-5:], 1)[0])
 46    return slope, counts
 47
 48
 49def correlation_dimension(z, seed=123, pairs=30000):
 50    rng = np.random.default_rng(seed)
 51    n = len(z)
 52    ii = rng.integers(0, n, pairs)
 53    jj = rng.integers(0, n, pairs)
 54    keep = ii != jj
 55    d = np.linalg.norm(z[ii[keep]] - z[jj[keep]], axis=1)
 56    d = d[d > 1e-15]
 57    if len(d) < 100:
 58        return float("nan")
 59    # Pairwise distance quantiles avoid the zero/finite-sample extremes.
 60    radii = np.geomspace(np.quantile(d, .01), np.quantile(d, .35), 14)
 61    c = np.array([(d < q).mean() for q in radii])
 62    good = (c > 3 / len(d)) & (c < .4)
 63    slope = float(np.polyfit(np.log(radii[good]), np.log(c[good]), 1)[0])
 64    return slope
 65
 66
 67def contraction_check(r=.8, ell=2, n=30):
 68    # Same driver, two fibers: exact difference is r^n times its initial value.
 69    x = .371234
 70    a, b = np.array([1.0, -.4]), np.array([-.3, .8])
 71    vals = []
 72    for t in range(n + 1):
 73        vals.append(np.linalg.norm(a - b))
 74        a = r * a + phi(x)
 75        b = r * b + phi(x)
 76        x = (ell * x) % 1
 77    slope = float(np.polyfit(np.arange(1, n + 1), np.log(vals[1:]), 1)[0])
 78    return slope, math.log(r), float(vals[-1] / vals[0]), r ** n
 79
 80
 81def main():
 82    np.random.seed(7)
 83    random.seed(7)
 84    ell = 2
 85    rs = [0.55, 0.65, 0.70, 1 / math.sqrt(2), 0.75, 0.85, 0.95]
 86    rows = []
 87    for k, r in enumerate(rs):
 88        z = generate(r, ell=ell, seed=100 + k)
 89        predicted = 1.0 + min(2.0, math.log(ell) / abs(math.log(r)))
 90        bd, boxes = box_dimension(z)
 91        cd = correlation_dimension(z, seed=200 + k)
 92        g = math.log(ell) + 2 * math.log(r)
 93        rows.append({
 94            "r": r, "ell_det": ell * r * r, "g": g,
 95            "predicted_dimension": predicted,
 96            "box_dimension": bd, "correlation_dimension": cd,
 97            "boxes": boxes,
 98        })
 99    cslope, clog, cfinal, cpred = contraction_check()
100    result = {
101        "seed": 7, "ell": ell, "dimension": 2,
102        "threshold_r": 1 / math.sqrt(ell),
103        "contraction": {"observed_log_slope": cslope, "predicted_log_slope": clog,
104                        "observed_ratio_n30": cfinal, "predicted_ratio_n30": cpred},
105        "sweep": rows,
106        "notes": "Volume exponent is exact for A=rI: log(ell)+log|det(A)|=log(ell*r^2).",
107    }
108    Path("toy_results.json").write_text(json.dumps(result, indent=2))
109    print(json.dumps(result, indent=2))
110
111
112if __name__ == "__main__":
113    main()