indicatrix 0.4.0

//! Next-event estimation (NEE), balance-heuristic MIS, and 1D/2D distribution
//! importance sampling ULP equivalence checks.
//!
//! Compares `transport_functions.wgsl`'s GPU implementations against the real CPU
//! functions in `optics::raytracer::scattering` and `renderer::env_map`.

use glam::Vec3;

use super::{SHADER_SRC, UlpAccumulator, UlpCheckResult};
use crate::{
    geometry::plane::GpuFacetPlane,
    optics::{
        polarization::StokesVector,
        raytracer::{
            EnvironmentSource, FacetFinish, balance_heuristic, build_plane_soa,
            scattering::{
                NeeContext, nee_contribution_frosted_exterior, nee_contribution_hg_scatter,
            },
        },
    },
    renderer::{
        env_map::{Distribution1D, Distribution2D, EnvironmentMap},
        env_map_gpu::HdrEnvGpuData,
        gpu::compute,
    },
};

// ---------------------------------------------------------------------------------
// 1. Balance heuristic
// ---------------------------------------------------------------------------------

#[repr(C)]
#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
pub struct BalanceHeuristicCase {
    pdf_a: f32,
    pdf_b: f32,
    _pad0: f32,
    _pad1: f32,
}

const _: () = assert!(size_of::<BalanceHeuristicCase>() == 16);

const BALANCE_HEURISTIC_ULP_BUDGET: u32 = 4;
const BALANCE_HEURISTIC_ABS_FLOOR: f32 = 1e-7;

fn build_balance_heuristic_cases() -> Vec<BalanceHeuristicCase> {
    let mut cases = Vec::new();
    let pdfs = [
        0.0f32, 1e-6, 1e-4, 0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 5.0, 10.0, 100.0, 1e4,
    ];
    for &a in &pdfs {
        for &b in &pdfs {
            cases.push(BalanceHeuristicCase {
                pdf_a: a,
                pdf_b: b,
                _pad0: 0.0,
                _pad1: 0.0,
            });
        }
    }
    let adversarial = [
        (-1.0f32, 1.0f32),
        (1.0, -1.0),
        (-0.5, -0.5),
        (0.0, 0.0),
        (1.0, 0.0),
        (0.0, 1.0),
        (1e-7, 1e-7),
        (1e7, 1e7),
        (1e-5, 1e5),
    ];
    for &(a, b) in &adversarial {
        cases.push(BalanceHeuristicCase {
            pdf_a: a,
            pdf_b: b,
            _pad0: 0.0,
            _pad1: 0.0,
        });
    }
    cases
}

/// Runs the Tier 2 ULP check for [`balance_heuristic`].
#[must_use]
pub fn run_balance_heuristic(
    ctx: &crate::renderer::gpu::GpuContext,
) -> UlpCheckResult<BalanceHeuristicCase> {
    let cases = build_balance_heuristic_cases();
    let total = cases.len();
    let in_buf = compute::upload(
        &ctx.device,
        "balance heuristic in",
        &cases,
        wgpu::BufferUsages::STORAGE,
    );
    let out_buf = compute::zeroed_buffer::<f32>(
        &ctx.device,
        "balance heuristic out",
        total,
        wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
    );
    let pipeline = compute::create_compute_pipeline(
        &ctx.device,
        "balance_heuristic_main",
        SHADER_SRC,
        "balance_heuristic_main",
    );
    let bind_group = compute::bind_buffers(
        &ctx.device,
        "balance heuristic bind group",
        &pipeline,
        &[(60, &in_buf), (61, &out_buf)],
    );
    let workgroups = (total as u32).div_ceil(64);
    compute::dispatch_and_wait(
        &ctx.device,
        &ctx.queue,
        &pipeline,
        &bind_group,
        (workgroups, 1, 1),
    );
    let gpu_out: Vec<f32> = compute::readback(&ctx.device, &ctx.queue, &out_buf, total);

    let mut acc = UlpAccumulator::new(
        "balance_heuristic",
        BALANCE_HEURISTIC_ULP_BUDGET,
        BALANCE_HEURISTIC_ABS_FLOOR,
    );
    for (idx, case) in cases.iter().enumerate() {
        let cpu = balance_heuristic(case.pdf_a, case.pdf_b);
        acc.record(case, "weight", cpu, gpu_out[idx]);
    }
    acc.finish()
}

// ---------------------------------------------------------------------------------
// Synthetic environment map for distribution tests
// ---------------------------------------------------------------------------------

fn synthetic_test_map() -> EnvironmentMap {
    let width = 8;
    let height = 4;
    let mut pixels = Vec::with_capacity(width * height);
    for y in 0..height {
        for x in 0..width {
            let r = (((x * 3 + y * 5) % 7) as f32).mul_add(0.4, 0.1);
            let g = (((x * 2 + y * 3) % 5) as f32).mul_add(0.3, 0.1);
            let b = (((x * 5 + y * 2) % 6) as f32).mul_add(0.5, 0.1);
            pixels.push([r, g, b]);
        }
    }
    EnvironmentMap::from_rgb(width, height, pixels).expect("valid synthetic env map")
}

// ---------------------------------------------------------------------------------
// 2. 1D distribution binary search (`dist1d_find_bucket`)
// ---------------------------------------------------------------------------------

#[repr(C)]
#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
pub struct Dist1dFindBucketCase {
    cdf_start: u32,
    count: u32,
    u: f32,
    _pad0: f32,
}

const _: () = assert!(size_of::<Dist1dFindBucketCase>() == 16);

fn build_dist1d_find_bucket_cases(map: &EnvironmentMap) -> (Vec<Dist1dFindBucketCase>, Vec<usize>) {
    let mut cases = Vec::new();
    let mut cpu_expected = Vec::new();
    let width = map.width() as u32;
    let height = map.height() as u32;
    let dist: &Distribution2D = map.distribution();

    let u_vals = [
        0.0f32, 0.001, 0.05, 0.12, 0.25, 0.38, 0.5, 0.63, 0.77, 0.89, 0.95, 0.999, 1.0,
    ];

    // Conditional distributions (rows)
    for (row_idx, row_dist) in dist.conditional().iter().enumerate() {
        let row_dist: &Distribution1D = row_dist;
        let cdf_start = (row_idx as u32) * (width + 1);
        for &u in &u_vals {
            cases.push(Dist1dFindBucketCase {
                cdf_start,
                count: width,
                u,
                _pad0: 0.0,
            });
            cpu_expected.push(row_dist.find_bucket(u));
        }
    }

    // Marginal distribution
    let marginal = dist.marginal();
    let marginal_cdf_start = height * (width + 1);
    for &u in &u_vals {
        cases.push(Dist1dFindBucketCase {
            cdf_start: marginal_cdf_start,
            count: height,
            u,
            _pad0: 0.0,
        });
        cpu_expected.push(marginal.find_bucket(u));
    }

    (cases, cpu_expected)
}

/// Runs the Tier 2 check for [`Distribution1D::find_bucket`].
#[must_use]
pub fn run_dist1d_find_bucket(
    ctx: &crate::renderer::gpu::GpuContext,
) -> UlpCheckResult<Dist1dFindBucketCase> {
    let map = synthetic_test_map();
    let gpu_data = HdrEnvGpuData::upload(&ctx.device, &map);
    let (cases, cpu_expected) = build_dist1d_find_bucket_cases(&map);
    let total = cases.len();

    let in_buf = compute::upload(
        &ctx.device,
        "dist1d find bucket in",
        &cases,
        wgpu::BufferUsages::STORAGE,
    );
    let out_buf = compute::zeroed_buffer::<u32>(
        &ctx.device,
        "dist1d find bucket out",
        total,
        wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
    );
    let pipeline = compute::create_compute_pipeline(
        &ctx.device,
        "dist1d_find_bucket_main",
        SHADER_SRC,
        "dist1d_find_bucket_main",
    );
    let bind_group = compute::bind_buffers(
        &ctx.device,
        "dist1d find bucket bind group",
        &pipeline,
        &[(63, &gpu_data.dist_cdf), (67, &in_buf), (68, &out_buf)],
    );
    let workgroups = (total as u32).div_ceil(64);
    compute::dispatch_and_wait(
        &ctx.device,
        &ctx.queue,
        &pipeline,
        &bind_group,
        (workgroups, 1, 1),
    );
    let gpu_out: Vec<u32> = compute::readback(&ctx.device, &ctx.queue, &out_buf, total);

    let mut acc = UlpAccumulator::new("dist1d_find_bucket", 0, 0.0);
    for (idx, case) in cases.iter().enumerate() {
        acc.record(
            case,
            "bucket",
            cpu_expected[idx] as f32,
            gpu_out[idx] as f32,
        );
    }
    acc.finish()
}

// ---------------------------------------------------------------------------------
// 3. 2D distribution continuous sampling (`dist2d_sample`)
// ---------------------------------------------------------------------------------

#[repr(C)]
#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
pub struct Dist2dSampleCase {
    u0: f32,
    u1: f32,
    _pad0: f32,
    _pad1: f32,
}

const _: () = assert!(size_of::<Dist2dSampleCase>() == 16);

const DIST2D_SAMPLE_ULP_BUDGET: u32 = 48;
const DIST2D_SAMPLE_ABS_FLOOR: f32 = 1e-5;

fn build_dist2d_sample_cases() -> Vec<Dist2dSampleCase> {
    let mut cases = Vec::new();
    let steps = 14;
    for i in 0..=steps {
        let u0 = (i as f32 / steps as f32).mul_add(0.998, 0.001);
        for j in 0..=steps {
            let u1 = (j as f32 / steps as f32).mul_add(0.998, 0.001);
            cases.push(Dist2dSampleCase {
                u0,
                u1,
                _pad0: 0.0,
                _pad1: 0.0,
            });
        }
    }
    let adversarial = [
        (0.0f32, 0.0f32),
        (0.0, 0.5),
        (0.5, 0.0),
        (0.999_999, 0.999_999),
        (1e-6, 1e-6),
        (0.5, 0.999_999),
        (0.999_999, 0.5),
        (0.25, 0.75),
        (0.75, 0.25),
    ];
    for &(u0, u1) in &adversarial {
        cases.push(Dist2dSampleCase {
            u0,
            u1,
            _pad0: 0.0,
            _pad1: 0.0,
        });
    }
    cases
}

/// Runs the Tier 2 check for [`EnvironmentMap::sample`].
#[must_use]
pub fn run_dist2d_sample(
    ctx: &crate::renderer::gpu::GpuContext,
) -> UlpCheckResult<Dist2dSampleCase> {
    let map = synthetic_test_map();
    let gpu_data = HdrEnvGpuData::upload(&ctx.device, &map);
    let cases = build_dist2d_sample_cases();
    let total = cases.len();

    let in_buf = compute::upload(
        &ctx.device,
        "dist2d sample in",
        &cases,
        wgpu::BufferUsages::STORAGE,
    );
    let out_buf = compute::zeroed_buffer::<f32>(
        &ctx.device,
        "dist2d sample out",
        total * 7,
        wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
    );
    let pipeline = compute::create_compute_pipeline(
        &ctx.device,
        "dist2d_sample_main",
        SHADER_SRC,
        "dist2d_sample_main",
    );
    let bind_group = compute::bind_buffers(
        &ctx.device,
        "dist2d sample bind group",
        &pipeline,
        &[
            (62, &gpu_data.dist_func),
            (63, &gpu_data.dist_cdf),
            (64, &gpu_data.dist_dims),
            (65, &gpu_data.texels),
            (66, &gpu_data.dims),
            (69, &in_buf),
            (70, &out_buf),
        ],
    );
    let workgroups = (total as u32).div_ceil(64);
    compute::dispatch_and_wait(
        &ctx.device,
        &ctx.queue,
        &pipeline,
        &bind_group,
        (workgroups, 1, 1),
    );
    let gpu_out: Vec<f32> = compute::readback(&ctx.device, &ctx.queue, &out_buf, total * 7);

    let mut acc = UlpAccumulator::new(
        "dist2d_sample",
        DIST2D_SAMPLE_ULP_BUDGET,
        DIST2D_SAMPLE_ABS_FLOOR,
    );
    for (idx, case) in cases.iter().enumerate() {
        let (cpu_dir, cpu_rgb, cpu_pdf) = map.sample(case.u0, case.u1);
        let base = idx * 7;
        for (c_idx, comp) in ["dir_x", "dir_y", "dir_z"].iter().enumerate() {
            acc.record(case, comp, cpu_dir[c_idx], gpu_out[base + c_idx]);
        }
        for (c_idx, comp) in ["rgb_r", "rgb_g", "rgb_b"].iter().enumerate() {
            acc.record(case, comp, cpu_rgb[c_idx], gpu_out[base + 3 + c_idx]);
        }
        let gpu_pdf = gpu_out[base + 6];
        if cpu_dir.y < SOUTH_POLE_CONDITIONING_Y {
            let gpu_dir = Vec3::new(gpu_out[base], gpu_out[base + 1], gpu_out[base + 2]);
            acc.record(
                case,
                "pdf*sin_theta (south pole, see SOUTH_POLE_CONDITIONING_Y)",
                cpu_pdf * sin_theta_of(cpu_dir),
                gpu_pdf * sin_theta_of(gpu_dir),
            );
        } else {
            acc.record(case, "pdf", cpu_pdf, gpu_pdf);
        }
    }
    acc.finish()
}

/// `dir.y` below which [`run_dist2d_sample`] compares `pdf * sin(theta)` (each side's
/// own `sin(theta)`, read off its own returned direction) instead of the raw pdf.
///
/// The solid-angle pdf is `pdf_uv / (2 PI^2 sin(theta))` with `theta = v * PI`, and
/// `v` is an `f32` in `[0, 1)`: next to the SOUTH pole (`v -> 1`) its spacing is `6e-8`,
/// so one rounding difference in `v` -- e.g. the GPU's division in
/// `dist1d_sample_continuous` rounding the other way than the CPU's -- moves the
/// reported pdf by `6e-8 / (1 - v)` relative. At `u1 = 0.999999` that is `1 - v =
/// 1.6e-6`, a 3.7 % swing (559043 ULP) from a single ULP of
/// `v`, on a direction whose own components are already exempted as near-zero. Both
/// sides are self-consistent (each pdf IS the density of its own direction), so the
/// ULP comparison is meaningless there; `pdf * sin(theta)` cancels the `1/sin(theta)`
/// amplification and compares the well-conditioned `pdf_uv / (2 PI^2)` instead. The
/// budget-48 threshold: `6e-8 / (1 - v) < 48 * 6e-8` needs `1 - v > 0.02`, i.e.
/// `theta < 176.4 degrees`, `y > -cos(0.02 PI) = -0.998`. The north pole needs no such
/// treatment: `f32` is dense near `v = 0`.
const SOUTH_POLE_CONDITIONING_Y: f32 = -0.998;

/// `sin(theta)` of a unit direction, `|(x, z)|` -- the same expression
/// `EnvironmentMap::pdf` uses.
fn sin_theta_of(dir: Vec3) -> f32 {
    dir.x.hypot(dir.z)
}

// ---------------------------------------------------------------------------------
// 4. 2D distribution pdf lookup (`dist2d_pdf`)
// ---------------------------------------------------------------------------------

#[repr(C)]
#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
pub struct Dist2dPdfCase {
    dir: [f32; 3],
    _pad0: f32,
}

const _: () = assert!(size_of::<Dist2dPdfCase>() == 16);

const DIST2D_PDF_ULP_BUDGET: u32 = 32;
const DIST2D_PDF_ABS_FLOOR: f32 = 1e-5;

/// Every direction here must map to a `(u, v)` strictly INSIDE a texel bucket of the
/// 8x4 [`synthetic_test_map`], never onto a bucket edge.
///
/// `EnvironmentMap::pdf` is piecewise constant, so it is discontinuous at every
/// `u * width` / `v * height` integer: a direction sitting exactly on such an edge
/// gets bucket `k` or `k-1` depending on the last ULP of `atan2`/`acos`, and the two
/// neighbouring texels of this map differ by up to 10x in luminance. The CPU's `libm`
/// `atan2(x, x)` rounds to exactly `PI/4` (so `u = 0.125`, `u * 8 == 1.0`, column 1);
/// the GPU's lands one ULP below (column 0) -- `cpu = 5.706e-2`
/// vs `gpu = 5.915e-3`, a 27-million-ULP "divergence" purely from this edge
/// ambiguity (a naive grid puts many cases on an edge: 89 of 154 candidate directions
/// sat on an edge to within `1e-5`). The grid
/// below therefore offsets both angles by an irrational-ish fraction of a step, and
/// the hand-picked directions avoid `phi` multiples of `PI/8` and `theta` multiples of
/// `PI/4` (the two near-pole ones carry a small `x` so `phi` is neither `0` nor `PI`,
/// both of which are column edges); the two poles stay (both sides return exactly
/// `0.0` there).
fn build_dist2d_pdf_cases() -> Vec<Dist2dPdfCase> {
    let mut cases = Vec::new();
    let hand_picked = [
        Vec3::Y,
        Vec3::NEG_Y,
        Vec3::new(0.007, 0.9999, 0.01).normalize(),
        Vec3::new(0.007, -0.9999, 0.01).normalize(),
        Vec3::new(0.3, 0.4, 0.5).normalize(),
        Vec3::new(-1.0, 2.0, -0.5).normalize(),
        Vec3::new(0.9, -0.1, -0.2).normalize(),
        Vec3::new(-0.2, -0.7, 0.6).normalize(),
    ];
    for &d in &hand_picked {
        cases.push(Dist2dPdfCase {
            dir: d.to_array(),
            _pad0: 0.0,
        });
    }
    for theta_i in 0..10 {
        let theta = (theta_i as f32 + 0.37) * std::f32::consts::PI / 10.0;
        let (sin_t, cos_t) = theta.sin_cos();
        for phi_i in 0..16 {
            let phi = (phi_i as f32 + 0.29) * 2.0 * std::f32::consts::PI / 16.0;
            let (sin_p, cos_p) = phi.sin_cos();
            let d = Vec3::new(sin_t * sin_p, cos_t, sin_t * cos_p);
            cases.push(Dist2dPdfCase {
                dir: d.to_array(),
                _pad0: 0.0,
            });
        }
    }
    cases
}

/// Runs the Tier 2 check for [`EnvironmentMap::pdf`].
#[must_use]
pub fn run_dist2d_pdf(ctx: &crate::renderer::gpu::GpuContext) -> UlpCheckResult<Dist2dPdfCase> {
    let map = synthetic_test_map();
    let gpu_data = HdrEnvGpuData::upload(&ctx.device, &map);
    let cases = build_dist2d_pdf_cases();
    let total = cases.len();

    let in_buf = compute::upload(
        &ctx.device,
        "dist2d pdf in",
        &cases,
        wgpu::BufferUsages::STORAGE,
    );
    let out_buf = compute::zeroed_buffer::<f32>(
        &ctx.device,
        "dist2d pdf out",
        total,
        wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
    );
    let pipeline = compute::create_compute_pipeline(
        &ctx.device,
        "dist2d_pdf_main",
        SHADER_SRC,
        "dist2d_pdf_main",
    );
    let bind_group = compute::bind_buffers(
        &ctx.device,
        "dist2d pdf bind group",
        &pipeline,
        &[
            (62, &gpu_data.dist_func),
            (64, &gpu_data.dist_dims),
            (71, &in_buf),
            (72, &out_buf),
        ],
    );
    let workgroups = (total as u32).div_ceil(64);
    compute::dispatch_and_wait(
        &ctx.device,
        &ctx.queue,
        &pipeline,
        &bind_group,
        (workgroups, 1, 1),
    );
    let gpu_out: Vec<f32> = compute::readback(&ctx.device, &ctx.queue, &out_buf, total);

    let mut acc = UlpAccumulator::new("dist2d_pdf", DIST2D_PDF_ULP_BUDGET, DIST2D_PDF_ABS_FLOOR);
    for (idx, case) in cases.iter().enumerate() {
        let cpu = map.pdf(Vec3::from_array(case.dir));
        acc.record(case, "pdf", cpu, gpu_out[idx]);
    }
    acc.finish()
}

// ---------------------------------------------------------------------------------
// 5. NEE contribution: frosted exterior
// ---------------------------------------------------------------------------------

#[repr(C)]
#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
pub struct NeeFrostedExteriorCase {
    ext_normal: [f32; 3],
    rng_seed: u32,
    bounce: u32,
    _pad0: f32,
    _pad1: f32,
    _pad2: f32,
    lambdas: [f32; 8],
    stokes: [[f32; 4]; 8],
    radiance_in: [f32; 8],
}

const _: () = assert!(size_of::<NeeFrostedExteriorCase>() == 224);

const NEE_FROSTED_ULP_BUDGET: u32 = 64;
const NEE_FROSTED_ABS_FLOOR: f32 = 1e-4;

fn build_nee_frosted_cases() -> Vec<NeeFrostedExteriorCase> {
    let mut cases = Vec::new();
    let normals = [
        Vec3::Y,
        Vec3::NEG_Y,
        Vec3::X,
        Vec3::Z,
        Vec3::new(0.5, 0.8, 0.3).normalize(),
        Vec3::new(-0.7, 0.2, -0.6).normalize(),
    ];
    let seeds = [12_345u32, 98_765, 314_159, 271_828];
    let bounces = [0u32, 1, 2, 4];
    let lambdas: [f32; 8] = std::array::from_fn(|k| (k as f32).mul_add(45.0, 400.0));
    let stokes: [[f32; 4]; 8] = std::array::from_fn(|k| {
        [
            (k as f32).mul_add(0.1, 1.0),
            (k as f32) * 0.05,
            -(k as f32) * 0.03,
            0.01,
        ]
    });
    let radiance_in: [f32; 8] = std::array::from_fn(|k| 0.1 * (k as f32));

    for &n in &normals {
        for &s in &seeds {
            for &b in &bounces {
                cases.push(NeeFrostedExteriorCase {
                    ext_normal: n.to_array(),
                    rng_seed: s,
                    bounce: b,
                    _pad0: 0.0,
                    _pad1: 0.0,
                    _pad2: 0.0,
                    lambdas,
                    stokes,
                    radiance_in,
                });
            }
        }
    }
    cases
}

/// Runs the Tier 2 check for [`nee_contribution_frosted_exterior`].
#[must_use]
pub fn run_nee_frosted_exterior(
    ctx: &crate::renderer::gpu::GpuContext,
) -> UlpCheckResult<NeeFrostedExteriorCase> {
    let map = synthetic_test_map();
    let gpu_data = HdrEnvGpuData::upload(&ctx.device, &map);
    let cases = build_nee_frosted_cases();
    let total = cases.len();

    let in_buf = compute::upload(
        &ctx.device,
        "nee frosted in",
        &cases,
        wgpu::BufferUsages::STORAGE,
    );
    let out_buf = compute::zeroed_buffer::<f32>(
        &ctx.device,
        "nee frosted out",
        total * 8,
        wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
    );
    let pipeline = compute::create_compute_pipeline(
        &ctx.device,
        "nee_frosted_exterior_main",
        SHADER_SRC,
        "nee_frosted_exterior_main",
    );
    let bind_group = compute::bind_buffers(
        &ctx.device,
        "nee frosted bind group",
        &pipeline,
        &[
            (62, &gpu_data.dist_func),
            (63, &gpu_data.dist_cdf),
            (64, &gpu_data.dist_dims),
            (65, &gpu_data.texels),
            (66, &gpu_data.dims),
            (73, &in_buf),
            (74, &out_buf),
        ],
    );
    let workgroups = (total as u32).div_ceil(64);
    compute::dispatch_and_wait(
        &ctx.device,
        &ctx.queue,
        &pipeline,
        &bind_group,
        (workgroups, 1, 1),
    );
    let gpu_out: Vec<f32> = compute::readback(&ctx.device, &ctx.queue, &out_buf, total * 8);

    let empty_planes: [GpuFacetPlane; 0] = [];
    let empty_soa = build_plane_soa(&empty_planes);
    let nee_ctx = NeeContext {
        environment: EnvironmentSource::HdrMap(&map),
        plane_soa: &empty_soa,
        enabled: true,
    };

    let mut acc = UlpAccumulator::new(
        "nee_frosted_exterior",
        NEE_FROSTED_ULP_BUDGET,
        NEE_FROSTED_ABS_FLOOR,
    );
    for (idx, case) in cases.iter().enumerate() {
        let mut cpu_rad = case.radiance_in;
        let stokes_cpu: [StokesVector; 8] = std::array::from_fn(|k| {
            StokesVector::new(
                case.stokes[k][0],
                case.stokes[k][1],
                case.stokes[k][2],
                case.stokes[k][3],
            )
        });
        nee_contribution_frosted_exterior(
            nee_ctx,
            &case.lambdas,
            Vec3::from_array(case.ext_normal),
            case.rng_seed,
            case.bounce,
            &stokes_cpu,
            &mut cpu_rad,
        );
        let base = idx * 8;
        for k in 0..8 {
            acc.record(case, "rad", cpu_rad[k], gpu_out[base + k]);
        }
    }
    acc.finish()
}

// ---------------------------------------------------------------------------------
// 6. NEE contribution: Henyey-Greenstein scattering
// ---------------------------------------------------------------------------------

#[repr(C)]
#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
pub struct NeeHgScatterCase {
    scatter_point: [f32; 3],
    n_inside_hero: f32,
    scatter_dir_in: [f32; 3],
    g: f32,
    rng_seed: u32,
    bounce: u32,
    /// `optics::materials::GemMaterial::scattering_sigma_s`, the SAME
    /// quantity `maybe_scatter_or_extinguish`'s survive branch uses.
    sigma_s: f32,
    /// `optics::materials::GemMaterial::absorption_path_scale`.
    absorption_path_scale: f32,
    /// `1` marks every one of the cube's six exit facets
    /// [`crate::optics::raytracer::FacetFinish::Frosted`] (see [`build_nee_hg_cases`]) so
    /// `nee_contribution_hg_scatter`'s frosted-exit skip fires regardless of which facet
    /// the sampled shadow ray actually hits; `0` leaves every facet
    /// [`crate::optics::raytracer::FacetFinish::Polished`] (the CPU function's default).
    frosted_exit: u32,
    _pad0: f32,
    _pad1: f32,
    _pad2: f32,
    /// `optics::absorption::channel_absorption_alphas_assigned`'s per-channel
    /// output, the same array `maybe_scatter_or_extinguish` and `nee_contribution_hg_scatter`
    /// both read.
    alphas: [f32; 8],
    lambdas: [f32; 8],
    stokes: [[f32; 4]; 8],
    radiance_in: [f32; 8],
}

const _: () = assert!(size_of::<NeeHgScatterCase>() == 288);

const NEE_HG_ULP_BUDGET: u32 = 64;
const NEE_HG_ABS_FLOOR: f32 = 1e-4;

fn build_nee_hg_cases() -> Vec<NeeHgScatterCase> {
    let mut cases = Vec::new();
    let points = [
        Vec3::ZERO,
        Vec3::new(0.2, -0.3, 0.1),
        Vec3::new(-0.4, 0.2, 0.3),
    ];
    let dirs_in = [
        Vec3::Z,
        Vec3::new(0.3, 0.9, 0.1).normalize(),
        Vec3::new(-0.7, -0.1, 0.5).normalize(),
    ];
    let gs = [-0.6f32, 0.0, 0.4, 0.8];
    let n_heroes = [1.33f32, 1.5, 1.77];
    let seeds = [42u32, 9999, 123_456];
    let bounces = [0u32, 1, 3];
    let lambdas: [f32; 8] = std::array::from_fn(|k| (k as f32).mul_add(45.0, 400.0));
    let stokes: [[f32; 4]; 8] = std::array::from_fn(|k| {
        [
            (k as f32).mul_add(0.1, 1.0),
            (k as f32) * 0.05,
            -(k as f32) * 0.03,
            0.01,
        ]
    });
    let radiance_in: [f32; 8] = std::array::from_fn(|k| 0.05 * (k as f32));

    // Findings 2b/2d: the medium-transmittance and frosted-exit-skip inputs the old
    // twin never modeled. Cycled by case index (not cross-producted with the six axes
    // above) so the case count stays the same order of magnitude while still covering
    // lossless/unit-scale/polished alongside real absorbing, scattering and frosted
    // combinations.
    let sigma_s_values = [0.0f32, 0.4, 1.2];
    let path_scale_values = [1.0f32, 1.6];
    let alpha_sets: [[f32; 8]; 3] = [
        [0.0; 8],
        std::array::from_fn(|k| (k as f32).mul_add(0.02, 0.05)),
        std::array::from_fn(|k| (k as f32).mul_add(-0.01, 0.15)),
    ];

    let mut i: usize = 0;
    for &p in &points {
        for &d in &dirs_in {
            for &g in &gs {
                for &n_hero in &n_heroes {
                    for &s in &seeds {
                        for &b in &bounces {
                            let sigma_s = sigma_s_values[i % sigma_s_values.len()];
                            let absorption_path_scale =
                                path_scale_values[i % path_scale_values.len()];
                            let alphas = alpha_sets[i % alpha_sets.len()];
                            // Every seventh case exercises the frosted-exit
                            // skip; the rest stay polished.
                            let frosted_exit = u32::from(i.is_multiple_of(7));
                            cases.push(NeeHgScatterCase {
                                scatter_point: p.to_array(),
                                n_inside_hero: n_hero,
                                scatter_dir_in: d.to_array(),
                                g,
                                rng_seed: s,
                                bounce: b,
                                sigma_s,
                                absorption_path_scale,
                                frosted_exit,
                                _pad0: 0.0,
                                _pad1: 0.0,
                                _pad2: 0.0,
                                alphas,
                                lambdas,
                                stokes,
                                radiance_in,
                            });
                            i += 1;
                        }
                    }
                }
            }
        }
    }
    cases
}

/// Runs the Tier 2 check for [`nee_contribution_hg_scatter`].
#[must_use]
pub fn run_nee_hg_scatter(
    ctx: &crate::renderer::gpu::GpuContext,
) -> UlpCheckResult<NeeHgScatterCase> {
    let map = synthetic_test_map();
    let gpu_data = HdrEnvGpuData::upload(&ctx.device, &map);
    let cases = build_nee_hg_cases();
    let total = cases.len();

    let cube_planes = [
        GpuFacetPlane::new(Vec3::X, -1.0),
        GpuFacetPlane::new(Vec3::NEG_X, -1.0),
        GpuFacetPlane::new(Vec3::Y, -1.0),
        GpuFacetPlane::new(Vec3::NEG_Y, -1.0),
        GpuFacetPlane::new(Vec3::Z, -1.0),
        GpuFacetPlane::new(Vec3::NEG_Z, -1.0),
    ];
    let planes_buf = compute::upload(
        &ctx.device,
        "tf_planes",
        &cube_planes,
        wgpu::BufferUsages::STORAGE,
    );

    let in_buf = compute::upload(
        &ctx.device,
        "nee hg in",
        &cases,
        wgpu::BufferUsages::STORAGE,
    );
    let out_buf = compute::zeroed_buffer::<f32>(
        &ctx.device,
        "nee hg out",
        total * 8,
        wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
    );
    let pipeline = compute::create_compute_pipeline(
        &ctx.device,
        "nee_hg_scatter_main",
        SHADER_SRC,
        "nee_hg_scatter_main",
    );
    let bind_group = compute::bind_buffers(
        &ctx.device,
        "nee hg bind group",
        &pipeline,
        &[
            (62, &gpu_data.dist_func),
            (63, &gpu_data.dist_cdf),
            (64, &gpu_data.dist_dims),
            (65, &gpu_data.texels),
            (66, &gpu_data.dims),
            (75, &planes_buf),
            (76, &in_buf),
            (77, &out_buf),
        ],
    );
    let workgroups = (total as u32).div_ceil(64);
    compute::dispatch_and_wait(
        &ctx.device,
        &ctx.queue,
        &pipeline,
        &bind_group,
        (workgroups, 1, 1),
    );
    let gpu_out: Vec<f32> = compute::readback(&ctx.device, &ctx.queue, &out_buf, total * 8);

    let plane_soa = build_plane_soa(&cube_planes);
    let nee_ctx = NeeContext {
        environment: EnvironmentSource::HdrMap(&map),
        plane_soa: &plane_soa,
        enabled: true,
    };

    let mut acc = UlpAccumulator::new("nee_hg_scatter", NEE_HG_ULP_BUDGET, NEE_HG_ABS_FLOOR);
    accumulate_nee_hg_cpu_results(&cases, nee_ctx, &gpu_out, &mut acc);
    acc.finish()
}

/// Runs the real CPU [`nee_contribution_hg_scatter`] for every case and records each
/// channel's CPU-vs-GPU comparison into `acc`. Split out of [`run_nee_hg_scatter`] to
/// keep that function under the house line-count limit.
fn accumulate_nee_hg_cpu_results(
    cases: &[NeeHgScatterCase],
    nee_ctx: NeeContext<'_>,
    gpu_out: &[f32],
    acc: &mut UlpAccumulator<NeeHgScatterCase>,
) {
    // Same six-facet cube `run_nee_hg_scatter`'s GPU side probes; marking
    // every facet `Frosted` (rather than trying to predict which one a given case's
    // shadow ray actually exits through) guarantees the skip fires whenever the
    // per-case `frosted_exit` flag is set, matching `frosted_exit: u32` driving the
    // WGSL twin's own skip (see `NeeHgScatterCase::frosted_exit`'s doc comment).
    let all_frosted = [FacetFinish::Frosted; 6];
    for (idx, case) in cases.iter().enumerate() {
        let mut cpu_rad = case.radiance_in;
        let stokes_cpu: [StokesVector; 8] = std::array::from_fn(|k| {
            StokesVector::new(
                case.stokes[k][0],
                case.stokes[k][1],
                case.stokes[k][2],
                case.stokes[k][3],
            )
        });
        let facet_finishes: &[FacetFinish] = if case.frosted_exit != 0 {
            &all_frosted
        } else {
            &[]
        };
        nee_contribution_hg_scatter(
            nee_ctx,
            &case.lambdas,
            case.n_inside_hero,
            Vec3::from_array(case.scatter_point),
            Vec3::from_array(case.scatter_dir_in),
            case.g,
            case.rng_seed,
            case.bounce,
            &stokes_cpu,
            &mut cpu_rad,
            &case.alphas,
            case.sigma_s,
            case.absorption_path_scale,
            facet_finishes,
        );
        let base = idx * 8;
        for k in 0..8 {
            acc.record(case, "rad", cpu_rad[k], gpu_out[base + k]);
        }
    }
}