React

d3gl ships first-class React bindings at @mapequation/d3gl/react. For maps there is <GeoMap>, which mounts the project-once map engine for you. For everything else — trees, generic 2D plots — there is a fully declarative API: a <Plot> element with <Layer> and <Points> children. You describe the chart as JSX; <Plot> builds the imperative plot engine under the hood, translates each child into an engine.layer / engine.points call (sibling order = paint order), and manages the lifecycle — no useEffect, no DOM, no engine plumbing in your component. Both examples below render through the same universal control bar as every other example on this site (backend switch, export, live perf); the components just receive backend and the canvas size as props.

World map

<GeoMap> mounts the map engine and hands it back via the onReady callback, where you add layers, enable zoom, and render once. It takes width, height, a d3-geo projection, and an optional backend ("webgl", "canvas", or "svg"). When the backend prop changes the component keeps the engine alive and calls map.setBackend() — no remount — so the layers and the current zoom/pan are preserved.

autowebglcanvassvg

Export PNG

fps 0frame 0 ms

WorldMapReact.tsx

WorldMapReact.tsx

import { geoNaturalEarth1 } from "d3-geo";
import { GeoMap } from "@mapequation/d3gl/react";
import { fitProjection } from "@mapequation/d3gl/geo";
import Example from "../../components/Example.js";
import { loadWorld } from "../shared/geo-data.js";

const OCEAN = "#d4e6f5";
const LAND = "#e3e6ea";

/**
 * The zoomable world map built with the React `<GeoMap>` wrapper. The universal
 * <Example> harness supplies `backend` + the measured `width`/`height` and an
 * export hook via `registerEngine` — this file is just the viz.
 */
export default function WorldMapReact() {
  return (
    <Example width={720} height={380}>
      {({ backend, width, height, registerEngine }) => (
        <GeoMap
          width={width}
          height={height}
          backend={backend}
          projection={fitProjection(geoNaturalEarth1(), { type: "Sphere" }, width, height)}
          onReady={(map) => {
            const world = loadWorld();
            map.layer("ocean", [world.sphere], { fill: OCEAN });
            map.layer("land", [world.land], { fill: LAND, stroke: "#9aa3ad", lineWidth: 0.5 });
            map.enableZoom([1, 50]);
            map.render();
            registerEngine(map);
          }}
        />
      )}
    </Example>
  );
}

geo-data.ts

geo-data.ts

import { geoGraticule } from "d3-geo";
import type { Feature, FeatureCollection, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon } from "geojson";
import { feature } from "topojson-client";
import land110m from "world-atlas/land-110m.json";

/** A synthetic grid cell with a continuous value and a categorical bioregion. */
export interface Cell {
  id: string;
  geometry: Polygon;
  /** Cell centroid [lon, lat]. */
  center: [number, number];
  /** Continuous field in [0, 1] (value field). */
  value: number;
  /** Categorical bioregion id in 0..7. */
  bioregion: number;
}

/** Base cell size in degrees; the example scales it by powers of two via a slider. */
export const BASE_STEP = 1;

function clamp01(x: number): number {
  return Math.max(0, Math.min(1, x));
}

/** Generate a global grid of `step`°×`step`° cells with smooth synthetic fields. */
export function makeCells(step: number = BASE_STEP): Cell[] {
  const cells: Cell[] = [];
  let col = 0;
  for (let lon = -180; lon < 180; lon += step, col++) {
    let row = 0;
    for (let lat = -90; lat < 90; lat += step, row++) {
      const lonR = (lon * Math.PI) / 180;
      const latR = (lat * Math.PI) / 180;
      const value = clamp01(0.5 + 0.5 * Math.sin(lonR * 2) * Math.cos(latR * 3));
      const field = (Math.sin(lon / 40) + Math.cos(lat / 30)) * 0.5 + 1; // ~[0,2]
      const bioregion = Math.min(7, Math.max(0, Math.floor((field / 2) * 8)));
      const geometry: Polygon = {
        type: "Polygon",
        coordinates: [
          [
            [lon, lat],
            [lon, lat + step],
            [lon + step, lat + step],
            [lon + step, lat],
            [lon, lat],
          ],
        ],
      };
      cells.push({
        id: `${col}-${row}`,
        geometry,
        center: [lon + step / 2, lat + step / 2],
        value,
        bioregion,
      });
    }
  }
  return cells;
}

/** A fine grid over the central third of the globe (lon ±60°, lat ±30°), 4° cells —
 *  a "dense" demo layer (used clipped to land). */
export function centreCells(): Cell[] {
  return makeCells(4).filter((c) => Math.abs(c.center[0]) <= 60 && Math.abs(c.center[1]) <= 30);
}

/** Wrap cells as a FeatureCollection for projection fitting. */
export function cellsToFeatureCollection(cells: readonly Cell[]): FeatureCollection {
  const features: Feature[] = cells.map((c) => ({
    type: "Feature",
    properties: { id: c.id },
    geometry: c.geometry,
  }));
  return { type: "FeatureCollection", features };
}

/** A GeoJSON object d3-geo can fill that isn't part of the strict GeoJSON spec. */
export type Sphere = { type: "Sphere" };

/** The land outline (Natural Earth 110m) plus a sphere to fill as ocean. */
export interface World {
  sphere: Sphere;
  land: MultiPolygon;
}

// Derive the topojson Topology type from feature()'s own signature so we don't
// take a direct dependency on the (transitive) topojson-specification types.
type Topology = Parameters<typeof feature>[0];

/**
 * Convert the bundled world-atlas TopoJSON into a land MultiPolygon and a sphere.
 */
export function loadWorld(): World {
  const topo = land110m as unknown as Topology;
  const fc = feature(topo, topo.objects.land!) as unknown as FeatureCollection<MultiPolygon>;
  return { sphere: { type: "Sphere" }, land: fc.features[0]!.geometry };
}

/** A few well-known cities to show point geometry rendered alongside the grid. */
export interface City {
  id: string;
  name: string;
  geometry: Point;
}

export function makeCities(): City[] {
  const places: [string, number, number][] = [
    ["London", -0.13, 51.51],
    ["New York", -74.01, 40.71],
    ["Tokyo", 139.69, 35.69],
    ["Sydney", 151.21, -33.87],
    ["Cape Town", 18.42, -33.92],
    ["Rio de Janeiro", -43.2, -22.91],
    ["Nairobi", 36.82, -1.29],
    ["Mumbai", 72.88, 19.08],
  ];
  return places.map(([name, lon, lat]) => ({
    id: name!,
    name: name!,
    geometry: { type: "Point", coordinates: [lon!, lat!] },
  }));
}

/** A 20° graticule as one MultiLineString feature. */
export function makeGraticule(): Feature<MultiLineString> {
  return { type: "Feature", properties: {}, geometry: geoGraticule().step([20, 20])() };
}

/** A great-circle-ish route as a LineString feature (London -> New York -> Tokyo). */
export function makeRoute(): Feature<LineString> {
  return {
    type: "Feature", properties: {},
    geometry: { type: "LineString", coordinates: [[-0.13, 51.51], [-74.01, 40.71], [139.69, 35.69]] },
  };
}

/** A cluster of locations as one MultiPoint feature. */
export function makeCluster(): Feature<MultiPoint> {
  return {
    type: "Feature", properties: {},
    geometry: { type: "MultiPoint", coordinates: [[18.42, -33.92], [151.21, -33.87], [-43.2, -22.91], [36.82, -1.29], [72.88, 19.08]] },
  };
}

/** A standalone Polygon feature (a box over the Sahara) to showcase polygon geometry.
 *  Wound CLOCKWISE so d3-geo fills the small box (not the sphere complement). */
export function makeDemoPolygon(): Feature<Polygon> {
  return {
    type: "Feature", properties: { name: "demo-region" },
    geometry: { type: "Polygon", coordinates: [[[0, 15], [0, 30], [30, 30], [30, 15], [0, 15]]] },
  };
}

/** A handful of major rivers as rough named polylines (the bundled world-atlas data
 *  has no rivers), shown on the GeoJSON-features map and used as streaming cluster
 *  centers. Coordinates are approximate [lon, lat] traces, mouth → source. */
export function makeMajorRivers(): Feature<LineString, { name: string }>[] {
  const rivers: [string, [number, number][]][] = [
    ["Amazon", [[-50.0, -0.7], [-55.5, -2.5], [-60.0, -3.1], [-67.9, -3.5], [-73.2, -4.5]]],
    ["Nile", [[31.3, 31.4], [32.9, 24.1], [32.5, 15.6], [32.5, 9.5], [31.6, 2.3]]],
    ["Mississippi", [[-89.2, 29.2], [-90.1, 32.3], [-90.2, 38.6], [-91.2, 43.5], [-95.0, 47.2]]],
    ["Yangtze", [[121.8, 31.4], [114.3, 30.6], [106.5, 29.6], [100.2, 26.9], [94.7, 33.5]]],
    ["Congo", [[12.4, -6.0], [16.2, -4.3], [20.0, -1.0], [25.2, 0.5], [27.2, 3.0]]],
    ["Volga", [[48.0, 46.3], [45.0, 48.7], [44.5, 51.6], [47.5, 54.3], [37.0, 57.3]]],
    ["Ganges", [[90.5, 22.5], [88.0, 24.5], [83.0, 25.4], [78.0, 26.5], [78.9, 30.1]]],
  ];
  return rivers.map(([name, coordinates]) => ({
    type: "Feature",
    properties: { name },
    geometry: { type: "LineString", coordinates },
  }));
}

// ---------------------------------------------------------------------------
// Streaming sources — async generators that emit batches of features lazily
// (only `batchSize` are materialized per tick, never the whole `total`), so a
// consumer can `await`-iterate and append them live. Points/cells are CLUSTERED
// around cities + major-river vertices (not uniform), which the world-map
// examples then clip to land. Used by the "streaming data" examples.
// ---------------------------------------------------------------------------

/** A solid default color all streamed features start with (the example's
 *  "randomize" button swaps in a new color for new + retained features). */
export const DEFAULT_STREAM_COLOR = "#e23b2f";

/** Per-feature properties: a stable id (continues across batches) + a color the
 *  example owns (constant by default; swapped by the randomize button). */
export interface StreamProps {
  id: number;
  color: string;
}
export type StreamPoint = Feature<Point, StreamProps>;
export type StreamPolygon = Feature<Polygon, StreamProps>;

export interface StreamOptions {
  /** Total features emitted before the generator completes. Default 10,000,000. */
  total?: number;
  /** Features per yielded batch. A function is re-read every batch, so the caller can
   *  resize adaptively (e.g. to fit a frame budget). Default 1000. */
  batchSize?: number | (() => number);
  /** Artificial delay between batches (ms), to mirror loading from a file/network.
   *  Even 0 yields a macrotask so the browser can paint between batches. Default 0. */
  delayMs?: number;
  /** Seed for the deterministic PRNG (reproducible streams). Default 1. */
  seed?: number;
  /** Cooperative cancellation: iteration stops once `signal.aborted` is true. */
  signal?: { aborted: boolean };
}

/** Small, fast, seedable PRNG (mulberry32) — deterministic so streams reproduce. */
function mulberry32(seed: number): () => number {
  let a = seed >>> 0;
  return () => {
    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;
  };
}

const clamp = (x: number, lo: number, hi: number): number => Math.max(lo, Math.min(hi, x));

/**
 * A clustered point process (Thomas-like): many weighted "hotspot" parents at a
 * range of scales — tight, dense cities; medium clumps along rivers; and a field
 * of random global hotspots with varied spread/weight. Offspring fall around a
 * parent with an exponential (heavy-tailed) radius, so the result is patchy and
 * multi-scale — far closer to real species-occurrence density than a smooth
 * gaussian. Clip-to-land then carves the continental outline.
 */
interface Parent {
  lon: number;
  lat: number;
  /** Mean offspring radius in degrees. */
  spread: number;
  /** Relative share of points drawn from this parent. */
  weight: number;
}

function buildParents(rng: () => number): Parent[] {
  const parents: Parent[] = [];
  for (const c of makeCities())
    parents.push({ lon: c.geometry.coordinates[0]!, lat: c.geometry.coordinates[1]!, spread: 1.5, weight: 3 });
  for (const river of makeMajorRivers())
    for (const p of river.geometry.coordinates) parents.push({ lon: p[0]!, lat: p[1]!, spread: 2.5, weight: 2 });
  // Random global hotspots: rng()*rng() biases toward small spreads/weights, so
  // most clumps are tight with a few diffuse ones — a wide range of scales.
  for (let i = 0; i < 240; i++) {
    parents.push({
      lon: rng() * 360 - 180,
      lat: rng() * 160 - 80,
      spread: 0.6 + 9 * rng() * rng(),
      weight: 0.2 + 3 * rng() * rng(),
    });
  }
  return parents;
}

function cumulativeWeights(parents: readonly Parent[]): number[] {
  const cum: number[] = [];
  let s = 0;
  for (const p of parents) {
    s += p.weight;
    cum.push(s);
  }
  return cum;
}

/** Pick a parent by weight (binary search over cumulative weights). */
function pickParent(rng: () => number, parents: readonly Parent[], cum: readonly number[]): Parent {
  const r = rng() * cum[cum.length - 1]!;
  let lo = 0;
  let hi = cum.length - 1;
  while (lo < hi) {
    const mid = (lo + hi) >> 1;
    if (cum[mid]! < r) lo = mid + 1;
    else hi = mid;
  }
  return parents[lo]!;
}

/** Offspring location around a parent: exponential radius, uniform direction. */
function clusteredLonLat(rng: () => number, p: Parent): [number, number] {
  const radius = -p.spread * Math.log(Math.max(1e-9, rng())); // exponential, mean = spread
  const ang = rng() * 2 * Math.PI;
  return [
    clamp(p.lon + radius * Math.cos(ang), -180, 180),
    clamp(p.lat + radius * Math.sin(ang), -90, 90),
  ];
}

/**
 * An irregular, star-convex polygon ring (3–10 vertices, varied per-vertex radius)
 * centered at [clon, clat] with overall extent ≤ ~`size`° — a rough species range.
 * Angles are evenly spaced with bounded jitter so they stay monotonic ⇒ the ring is
 * simple (non-self-intersecting) and closed. Longitude offsets are widened by
 * 1/cos(lat) so ranges don't look squished toward the poles.
 *
 * WINDING: built CLOCKWISE in [lon, lat] (note the NEGATIVE angle). d3-geo fills on
 * the sphere, and a small exterior ring wound counter-clockwise is treated as its
 * complement → it fills the whole map. See `AGENTS.md` and `geo/project.ts`.
 */
function randomRangeRing(rng: () => number, clon: number, clat: number, size: number): [number, number][] {
  const verts = 3 + Math.floor(rng() * 8); // 3..10
  // Strongly heavy-tailed size: the vast majority of ranges are TINY and only ~3% are
  // visibly large. At high counts the translucent fill then reads as a density gradient
  // (clustered richness hotspots) instead of saturating the whole map red.
  const base =
    rng() < 0.03
      ? size * (0.1 + 0.15 * rng()) // ~3% larger ranges: 0.10..0.25 * size
      : size * (0.02 + 0.07 * rng() * rng()); // most tiny: 0.02..0.09 * size, biased small (rng²)
  const latScale = 1 / Math.max(0.25, Math.cos((clat * Math.PI) / 180));
  const ring: [number, number][] = [];
  for (let i = 0; i < verts; i++) {
    const ang = -((i + 0.5 * rng()) / verts) * 2 * Math.PI; // NEGATIVE ⇒ clockwise ⇒ fills interior
    const r = base * (0.5 + rng()); // per-vertex radius variation
    ring.push([
      clamp(clon + r * Math.cos(ang) * latScale, -180, 180),
      clamp(clat + r * Math.sin(ang), -90, 90),
    ]);
  }
  ring.push(ring[0]!); // close the ring
  return ring;
}

const tick = (ms: number): Promise<void> => new Promise((r) => setTimeout(r, ms));

/** Stream world points as GeoJSON `Feature<Point>` batches, clustered around
 *  many multi-scale hotspots. All start with DEFAULT_STREAM_COLOR. */
export async function* makeStreamingPoints(opts: StreamOptions = {}): AsyncGenerator<StreamPoint[]> {
  const { total = 10_000_000, batchSize = 1000, delayMs = 0, seed = 1, signal } = opts;
  const sizeOf = (): number => Math.max(1, Math.floor(typeof batchSize === "function" ? batchSize() : batchSize));
  const rng = mulberry32(seed);
  const parents = buildParents(rng);
  const cum = cumulativeWeights(parents);
  let id = 0;
  while (id < total) {
    if (signal?.aborted) return;
    const n = Math.min(sizeOf(), total - id);
    const batch: StreamPoint[] = new Array(n);
    for (let k = 0; k < n; k++) {
      const [lon, lat] = clusteredLonLat(rng, pickParent(rng, parents, cum));
      batch[k] = {
        type: "Feature",
        properties: { id: id++, color: DEFAULT_STREAM_COLOR },
        geometry: { type: "Point", coordinates: [lon, lat] },
      };
    }
    yield batch;
    await tick(delayMs);
  }
}

/** Stream irregular polygon "ranges" clustered around hotspots — each a 3–10
 *  vertex star-convex polygon of varied size ≤ ~`size`°, to mimic species ranges.
 *  All start with DEFAULT_STREAM_COLOR; the example renders them very transparent
 *  so overlapping ranges build up richness. */
export async function* makeStreamingPolygons(
  opts: StreamOptions & { size?: number } = {},
): AsyncGenerator<StreamPolygon[]> {
  const { total = 10_000_000, batchSize = 1000, delayMs = 0, seed = 1, size = 16, signal } = opts;
  const sizeOf = (): number => Math.max(1, Math.floor(typeof batchSize === "function" ? batchSize() : batchSize));
  const rng = mulberry32(seed);
  const parents = buildParents(rng);
  const cum = cumulativeWeights(parents);
  let id = 0;
  while (id < total) {
    if (signal?.aborted) return;
    const n = Math.min(sizeOf(), total - id);
    const batch: StreamPolygon[] = new Array(n);
    for (let k = 0; k < n; k++) {
      const [clon, clat] = clusteredLonLat(rng, pickParent(rng, parents, cum));
      batch[k] = {
        type: "Feature",
        properties: { id: id++, color: DEFAULT_STREAM_COLOR },
        geometry: { type: "Polygon", coordinates: [randomRangeRing(rng, clon, clat, size)] },
      };
    }
    yield batch;
    await tick(delayMs);
  }
}

Phylogenetic tree

The same 64-tip rectangular phylogram as the vanilla simple-tree example (d3.link(curveStepBefore) branches, schemeCategory10 tip nodes), written declaratively with <Plot>. The whole chart is JSX: a <Layer> draws the branches with a d3-shape link generator straight into the d3gl context, a <Points> places the coloured tip nodes, and zoom={[0.5, 40]} enables scroll-to-zoom / drag-to-pan. The layout is recomputed for the harness-measured width/height. There is no useEffect and no engine handling in the component — <Plot> builds the plot engine, applies the children, and hands it back through onReady so the harness can drive export and backend switching (which preserves the current zoom/pan).

autowebglcanvassvg

Export PNG

fps 0frame 0 ms

PhyloTreeReact.tsx

PhyloTreeReact.tsx

import { schemeCategory10 } from "d3-scale-chromatic";
import { link as d3link, curveStepBefore } from "d3-shape";
import type { HierarchyPointNode, HierarchyPointLink } from "d3-hierarchy";
import { Plot, Layer, Points } from "@mapequation/d3gl/react";
import Example from "../../components/Example.js";
import { makeTree, type TreeNode } from "../shared/tree.js";
import { layoutRectangular, nodeXY } from "../shared/layout.js";

type PNode = HierarchyPointNode<TreeNode>;
type PLink = HierarchyPointLink<TreeNode>;

// Rectangular step links: a d3-shape link generator drawing straight into the d3gl context.
const gen = d3link<PLink, PNode>(curveStepBefore).x((d) => d.y).y((d) => d.x);

/**
 * The 64-tip rectangular phylogram, written declaratively in React with the
 * `<Plot>` / `<Layer>` / `<Points>` components from `@mapequation/d3gl/react`.
 * The layout is recomputed for the harness-measured `width`/`height`; everything
 * else is pure JSX — no `useEffect`, no DOM, no engine plumbing. `onReady` hands
 * the engine to the harness so it can drive export and backend switching.
 */
export default function PhyloTreeReact() {
  return (
    <Example width={720} height={460}>
      {({ backend, width, height, registerEngine }) => {
        const root = layoutRectangular(makeTree(64), width, height, "linear");
        return (
          <Plot
            width={width}
            height={height}
            backend={backend}
            zoom={[0.5, 40]}
            onReady={registerEngine}
          >
            <Layer
              name="links"
              data={root.links()}
              draw={(ctx, l) => {
                gen.context(ctx);
                gen(l);
              }}
              stroke="#555"
              lineWidth={0.8}
            />
            <Points
              name="nodes"
              data={root.leaves()}
              x={(n) => nodeXY(n, "rectangular")[0]}
              y={(n) => nodeXY(n, "rectangular")[1]}
              radius={2.6}
              fill={(n) => schemeCategory10[n.data.group % 10] ?? "#888"}
            />
          </Plot>
        );
      }}
    </Example>
  );
}

tree.ts

tree.ts

import type { RegionSet } from "./parsimony.js";

export interface TreeNode {
  name: string;
  group: number;
  /** Branch length = parent.time - this.time (time elapsed along the branch). */
  length: number;
  /** Age before present, in [0, 1]. Tips are at 0 (the present); the root is oldest. */
  time: number;
  children?: TreeNode[];
  /** Reconstructed ancestral range *set* (membership): species' present regions at tips,
   *  the most-parsimonious ancestral set at internal nodes. Set by calcMaximumParsimony(). */
  ranges?: RegionSet;
  /** Occurrence-count *distribution*: leaf counts summed up the tree, sorted by descending
   *  count. Sizes the pie wedges. Set by parsimony.aggregateClusters(). */
  clusters?: RegionSet;
  /** Number of subtended terminals (Fig. 3 branch-thickness metric). Set by
   *  parsimony.aggregateSpeciesCount(). */
  speciesCount?: number;
}

/** Deterministic LCG so trees are stable across renders (no Math.random in render path). */
function lcg(seed: number) { let s = seed >>> 0; return () => (s = (s * 1664525 + 1013904223) >>> 0) / 0xffffffff; }

/**
 * Build a roughly-balanced bifurcating tree with exactly `tips` leaves by
 * recursively splitting the tip count in two (with a little jitter so it looks
 * natural rather than perfectly symmetric). Branch lengths are random — a
 * phylogram, not a cladogram. Deterministic for a given seed.
 */
export function makeTree(tips: number, seed = 42): TreeNode {
  const rnd = lcg(seed);
  let leaf = 0;
  let groupCounter = 0;
  // Ultrametric ("dated") tree: every tip is at time 0 (the present); internal
  // nodes carry a divergence age younger than their parent's. `build` takes this
  // node's age and the parent's age (to compute the branch length).
  const build = (n: number, age: number, parentAge: number, group: number): TreeNode => {
    const length = parentAge - age;
    if (n <= 1) {
      leaf++;
      return { name: `tip_${leaf}`, group, length, time: 0 };
    }
    const jitter = (rnd() - 0.5) * n * 0.4;
    const left = Math.max(1, Math.min(n - 1, Math.round(n / 2 + jitter)));
    // Occasionally start a new colour group so tips cluster by clade.
    const g1 = rnd() < 0.22 ? ++groupCounter % 10 : group;
    const g2 = rnd() < 0.22 ? ++groupCounter % 10 : group;
    // Children diverge at younger ages (closer to the present); leaves sit at 0.
    const ageL = left <= 1 ? 0 : age * (0.3 + rnd() * 0.5);
    const ageR = n - left <= 1 ? 0 : age * (0.3 + rnd() * 0.5);
    return {
      name: "node", group, length, time: age,
      children: [build(left, ageL, age, g1), build(n - left, ageR, age, g2)],
    };
  };
  return build(Math.max(2, tips), 1, 1, 0);
}

layout.ts

layout.ts

import { hierarchy, cluster, type HierarchyPointNode } from "d3-hierarchy";
import { pointRadial } from "d3-shape";
import { scaleLinear, scaleSymlog } from "d3-scale";
import type { TreeNode } from "./tree.js";

export type LayoutMode = "rectangular" | "radial";
export type TimeScaleKind = "linear" | "log";

/**
 * Map a node's age (time before present; tips = 0, root = maxAge) to a position.
 * `present` is the coordinate for age 0, `root` the coordinate for the oldest node.
 * `scaleSymlog` is used for "log" because it is defined at 0 (unlike scaleLog).
 */
function timePosition(kind: TimeScaleKind, maxAge: number, present: number, root: number): (age: number) => number {
  if (kind === "log") {
    const s = scaleSymlog().domain([0, maxAge]).range([present, root]);
    return (age) => s(age);
  }
  const s = scaleLinear().domain([0, maxAge]).range([present, root]);
  return (age) => s(age);
}

/**
 * Rectangular dated phylogram. Uses d3's HierarchyPointNode coordinates directly:
 * `x` = vertical leaf-spacing axis, `y` = horizontal time axis. Tips (age 0) align at
 * the right (the present); the root (oldest) sits at the left.
 */
export function layoutRectangular(
  root: TreeNode,
  width: number,
  height: number,
  time: TimeScaleKind,
  pad = 40,
): HierarchyPointNode<TreeNode> {
  const h = cluster<TreeNode>().size([height - 2 * pad, 1])(hierarchy(root, (d) => d.children));
  const maxAge = h.data.time || 1;
  const pos = timePosition(time, maxAge, width - pad, pad); // age 0 → right, age max → left
  h.each((n) => {
    n.x += pad;
    n.y = pos(n.data.time);
  });
  return h;
}

/**
 * Radial dated phylogram. d3 convention: `x` = angle (0..2π), `y` = radius (= time).
 * Tips (age 0) sit on the outer rim, the root at the centre. Cartesian positions come
 * from `d3.pointRadial(x, y)` and are origin-centred — the caller centres the view with
 * a `translate(CX, CY)` transform (which also lets `d3.linkRadial()` work unmodified).
 */
export function layoutRadial(
  root: TreeNode,
  width: number,
  height: number,
  time: TimeScaleKind,
  pad = 30,
  angleExtent = 2 * Math.PI,
  angleStart = 0,
): HierarchyPointNode<TreeNode> {
  const h = cluster<TreeNode>().size([angleExtent, 1])(hierarchy(root, (d) => d.children));
  const maxAge = h.data.time || 1;
  // A partial fan (e.g. a half-circle "sunset") can use a larger radius: the leaves span only
  // part of the disc, so it is bounded by half the width and the full height rather than the
  // inscribed circle.
  const half = angleExtent < 2 * Math.PI - 1e-9;
  const R = (half ? Math.min(width / 2, height) : Math.min(width, height) / 2) - pad;
  const pos = timePosition(time, maxAge, R, 0); // age 0 → R (rim), age max → 0 (centre)
  h.each((n) => {
    n.x += angleStart; // shift the fan's angular origin (e.g. centre a half-fan on "north")
    n.y = pos(n.data.time);
  });
  return h;
}

/** Cartesian world coordinates for a node, for points / labels / hit-testing. */
export function nodeXY(n: HierarchyPointNode<TreeNode>, mode: LayoutMode): [number, number] {
  return mode === "radial" ? pointRadial(n.x, n.y) : [n.y, n.x];
}

parsimony.ts

parsimony.ts

import type { TreeNode } from "./tree.js";

/**
 * Fitch maximum-parsimony reconstruction of ancestral bioregion ranges (Fitch 1971),
 * plus occurrence-count aggregation. A region is `{clusterId, count}` — a bioregion id
 * and a count of species occurrences in it. Two distinct per-node products:
 *
 *  - `node.ranges`  — the reconstructed ancestral *range set* (membership): at a tip the
 *                     species' present regions, at an internal node the most-parsimonious
 *                     ancestral set. Counts here are presence (1); membership is the point.
 *  - `node.clusters`— the occurrence-count *distribution*: leaf counts summed up the tree,
 *                     sorted by descending count. Used to size the pie wedges.
 *
 * Fitch's two-phase algorithm. Preliminary phase (post-order): a leaf's set is its
 * regions; an internal node's set is the intersection of its children's sets, or — when
 * that is empty — their union (flagged `byUnion`). Final phase (pre-order), applying the
 * classic rules to each non-root node given its already-final parent set:
 *
 *   I.   If the preliminary set contains all regions of the parent's final set, go to II,
 *        else go to III.
 *   II.  (diminished ambiguity) Drop every region not in the parent's final set — i.e. keep
 *        the intersection with the parent. Done.
 *   III. If the preliminary set was formed by a union of its children's sets, go to IV,
 *        else go to V.
 *   IV.  (expanded ambiguity) Add every parent-final region not already present — i.e. the
 *        union with the parent. Done.
 *   V.   (encompassing ambiguity) Add any region not already present that is in BOTH the
 *        parent's final set AND at least one child's preliminary set. Done.
 *
 * NOTE: this fixes a bug in the bioregions1 reference, whose final phase tests
 * `node.byUnion` (always undefined — the flag lives on `node.clusters.byUnion`), so Rule IV
 * never fires and Rule V adds nothing. That produces incorrect sets for Fitch figs 2d/2f
 * (e.g. fig 2d node 0 wrongly collapses to {A,C} instead of the correct {A,C,G}). We follow
 * the written rules above, verified against the most-parsimonious reconstructions.
 */
export interface Region { clusterId: number; count: number; }
export interface RegionSet { totCount: number; clusters: Region[]; byUnion?: boolean; }
export type ClustersPerSpecies = Record<string, RegionSet>;

const idSet = (s: Region[]): Set<number> => new Set(s.map((r) => r.clusterId));
/** Regions of `a` also present in `b` (by clusterId), preserving `a`'s order. */
function intersectBy(a: Region[], b: Region[]): Region[] { const bi = idSet(b); return a.filter((r) => bi.has(r.clusterId)); }
/** Regions of `a`, then regions of `b` not already in `a`. */
function unionBy(a: Region[], b: Region[]): Region[] { const ai = idSet(a); return a.concat(b.filter((r) => !ai.has(r.clusterId))); }
const set = (clusters: Region[], byUnion?: boolean): RegionSet => ({ totCount: clusters.length, clusters, byUnion });

function visitPostOrder(node: TreeNode, cb: (n: TreeNode) => void): void {
  node.children?.forEach((c) => visitPostOrder(c, cb));
  cb(node);
}
function visitPreOrder(node: TreeNode, cb: (n: TreeNode, parent: TreeNode | null) => void, parent: TreeNode | null = null): void {
  cb(node, parent);
  node.children?.forEach((c) => visitPreOrder(c, cb, node));
}

/** Bottom-up pass: leaves presence-count their regions; each internal node takes the
 *  intersection of its children's sets, falling back to the union (flagged) if empty. */
export function calcMaximumParsimonyPreliminaryPhase(root: TreeNode, clustersPerSpecies: ClustersPerSpecies): TreeNode {
  visitPostOrder(root, (node) => {
    if (!node.children) {
      const cl = clustersPerSpecies[node.name];
      node.ranges = set(cl ? cl.clusters.map((r) => ({ clusterId: r.clusterId, count: 1 })) : []);
      return;
    }
    const childSets = node.children.map((c) => c.ranges!.clusters);
    let anc = childSets.reduce((acc, s) => intersectBy(acc, s));
    const byUnion = anc.length === 0;
    if (byUnion) anc = childSets.reduce((acc, s) => unionBy(acc, s), [] as Region[]);
    node.ranges = set(anc, byUnion);
  });
  return root;
}

/** Top-down pass applying Fitch's rules I–V (see the module doc). Root and leaves are
 *  already final. Children are visited after their parent, so a node's Rule-V lookup of
 *  its children's *preliminary* sets reads them before they are overwritten. */
export function calcMaximumParsimonyFinalPhase(root: TreeNode): TreeNode {
  visitPreOrder(root, (node, parent) => {
    if (!parent || !node.children) return;
    const prelim = node.ranges!.clusters;
    const pfinal = parent.ranges!.clusters;
    if (intersectBy(prelim, pfinal).length === pfinal.length) {
      node.ranges = set(intersectBy(prelim, pfinal)); // I → II (diminished)
    } else if (node.ranges!.byUnion) {
      node.ranges = set(unionBy(prelim, pfinal)); // III → IV (expanded)
    } else {
      // III → V (encompassing): add parent regions present in ≥1 child's preliminary set.
      const childrenUnion = node.children.map((c) => c.ranges!.clusters).reduce((a, b) => unionBy(a, b), [] as Region[]);
      node.ranges = set(unionBy(prelim, intersectBy(pfinal, childrenUnion)));
    }
  });
  return root;
}

export function calcMaximumParsimony(root: TreeNode, clustersPerSpecies: ClustersPerSpecies): TreeNode {
  calcMaximumParsimonyPreliminaryPhase(root, clustersPerSpecies);
  calcMaximumParsimonyFinalPhase(root);
  return root;
}

/**
 * Sum occurrence counts up the tree into `node.clusters`, sorted by descending count (ties
 * broken by clusterId). A leaf's counts come from the data; an internal node's count for a
 * region is the sum over its descendants. Mirrors bioregions1 `_aggregateClusters`, minus
 * the fraction-threshold "rest" grouping.
 */
export function aggregateClusters(root: TreeNode, clustersPerSpecies: ClustersPerSpecies): TreeNode {
  const sorted = (clusters: Region[], totCount: number): RegionSet => ({
    totCount,
    clusters: [...clusters].sort((a, b) => b.count - a.count || a.clusterId - b.clusterId),
  });
  visitPostOrder(root, (node) => {
    if (!node.children) {
      const cl = clustersPerSpecies[node.name];
      const clusters = cl ? cl.clusters.map((r) => ({ clusterId: r.clusterId, count: r.count })) : [];
      node.clusters = sorted(clusters, clusters.reduce((s, r) => s + r.count, 0));
      return;
    }
    const agg = new Map<number, number>();
    let totCount = 0;
    for (const child of node.children) {
      for (const r of child.clusters!.clusters) {
        agg.set(r.clusterId, (agg.get(r.clusterId) ?? 0) + r.count);
        totCount += r.count;
      }
    }
    node.clusters = sorted([...agg].map(([clusterId, count]) => ({ clusterId, count })), totCount);
  });
  return root;
}

/**
 * Post-order aggregation of `speciesCount` = number of subtended terminals (the Fig. 3
 * branch-thickness metric). With no `presence` map every leaf counts as one; with one,
 * a leaf counts only if present (absent species contribute 0).
 */
export function aggregateSpeciesCount(root: TreeNode, presence?: Record<string, number>): TreeNode {
  visitPostOrder(root, (node) => {
    if (!node.children) {
      node.speciesCount = presence ? (presence[node.name] ? 1 : 0) : 1;
      return;
    }
    node.speciesCount = node.children.reduce((sum, c) => sum + (c.speciesCount ?? 0), 0);
  });
  return root;
}