14.2 光照的批处理与代价

windowed(dd, rr) = fdiv(rr − dd, rr + dd)     dd = d², rr = range²

#!/usr/bin/env python3
# file: lantern/sim/cg_shade_14_2_light_batching.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_14_2_light_batching.py
"""14.2 光照的批处理与代价：光源数 1→64 的代价曲线，朴素全量 vs 分桶剔除。

**先换衰减形式**：14.1 的 `r²/(r²+d²)` 在数学上处处非零（本场景下最小也有几百个
Q16.16 单位），于是**每个光源都影响每个像素**，剔除无从谈起。本节改用**有限支撑**的
窗口衰减 `(r²−d²)/(r²+d²)`：`d=0` 时 1、`d=r` 时 0、`d>r` 时被 `dd <= rr` 挡掉——
支撑集恰好是射程球，这才让「哪些光源影响哪个像素」成为一条可判的整数谓词。

两级剔除：

  距离阈值（逐像素）  `dd = Σd² <= rr`，纯整数、无开方；命中才做昂贵的归一化与除法
  网格分桶（逐瓦片）  把光源的**屏幕圆盘**与 16×16 的瓦片矩形求交，得到每瓦片的光源表

分桶是**保守**的（瓦片表是超集），所以分桶渲染必须与朴素渲染**逐字节相同**（sha256）。
代价模型 = 求值次数 × 每光源运算量。

产物：out/cg_shade_14_2_light_batching.{json,csv}、out/light_batch_{naive,bucket}_64.png
"""

import hashlib
import json
import os
import sys
import time

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib  # noqa: E402
import cgshade  # noqa: E402
import cg_shade_14_lightlib as L  # noqa: E402

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
SUB1 = cglib.SUB1
SPEC = os.path.join(ROOT, "config", "cg_shade_14_2_lights.json")

# 批处理原语（Ops / gen_lights / screen_disc / build_tiles / camera_basis）已收进
# 本章公共库 cg_shade_14_lightlib，本节只做编排与度量，不重复实现定点数学。


def main():
    cfg = json.load(open(SPEC, encoding="utf-8"))
    spec = cglib.load_spec(os.path.join(ROOT, "config", cfg["sceneSpec"]))
    w, h, bg = spec["width"], spec["height"], spec["backgroundQ8"]
    ambient = cfg["ambientQ16"]
    tile = cfg["tilePixels"]

    verts, normals, albs, tris = L.build_scene_lit(spec)
    cm = L.camera_mvp(spec)
    ltris, dropped = L.build_lit_tris(verts, normals, albs, tris, cm, w, h)
    basis = L.camera_basis(spec)

    all_lights = L.gen_lights(max(cfg["lightCounts"]), spec, cfg)
    out = os.path.join(ROOT, "out")
    os.makedirs(out, exist_ok=True)

    rows, curves = [], []
    evals_max, table_max = None, None
    for n in cfg["lightCounts"]:
        lights = all_lights[:n]
        tiles, nx, ny = L.build_tiles(lights, cm, w, h, basis, tile)

        def naive_cand(px, py, n=n):
            return range(n)

        def bucket_cand(px, py, tl=tiles, t=tile):
            return tl[py // t][px // t]

        rec = set() if n == max(cfg["lightCounts"]) else None
        t0 = time.perf_counter()
        c_n, d_n, cov_n, _cl, loops_n, full_n = L.rasterize_lights(
            ltris, w, h, bg, ambient, lights, naive_cand, L.windowed_attenuation, evals=rec)
        ms_n = (time.perf_counter() - t0) * 1000

        t0 = time.perf_counter()
        c_b, d_b, cov_b, _cl, loops_b, full_b = L.rasterize_lights(
            ltris, w, h, bg, ambient, lights, bucket_cand, L.windowed_attenuation)
        ms_b = (time.perf_counter() - t0) * 1000

        sha_n = cgshade.digest(np.clip(c_n, 0, 255).astype(np.uint8))
        sha_b = cgshade.digest(np.clip(c_b, 0, 255).astype(np.uint8))
        covered = int(cov_n.sum())
        tile_lights = sum(len(tiles[ty][tx]) for ty in range(ny) for tx in range(nx))

        curves.append({
            "n": n, "covered": covered, "tiles": nx * ny, "tile_light_pairs": tile_lights,
            "naive_loops": loops_n, "naive_evals": full_n,
            "bucket_loops": loops_b, "bucket_evals": full_b,
            "eval_ratio": round(full_b / max(1, full_n), 6),
            "naive_ms": round(ms_n, 1), "bucket_ms": round(ms_b, 1),
            "speedup_ms": round(ms_n / max(1e-6, ms_b), 2),
            "identical": sha_n == sha_b, "sha256": sha_n,
        })
        rows.append([n, loops_n, full_n, loops_b, full_b, f"{full_b / max(1, full_n):.6f}",
                     f"{ms_n:.1f}", f"{ms_b:.1f}", f"{ms_n / max(1e-6, ms_b):.2f}",
                     int(sha_n == sha_b)])

        if rec is not None:
            evals_max = rec
            table_max = {(tx, ty, li) for ty in range(ny) for tx in range(nx)
                         for li in tiles[ty][tx]}

        if n == max(cfg["lightCounts"]):
            cglib.save_png(os.path.join(out, "light_batch_naive_64.png"), c_n)
            cglib.save_png(os.path.join(out, "light_batch_bucket_64.png"), c_b)

    # 每光源运算量：真跑一次「距离测试」与一次「完整求值」
    lt = all_lights[0]
    o_test = L.Ops()
    dx, dy, dz = 1000, -2000, 3000
    a = o_test.s(dx, 0)
    b = o_test.s(dy, 0)
    c = o_test.s(dz, 0)
    aa, bb, cc = o_test.m(a, a), o_test.m(b, b), o_test.m(c, c)
    _dd = o_test.s(o_test.s(aa, -bb), -cc)   # 3 减 + 3 乘 + 2 加 = 8 次定点运算
    o_eval = L.Ops()
    world = (0, ONE, 0)
    n = cglib.v_norm((0, ONE, 0))
    to_light = o_eval.s(lt["posQ16"][0], world[0]), o_eval.s(lt["posQ16"][1], world[1]), \
        o_eval.s(lt["posQ16"][2], world[2])
    ndl = o_eval.dot(n, o_eval.norm(to_light))
    if ndl < 0:
        ndl = 0
    _k = o_eval.m(o_eval.m(ndl, o_eval.d(lt["rr"], lt["rr"])), lt["intensityQ16"])

    # 保守性不变量（直接核对）：朴素全量**实际求值**过的 (像素, 光源) 对，其所属瓦片必须
    # 在分桶表里——即 `表 ⊇ 实际求值集`。这是分桶结果与朴素逐字节相同的**充分条件**，
    # 且不依赖任何「射程内像素」的中间缓冲。
    nx, ny = (w + tile - 1) // tile, (h + tile - 1) // tile
    need = {(px // tile, py // tile, li) for (px, py, li) in evals_max}
    missing = sorted(need - table_max)
    violations = len(missing)
    disc_rows = [{"eval_pairs": len(need), "table_pairs": len(table_max),
                  "missing": len(missing), "sample_missing": missing[:8]}]

    # 供 P 层复核的采样：光源屏幕圆盘 + 窗口衰减曲线
    disc_samples = []
    for li in range(min(8, len(all_lights))):
        d = L.screen_disc(cm, all_lights[li]["posQ16"], all_lights[li]["rangeQ16"], w, h, basis)
        disc_samples.append({"i": li, "sx": None if d is None else d[0],
                             "sy": None if d is None else d[1],
                             "rad": None if d is None else d[2]})
    rr0 = all_lights[0]["rr"]
    att_samples = [{"dd": k * (rr0 // 8), "rr": rr0,
                    "att": L.windowed_attenuation(k * (rr0 // 8), rr0)} for k in range(17)]

    last = curves[-1]
    test_ops = o_test.mul + o_test.div + o_test.sqrt + o_test.addsub
    eval_ops = o_eval.mul + o_eval.div + o_eval.sqrt + o_eval.addsub
    naive_ops = last["naive_loops"] * test_ops + last["naive_evals"] * eval_ops
    bucket_ops = last["bucket_loops"] * test_ops + last["bucket_evals"] * eval_ops

    cglib.write_csv(
        os.path.join(out, "cg_shade_14_2_light_batching.csv"),
        ["n", "naive_loops", "naive_evals", "bucket_loops", "bucket_evals", "eval_ratio",
         "naive_ms", "bucket_ms", "speedup", "identical"],
        rows)

    result = {
        "ambientQ16": ambient, "canvas": [w, h], "tile": tile,
        "tris_drawn": len(ltris), "dropped": dropped,
        "max_lights": max(cfg["lightCounts"]), "light_rangeQ16": cfg["lightRangeQ16"],
        "light_positions": [list(lt["posQ16"]) for lt in all_lights],
        "camera": spec["camera"], "width": w, "height": h,
        "table_invariant": {"violations": violations, "rows": disc_rows},
        "disc_samples": disc_samples, "att_samples": att_samples,
        "rangeQ16": cfg["lightRangeQ16"], "rr": rr0, "lightSeed": cfg["lightSeed"],
        "heightBaseQ16": cfg["heightBaseQ16"], "heightSpanQ16": cfg["heightSpanQ16"],
        "ops_per_eval": {"mul": o_eval.mul, "div": o_eval.div, "sqrt": o_eval.sqrt,
                         "addsub": o_eval.addsub, "total": eval_ops},
        "ops_per_test": test_ops,
        "cost_model_64": {"naive_ops": naive_ops, "bucket_ops": bucket_ops,
                          "ratio": round(bucket_ops / max(1, naive_ops), 6)},
        "curves": curves,
    }
    with open(os.path.join(out, "cg_shade_14_2_light_batching.json"), "w", encoding="utf-8") as fh:
        json.dump(result, fh, ensure_ascii=False, indent=2)

    print(f"几何            : {len(verts)} 顶点 / {len(tris)} 三角形  覆盖 {curves[0]['covered']}")
    print(f"瓦片            : {tile}x{tile} 像素  共 {curves[0]['tiles']} 块")
    print(f"{'N':>3} {'朴素循环':>10} {'朴素求值':>9} {'分桶循环':>10} {'分桶求值':>9} "
          f"{'比例':>7} {'朴素ms':>7} {'分桶ms':>7} {'加速':>6} {'逐字节':>6}")
    for r in curves:
        print(f"{r['n']:>3} {r['naive_loops']:>10} {r['naive_evals']:>9} "
              f"{r['bucket_loops']:>10} {r['bucket_evals']:>9} {r['eval_ratio']:>7.3f} "
              f"{r['naive_ms']:>7.1f} {r['bucket_ms']:>7.1f} {r['speedup_ms']:>6.2f} "
              f"{str(bool(r['identical'])):>6}")
    print(f"每光源运算量    : 距离测试 {test_ops} ops  完整求值 {eval_ops} ops "
          f"(mul={o_eval.mul} div={o_eval.div} sqrt={o_eval.sqrt} addsub={o_eval.addsub})")
    print(f"代价模型 (N=64) : 朴素 {naive_ops} ops vs 分桶 {bucket_ops} ops  "
          f"比例 {result['cost_model_64']['ratio']:.4f}")
    print(f"保守性不变量    : 实际求值 {len(need)} 对，分桶表 {len(table_max)} 对，"
          f"表未覆盖的 = {violations}")
    return 0


if __name__ == "__main__":
    sys.exit(main())

// snippet: core/test/cg_shade_14_2_light_batching.test.ts（节选：圆盘与代价曲线）
test("光源屏幕圆盘的 (sx, sy, rad) 与 S 侧逐项相同", () => {
  for (const d of stats.disc_samples) {
    const got = screenDisc(camera, V(d.pos), BigInt(d.rangeQ16), stats.width, stats.height, basis);
    assert.equal(got.sx, d.sx, `sx @${d.i}`);
    assert.equal(got.sy, d.sy, `sy @${d.i}`);
    assert.equal(got.rad, d.rad, `rad @${d.i}`);
  }
});

test("代价曲线：分桶完整求值恒等于朴素，且逐字节相同", () => {
  for (const row of stats.curves) {
    assert.equal(row.bucket_evals, row.naive_evals, `evals @N=${row.n}`);
    assert.equal(row.eval_ratio, 1, `eval_ratio @N=${row.n}`);
    assert.equal(row.identical, true, `identical @N=${row.n}`);
  }
  const last = stats.curves[stats.curves.length - 1];
  assert.equal(last.n, 64);
  assert.ok(last.bucket_loops < last.naive_loops, "分桶必须减少距离测试循环");
});

✔ 光源屏幕圆盘的 (sx, sy, rad) 与 S 侧逐项相同
✔ 窗口衰减：d=0 为 1、d=r 为 0、d>r 为负
✔ 代价曲线：分桶完整求值恒等于朴素，且逐字节相同
✔ 代价模型运算量：距离测试 8 ops、完整求值 20 ops，与 S 侧一致
ℹ tests 4  ℹ pass 4  ℹ fail 0