// ---------------------------------------------------------------------------
// ldk — a tiny, dependency-free plotting + statistics library for these decks.
//
// Everything runs in the browser as plain JavaScript building SVG nodes, so a
// slider move costs a fraction of a millisecond instead of a round trip
// through a WebAssembly R runtime.
// ---------------------------------------------------------------------------
ldk = {
  const NS = "http://www.w3.org/2000/svg";
  let uid = 0;

  // Deck palette, matching ldkslides.css
  const C = {
    green:  "#74b83d",
    blue:   "#4072c2",
    pink:   "#ED017D",
    orange: "#d18b2a",
    grey:   "#6E6E6E",
    ink:    "#222222",
    axis:   "#9aa0a6",
    faint:  "#efefef"
  };

  function el(tag, attrs, text) {
    const e = document.createElementNS(NS, tag);
    if (attrs) for (const k in attrs) if (attrs[k] != null) e.setAttribute(k, attrs[k]);
    if (text != null)
      // SVG collapses runs of whitespace, which would eat the wide gaps that
      // separate the readout columns; non-breaking spaces survive. (xml:space
      // is a namespaced attribute and does not survive setAttribute here.)
      e.textContent = tag === "text"
        ? String(text).replace(/ {2,}/g, m => " ".repeat(m.length))
        : text;
    return e;
  }

  function alpha(hex, a) {
    const n = parseInt(hex.slice(1), 16);
    return `rgba(${(n >> 16) & 255},${(n >> 8) & 255},${n & 255},${a})`;
  }

  // ------------------------------------------------------------------ numerics

  const LGC = [676.5203681218851, -1259.1392167224028, 771.32342877765313,
               -176.61502916214059, 12.507343278686905, -0.13857109526572012,
               9.9843695780195716e-6, 1.5056327351493116e-7];

  function lgamma(z) {
    if (z < 0.5) return Math.log(Math.PI / Math.sin(Math.PI * z)) - lgamma(1 - z);
    z -= 1;
    let x = 0.99999999999980993;
    for (let i = 0; i < 8; i++) x += LGC[i] / (z + i + 1);
    const t = z + 7.5;
    return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
  }

  const lchoose = (n, k) => lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1);

  // Chebyshev erfc (Numerical Recipes); plenty of accuracy for plotting.
  const ERFC_COF = [
    -1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2,
    -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4,
     4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6,
     1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8,
     6.529054439e-9, 5.059343495e-9, -9.91364156e-10, -2.27365122e-10,
     9.6467911e-11, 2.394038e-12, -6.886027e-12, 8.94487e-13,
     3.13092e-13, -1.12708e-13, 3.81e-16, 7.106e-15];

  function erfc(x) {
    const z = Math.abs(x), t = 2 / (2 + z), ty = 4 * t - 2;
    let d = 0, dd = 0;
    for (let j = ERFC_COF.length - 1; j > 0; j--) {
      const tmp = d;
      d = ty * d - dd + ERFC_COF[j];
      dd = tmp;
    }
    const ans = t * Math.exp(-z * z + 0.5 * (ERFC_COF[0] + ty * d) - dd);
    return x >= 0 ? ans : 2 - ans;
  }

  const dnorm = (x, m = 0, s = 1) => {
    const z = (x - m) / s;
    return Math.exp(-0.5 * z * z) / (s * Math.sqrt(2 * Math.PI));
  };
  const pnorm = (x, m = 0, s = 1) => 0.5 * erfc(-(x - m) / (s * Math.SQRT2));

  // Acklam's inverse normal CDF (relative error < 1.2e-9)
  const QN_A = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
                 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00];
  const QN_B = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
                 6.680131188771972e+01, -1.328068155288572e+01];
  const QN_C = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
                -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00];
  const QN_D = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
                3.754408661907416e+00];

  function qnorm(p, m = 0, s = 1) {
    if (p <= 0) return -Infinity;
    if (p >= 1) return Infinity;
    const lo = 0.02425, hi = 1 - lo;
    let x, q, r;
    if (p < lo) {
      q = Math.sqrt(-2 * Math.log(p));
      x = (((((QN_C[0] * q + QN_C[1]) * q + QN_C[2]) * q + QN_C[3]) * q + QN_C[4]) * q + QN_C[5]) /
          ((((QN_D[0] * q + QN_D[1]) * q + QN_D[2]) * q + QN_D[3]) * q + 1);
    } else if (p <= hi) {
      q = p - 0.5; r = q * q;
      x = (((((QN_A[0] * r + QN_A[1]) * r + QN_A[2]) * r + QN_A[3]) * r + QN_A[4]) * r + QN_A[5]) * q /
          (((((QN_B[0] * r + QN_B[1]) * r + QN_B[2]) * r + QN_B[3]) * r + QN_B[4]) * r + 1);
    } else {
      q = Math.sqrt(-2 * Math.log(1 - p));
      x = -(((((QN_C[0] * q + QN_C[1]) * q + QN_C[2]) * q + QN_C[3]) * q + QN_C[4]) * q + QN_C[5]) /
           ((((QN_D[0] * q + QN_D[1]) * q + QN_D[2]) * q + QN_D[3]) * q + 1);
    }
    return m + s * x;
  }

  function dbinom(k, n, p) {
    if (k < 0 || k > n) return 0;
    if (p <= 0) return k === 0 ? 1 : 0;
    if (p >= 1) return k === n ? 1 : 0;
    return Math.exp(lchoose(n, k) + k * Math.log(p) + (n - k) * Math.log1p(-p));
  }
  const dpois = (k, lam) => k < 0 ? 0 : Math.exp(-lam + k * Math.log(lam) - lgamma(k + 1));
  const dexp  = (x, rate) => x < 0 ? 0 : rate * Math.exp(-rate * x);
  const pexp  = (x, rate) => x <= 0 ? 0 : 1 - Math.exp(-rate * x);
  const dunif = (x, a, b) => (x >= a && x <= b) ? 1 / (b - a) : 0;
  const punif = (x, a, b) => x <= a ? 0 : x >= b ? 1 : (x - a) / (b - a);
  const dlnorm = (x, ml = 0, sl = 1) => x <= 0 ? 0 : dnorm(Math.log(x), ml, sl) / x;
  const dchisq1 = (x) => x <= 0 ? 0 : dnorm(Math.sqrt(x)) / Math.sqrt(x);

  // ---- regularized lower incomplete gamma P(a,x): series + continued fraction
  function gammaP(a, x) {
    if (!(x >= 0) || !(a > 0)) return NaN;
    if (x === 0) return 0;
    const pref = Math.exp(-x + a * Math.log(x) - lgamma(a));
    if (x < a + 1) {                       // series
      let ap = a, sum = 1 / a, del = sum;
      for (let i = 0; i < 500; i++) {
        ap++; del *= x / ap; sum += del;
        if (Math.abs(del) < Math.abs(sum) * 1e-16) break;
      }
      return sum * pref;
    }
    const TINY = 1e-300;                   // continued fraction for Q = 1 - P
    let b = x + 1 - a, c = 1 / TINY, d = 1 / b, h = d;
    for (let i = 1; i <= 500; i++) {
      const an = -i * (i - a);
      b += 2;
      d = an * d + b; if (Math.abs(d) < TINY) d = TINY;
      c = b + an / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d;
      const del = d * c;
      h *= del;
      if (Math.abs(del - 1) < 1e-16) break;
    }
    return 1 - pref * h;
  }

  // ---- regularized incomplete beta I_x(a,b)
  function betacf(a, b, x) {
    const TINY = 1e-300, qab = a + b, qap = a + 1, qam = a - 1;
    let c = 1, d = 1 - qab * x / qap;
    if (Math.abs(d) < TINY) d = TINY;
    d = 1 / d;
    let h = d;
    for (let m = 1; m <= 400; m++) {
      const m2 = 2 * m;
      let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY;
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d; h *= d * c;
      aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY;
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d;
      const del = d * c;
      h *= del;
      if (Math.abs(del - 1) < 3e-16) break;
    }
    return h;
  }
  function betai(a, b, x) {
    if (!(x > 0)) return 0;
    if (x >= 1) return 1;
    const bt = Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
                        a * Math.log(x) + b * Math.log1p(-x));
    return x < (a + 1) / (a + b + 2)
      ? bt * betacf(a, b, x) / a
      : 1 - bt * betacf(b, a, 1 - x) / b;
  }

  // ---- chi-squared and Student t
  const dchisq = (x, k) => x <= 0 ? 0 :
    Math.exp((k / 2 - 1) * Math.log(x) - x / 2 - lgamma(k / 2) - (k / 2) * Math.LN2);
  const pchisq = (x, k) => x <= 0 ? 0 : gammaP(k / 2, x / 2);

  const dt = (x, nu) => Math.exp(lgamma((nu + 1) / 2) - lgamma(nu / 2)
    - 0.5 * Math.log(nu * Math.PI) - (nu + 1) / 2 * Math.log1p(x * x / nu));
  const pt = (x, nu) => {
    const half = 0.5 * betai(nu / 2, 0.5, nu / (nu + x * x));
    return x > 0 ? 1 - half : half;
  };

  // quantiles by bisection — a handful of calls per redraw, so speed is moot
  function invert(cdf, p, lo, hi) {
    for (let i = 0; i < 90; i++) {
      const mid = (lo + hi) / 2;
      if (cdf(mid) < p) lo = mid; else hi = mid;
    }
    return (lo + hi) / 2;
  }
  const qchisq = (p, k) => invert(x => pchisq(x, k), p, 0, Math.max(60, k * 25));
  const qt = (p, nu) => invert(x => pt(x, nu), p, -1e4, 1e4);

  // Cumulative sums of a pmf enumerated over 0..n — cheaper and more accurate
  // here than a general incomplete-beta / incomplete-gamma implementation.
  function cdfFromPmf(pmf) {
    const out = new Array(pmf.length);
    let s = 0;
    for (let i = 0; i < pmf.length; i++) { s += pmf[i]; out[i] = Math.min(s, 1); }
    return out;
  }

  // ----------------------------------------------------------------- sequences

  function seq(a, b, n) {
    const out = new Array(n);
    for (let i = 0; i < n; i++) out[i] = a + (b - a) * i / (n - 1);
    return out;
  }
  const cumsum = (v) => { let s = 0; return v.map(x => (s += x)); };
  const sum = (v) => v.reduce((a, b) => a + b, 0);
  const mean = (v) => sum(v) / v.length;

  function sd(v) {
    const m = mean(v);
    return Math.sqrt(v.reduce((a, x) => a + (x - m) * (x - m), 0) / (v.length - 1));
  }

  function quantileSorted(s, p) {
    const h = (s.length - 1) * p, i = Math.floor(h);
    const j = Math.min(i + 1, s.length - 1);
    return s[i] + (h - i) * (s[j] - s[i]);
  }

  function hist(v, lo, hi, nb) {
    const cnt = new Float64Array(nb), w = (hi - lo) / nb;
    for (let i = 0; i < v.length; i++) {
      const k = Math.floor((v[i] - lo) / w);
      if (k >= 0 && k < nb) cnt[k]++;
      else if (v[i] === hi) cnt[nb - 1]++;
    }
    const dens = new Array(nb);
    for (let i = 0; i < nb; i++) dens[i] = cnt[i] / (v.length * w);
    return { lo, hi, w, nb, dens, count: cnt };
  }

  // --------------------------------------------------------------------- rngs

  // mulberry32: small, fast, seedable — the demos stay reproducible.
  function rng(seed) {
    let a = seed >>> 0;
    const u = () => {
      a |= 0; a = (a + 0x6D2B79F5) | 0;
      let t = Math.imul(a ^ (a >>> 15), 1 | a);
      t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
      return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
    };
    let spare = null;
    u.normal = () => {                       // Marsaglia polar
      if (spare !== null) { const v = spare; spare = null; return v; }
      let x, y, s;
      do { x = 2 * u() - 1; y = 2 * u() - 1; s = x * x + y * y; } while (s >= 1 || s === 0);
      const f = Math.sqrt(-2 * Math.log(s) / s);
      spare = y * f;
      return x * f;
    };
    u.exp = () => -Math.log(1 - u());
    u.cauchy = () => Math.tan(Math.PI * (u() - 0.5));
    return u;
  }

  // Weighted sampling without replacement (Efraimidis–Spirakis exponential race)
  function sampleWeighted(n, w, u) {
    const N = w.length;
    const keys = new Array(N);
    for (let i = 0; i < N; i++) keys[i] = w[i] > 0 ? -Math.log(1 - u()) / w[i] : Infinity;
    const idx = Array.from({ length: N }, (_, i) => i);
    idx.sort((a, b) => keys[a] - keys[b]);
    return idx.slice(0, n);
  }

  // Sampling *with* replacement — what the bootstrap needs, and what
  // sampleWeighted deliberately is not.
  function sampleReplace(n, v, u) {
    const out = new Array(n);
    for (let i = 0; i < n; i++) out[i] = v[Math.floor(u() * v.length)];
    return out;
  }

  // Beta density — the conjugate posterior for a proportion.
  const dbeta = (x, a, b) => (x <= 0 || x >= 1) ? 0 :
    Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
             (a - 1) * Math.log(x) + (b - 1) * Math.log(1 - x));

  // Beta quantile, by inverting the regularized incomplete beta above. Used for
  // the Clopper-Pearson interval, whose endpoints are exactly beta quantiles.
  const qbeta = (p, a, b) => invert(x => betai(a, b, x), p, 0, 1);

  // Gini coefficient of a sorted, non-negative sample: twice the area between
  // the Lorenz curve and the diagonal. No closed-form standard error exists,
  // which is exactly why the bootstrap earns its place.
  function gini(sorted) {
    const n = sorted.length;
    let num = 0, tot = 0;
    for (let i = 0; i < n; i++) { num += (i + 1) * sorted[i]; tot += sorted[i]; }
    return tot > 0 ? (2 * num) / (n * tot) - (n + 1) / n : 0;
  }

  // ------------------------------------------------------- regression + series

  // Simple OLS of y on x, with both the classical and the heteroskedasticity-
  // robust (HC1) standard error. Deck 08 leans on the pair: the estimate and its
  // standard error fail independently, and the standard error fails first.
  function ols(x, y) {
    const n = x.length;
    const mx = mean(x), my = mean(y);
    let sxx = 0, sxy = 0;
    for (let i = 0; i < n; i++) { const dx = x[i] - mx; sxx += dx * dx; sxy += dx * (y[i] - my); }
    const b = sxx > 0 ? sxy / sxx : 0;
    const a = my - b * mx;
    let rss = 0, meat = 0, tss = 0;
    for (let i = 0; i < n; i++) {
      const e = y[i] - a - b * x[i], dx = x[i] - mx;
      rss += e * e; meat += dx * dx * e * e; tss += (y[i] - my) * (y[i] - my);
    }
    const se = sxx > 0 && n > 2 ? Math.sqrt(rss / (n - 2) / sxx) : NaN;
    const seR = sxx > 0 && n > 2
      ? Math.sqrt(meat * n / (n - 2)) / sxx                    // HC1
      : NaN;
    return { a, b, se, seR, t: b / se, tR: b / seR,
             r2: tss > 0 ? 1 - rss / tss : 0, n };
  }

  const corr = (x, y) => {
    const mx = mean(x), my = mean(y);
    let sxy = 0, sxx = 0, syy = 0;
    for (let i = 0; i < x.length; i++) {
      const dx = x[i] - mx, dy = y[i] - my;
      sxy += dx * dy; sxx += dx * dx; syy += dy * dy;
    }
    return sxx > 0 && syy > 0 ? sxy / Math.sqrt(sxx * syy) : 0;
  };

  // Stationary AR(1), started from its stationary distribution so the first
  // observations are not systematically closer to zero than the rest.
  function ar1(n, rho, u, sd = 1) {
    const out = new Float64Array(n);
    let v = u.normal() * sd / Math.sqrt(Math.max(1 - rho * rho, 1e-9));
    for (let i = 0; i < n; i++) { v = rho * v + sd * u.normal(); out[i] = v; }
    return out;
  }

  const randomWalk = (n, u, sd = 1) => {
    const out = new Float64Array(n);
    let v = 0;
    for (let i = 0; i < n; i++) { v += sd * u.normal(); out[i] = v; }
    return out;
  };

  // G clusters of m units with intra-cluster correlation rho: a shared group
  // draw plus an idiosyncratic one, weighted so the total variance stays 1.
  function clustered(G, m, rho, u) {
    const out = new Float64Array(G * m);
    const sg = Math.sqrt(Math.max(rho, 0)), si = Math.sqrt(Math.max(1 - rho, 0));
    for (let g = 0; g < G; g++) {
      const shared = sg * u.normal();
      for (let i = 0; i < m; i++) out[g * m + i] = shared + si * u.normal();
    }
    return out;
  }

  // ---------------------------------------------------------------- formatting

  const fmt = (v, d = 2) => (Number.isFinite(v) ? v.toFixed(d) : "—");
  const signed = (v, d = 2) => (Number.isFinite(v) ? (v >= 0 ? "+" : "") + v.toFixed(d) : "—");

  function smart(v) {
    if (!Number.isFinite(v)) return "—";
    const a = Math.abs(v);
    if (a >= 1e5 || (a > 0 && a < 1e-3)) return v.toExponential(1);
    return v.toFixed(a >= 100 ? 0 : a >= 10 ? 1 : 2);
  }

  function nice(lo, hi, n) {
    const span = (hi - lo) || Math.abs(hi) || 1;
    const raw = span / n;
    const mag = Math.pow(10, Math.floor(Math.log10(raw)));
    const norm = raw / mag;
    const step = mag * (norm < 1.5 ? 1 : norm < 3 ? 2 : norm < 7 ? 5 : 10);
    const out = [];
    for (let v = Math.ceil(lo / step - 1e-9) * step; v <= hi + 1e-9 * span; v += step)
      out.push(Math.abs(v) < step * 1e-9 ? 0 : v);
    return { ticks: out, step };
  }

  function logTicks(lo, hi) {
    const decades = Math.log10(hi / lo);
    const mults = decades > 2.2 ? [1] : decades > 1.1 ? [1, 3] : [1, 2, 5];
    const out = [];
    for (let e = Math.floor(Math.log10(lo)); e <= Math.ceil(Math.log10(hi)); e++)
      for (const k of mults) {
        const v = k * Math.pow(10, e);
        if (v >= lo * 0.9999 && v <= hi * 1.0001) out.push(v);
      }
    return out;
  }

  const decimalsFor = (step) =>
    Math.max(0, Math.min(6, -Math.floor(Math.log10(Math.abs(step)) + 1e-9)));

  // -------------------------------------------------------------------- charts

  function chart(W, H, extra) {
    const svg = el("svg", Object.assign({
      viewBox: `0 0 ${W} ${H}`,
      style: `width:100%;height:auto;display:block;font-family:inherit;color:${C.ink}`
    }, extra || {}));
    svg.__W = W;
    svg.__H = H;
    return svg;
  }

  function panel(svg, o) {
    const L = o.left ?? 0, T = o.top ?? 0;
    const W = o.w ?? svg.__W, H = o.h ?? svg.__H;
    const m = Object.assign({ l: 74, r: 22, t: 52, b: 54 }, o.margin || {});
    const fs = o.fontSize ?? 15;
    const iw = W - m.l - m.r, ih = H - m.t - m.b;
    const px0 = L + m.l, py0 = T + m.t;

    let [xa, xb] = o.xdom, [ya, yb] = o.ydom;
    if (!(xb > xa)) xb = xa + 1;
    if (!(yb > ya)) yb = ya + 1;

    const lg = o.xlog ? Math.log : (v => v);
    const xA = lg(xa), xB = lg(xb);
    const x = v => px0 + (lg(v) - xA) / (xB - xA) * iw;
    const y = v => py0 + ih - (v - ya) / (yb - ya) * ih;

    const gAxis = el("g");
    const id = `ldkclip${++uid}`;
    const defs = el("defs");
    const cp = el("clipPath", { id });
    cp.append(el("rect", { x: px0 - 3, y: py0 - 8, width: iw + 6, height: ih + 11 }));
    defs.append(cp);
    const gc = el("g", { "clip-path": `url(#${id})` });
    const gt = el("g");
    svg.append(defs, gAxis, gc, gt);

    // ticks
    let xt, xfmt, yt, yfmt;
    if (o.xticks) {
      xt = o.xticks;
      xfmt = o.xfmt || (v => String(v));
    } else if (o.xlog) {
      xt = logTicks(xa, xb);
      xfmt = o.xfmt || (v => v >= 1 ? String(Math.round(v)) : String(v));
    } else {
      const n = nice(xa, xb, o.nx ?? 6);
      xt = n.ticks;
      xfmt = o.xfmt || (v => v.toFixed(decimalsFor(n.step)));
    }
    if (o.yticks) {
      yt = o.yticks;
      yfmt = o.yfmt || (v => String(v));
    } else {
      const n = nice(ya, yb, o.ny ?? 5);
      yt = n.ticks;
      yfmt = o.yfmt || (v => v.toFixed(decimalsFor(n.step)));
    }

    if (o.axes === false) { xt = []; yt = []; }   // bare drawing surface

    if (o.axes !== false) {
    gAxis.append(el("line", { x1: px0, y1: py0 + ih, x2: px0 + iw, y2: py0 + ih,
                              stroke: C.axis, "stroke-width": 1.3 }));
    gAxis.append(el("line", { x1: px0, y1: py0, x2: px0, y2: py0 + ih,
                              stroke: C.axis, "stroke-width": 1.3 }));
    }
    for (const t of xt) {
      const px = x(t);
      if (px < px0 - 0.5 || px > px0 + iw + 0.5) continue;
      gAxis.append(el("line", { x1: px, y1: py0 + ih, x2: px, y2: py0 + ih + 5,
                                stroke: C.axis, "stroke-width": 1.3 }));
      gAxis.append(el("text", { x: px, y: py0 + ih + 6 + fs, "text-anchor": "middle",
                                "font-size": fs, fill: C.grey }, xfmt(t)));
    }
    for (const t of yt) {
      const py = y(t);
      if (py < py0 - 0.5 || py > py0 + ih + 0.5) continue;
      gAxis.append(el("line", { x1: px0 - 5, y1: py, x2: px0, y2: py,
                                stroke: C.axis, "stroke-width": 1.3 }));
      gAxis.append(el("text", { x: px0 - 10, y: py + fs * 0.35, "text-anchor": "end",
                                "font-size": fs, fill: C.grey }, yfmt(t)));
    }
    if (o.xlab)
      gAxis.append(el("text", { x: px0 + iw / 2, y: T + H - 8, "text-anchor": "middle",
                                "font-size": fs + 2, fill: C.ink }, o.xlab));
    if (o.ylab) {
      const ty = py0 + ih / 2, tx = L + 20;
      gAxis.append(el("text", { x: tx, y: ty, "text-anchor": "middle", "font-size": fs + 2,
                                fill: C.ink, transform: `rotate(-90 ${tx} ${ty})` }, o.ylab));
    }

    const P = {
      x, y, svg, gAxis, gc, gt, fs,
      xdom: [xa, xb], ydom: [ya, yb],
      area: { x0: px0, y0: py0, w: iw, h: ih },

      add(node, top) { (top ? gt : gc).append(node); return node; },

      path(pts, a) {
        let d = "";
        for (let i = 0; i < pts.length; i++)
          d += (i ? "L" : "M") + x(pts[i][0]).toFixed(2) + "," + y(pts[i][1]).toFixed(2);
        return P.add(el("path", Object.assign(
          { d, fill: "none", "stroke-linejoin": "round", "stroke-linecap": "round" }, a)));
      },

      fill(pts, base, a) {
        let d = "";
        for (let i = 0; i < pts.length; i++)
          d += (i ? "L" : "M") + x(pts[i][0]).toFixed(2) + "," + y(pts[i][1]).toFixed(2);
        d += "L" + x(pts[pts.length - 1][0]).toFixed(2) + "," + y(base).toFixed(2);
        d += "L" + x(pts[0][0]).toFixed(2) + "," + y(base).toFixed(2) + "Z";
        return P.add(el("path", Object.assign({ d, stroke: "none" }, a)));
      },

      seg(x1, y1, x2, y2, a) {
        return P.add(el("line", Object.assign(
          { x1: x(x1), y1: y(y1), x2: x(x2), y2: y(y2) }, a)));
      },
      vline(v, a) { return P.seg(v, ya, v, yb, a); },
      hline(v, a) { return P.seg(xa, v, xb, v, a); },

      bar(x1, x2, y1, y2, a) {
        const X1 = x(x1), X2 = x(x2), Y1 = y(y1), Y2 = y(y2);
        return P.add(el("rect", Object.assign({
          x: Math.min(X1, X2), y: Math.min(Y1, Y2),
          width: Math.max(Math.abs(X2 - X1), 0.5),
          height: Math.max(Math.abs(Y2 - Y1), 0)
        }, a)));
      },

      dot(vx, vy, r, a) {
        return P.add(el("circle", Object.assign({ cx: x(vx), cy: y(vy), r }, a)));
      },

      // mtext-style annotation above the plot region
      header(s, a = {}) {
        const adj = a.adj ?? 0;
        const py = py0 - 14 - (a.line ?? 0) * (fs + 7);
        return P.add(el("text", {
          x: adj === 1 ? px0 + iw : adj === 0.5 ? px0 + iw / 2 : px0,
          y: py,
          "text-anchor": adj === 1 ? "end" : adj === 0.5 ? "middle" : "start",
          "font-size": a.size ?? fs + 2,
          "font-weight": a.weight ?? "bold",
          fill: a.fill ?? C.ink
        }, s), true);
      },

      legend(items, opt = {}) {
        const size = opt.size ?? fs;
        const pad = opt.pad ?? 8;
        const rowH = size + 8, swW = 30;
        const g = el("g");
        let maxw = 0;
        items.forEach((it, i) => {
          const ty = i * rowH + size;
          if (it.swatch)
            g.append(el("rect", { x: 0, y: ty - size * 0.82, width: swW,
                                  height: size * 0.85, fill: it.swatch, stroke: "none" }));
          else
            g.append(el("line", { x1: 0, y1: ty - size * 0.35, x2: swW, y2: ty - size * 0.35,
                                  stroke: it.color, "stroke-width": it.lw ?? 3.5,
                                  "stroke-dasharray": it.dash ?? null,
                                  "stroke-linecap": "round" }));
          g.append(el("text", { x: swW + 9, y: ty, "font-size": size,
                                fill: opt.fill ?? C.ink }, it.label));
          maxw = Math.max(maxw, swW + 9 + it.label.length * size * 0.54);
        });
        // opt-in backing panel, for legends that must sit over the data
        if (opt.bg) {
          const q = 8;
          g.insertBefore(el("rect", {
            x: -q, y: -q + 3, width: maxw + 2 * q,
            height: items.length * rowH + 2 * q - 6,
            fill: "white", "fill-opacity": opt.bg === true ? 0.85 : opt.bg, rx: 5
          }), g.firstChild);
        }
        const corner = opt.corner ?? "topright";
        const gx = corner.includes("right") ? px0 + iw - maxw - pad : px0 + pad;
        const gy = corner.includes("top") ? py0 + pad
                                          : py0 + ih - items.length * rowH - pad;
        g.setAttribute("transform", `translate(${gx},${gy})`);
        gt.append(g);
        return g;
      }
    };
    return P;
  }

  // Convenience: one chart with exactly one panel filling it.
  function single(W, H, o) {
    const svg = chart(W, H);
    return { svg, p: panel(svg, o) };
  }

  return {
    C, el, alpha, chart, panel, single,
    lgamma, lchoose, erfc, dnorm, pnorm, qnorm, dbinom, dpois, dexp, pexp,
    dunif, punif, dlnorm, dchisq1, cdfFromPmf,
    gammaP, betai, dbeta, qbeta, dchisq, pchisq, qchisq, dt, pt, qt, invert,
    gini,
    seq, cumsum, sum, mean, sd, quantileSorted, hist, rng, sampleWeighted,
    sampleReplace, ols, corr, ar1, randomWalk, clustered,
    fmt, signed, smart, nice, logTicks
  };
}