ggplot-style WebGPU Chirp Plot with GPU Autoscale and Infinite Zoom

@compute @workgroup_size(64)
fn cs_points(@builtin(global_invocation_id) gid : vec3<u32>) {
  let i = gid.x;
  if (i >= cu.pointCount) { return; }
  let n = f32(cu.pointCount);

  // Delta (time relative to xMin) — small, full f32 precision.
  let delta = f32(i) / max(n - 1.0, 1.0) * (cu.xMax - cu.xMin);

  // Decompose phase = phase0 + dPhase to avoid f32 quantization
  // at deep zoom (see "f32 precision via phase decomposition" below).
  let phase0 = 6.28318530718 * (cu.chirpF0 * cu.xMin + 0.5 * cu.chirpK * cu.xMin * cu.xMin);
  let fInst  = cu.chirpF0 + cu.chirpK * cu.xMin;
  let dPhase = 6.28318530718 * (fInst * delta + 0.5 * cu.chirpK * delta * delta);

  // sin(a+b) = sin(a)cos(b) + cos(a)sin(b) — keeps dPhase in full
  // precision while absorbing phase0's f32 error as a constant offset.
  let y = cu.amplitude * (sin(dPhase) * cos(phase0) + cos(dPhase) * sin(phase0));

  // Store delta (not absolute t) — the vertex shader maps delta
  // directly to clip space, avoiding the f32 precision loss of
  // computing xMin + delta at deep zoom.
  points[i] = vec2<f32>(delta, y);
}

fn gridLineDistPx(base : f32, delta : f32, step : f32, pxPerWorld : f32) -> f32 {
  let baseFrac = fract(base / step);
  let p = baseFrac + delta / step;
  let frac_p = fract(p);
  let dWorld = min(frac_p, 1.0 - frac_p);
  return dWorld * pxPerWorld;
}

const overlay = document.createElement('canvas');
overlay.style.position = 'absolute';
overlay.style.left = '0';
overlay.style.top = '0';
overlay.style.width = '100%';
overlay.style.height = '100%';
overlay.style.pointerEvents = 'none';
canvas.parentElement.appendChild(overlay);
const octx = overlay.getContext('2d');

readbackBuffer.mapAsync(GPUMapMode.READ).then(() => {
  const arr = new Float32Array(readbackBuffer.getMappedRange(0, 8));
  const min = arr[0], max = arr[1];
  readbackBuffer.unmap();
  autoscalePending = false;
  if (min <= max && isFinite(min) && isFinite(max)) {
    const range = max - min;
    const margin = range > 0 ? range * 0.1 : Math.abs(max) * 0.1 + 0.1;
    yMin = min - margin;
    yMax = max + margin;
  }
});

const worldX = cssXToWorldX(cssX);
const factor = e.deltaY > 0 ? 1.15 : 1 / 1.15;
xMin = worldX - (worldX - xMin) * factor;
xMax = worldX + (xMax - worldX) * factor;

<div style="position:relative;width:100%;height:450px;">
  <canvas class="chirp-plot-canvas"
          style="width:100%;height:100%;display:block;"></canvas>
</div>
<script src="/scripts/chirp-plot-webgpu.js" defer></script>

// SPDX-FileCopyrightText: 2026 Uli Köhler <gitlab@techoverflow.net>
// SPDX-License-Identifier: CC0-1.0
//
// ggplot-style WebGPU line plot of a chirp signal with GPU-side
// function evaluation, MSAA anti-aliasing, dynamic grid, autoscale,
// and infinite zoom/pan.
//
// The chirp is y = A·sin(2π·(f₀·t + ½·k·t²)) — a sine whose frequency
// increases linearly with t. The function is evaluated on the GPU at
// one sample per pixel (× oversampling), so zooming in reveals more
// detail without increasing the sample count.
//
// Auto-initializes every <canvas class="chirp-plot-canvas"> on the page.
(function () {
  'use strict';

  // ─────────────────────────────────────────────────────────────────────
  // Configuration
  // ─────────────────────────────────────────────────────────────────────

  // Chirp parameters: y = A·sin(2π·(f₀·t + ½·k·t²))
  // Frequency at time t: f(t) = f₀ + k·t
  // f₀=1, k=10 → frequency goes from 1 Hz to 101 Hz over t∈[0,10].
  const CHIRP_F0 = 1.0;
  const CHIRP_K = 10.0;
  const AMPLITUDE = 1.0;

  // Initial X viewport
  const X_MIN_INIT = 0.0;
  const X_MAX_INIT = 10.0;

  // Sampling: one point per pixel × oversampling factor.
  // This gives infinite zoom — the sample density is always proportional
  // to the canvas width, not the viewport range.
  const OVERSAMPLE = 2;
  const MAX_POINTS = 16384;

  // MSAA sample counts to probe (WebGPU guarantees only 1× and 4×).
  const MSAA_SAMPLE_COUNTS = [1, 2, 4, 8, 16];

  // Plot margins in CSS pixels (for axis labels).
  const MARGIN_LEFT = 64;
  const MARGIN_RIGHT = 16;
  const MARGIN_TOP = 16;
  const MARGIN_BOTTOM = 48;

  // ─────────────────────────────────────────────────────────────────────
  // WGSL shader
  // ─────────────────────────────────────────────────────────────────────
  // Three pipelines share this module:
  //   1. compute  — fills point buffer with chirp values
  //   2. line     — renders the chirp as an MSAA line strip
  //   3. grid     — renders the plot-area background with grid lines
  //
  // The compute shader evaluates the chirp at evenly-spaced X positions
  // across the current viewport. The point count is proportional to the
  // canvas width (not the viewport range), so zooming in reveals more
  // detail without increasing the sample count.

  const SHADER = /* wgsl */ `
  // ── Compute: fill points[] with chirp values ────────────────────────
  struct ComputeUniforms {
    pointCount : u32,
    _pad0      : u32,
    xMin       : f32,
    xMax       : f32,
    chirpF0    : f32,
    chirpK     : f32,
    amplitude  : f32,
    _pad1      : f32,
  };

  @group(0) @binding(0) var<uniform>  cu : ComputeUniforms;
  @group(0) @binding(1) var<storage, read_write> points : array<vec2<f32>>;

  @compute @workgroup_size(64)
  fn cs_points(@builtin(global_invocation_id) gid : vec3<u32>) {
    let i = gid.x;
    if (i >= cu.pointCount) { return; }
    let n = f32(cu.pointCount);

    // Delta (time relative to xMin) — small, full f32 precision.
    let delta = f32(i) / max(n - 1.0, 1.0) * (cu.xMax - cu.xMin);

    // Decompose phase = phase0 + dPhase to avoid f32 quantization at
    // deep zoom.  When |xMin| is large and |xMax - xMin| is small,
    // computing t = mix(xMin, xMax, ...) in f32 loses precision: the
    // step between consecutive samples can be smaller than the f32 ulp
    // at xMin, so multiple samples get the same t and the same phase,
    // producing a visibly quantized sine.
    //
    // phase0 is the constant phase at xMin (large, same for all samples).
    // dPhase is the per-sample increment, computed from small numbers
    // (delta, instantaneous frequency at xMin) in full f32 precision.
    //
    // sin(phase0 + dPhase) via the addition formula:
    //   sin(dPhase)*cos(phase0) + cos(dPhase)*sin(phase0)
    // sin/cos(phase0) lose precision but the error is identical for
    // every sample — a constant phase offset, not quantization.
    // sin/cos(dPhase) are precise (small argument).
    let phase0 = 6.28318530718 * (cu.chirpF0 * cu.xMin + 0.5 * cu.chirpK * cu.xMin * cu.xMin);
    let fInst  = cu.chirpF0 + cu.chirpK * cu.xMin;
    let dPhase = 6.28318530718 * (fInst * delta + 0.5 * cu.chirpK * delta * delta);

    let y = cu.amplitude * (sin(dPhase) * cos(phase0) + cos(dPhase) * sin(phase0));

    // Store delta (not absolute t) in the point buffer.  The line vertex
    // shader maps delta directly to clip space without computing
    // xMin + delta, which would lose f32 precision at deep zoom (the
    // per-sample delta can be smaller than the f32 ulp at xMin).
    points[i] = vec2<f32>(delta, y);
  }

  // ── Compute: reduce min/max Y for autoscale ─────────────────────────
  // Each workgroup computes a local min/max over 64 points, then
  // workgroup 0 reduces all workgroup results into a single vec2<f32>
  // (minY, maxY) stored in the reduce buffer.
  //
  // This is a two-pass reduce for simplicity: the first pass writes
  // per-workgroup min/max to a workgroup buffer; the second pass (only
  // workgroup 0) reads all workgroup results and writes the final
  // min/max. For MAX_POINTS=16384, that's at most 256 workgroups, so
  // the second pass is a single workgroup of 64 threads reading 256
  // values — well within the 256-element workgroup uniform limit.

  struct ReduceUniforms {
    pointCount : u32,
    workgroupCount : u32,
    _pad0 : u32,
    _pad1 : u32,
  };

  @group(0) @binding(0) var<uniform>  rdu : ReduceUniforms;
  @group(0) @binding(1) var<storage, read> rdpoints : array<vec2<f32>>;
  @group(0) @binding(2) var<storage, read_write> wgMinMax : array<vec2<f32>>;
  @group(0) @binding(3) var<storage, read_write> finalMinMax : vec2<f32>;

  var<workgroup> wgMin : f32;
  var<workgroup> wgMax : f32;

  @compute @workgroup_size(64)
  fn cs_reduce(@builtin(workgroup_id) wid : vec3<u32>,
               @builtin(local_invocation_id) lid : vec3<u32>,
               @builtin(local_invocation_index) lidx : u32) {
    // Pass 1: each workgroup computes its local min/max.
    if (lidx == 0) { wgMin = 1e30; wgMax = -1e30; }
    workgroupBarrier();

    let i = wid.x * 64 + lidx;
    if (i < rdu.pointCount) {
      let y = rdpoints[i].y;
      // Atomic-free reduction: use workgroup shared variables with
      // a simple serial reduction within the workgroup.
      // For 64 threads this is fast enough.
    }

    // Each thread loads its Y value (if valid) into a local variable.
    let myY = select(rdpoints[i].y, 0.0, i >= rdu.pointCount);
    let valid = i < rdu.pointCount;

    // Serial reduction within workgroup (64 threads → simple loop).
    // We use a single-thread accumulation via workgroupBarrier + shared var.
    if (lidx == 0) {
      for (var j = 0u; j < 64u; j++) {
        let idx = wid.x * 64 + j;
        if (idx < rdu.pointCount) {
          let yv = rdpoints[idx].y;
          if (yv < wgMin) { wgMin = yv; }
          if (yv > wgMax) { wgMax = yv; }
        }
      }
      wgMinMax[wid.x] = vec2<f32>(wgMin, wgMax);
    }

    workgroupBarrier();

    // Pass 2: workgroup 0 reduces all workgroup results.
    if (wid.x == 0 && lidx == 0) {
      var gMin = 1e30;
      var gMax = -1e30;
      for (var w = 0u; w < rdu.workgroupCount; w++) {
        let mm = wgMinMax[w];
        if (mm.x < gMin) { gMin = mm.x; }
        if (mm.y > gMax) { gMax = mm.y; }
      }
      finalMinMax = vec2<f32>(gMin, gMax);
    }
  }

  // ── Shared render uniforms (line + grid pipelines) ──────────────────
  // Maps world coordinates → canvas pixels → clip space.
  struct RenderUniforms {
    xMin       : f32,
    xMax       : f32,
    yMin       : f32,
    yMax       : f32,
    plotX      : f32,  // plot area origin X in canvas pixels
    plotY      : f32,  // plot area origin Y in canvas pixels
    plotW      : f32,  // plot area width in canvas pixels
    plotH      : f32,  // plot area height in canvas pixels
    canvasW    : f32,
    canvasH    : f32,
    xMajorStep : f32,
    yMajorStep : f32,
  };

  @group(0) @binding(0) var<uniform> ru : RenderUniforms;

  // World → clip space helper (used by grid vertex shader only).
  fn worldToClip(worldX : f32, worldY : f32) -> vec4<f32> {
    let plotPx = (worldX - ru.xMin) / (ru.xMax - ru.xMin) * ru.plotW;
    let plotPy = (ru.yMax - worldY) / (ru.yMax - ru.yMin) * ru.plotH;
    let canvasPx = ru.plotX + plotPx;
    let canvasPy = ru.plotY + plotPy;
    return vec4<f32>(
      (canvasPx / ru.canvasW) * 2.0 - 1.0,
      1.0 - (canvasPy / ru.canvasH) * 2.0,
      0.0, 1.0
    );
  }

  // Delta → clip space (used by line vertex shader).
  // delta is the time relative to xMin, stored in the point buffer
  // instead of absolute t to avoid f32 precision loss at deep zoom.
  // Maps delta directly to plot pixels without computing xMin + delta.
  fn deltaToClip(delta : f32, worldY : f32) -> vec4<f32> {
    let plotPx = delta / (ru.xMax - ru.xMin) * ru.plotW;
    let plotPy = (ru.yMax - worldY) / (ru.yMax - ru.yMin) * ru.plotH;
    let canvasPx = ru.plotX + plotPx;
    let canvasPy = ru.plotY + plotPy;
    return vec4<f32>(
      (canvasPx / ru.canvasW) * 2.0 - 1.0,
      1.0 - (canvasPy / ru.canvasH) * 2.0,
      0.0, 1.0
    );
  }

  // ── Line strip pipeline ─────────────────────────────────────────────
  @group(0) @binding(1) var<storage, read> rpoints : array<vec2<f32>>;

  @vertex
  fn vs_line(@builtin(vertex_index) vi : u32) -> @builtin(position) vec4<f32> {
    let p = rpoints[vi];
    // p.x is delta (time relative to xMin), not absolute t.
    return deltaToClip(p.x, p.y);
  }

  @fragment
  fn fs_line() -> @location(0) vec4<f32> {
    return vec4<f32>(0.0, 0.45, 0.75, 1.0);  // ggplot-like blue
  }

  // ── Grid pipeline ───────────────────────────────────────────────────
  // Renders a fullscreen quad. The fragment shader computes:
  //   1. Whether the pixel is inside the plot area or in the margin.
  //   2. If inside, anti-aliased major/minor grid lines using fwidth().
  //
  // The grid lines are computed in world coordinates — the fragment
  // shader converts pixel position back to world X/Y, then uses fract()
  // and fwidth() to draw 1px-wide anti-aliased lines at each tick step.
  // This is the same technique as the mexhat floor grid, but in 2D.

  struct GridVSOut {
    @builtin(position) clipPos : vec4<f32>,
    @location(0)       canvasPx : vec2<f32>,  // pixel coords [0,cw]×[0,ch]
  };

  @vertex
  fn vs_grid(@builtin(vertex_index) vi : u32) -> GridVSOut {
    var pos = array<vec2<f32>, 6>(
      vec2<f32>(-1.0, -1.0),
      vec2<f32>( 1.0, -1.0),
      vec2<f32>( 1.0,  1.0),
      vec2<f32>(-1.0, -1.0),
      vec2<f32>( 1.0,  1.0),
      vec2<f32>(-1.0,  1.0),
    );
    var out : GridVSOut;
    out.clipPos = vec4<f32>(pos[vi], 0.0, 1.0);
    // Convert clip space → canvas pixels for the fragment shader.
    out.canvasPx = vec2<f32>(
      (pos[vi].x * 0.5 + 0.5) * ru.canvasW,
      (1.0 - pos[vi].y * 0.5 - 0.5) * ru.canvasH,
    );
    return out;
  }

  // Distance (in pixels) from base+delta to the nearest grid line at
  // multiples of step.  base is the viewport origin (xMin or yMin),
  // delta is the per-pixel offset from that origin (small, full f32
  // precision).  Decomposes fract((base+delta)/step) into
  // fract(fract(base/step) + delta/step) so the large base value
  // doesn't quantize the per-pixel delta.
  fn gridLineDistPx(base : f32, delta : f32, step : f32, pxPerWorld : f32) -> f32 {
    let baseFrac = fract(base / step);
    let p = baseFrac + delta / step;
    let frac_p = fract(p);
    let dWorld = min(frac_p, 1.0 - frac_p);
    return dWorld * pxPerWorld;
  }

  // Grid line intensity from a pixel distance: 1.0 on the line,
  // anti-aliased to 0 over ~1px. width is the line half-width in px.
  fn gridLineIntensity(distPx : f32, width : f32) -> f32 {
    return 1.0 - smoothstep(0.0, max(width, 0.0001), distPx);
  }

  @fragment
  fn fs_grid(in : GridVSOut) -> @location(0) vec4<f32> {
    // Margin color (outside plot area).
    let marginColor = vec3<f32>(0.94, 0.94, 0.94);
    // Plot background (ggplot-style off-white).
    let plotBg = vec3<f32>(0.97, 0.97, 0.97);

    // Check if inside plot area.
    let px = in.canvasPx.x;
    let py = in.canvasPx.y;
    let inside = px >= ru.plotX && px < ru.plotX + ru.plotW
              && py >= ru.plotY && py < ru.plotY + ru.plotH;

    if (!inside) {
      return vec4<f32>(marginColor, 1.0);
    }

    // Delta from viewport origin (small, full f32 precision).
    // Avoids computing xMin + delta which loses f32 precision at deep zoom.
    let deltaX = (px - ru.plotX) / ru.plotW * (ru.xMax - ru.xMin);
    let deltaY = (py - ru.plotY) / ru.plotH * (ru.yMax - ru.yMin);

    // World → pixel scale for each axis.
    let pxPerWorldX = ru.plotW / (ru.xMax - ru.xMin);
    let pxPerWorldY = ru.plotH / (ru.yMax - ru.yMin);

    // Distance to the nearest minor/major grid line (horizontal or
    // vertical) in pixels. min() across axes gives the closest line
    // of either orientation.
    let minorD = min(gridLineDistPx(ru.xMin, deltaX, ru.xMajorStep * 0.5, pxPerWorldX),
                     gridLineDistPx(ru.yMin, deltaY, ru.yMajorStep * 0.5, pxPerWorldY));
    let majorD = min(gridLineDistPx(ru.xMin, deltaX, ru.xMajorStep, pxPerWorldX),
                     gridLineDistPx(ru.yMin, deltaY, ru.yMajorStep, pxPerWorldY));

    // ~1px wide lines, anti-aliased.
    let xMinor = gridLineIntensity(minorD, 1.0);
    let xMajor = gridLineIntensity(majorD, 1.0);

    let minorColor = vec3<f32>(0.88, 0.88, 0.88);
    let majorColor = vec3<f32>(0.75, 0.75, 0.78);

    // Blend: minor first, then major on top.
    var color = plotBg;
    color = mix(color, minorColor, xMinor * 0.5);
    color = mix(color, majorColor, xMajor * 0.7);

    return vec4<f32>(color, 1.0);
  }
`;

  // ─────────────────────────────────────────────────────────────────────
  // Nice tick step algorithm (ggplot-style)
  // ─────────────────────────────────────────────────────────────────────
  // Given a data range [min, max] and a target number of ticks,
  // returns a "nice" step size (1, 2, 5 × 10^n).

  function niceStep(range, targetCount) {
    const raw = range / targetCount;
    const mag = Math.pow(10, Math.floor(Math.log10(raw)));
    const norm = raw / mag;
    let step;
    if (norm <= 1.5) step = 1;
    else if (norm <= 3) step = 2;
    else if (norm <= 7) step = 5;
    else step = 10;
    return step * mag;
  }

  // ─────────────────────────────────────────────────────────────────────
  // Per-canvas initialization
  // ─────────────────────────────────────────────────────────────────────

  async function initCanvas(canvas) {
    if (!navigator.gpu) {
      canvas.replaceWith(document.createTextNode('WebGPU is not supported in this browser.'));
      return;
    }
    const adapter = await navigator.gpu.requestAdapter();
    if (!adapter) {
      canvas.replaceWith(document.createTextNode('No WebGPU adapter available.'));
      return;
    }
    const device = await adapter.requestDevice();
    const ctx = canvas.getContext('webgpu');
    const format = navigator.gpu.getPreferredCanvasFormat();
    ctx.configure({ device, format, alphaMode: 'premultiplied' });

    const module = device.createShaderModule({ code: SHADER });

    // --- Bind group layouts ---
    const computeLayout = device.createBindGroupLayout({
      entries: [
        { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform' } },
        { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
      ],
    });

    const renderLayout = device.createBindGroupLayout({
      entries: [
        { binding: 0, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: 'uniform' } },
        { binding: 1, visibility: GPUShaderStage.VERTEX, buffer: { type: 'read-only-storage' } },
      ],
    });

    // --- Compute pipeline (chirp point generation) ---
    const computePipeline = device.createComputePipeline({
      layout: device.createPipelineLayout({ bindGroupLayouts: [computeLayout] }),
      compute: { module, entryPoint: 'cs_points' },
    });

    // --- Reduce pipeline (autoscale min/max) ---
    // Computes min/max Y of the chirp points on the GPU. The result is
    // read back to the CPU asynchronously to adjust the Y viewport.
    const reduceLayout = device.createBindGroupLayout({
      entries: [
        { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform' } },
        { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
        { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
        { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
      ],
    });
    const reducePipeline = device.createComputePipeline({
      layout: device.createPipelineLayout({ bindGroupLayouts: [reduceLayout] }),
      compute: { module, entryPoint: 'cs_reduce' },
    });

    // Reduce buffers.
    // wgMinMax: one vec2 per workgroup (max 256 workgroups for 16384 points).
    // finalMinMax: single vec2 (minY, maxY) — the reduce result.
    // readbackBuffer: copy target for finalMinMax, mapped for CPU read.
    const maxWorkgroups = Math.ceil(MAX_POINTS / 64);
    const reduceUniforms = device.createBuffer({
      size: 16,
      usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
    });
    const wgMinMaxBuffer = device.createBuffer({
      size: maxWorkgroups * 8,
      usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
    });
    const finalMinMaxBuffer = device.createBuffer({
      size: 8,
      usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
    });
    const readbackBuffer = device.createBuffer({
      size: 8,
      usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
    });
    // reduceBindGroup is created after pointBuffer (see below).

    // --- MSAA probe + render pipelines ---
    // Not all sample counts are supported by every adapter. The WebGPU
    // spec only guarantees 1× and 4×; 2×, 8×, and 16× are
    // adapter-dependent. We probe by creating the pipeline with
    // createRenderPipelineAsync — if the sample count is unsupported,
    // the Promise rejects and we skip that count.
    const supportedSC = [];
    const linePipelines = [];
    const gridPipelines = [];
    const lineLayout = device.createPipelineLayout({ bindGroupLayouts: [renderLayout] });
    const gridLayout = device.createPipelineLayout({ bindGroupLayouts: [renderLayout] });
    for (const sc of MSAA_SAMPLE_COUNTS) {
      try {
        const [linePipe, gridPipe] = await Promise.all([
          device.createRenderPipelineAsync({
            layout: lineLayout,
            vertex:   { module, entryPoint: 'vs_line' },
            fragment: { module, entryPoint: 'fs_line', targets: [{ format }] },
            primitive: { topology: 'line-strip' },
            multisample: { count: sc },
          }),
          device.createRenderPipelineAsync({
            layout: gridLayout,
            vertex:   { module, entryPoint: 'vs_grid' },
            fragment: { module, entryPoint: 'fs_grid', targets: [{ format }] },
            primitive: { topology: 'triangle-list' },
            multisample: { count: sc },
          }),
        ]);
        supportedSC.push(sc);
        linePipelines.push(linePipe);
        gridPipelines.push(gridPipe);
      } catch {
        // Unsupported sample count — skip.
      }
    }
    // Use 4× MSAA if available, otherwise the highest supported count
    // below 4 (avoids the steeper cost of 8×/16× on high-DPI displays).
    const TARGET_SC = 4;
    let msaaIndex = supportedSC.length - 1;
    const targetIdx = supportedSC.indexOf(TARGET_SC);
    if (targetIdx !== -1) msaaIndex = targetIdx;

    // --- Buffers ---
    const computeUniforms = device.createBuffer({
      size: 32,
      usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
    });

    const renderUniforms = device.createBuffer({
      size: 48,
      usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
    });

    const pointBuffer = device.createBuffer({
      size: MAX_POINTS * 8,
      usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
    });

    const computeBindGroup = device.createBindGroup({
      layout: computeLayout,
      entries: [
        { binding: 0, resource: { buffer: computeUniforms } },
        { binding: 1, resource: { buffer: pointBuffer } },
      ],
    });

    const lineBindGroup = device.createBindGroup({
      layout: renderLayout,
      entries: [
        { binding: 0, resource: { buffer: renderUniforms } },
        { binding: 1, resource: { buffer: pointBuffer } },
      ],
    });

    const gridBindGroup = device.createBindGroup({
      layout: renderLayout,
      entries: [
        { binding: 0, resource: { buffer: renderUniforms } },
        { binding: 1, resource: { buffer: pointBuffer } },
      ],
    });

    // Reduce bind group (needs pointBuffer, so created here after it).
    const reduceBindGroup = device.createBindGroup({
      layout: reduceLayout,
      entries: [
        { binding: 0, resource: { buffer: reduceUniforms } },
        { binding: 1, resource: { buffer: pointBuffer } },
        { binding: 2, resource: { buffer: wgMinMaxBuffer } },
        { binding: 3, resource: { buffer: finalMinMaxBuffer } },
      ],
    });

    // --- MSAA texture ---
    let msaaTexture = null;
    let msaaSC = 0, msaaW = 0, msaaH = 0;
    function ensureMsaa(w, h, sc) {
      if (sc === 1) {
        if (msaaTexture) { msaaTexture.destroy(); msaaTexture = null; }
        msaaSC = 1; msaaW = w; msaaH = h;
        return;
      }
      if (msaaTexture && msaaSC === sc && msaaW === w && msaaH === h) return;
      if (msaaTexture) msaaTexture.destroy();
      msaaTexture = device.createTexture({
        size: [w, h, 1], format,
        usage: GPUTextureUsage.RENDER_ATTACHMENT,
        sampleCount: sc,
      });
      msaaSC = sc; msaaW = w; msaaH = h;
    }

    // --- 2D canvas overlay for axis labels ---
    // A sibling <canvas> positioned on top of the WebGPU canvas.
    // We draw axis tick labels and plot border using the Canvas 2D API,
    // which is simpler than rendering text in WebGPU and sharp at any DPR.
    const overlay = document.createElement('canvas');
    overlay.style.position = 'absolute';
    overlay.style.left = '0';
    overlay.style.top = '0';
    overlay.style.width = '100%';
    overlay.style.height = '100%';
    overlay.style.pointerEvents = 'none';
    // The WebGPU canvas must be position:relative for the overlay to anchor.
    canvas.style.position = canvas.style.position || 'relative';
    canvas.parentElement.appendChild(overlay);
    const octx = overlay.getContext('2d');

    // --- Viewport state ---
    let xMin = X_MIN_INIT, xMax = X_MAX_INIT;
    let yMin = -AMPLITUDE * 1.1, yMax = AMPLITUDE * 1.1;

    // --- Interaction: scroll to zoom, drag to pan ---
    // Scroll: zoom toward cursor position (the world point under the
    //   cursor stays fixed). Zoom only affects X; Y is controlled by
    //   autoscale (or by the user if they drag vertically).
    // Drag: pan in both X and Y. Dragging vertically disables autoscale
    //   (user is manually controlling Y). Double-click re-enables autoscale.
    let dragging = false, lastX = 0, lastY = 0;

    // Convert canvas CSS pixel X to world X.
    function cssXToWorldX(cssX) {
      const dpr = Math.min(window.devicePixelRatio || 1, 2);
      const ml = MARGIN_LEFT * dpr;
      const mr = MARGIN_RIGHT * dpr;
      const pw = canvas.width - ml - mr;
      const px = cssX * dpr - ml;
      return xMin + px / pw * (xMax - xMin);
    }

    canvas.style.cursor = 'crosshair';
    canvas.style.touchAction = 'none';

    canvas.addEventListener('pointerdown', (e) => {
      dragging = true; lastX = e.clientX; lastY = e.clientY;
      canvas.setPointerCapture(e.pointerId);
      canvas.style.cursor = 'grabbing';
    });
    canvas.addEventListener('pointerup', () => {
      dragging = false; canvas.style.cursor = 'crosshair';
    });
    canvas.addEventListener('pointermove', (e) => {
      if (!dragging) return;
      const dpr = Math.min(window.devicePixelRatio || 1, 2);
      const ml = MARGIN_LEFT * dpr, mr = MARGIN_RIGHT * dpr;
      const mt = MARGIN_TOP * dpr, mb = MARGIN_BOTTOM * dpr;
      const pw = canvas.width - ml - mr;
      const ph = canvas.height - mt - mb;
      const dx = (e.clientX - lastX) * dpr;
      const dy = (e.clientY - lastY) * dpr;
      lastX = e.clientX; lastY = e.clientY;

      // Pan X: convert pixel delta to world delta.
      const xRange = xMax - xMin;
      const xDelta = -dx / pw * xRange;
      xMin += xDelta; xMax += xDelta;

      // Pan Y: convert pixel delta to world delta (inverted — drag down = up).
      const yRange = yMax - yMin;
      const yDelta = dy / ph * yRange;
      yMin += yDelta; yMax += yDelta;

      // User is manually controlling Y → disable autoscale.
      autoscaleEnabled = false;
    });

    canvas.addEventListener('wheel', (e) => {
      e.preventDefault();
      const dpr = Math.min(window.devicePixelRatio || 1, 2);
      const ml = MARGIN_LEFT * dpr, mr = MARGIN_RIGHT * dpr;
      const pw = canvas.width - ml - mr;
      const rect = canvas.getBoundingClientRect();
      const cssX = e.clientX - rect.left;
      const worldX = cssXToWorldX(cssX);
      const factor = e.deltaY > 0 ? 1.15 : 1 / 1.15;
      // Zoom toward cursor: keep worldX fixed.
      xMin = worldX - (worldX - xMin) * factor;
      xMax = worldX + (xMax - worldX) * factor;
      // Re-enable autoscale on zoom (Y adapts to new X range).
      autoscaleEnabled = true;
    }, { passive: false });

    // Double-click: reset to initial view + re-enable autoscale.
    canvas.addEventListener('dblclick', () => {
      xMin = X_MIN_INIT; xMax = X_MAX_INIT;
      yMin = -AMPLITUDE * 1.1; yMax = AMPLITUDE * 1.1;
      autoscaleEnabled = true;
    });

    // --- Autoscale state ---
    // The GPU reduce runs asynchronously. We submit the reduce compute
    // + copy, then map the readback buffer. When the map completes, we
    // update yMin/yMax with a 10% margin. A new reduce is submitted
    // only when the previous readback completes (no queuing).
    let autoscalePending = false;
    let autoscaleEnabled = true;

    function submitAutoscale(pointCount) {
      if (autoscalePending || !autoscaleEnabled) return;
      autoscalePending = true;

      const wgCount = Math.ceil(pointCount / 64);
      const rduData = new ArrayBuffer(16);
      const rduView = new DataView(rduData);
      rduView.setUint32(0, pointCount, true);
      rduView.setUint32(4, wgCount, true);
      rduView.setUint32(8, 0, true);
      rduView.setUint32(12, 0, true);
      device.queue.writeBuffer(reduceUniforms, 0, rduData);

      const enc = device.createCommandEncoder();
      const rp = enc.beginComputePass();
      rp.setPipeline(reducePipeline);
      rp.setBindGroup(0, reduceBindGroup);
      rp.dispatchWorkgroups(wgCount);
      rp.end();
      enc.copyBufferToBuffer(finalMinMaxBuffer, 0, readbackBuffer, 0, 8);
      device.queue.submit([enc.finish()]);

      // Map for read. When complete, update yMin/yMax.
      readbackBuffer.mapAsync(GPUMapMode.READ).then(() => {
        const arr = new Float32Array(readbackBuffer.getMappedRange(0, 8));
        const min = arr[0], max = arr[1];
        readbackBuffer.unmap();
        autoscalePending = false;
        if (min <= max && isFinite(min) && isFinite(max)) {
          const range = max - min;
          const margin = range > 0 ? range * 0.1 : Math.abs(max) * 0.1 + 0.1;
          yMin = min - margin;
          yMax = max + margin;
        }
      }).catch(() => { autoscalePending = false; });
    }

    // --- Axis label drawing (Canvas 2D overlay) ---
    // Draws tick labels and plot border on the 2D overlay canvas.
    // Called every frame after the WebGPU render.
    function drawAxes(cssW, cssH, dpr) {
      const cw = Math.round(cssW * dpr);
      const ch = Math.round(cssH * dpr);
      if (overlay.width !== cw || overlay.height !== ch) {
        overlay.width = cw; overlay.height = ch;
      }
      octx.clearRect(0, 0, cw, ch);
      octx.save();
      octx.scale(dpr, dpr);

      const ml = MARGIN_LEFT, mr = MARGIN_RIGHT, mt = MARGIN_TOP, mb = MARGIN_BOTTOM;
      const px = ml, py = mt;
      const pw = cssW - ml - mr, ph = cssH - mt - mb;

      // Plot axes (ggplot style: only bottom and left axis lines,
      // no top/right border, with small outward-pointing tick marks).
      octx.strokeStyle = '#333';
      octx.lineWidth = 1;
      octx.beginPath();
      // Left axis.
      octx.moveTo(px, py);
      octx.lineTo(px, py + ph);
      // Bottom axis.
      octx.lineTo(px + pw, py + ph);
      octx.stroke();

      // Tick marks.
      octx.font = '11px sans-serif';
      octx.fillStyle = '#333';

      // X axis tick labels (bottom) — relative time from xMin,
      // auto-scaled to s / ms / µs / ns based on the step size.
      const xStep = niceStep(xMax - xMin, 8);
      const xStart = Math.ceil(xMin / xStep) * xStep;
      octx.textAlign = 'center';
      octx.textBaseline = 'top';
      for (let x = xStart; x <= xMax + xStep * 0.001; x += xStep) {
        const sx = px + (x - xMin) / (xMax - xMin) * pw;
        // Clip to plot area to avoid labels spilling outside.
        if (sx < px - 1 || sx > px + pw + 1) continue;
        // Outward tick mark.
        octx.beginPath();
        octx.moveTo(sx, py + ph);
        octx.lineTo(sx, py + ph + 4);
        octx.stroke();
        octx.fillText(formatTimeTick(x - xMin, xStep), sx, py + ph + 8);
      }

      // Y axis tick labels (left).
      const yStep = niceStep(yMax - yMin, 6);
      const yStart = Math.ceil(yMin / yStep) * yStep;
      octx.textAlign = 'right';
      octx.textBaseline = 'middle';
      for (let y = yStart; y <= yMax + yStep * 0.001; y += yStep) {
        const sy = py + (yMax - y) / (yMax - yMin) * ph;
        if (sy < py - 1 || sy > py + ph + 1) continue;
        // Outward tick mark.
        octx.beginPath();
        octx.moveTo(px, sy);
        octx.lineTo(px - 4, sy);
        octx.stroke();
        octx.fillText(formatTick(y), px - 8, sy);
      }

      // Axis titles.
      octx.font = '13px sans-serif';
      octx.fillStyle = '#222';
      // X title (centered below ticks).
      octx.textAlign = 'center';
      octx.textBaseline = 'bottom';
      octx.fillText('t (relative)', ml + pw / 2, cssH - 4);
      // Y title (rotated, left of ticks).
      octx.save();
      octx.translate(12, mt + ph / 2);
      octx.rotate(-Math.PI / 2);
      octx.textAlign = 'center';
      octx.textBaseline = 'top';
      octx.fillText('y', 0, 0);
      octx.restore();

      octx.restore();
    }

    // Format a tick value: show enough decimals for the step size.
    function formatTick(v) {
      const abs = Math.abs(v);
      if (abs === 0) return '0';
      if (abs >= 1000 || abs < 1e-3) return v.toExponential(1);
      // Fixed notation with up to 3 decimals, trailing zeros stripped.
      let s = v.toFixed(3);
      // Strip trailing zeros and a trailing decimal point.
      s = s.replace(/0+$/, '').replace(/\.$/, '');
      return s;
    }

    // Format a time tick (in seconds) with automatic unit selection.
    // Picks the largest unit where the step >= 1 in that unit, so
    // labels stay short (e.g. "1 ms" instead of "0.001 s").
    function formatTimeTick(delta, step) {
      let unit, scale;
      if (step >= 1)        { unit = 's';  scale = 1; }
      else if (step >= 1e-3) { unit = 'ms'; scale = 1e3; }
      else if (step >= 1e-6) { unit = 'µs'; scale = 1e6; }
      else                   { unit = 'ns'; scale = 1e9; }

      const stepInUnit = step * scale;
      let decimals = 0;
      if (stepInUnit < 1)   decimals = 1;
      if (stepInUnit < 0.1) decimals = 2;
      if (stepInUnit < 0.01) decimals = 3;

      let s = (delta * scale).toFixed(decimals);
      if (decimals > 0) s = s.replace(/0+$/, '').replace(/\.$/, '');
      return s + ' ' + unit;
    }

    // --- Render loop ---
    function frame() {
      const dpr = Math.min(window.devicePixelRatio || 1, 2);
      const cw = Math.max(1, Math.floor(canvas.clientWidth * dpr));
      const ch = Math.max(1, Math.floor(canvas.clientHeight * dpr));
      if (canvas.width !== cw || canvas.height !== ch) {
        canvas.width = cw; canvas.height = ch;
      }

      // Plot area in canvas pixels (account for DPR-scaled margins).
      const ml = MARGIN_LEFT * dpr;
      const mr = MARGIN_RIGHT * dpr;
      const mt = MARGIN_TOP * dpr;
      const mb = MARGIN_BOTTOM * dpr;
      const plotX = ml;
      const plotY = mt;
      const plotW = Math.max(1, cw - ml - mr);
      const plotH = Math.max(1, ch - mt - mb);

      // Point count: proportional to canvas width (infinite zoom).
      const pointCount = Math.min(MAX_POINTS, Math.max(2, Math.floor(plotW * OVERSAMPLE)));
      const sc = supportedSC[msaaIndex];

      // Compute nice tick steps.
      const xRange = xMax - xMin;
      const yRange = yMax - yMin;
      const xMajorStep = niceStep(xRange, 8);
      const yMajorStep = niceStep(yRange, 6);

      // Write compute uniforms.
      const cuData = new ArrayBuffer(32);
      const cuView = new DataView(cuData);
      cuView.setUint32(0, pointCount, true);
      cuView.setUint32(4, 0, true);
      cuView.setFloat32(8, xMin, true);
      cuView.setFloat32(12, xMax, true);
      cuView.setFloat32(16, CHIRP_F0, true);
      cuView.setFloat32(20, CHIRP_K, true);
      cuView.setFloat32(24, AMPLITUDE, true);
      cuView.setFloat32(28, 0, true);
      device.queue.writeBuffer(computeUniforms, 0, cuData);

      // Write render uniforms.
      const ruData = new Float32Array(12);
      ruData[0] = xMin; ruData[1] = xMax;
      ruData[2] = yMin; ruData[3] = yMax;
      ruData[4] = plotX; ruData[5] = plotY;
      ruData[6] = plotW; ruData[7] = plotH;
      ruData[8] = cw; ruData[9] = ch;
      ruData[10] = xMajorStep; ruData[11] = yMajorStep;
      device.queue.writeBuffer(renderUniforms, 0, ruData);

      ensureMsaa(cw, ch, sc);

      const encoder = device.createCommandEncoder();

      // Compute pass — generate chirp points on the GPU.
      const cp = encoder.beginComputePass();
      cp.setPipeline(computePipeline);
      cp.setBindGroup(0, computeBindGroup);
      cp.dispatchWorkgroups(Math.ceil(pointCount / 64));
      cp.end();

      // Render pass — grid (background) + line strip (chirp).
      const canvasView = ctx.getCurrentTexture().createView();
      const colorAttachment = {
        clearValue: { r: 0.94, g: 0.94, b: 0.94, a: 1 },
        loadOp: 'clear',
      };
      if (sc > 1) {
        colorAttachment.view = msaaTexture.createView();
        colorAttachment.resolveTarget = canvasView;
        colorAttachment.storeOp = 'discard';
      } else {
        colorAttachment.view = canvasView;
        colorAttachment.storeOp = 'store';
      }

      const rp = encoder.beginRenderPass({ colorAttachments: [colorAttachment] });

      // 1) Grid (background)
      rp.setPipeline(gridPipelines[msaaIndex]);
      rp.setBindGroup(0, gridBindGroup);
      rp.draw(6);

      // 2) Line strip (chirp)
      rp.setPipeline(linePipelines[msaaIndex]);
      rp.setBindGroup(0, lineBindGroup);
      rp.draw(pointCount);

      rp.end();
      device.queue.submit([encoder.finish()]);

      // Submit autoscale reduce (async — updates yMin/yMax when readback completes).
      submitAutoscale(pointCount);

      // Draw axis labels on the 2D overlay.
      drawAxes(canvas.clientWidth, canvas.clientHeight, dpr);

      requestAnimationFrame(frame);
    }
    requestAnimationFrame(frame);
  }

  // ─────────────────────────────────────────────────────────────────────
  // Auto-initialization
  // ─────────────────────────────────────────────────────────────────────
  function initAll() {
    document.querySelectorAll('canvas.chirp-plot-canvas').forEach(initCanvas);
  }
  if (document.readyState === 'loading') {
    document.addEventListener('DOMContentLoaded', initAll);
  } else {
    initAll();
  }
})();