- 0 %
+ 0.00 %
- 0 %
+ 0.00 %
- 0 %
+ 0.00 %
- 0
+ 0.00
ms
- 0
+ 0.00
ms
- 0 %
+ 0.00 %
- 0
+ 0.00
ms
@@ -1300,7 +1300,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -1311,7 +1311,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
- 0
+ 0.00
ms
@@ -1448,7 +1448,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -1459,7 +1459,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
- 0
+ 0.00
ms
@@ -1596,7 +1596,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -1607,7 +1607,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
- 0
+ 0.00
ms
@@ -1744,7 +1744,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -1755,7 +1755,7 @@ exports[`SubtestsResultsView Component Tests should render the subtests results
- 0.963 %
+ 0.96 %
- 1.135 %
+ 1.14 %
- 0.963 %
+ 0.96 %
- 1.135 %
+ 1.14 %
- 0.963 %
+ 0.96 %
- 1.135 %
+ 1.14 %
- 0.963 %
+ 0.96 %
- 1.135 %
+ 1.14 %
- 0.963 %
+ 0.96 %
- 1.135 %
+ 1.14 %
- 0.963 %
+ 0.96 %
- 1.135 %
+ 1.14 %
- 0
+ 0.00
ms
- 0
+ 0.00
ms
- 0 %
+ 0.00 %
- 0
+ 0.00
ms
@@ -8017,7 +8017,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -8028,7 +8028,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
- 0
+ 0.00
ms
@@ -8165,7 +8165,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -8176,7 +8176,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
- 0
+ 0.00
ms
@@ -8313,7 +8313,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -8324,7 +8324,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
- 0
+ 0.00
ms
@@ -8461,7 +8461,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
class="mann-witney-browser-name cell"
role="cell"
>
- 0
+ 0.00
ms
@@ -8472,7 +8472,7 @@ exports[`SubtestsViewCompareOverTime Component Tests should render the subtests
{
+ it('returns null when either group has fewer than 2 observations', () => {
+ // A single run per side (common for some subtests) has no resampling
+ // variability and the BCa jackknife is undefined — must not yield NaNs.
+ expect(bootstrapMedianDiffCI([8.53], [8.73])).toBeNull();
+ expect(bootstrapMedianDiffCI([], [1, 2, 3])).toBeNull();
+ expect(bootstrapMedianDiffCI([1, 2, 3], [])).toBeNull();
+ expect(bootstrapMedianDiffCI([5], [6])).toBeNull();
+ });
+
+ it('returns a finite interval bracketing the observed difference', () => {
+ const base = [100, 102, 98, 101, 99, 103, 97, 100, 101, 99];
+ const newData = [106, 108, 104, 107, 105, 109, 103, 106, 107, 105];
+ const ci = bootstrapMedianDiffCI(base, newData);
+ expect(ci).not.toBeNull();
+ const { medianDiff, ciLow, ciHigh } = ci!;
+ expect(Number.isFinite(ciLow)).toBe(true);
+ expect(Number.isFinite(ciHigh)).toBe(true);
+ expect(ciLow).toBeLessThanOrEqual(medianDiff);
+ expect(ciHigh).toBeGreaterThanOrEqual(medianDiff);
+ });
+
+ it('is deterministic for a fixed seed', () => {
+ const base = [10, 12, 11, 13, 9, 14, 10, 11];
+ const newData = [14, 16, 15, 17, 13, 18, 14, 15];
+ expect(bootstrapMedianDiffCI(base, newData)).toEqual(
+ bootstrapMedianDiffCI(base, newData),
+ );
+ });
+});
diff --git a/src/common/testVersions/mannWhitney.tsx b/src/common/testVersions/mannWhitney.tsx
index 536e514d8..b46e4011d 100644
--- a/src/common/testVersions/mannWhitney.tsx
+++ b/src/common/testVersions/mannWhitney.tsx
@@ -411,8 +411,10 @@ export const mannWhitneyStrategy = {
const baseRuns = mwResult.base_runs ?? [];
const newRuns = mwResult.new_runs ?? [];
+ // A BCa CI needs at least 2 runs per side (see bootstrapMedianDiffCI);
+ // with fewer it returns null and we render no interval rather than NaNs.
const ci =
- baseRuns.length > 0 && newRuns.length > 0
+ baseRuns.length >= 2 && newRuns.length >= 2
? bootstrapMedianDiffCI(baseRuns, newRuns)
: null;
const rawUnit =
diff --git a/src/components/CompareResults/CommonGraph.tsx b/src/components/CompareResults/CommonGraph.tsx
index f539ad3d3..21d4a8699 100644
--- a/src/components/CompareResults/CommonGraph.tsx
+++ b/src/components/CompareResults/CommonGraph.tsx
@@ -13,12 +13,12 @@ import { useAppSelector } from '../../hooks/app';
import { Colors } from '../../styles/Colors';
import { getDisplayScale } from '../../utils/format';
import {
- areaFracs,
- assignLetters,
fftkde,
- fitModesFromKde,
+ fitKdePeakModes,
improvedSheatherJones,
+ matchModeLetters,
silvermansRule,
+ type GmmComponent,
} from '../../utils/kde.js';
// This computes the min, max from a list of numbers.
@@ -54,64 +54,84 @@ function computeMax(a?: number, b?: number) {
// it down to see structure.
const LARGE_BW_RATIO = 0.5;
const KDE_GRID_POINTS = 1024;
-const LABEL_ROW_PX = 16; // vertical space per stagger level
const KDE_TOP_BASE = 28;
-const KDE_HEIGHT = 155;
-const SCATTER_TOP_BASE = 250;
+const KDE_HEIGHT = 230;
+const SCATTER_TOP_BASE = 340;
const SCATTER_HEIGHT = 50;
-const CHART_HEIGHT_BASE = 340;
+const LEGEND_TOP = SCATTER_TOP_BASE - 18;
+const CHART_HEIGHT_BASE = 440;
-// Valley-depth threshold bounds for the mode-detection slider.
+// One color per mode letter (A–E). Same letter = same color across Base/New.
+//
+// Visual grammar:
+// color — mode identity (A=blue, B=amber, …)
+// markLine type — solid thick = Base (reference), dashed thinner = New
+// horizontal span — shows the mode's x-extent near the KDE baseline
+// vertical tick — thin guide to the exact peak position (no label)
+// Darkened from a brighter starting palette so label text (11px, on a
+// near-white background) clears WCAG comfortably — every entry here is
+// 8:1+ against white, well past the 4.5:1 AA minimum for normal text.
+// (The original palette ranged 3.6:1-6.9:1; some letters failed AA outright.)
+const MODE_FILL_COLORS = [
+ '#0D47A1', // A – blue
+ '#7A4009', // B – amber
+ '#0F4D29', // C – green
+ '#8B2318', // D – red
+ '#5A2D77', // E – purple
+];
+
+// Mode-sensitivity slider bounds.
const VT_MIN = 0.1;
const VT_MAX = 0.99;
const VT_STEP = 0.01;
+// Map the slider (0.1–0.99) to a valley-depth ratio for KDE peak merging.
+// valleyRatio = fraction of the lower peak's height the valley must DROP BELOW
+// for two peaks to be treated as separate modes.
+// vt→0 (left): high ratio (0.95) → most valleys merge → few modes
+// vt→1 (right): low ratio (0.05) → only very shallow valleys merge → more modes
+// midpoint 0.5: ratio ≈ 0.50 (valley must be <50% of lower peak)
+function sensitivityToValleyRatio(vt: number): number {
+ return Math.max(0.05, Math.min(0.95, 1 - vt));
+}
+
// Per-series mode summary, suitable both for chart overlays and the blurb.
type ModeInfo = {
peakLocs: number[];
+ boundaries: number[];
fracs: number[];
letters: string[];
+ components: GmmComponent[];
};
-function computeModeInfo(x: number[], y: number[], vt: number): ModeInfo {
- if (!x.length || !y.length) {
- return { peakLocs: [], fracs: [], letters: [] };
- }
- const { peakLocs, boundaries } = fitModesFromKde(x, y, vt);
- if (!peakLocs.length) {
- return { peakLocs: [], fracs: [], letters: [] };
- }
- const fracs = areaFracs(x, y, boundaries);
- const letters = assignLetters(peakLocs);
- return { peakLocs, fracs, letters };
+const EMPTY_MODE_INFO: ModeInfo = {
+ peakLocs: [],
+ boundaries: [],
+ fracs: [],
+ letters: [],
+ components: [],
+};
+
+// Detect modes from the KDE curve: find local maxima, then merge adjacent
+// peaks whose separating valley is shallower than valleyRatio * lowerPeak.
+// This gives exactly one mode per visual bump — the KDE bandwidth controls
+// resolution, so GMM over-fitting (multiple tight Gaussians per peak) can't
+// happen.
+function computeModeInfo(
+ kde: { x: ArrayLike
; y: ArrayLike } | null,
+ values: number[],
+ valleyRatio: number,
+): Omit | null {
+ if (!kde || values.length < 2) return null;
+ const result = fitKdePeakModes(kde.x, kde.y, values, { valleyRatio });
+ if (!result.peakLocs.length) return null;
+ return result;
}
// Stagger levels (0, 1, 2 …) for peak labels: peaks closer than ~13% of the
// x-span get bumped to different levels so their labels don't overlap. Ported
// from kde-widget.js's allPeaks.level pass; we use a fixed 13% threshold
// because the chart's pixel width isn't known inside useMemo.
-type PeakRef = {
- loc: number;
- seriesIdx: number;
- peakIdx: number;
- level: number;
-};
-
-function assignStaggerLevels(peaks: PeakRef[], xSpan: number): void {
- peaks.sort((a, b) => a.loc - b.loc);
- const threshold = xSpan * 0.2;
- for (let idx = 0; idx < peaks.length; idx++) {
- const used = new Set();
- for (let k = 0; k < idx; k++) {
- if (Math.abs(peaks[k].loc - peaks[idx].loc) < threshold) {
- used.add(peaks[k].level);
- }
- }
- let level = 0;
- while (used.has(level)) level++;
- peaks[idx].level = level;
- }
-}
function quantileSorted(sorted: number[], q: number): number {
const pos = (sorted.length - 1) * q;
@@ -158,6 +178,35 @@ function safeKde(values: number[], bw?: number) {
}
}
+// Measure a label's rendered width in pixels so it can be shifted left by
+// exactly that amount, keeping its right edge (rather than its left edge)
+// anchored to the peak — this is what prevents right-edge clipping near the
+// chart boundary.
+let measureCanvas: HTMLCanvasElement | null = null;
+function measureTextWidth(
+ text: string,
+ fontSize: number,
+ fontWeight: string,
+): number {
+ measureCanvas ??= document.createElement('canvas');
+ const ctx = measureCanvas.getContext('2d');
+ if (!ctx) return text.length * fontSize * 0.6;
+ ctx.font = `${fontWeight} ${fontSize}px sans-serif`;
+ return ctx.measureText(text).width;
+}
+
+// Scale a hex color's channels toward black. Used to give the New label a
+// visibly distinct (and only-higher-contrast) text color from Base's, rather
+// than relying on font-weight alone to tell them apart.
+function darkenHex(hex: string, factor: number): string {
+ const num = parseInt(hex.slice(1), 16);
+ const channel = (shift: number) =>
+ Math.round(((num >> shift) & 0xff) * factor)
+ .toString(16)
+ .padStart(2, '0');
+ return `#${channel(16)}${channel(8)}${channel(0)}`;
+}
+
// Linearly resample a uniform-grid KDE curve onto an arbitrary target x array.
// Outside the source range we return 0: each KDE's grid is padded so its
// density has already tapered to ≈0 at the edges.
@@ -334,44 +383,37 @@ function CommonGraph({
newScatterData,
min,
max,
+ sharedBw,
};
}, [baseValues, newValues, isSubtest, rawBandwidths, bwMultiplier]);
- // Mode detection (peaks, area fractions, label assignment, stagger levels)
- // lives in its own memo so it only re-runs when the threshold or the
- // underlying curves change — not on theme switch, scatter strip toggle, or
+ // Mode detection (Gaussian-mixture fit, label assignment, stagger levels)
+ // lives in its own memo so it only re-runs when the sensitivity slider or the
+ // underlying samples change — not on theme switch, scatter strip toggle, or
// unit changes. Uses localVt (the live slider position) so mode lines track
- // the thumb in real time. fitModesFromKde is cheap (array ops on a
- // pre-computed 1024-point grid), so running it on every drag pixel is fine.
+ // the thumb in real time.
const modes = useMemo(() => {
- const { bKde, nKde, sharedX, baseY, newY, min, max } = analysis;
-
- const baseModes = bKde
- ? computeModeInfo(sharedX, baseY, localVt)
- : { peakLocs: [], fracs: [], letters: [] };
- const newModes = nKde
- ? computeModeInfo(sharedX, newY, localVt)
- : { peakLocs: [], fracs: [], letters: [] };
-
- // Assign vertical stagger levels across all peaks so labels don't collide.
- const allPeaks: PeakRef[] = [];
- baseModes.peakLocs.forEach((loc, peakIdx) =>
- allPeaks.push({ loc, seriesIdx: 0, peakIdx, level: 0 }),
- );
- newModes.peakLocs.forEach((loc, peakIdx) =>
- allPeaks.push({ loc, seriesIdx: 1, peakIdx, level: 0 }),
+ const { bKde, nKde, min, max } = analysis;
+ const valleyRatio = sensitivityToValleyRatio(localVt);
+
+ const baseRaw = computeModeInfo(bKde, baseValues, valleyRatio);
+ const newRaw = computeModeInfo(nKde, newValues, valleyRatio);
+
+ const { baseLetters, newLetters } = matchModeLetters(
+ baseRaw?.peakLocs ?? [],
+ baseRaw?.fracs ?? [],
+ newRaw?.peakLocs ?? [],
+ newRaw?.fracs ?? [],
);
- const xSpan = max - min;
- if (xSpan > 0) assignStaggerLevels(allPeaks, xSpan);
- const levelLookup = new Map();
- for (const p of allPeaks) {
- levelLookup.set(`${p.seriesIdx}-${p.peakIdx}`, p.level);
- }
+ const baseModes: ModeInfo = baseRaw
+ ? { ...baseRaw, letters: baseLetters }
+ : EMPTY_MODE_INFO;
+ const newModes: ModeInfo = newRaw
+ ? { ...newRaw, letters: newLetters }
+ : EMPTY_MODE_INFO;
- const maxLevel =
- levelLookup.size > 0 ? Math.max(...levelLookup.values()) : 0;
- return { baseModes, newModes, levelLookup, maxLevel };
- }, [analysis, localVt]);
+ return { baseModes, newModes };
+ }, [analysis, localVt, baseValues, newValues]);
const option: EChartsOption = useMemo(() => {
const textColor =
@@ -383,19 +425,19 @@ function CommonGraph({
newScatterData,
min,
max,
+ sharedBw,
} = analysis;
- const { baseModes, newModes, levelLookup, maxLevel } = modes;
- const extraTop = maxLevel * LABEL_ROW_PX;
+ const { baseModes, newModes } = modes;
const kdeGrid = {
left: 70,
right: 70,
- top: KDE_TOP_BASE + extraTop,
+ top: KDE_TOP_BASE,
height: KDE_HEIGHT,
};
const scatterGrid = {
left: 70,
right: 70,
- top: SCATTER_TOP_BASE + extraTop,
+ top: SCATTER_TOP_BASE,
height: SCATTER_HEIGHT,
};
@@ -406,49 +448,222 @@ function CommonGraph({
const totalCount = baseValues.length + newValues.length;
const symbolSize = totalCount < 20 ? 14 : 10;
const tickFormatter = (value: number) => (value / scale).toFixed(decimals);
+ // Scale density y-values to match the display unit. KDE/GMM densities are
+ // computed in raw-unit space (e.g. per uWh); when the x-axis is shown in a
+ // larger unit (e.g. mWh, scale=1000), the probability density must be
+ // multiplied by scale so ∫f(x)dx = 1 still holds visually.
+ const scaleDensity = ([x, y]: [number, number]): [number, number] => [
+ x,
+ y * scale,
+ ];
+ const scaledBaseRunsDensity = baseRunsDensity.map(scaleDensity);
+ const scaledNewRunsDensity = newRunsDensity.map(scaleDensity);
- // Build the per-peak markLine overlays. Each is a dataless line series
- // whose markLine renders on its own. They share names with their parent
- // series ('Base' / 'New') so the legend's toggle action cascades to the
- // overlays automatically — clicking 'Base' hides every series named
- // 'Base', including the overlays. The legend itself only renders the two
- // entries listed in `legend.data` regardless of how many series share
- // those names.
+ // Build mode overlays: one horizontal span line per mode per series.
+ //
+ // Each span runs from the mode's left boundary to its right boundary, drawn
+ // near the KDE baseline. Base = thick solid, New = thinner dashed. Same
+ // color = same mode letter (matched across Base/New). A thin vertical tick
+ // marks the exact peak position; the label on the span carries the detail.
+ const maxDensity = Math.max(
+ ...scaledBaseRunsDensity.map(([, y]) => y),
+ ...scaledNewRunsDensity.map(([, y]) => y),
+ 0.001,
+ );
+ // Highest letter index across both series determines the topmost span row.
+ // We need to set the KDE y-axis max explicitly so ECharts doesn't clip spans
+ // (markLines don't participate in auto-range computation).
+ const maxLetterIdx = Math.max(
+ 0,
+ ...[...baseModes.letters, ...newModes.letters].map(
+ (l) => l.charCodeAt(0) - 65,
+ ),
+ );
+ const kdeYAxisMax = showModes
+ ? maxDensity * (1.08 + maxLetterIdx * 0.45 + 0.45)
+ : undefined;
const modeOverlays: EChartsOption['series'] = [];
- function pushOverlays(
- seriesIdx: 0 | 1,
+
+ // Peak-label placement: by default a label sits just to the left of the
+ // leftmost of its own tick and its matched Base/New counterpart's tick
+ // (so two nearby peaks share one consistent reference instead of each
+ // pushing left by its own, possibly very different, width). If that
+ // would run the label past the left axis, flip it to sit just to the
+ // right of its own tick instead.
+ const LABEL_GAP_PX = 6;
+ const chartWidthPx = chartInstanceRef.current?.getWidth?.() ?? 640;
+ const plotWidthPx = Math.max(1, chartWidthPx - 70 - 70);
+ const pxPerUnit = plotWidthPx / ((max - min) || 1);
+ const baseLocByLetter = new Map(
+ baseModes.letters.map((l, i) => [l, baseModes.peakLocs[i]]),
+ );
+ const newLocByLetter = new Map(
+ newModes.letters.map((l, i) => [l, newModes.peakLocs[i]]),
+ );
+
+ function pushModeOverlays(
seriesName: 'Base' | 'New',
- modes: ModeInfo,
- color: string,
+ modeInfo: ModeInfo,
+ seriesValues: number[],
+ xStart: number,
+ xEnd: number,
) {
- modes.peakLocs.forEach((loc, peakIdx) => {
- const level = levelLookup.get(`${seriesIdx}-${peakIdx}`) ?? 0;
+ if (!modeInfo.peakLocs.length) return;
+ const bounds = [xStart, ...modeInfo.boundaries, xEnd];
+ const isBase = seriesName === 'Base';
+ // Horizontal spans must not extend past the series' own actual data —
+ // xStart/xEnd come from the (padded) shared KDE grid, so clamp to the
+ // real leftmost/rightmost sample value.
+ const dataMin = Math.min(...seriesValues);
+ const dataMax = Math.max(...seriesValues);
+
+ modeInfo.peakLocs.forEach((loc, peakIdx) => {
+ const regionStart = Math.max(bounds[peakIdx], dataMin);
+ const regionEnd = Math.min(bounds[peakIdx + 1], dataMax);
+ const letterIdx = modeInfo.letters[peakIdx].charCodeAt(0) - 65;
+ // Each mode letter occupies its own row (stagger by letterIdx).
+ // Within each row, Base sits above New with a visible gap.
+ // Row height 0.16 keeps adjacent mode rows clearly separated.
+ const spanY = maxDensity * (1.08 + letterIdx * 0.45 + (isBase ? 0.20 : 0));
+ const modeColor = MODE_FILL_COLORS[letterIdx % MODE_FILL_COLORS.length];
+ const letter = modeInfo.letters[peakIdx];
+ const frac = modeInfo.fracs[peakIdx];
+ const valueStr = (loc / scale).toFixed(decimals);
+ const fracPct = Math.round(frac * 100);
+ const labelText = `${seriesName} ${letter}: ${valueStr} ${displayUnit} (${fracPct}%)`;
+
+ // Horizontal span showing the mode's x-extent (unlabeled — the
+ // combined label lives on the vertical tick below).
+ (modeOverlays as unknown[]).push({
+ name: seriesName,
+ type: 'line',
+ xAxisIndex: 0,
+ yAxisIndex: 0,
+ data: [],
+ z: 3,
+ markLine: {
+ silent: true,
+ symbol: 'none',
+ data: [
+ [{ coord: [regionStart, spanY] }, { coord: [regionEnd, spanY] }],
+ ],
+ lineStyle: {
+ color: modeColor,
+ type: isBase ? 'solid' : 'dashed',
+ width: isBase ? 2 : 1.5,
+ opacity: isBase ? 1 : 0.6,
+ },
+ label: { show: false },
+ },
+ });
+
+ // Vertical tick at the peak position, carrying the single combined
+ // label (series, letter, value, fraction) above its top.
+ const tickTopY = spanY + maxDensity * 0.001;
+ const labelFontSize = 11;
+ // Both bold — thin (normal-weight) text at 11px reads poorly even
+ // with good color contrast. Base and New are told apart by color
+ // (New is a darkened variant of the letter color), not weight.
+ const labelFontWeight = 'bold';
+ const labelColor = isBase ? modeColor : darkenHex(modeColor, 0.6);
+ const labelWidthPx = measureTextWidth(
+ labelText,
+ labelFontSize,
+ labelFontWeight,
+ );
+ const counterpartLoc = isBase
+ ? newLocByLetter.get(letter)
+ : baseLocByLetter.get(letter);
+ const anchorX = counterpartLoc !== undefined
+ ? Math.min(loc, counterpartLoc)
+ : loc;
+ const anchorPxFromLeft = (anchorX - min) * pxPerUnit;
+ const wouldCrossLeft =
+ anchorPxFromLeft - LABEL_GAP_PX - labelWidthPx < 0;
+ const deltaPx = (anchorX - loc) * pxPerUnit;
+ const labelAlign: 'left' | 'right' = wouldCrossLeft ? 'left' : 'right';
+ const labelOffset: [number, number] = wouldCrossLeft
+ ? [LABEL_GAP_PX, 4]
+ : [deltaPx - LABEL_GAP_PX, 4];
(modeOverlays as unknown[]).push({
name: seriesName,
type: 'line',
xAxisIndex: 0,
yAxisIndex: 0,
data: [],
+ z: 2,
markLine: {
silent: true,
symbol: 'none',
- data: [{ xAxis: loc }],
- lineStyle: { color, type: 'solid', width: 1.5 },
+ data: [[{ coord: [loc, 0] }, { coord: [loc, tickTopY] }]],
+ lineStyle: {
+ color: modeColor,
+ type: isBase ? 'solid' : 'dashed',
+ width: 1,
+ opacity: 0.5,
+ },
label: {
- formatter:
- `${seriesName} ${modes.letters[peakIdx]}: ` +
- `${tickFormatter(loc)} (${Math.round(modes.fracs[peakIdx] * 100)}%)`,
- distance: [0, level * 16],
- color,
- fontSize: 12,
+ show: true,
+ position: 'end',
+ align: labelAlign,
+ offset: labelOffset,
+ formatter: labelText,
+ color: labelColor,
+ opacity: 1,
+ fontSize: labelFontSize,
+ fontWeight: labelFontWeight,
+ backgroundColor: 'rgba(255,255,255,0.85)',
+ padding: [1, 3],
+ borderRadius: 2,
},
},
});
});
}
+
if (showModes) {
- pushOverlays(0, 'Base', baseModes, Colors.ChartBase);
- pushOverlays(1, 'New', newModes, Colors.ChartNew);
+ const xStart = analysis.sharedX[0] ?? min;
+ const xEnd = analysis.sharedX[analysis.sharedX.length - 1] ?? max;
+ pushModeOverlays('Base', baseModes, baseValues, xStart, xEnd);
+ pushModeOverlays('New', newModes, newValues, xStart, xEnd);
+
+ // Shift arrows: for each matched mode letter present in both series,
+ // draw a horizontal arrow from the Base peak to the New peak.
+ // Green = value decreased (good for lower-is-better), red = increased.
+ const basePeakByLetter = new Map(
+ baseModes.peakLocs.map((loc, i) => [baseModes.letters[i], loc]),
+ );
+ const newPeakByLetter = new Map(
+ newModes.peakLocs.map((loc, i) => [newModes.letters[i], loc]),
+ );
+ for (const [letter, baseLoc] of basePeakByLetter) {
+ const newLoc = newPeakByLetter.get(letter);
+ if (newLoc === undefined || newLoc === baseLoc) continue;
+ // Skip shifts smaller than the KDE bandwidth — below that scale the
+ // two peaks aren't distinguishable from smoothing noise, so drawing
+ // an arrow would overstate a difference that isn't significant.
+ if (sharedBw && Math.abs(newLoc - baseLoc) < sharedBw) continue;
+ const letterIdx = letter.charCodeAt(0) - 65;
+ // Place arrow between the New span row and the Base span row.
+ const arrowY = maxDensity * (1.08 + letterIdx * 0.45 + 0.10);
+ const arrowColor = newLoc < baseLoc ? '#1E8A4A' : '#C0392B';
+ (modeOverlays as unknown[]).push({
+ name: 'Base',
+ type: 'line',
+ xAxisIndex: 0,
+ yAxisIndex: 0,
+ data: [],
+ z: 4,
+ markLine: {
+ silent: true,
+ symbol: ['none', 'arrow'],
+ symbolSize: 8,
+ data: [[{ coord: [baseLoc, arrowY] }, { coord: [newLoc, arrowY] }]],
+ lineStyle: { color: arrowColor, width: 1.5 },
+ label: { show: false },
+ },
+ });
+ }
}
return {
@@ -486,10 +701,16 @@ function CommonGraph({
gridIndex: 0,
type: 'value',
min: 0,
+ max: kdeYAxisMax,
splitLine: { show: true, lineStyle: { color: '#eee' } },
axisLine: { show: true, lineStyle: { color: '#999' } },
axisTick: { show: false },
- axisLabel: { show: true, color: textColor, fontSize: 12 },
+ axisLabel: {
+ show: true,
+ color: textColor,
+ fontSize: 12,
+ formatter: (v: number) => (v === 0 ? '0' : v.toPrecision(2)),
+ },
},
{
gridIndex: 1,
@@ -570,7 +791,7 @@ function CommonGraph({
data: ['Base', 'New'],
// Sit below the centered x-axis unit label, between the KDE grid and
// the scatter strip, with a small gap above and below.
- top: 232,
+ top: LEGEND_TOP,
left: 'center',
itemHeight: 10,
itemWidth: 30,
@@ -582,7 +803,8 @@ function CommonGraph({
triggerLineEvent: true,
xAxisIndex: 0,
yAxisIndex: 0,
- data: baseRunsDensity,
+ z: 2,
+ data: scaledBaseRunsDensity,
showSymbol: false,
lineStyle: { width: 3, color: Colors.ChartBase },
itemStyle: { color: Colors.ChartBase },
@@ -594,7 +816,8 @@ function CommonGraph({
triggerLineEvent: true,
xAxisIndex: 0,
yAxisIndex: 0,
- data: newRunsDensity,
+ z: 2,
+ data: scaledNewRunsDensity,
showSymbol: false,
lineStyle: { width: 3, color: Colors.ChartNew },
itemStyle: { color: Colors.ChartNew },
@@ -700,53 +923,6 @@ function CommonGraph({
mb: 0.5,
}}
>
-
- Valley depth threshold
-
-
-
- :
-
- {/*
- MUI Slider exposes two events: `onChange` fires continuously during
- drag (we send it to local state for a smooth thumb), and
- `onChangeCommitted` fires once when the user releases (we push the
- final value up to the parent then). This is the moral equivalent of
- a debounce — the expensive consumer (`computeModeInfo`) runs once
- per drag instead of on every pixel of movement.
- */}
- setLocalVt(value)}
- onChangeCommitted={(_, value) => onVtChange(value)}
- aria-label='Valley depth threshold'
- sx={{ maxWidth: 240 }}
- />
-
- {Math.round(localVt * 100)}%
-
onShowModesChange(checked)}
/>
}
- label='Show modes'
- sx={{ ml: 1, '& .MuiFormControlLabel-label': { fontSize: 14 } }}
+ label='Modal analysis'
+ sx={{ '& .MuiFormControlLabel-label': { fontSize: 14 } }}
/>
+ {/*
+ The sensitivity slider only makes sense while modal analysis is on —
+ with it off the chart is just the KDE, so we hide the slider entirely
+ rather than showing a disabled control.
+ */}
+ {showModes && (
+ <>
+
+ Mode sensitivity
+
+
+
+ :
+
+ {/*
+ MUI Slider exposes two events: `onChange` fires continuously during
+ drag (we send it to local state for a smooth thumb), and
+ `onChangeCommitted` fires once when the user releases (we push the
+ final value up to the parent then). This is the moral equivalent of
+ a debounce — the expensive consumer (`computeModeInfo`) runs once
+ per drag instead of on every pixel of movement.
+ */}
+ setLocalVt(value)}
+ onChangeCommitted={(_, value) => onVtChange(value)}
+ aria-label='Mode sensitivity'
+ sx={{ maxWidth: 240 }}
+ />
+
+ {Math.round(localVt * 100)}%
+
+ >
+ )}
@@ -785,6 +1016,7 @@ function CommonGraph({
/>
)}
+
>
);
}
diff --git a/src/components/CompareResults/MannWhitneyCompareMetrics.tsx b/src/components/CompareResults/MannWhitneyCompareMetrics.tsx
index 1cead497e..6eb90d2c7 100644
--- a/src/components/CompareResults/MannWhitneyCompareMetrics.tsx
+++ b/src/components/CompareResults/MannWhitneyCompareMetrics.tsx
@@ -1,6 +1,7 @@
import { Box } from '@mui/material';
import { MannWhitneyResultsItem } from '../../types/state';
+import { getDisplayScale } from '../../utils/format';
import { getModeInterpretation } from '../../utils/helpers';
const METRIC_HEADERS = ['Metric', 'Base', 'New', 'Interpretation'];
@@ -46,6 +47,17 @@ export const MannWhitneyCompareMetrics = ({
min: null,
max: null,
};
+ const rawUnit =
+ result.base_measurement_unit ?? result.new_measurement_unit ?? 'ms';
+ const metricValues = [
+ baseMean, newMean, baseMedian, newMedian,
+ baseStandardDev, newStandardDev, baseMin, newMin, baseMax, newMax,
+ ].filter((v): v is number => v != null);
+ const { scale, displayUnit, decimals } = getDisplayScale(metricValues, rawUnit);
+ const fmtMetric = (v: number | null | undefined) =>
+ v != null ? (v / scale).toFixed(decimals) : 'N/A';
+ const unitLabel = displayUnit ? ` (${displayUnit})` : '';
+
const baseShapiroWilkPVal = result.shapiro_wilk_test_base?.pvalue ?? 'N/A';
const newShapiroWilkPVal = result.shapiro_wilk_test_new?.pvalue ?? 'N/A';
const baseShapiroWilkInterpretation =
@@ -85,12 +97,12 @@ export const MannWhitneyCompareMetrics = ({
>
- {METRIC_HEADERS.map((header) => (
+ {METRIC_HEADERS.map((header, i) => (
|
- {header}
+ {i === 0 ? `${header}${unitLabel}` : header}
|
))}
@@ -98,14 +110,14 @@ export const MannWhitneyCompareMetrics = ({
| Mean |
- {baseMean?.toFixed(2) ?? 'N/A'} |
- {newMean?.toFixed(2) ?? 'N/A'} |
+ {fmtMetric(baseMean)} |
+ {fmtMetric(newMean)} |
|
| Median |
- {baseMedian?.toFixed(2) ?? 'N/A'} |
- {newMedian?.toFixed(2) ?? 'N/A'} |
+ {fmtMetric(baseMedian)} |
+ {fmtMetric(newMedian)} |
|
@@ -116,20 +128,20 @@ export const MannWhitneyCompareMetrics = ({
| Standard Deviation |
- {baseStandardDev?.toFixed(2) ?? 'N/A'} |
- {newStandardDev?.toFixed(2) ?? 'N/A'} |
+ {fmtMetric(baseStandardDev)} |
+ {fmtMetric(newStandardDev)} |
|
| Min |
- {baseMin?.toFixed(2) ?? 'N/A'} |
- {newMin?.toFixed(2) ?? 'N/A'} |
+ {fmtMetric(baseMin)} |
+ {fmtMetric(newMin)} |
|
| Max |
- {baseMax?.toFixed(2) ?? 'N/A'} |
- {newMax?.toFixed(2) ?? 'N/A'} |
+ {fmtMetric(baseMax)} |
+ {fmtMetric(newMax)} |
|
diff --git a/src/utils/bootstrap-ci.ts b/src/utils/bootstrap-ci.ts
index eb6155b37..86c7cb209 100644
--- a/src/utils/bootstrap-ci.ts
+++ b/src/utils/bootstrap-ci.ts
@@ -136,6 +136,11 @@ export type BootstrapCI = {
* - Skewness: the interval may need to be asymmetric (a, the acceleration,
* estimated via leave-one-out jackknife).
*
+ * Returns null when either group has fewer than 2 observations: BCa relies on
+ * a leave-one-out jackknife for the acceleration term, which is undefined for a
+ * single observation (leaving it out yields an empty sample), and a one-point
+ * sample carries no resampling variability to form an interval from.
+ *
* @param base - baseline sample values
* @param newData - new/comparison sample values
* @param nIter - bootstrap resamples; 9999 is standard for BCa
@@ -148,7 +153,10 @@ export function bootstrapMedianDiffCI(
nIter: number = 9999,
alpha: number = 0.05,
seed: number = 42,
-): BootstrapCI {
+): BootstrapCI | null {
+ if (base.length < 2 || newData.length < 2) {
+ return null;
+ }
const rng = mulberry32(seed);
const baseArr = new Float64Array(base);
const newArr = new Float64Array(newData);
diff --git a/src/utils/format.ts b/src/utils/format.ts
index b0e2aec2f..f2e404fe4 100644
--- a/src/utils/format.ts
+++ b/src/utils/format.ts
@@ -2,7 +2,10 @@ const dateFormatter = new Intl.DateTimeFormat('en-US', {
dateStyle: 'medium',
});
-const numberFormatter = new Intl.NumberFormat('en-US');
+const numberFormatter = new Intl.NumberFormat('en-US', {
+ minimumFractionDigits: 2,
+ maximumFractionDigits: 2,
+});
export function formatDate(date: Date) {
return dateFormatter.format(date);
@@ -40,9 +43,14 @@ export function getDisplayScale(
if (rawUnit === 'ms') {
if (maxAbs >= 60000)
return { scale: 60000, displayUnit: 'min', decimals: 2 };
- if (maxAbs >= 1000) return { scale: 1000, displayUnit: 's', decimals: 2 };
+ if (maxAbs >= 10000) return { scale: 1000, displayUnit: 's', decimals: 2 };
return { scale: 1, displayUnit: 'ms', decimals: 2 };
}
+ if (rawUnit === 'uWh') {
+ if (maxAbs >= 1e6) return { scale: 1e6, displayUnit: 'Wh', decimals: 2 };
+ if (maxAbs >= 1e3) return { scale: 1e3, displayUnit: 'mWh', decimals: 2 };
+ return { scale: 1, displayUnit: 'uWh', decimals: 1 };
+ }
return { scale: 1, displayUnit: rawUnit, decimals: 1 };
}
diff --git a/src/utils/kde.d.ts b/src/utils/kde.d.ts
index 539097dc4..4949b236e 100644
--- a/src/utils/kde.d.ts
+++ b/src/utils/kde.d.ts
@@ -67,7 +67,9 @@ export declare function fitModesFromKde(
x: ArrayLike,
y: ArrayLike,
vt: number,
+ n?: number,
mpf?: number,
+ minSamples?: number,
mdf?: number,
): FitModesFromKdeResult;
export declare function areaFracs(
@@ -75,7 +77,55 @@ export declare function areaFracs(
y: ArrayLike,
boundaries: number[],
): number[];
+export type AdaptiveKDEResult = {
+ x: number[];
+ y: number[];
+ bandwidth: number;
+ localBandwidths: number[];
+};
+export declare function adaptiveKde(
+ data: number[],
+ opts?: {
+ numGridPoints?: number;
+ alpha?: number;
+ pilotBandwidth?: number;
+ bwMultiplier?: number;
+ },
+): AdaptiveKDEResult;
+export type GmmComponent = { mu: number; sigma: number; weight: number };
+export type GmmModeResult = {
+ peakLocs: number[];
+ boundaries: number[];
+ fracs: number[];
+ components: GmmComponent[];
+ nModes: number;
+};
+export declare function fitGmmModes(
+ data: number[],
+ opts?: {
+ penaltyScale?: number;
+ maxComponents?: number;
+ varFloorFrac?: number;
+ minSamples?: number;
+ },
+): GmmModeResult;
+export declare function gmmDensity(
+ components: GmmComponent[],
+ x: ArrayLike,
+): number[];
+export declare function fitKdePeakModes(
+ kdeX: ArrayLike,
+ kdeY: ArrayLike,
+ data: number[],
+ opts?: { valleyRatio?: number; minFrac?: number },
+): { peakLocs: number[]; boundaries: number[]; fracs: number[]; components: GmmComponent[] };
export declare function assignLetters(locs: number[]): string[];
+export declare function matchModeLetters(
+ baseLocs: number[],
+ baseFracs: number[],
+ newLocs: number[],
+ newFracs: number[],
+): { baseLetters: string[]; newLetters: string[] };
export type MatchModesResult = {
pairs: Array<[number, number]>;
ub: number[];
diff --git a/src/utils/kde.js b/src/utils/kde.js
index f91137ab5..eee6ffeb4 100644
--- a/src/utils/kde.js
+++ b/src/utils/kde.js
@@ -255,7 +255,11 @@ function gaussianKernel1D(x, bw) {
function gaussianPracticalSupport(bw, atol = 10e-5) {
const inner = atol * Math.sqrt(2 * Math.PI) * bw;
if (inner >= 1) return 3 * bw;
- return bw * Math.sqrt(-2 * Math.log(inner)) + 1e-3;
+ // Enforce at least 3σ: the atol formula finds where the kernel *value* hits
+ // atol, but for large bw the kernel peak is low so that threshold is reached
+ // at only ~1-2σ, truncating >20% of the probability mass and causing
+ // convolution ringing that looks like aliasing on the chart.
+ return Math.max(bw * Math.sqrt(-2 * Math.log(inner)) + 1e-3, 3 * bw);
}
// ---------------------------------------------------------------------------
// ISJ fixed-point function — port of KDEpy's _fixed_point
@@ -549,6 +553,100 @@ function fftkdeReflection(data, bw, weights, numGridPoints) {
return { x, y, bandwidth };
}
// ---------------------------------------------------------------------------
+// Adaptive (variable-bandwidth) KDE — Abramson's sample-point estimator
+//
+// Abramson, I. S. (1982): On bandwidth variation in kernel estimates — a
+// square root law. Ann. Stat. 10(4), 1217-1223.
+//
+// Why adaptive?
+// ------------
+// A single global bandwidth can't serve a distribution whose modes have very
+// different spreads. Performance data routinely does: a tight "fast path"
+// cluster sitting next to a diffuse "slow path" spread. A bandwidth wide
+// enough to coalesce the slow samples into a visible bump blurs the fast peak;
+// one narrow enough to keep the fast peak sharp shatters the slow samples into
+// isolated spikes that read as noise (and get filtered away), so the slow mode
+// vanishes entirely.
+//
+// The fix: give every sample its own bandwidth. Run a fixed-bandwidth pilot
+// estimate, then scale each sample's kernel by the local sparsity —
+// h_i = h0 · (g / f̃(x_i))^alpha, g = geometric mean of the pilot densities.
+// Samples in sparse regions (the slow path) get wider kernels and merge into a
+// single bump; samples in the dense fast cluster get narrower kernels and stay
+// sharp. alpha = 0.5 is Abramson's square-root law; larger alpha adapts more
+// aggressively. The slow mode then forms a real peak with a real valley, so
+// the existing mode-bucketing recovers it without over-smoothing the fast one.
+//
+// Cost is O(n·G) (direct evaluation, no FFT) because each sample carries a
+// different kernel width — fine for the small n of performance runs.
+//
+// @param data - raw sample values (at least 2)
+// @param opts - { numGridPoints=1024, alpha=0.5, pilotBandwidth?, bwMultiplier=1 }
+// pilotBandwidth overrides the internal ISJ pilot (e.g. to share
+// one pilot across base/new for comparability); bwMultiplier
+// scales the pilot (the smoothing slider).
+// @returns { x, y, bandwidth, localBandwidths } — bandwidth is the pilot h0
+// ---------------------------------------------------------------------------
+export function adaptiveKde(data, opts = {}) {
+ const numGridPoints = opts.numGridPoints ?? 1024;
+ const alpha = opts.alpha ?? 0.5;
+ const bwMultiplier = opts.bwMultiplier ?? 1;
+ const n = data.length;
+ if (n < 2) throw new Error('adaptiveKde requires at least 2 data points.');
+ // 1. Pilot fixed bandwidth — ISJ, falling back to Silverman's rule.
+ let h0 = opts.pilotBandwidth;
+ if (!(h0 > 0)) {
+ try {
+ h0 = improvedSheatherJones(data);
+ } catch {
+ h0 = silvermansRule(data);
+ }
+ if (!(h0 > 0)) h0 = silvermansRule(data);
+ }
+ h0 *= bwMultiplier;
+ // 2. Pilot density at each sample (direct fixed-bandwidth Gaussian KDE).
+ const pilot = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ let s = 0;
+ for (let j = 0; j < n; j++) s += gaussianKernel1D(data[i] - data[j], h0);
+ pilot[i] = s / n;
+ }
+ // 3. Per-sample bandwidth from the local sparsity (geometric-mean normalised).
+ let logSum = 0;
+ for (let i = 0; i < n; i++) {
+ logSum += Math.log(pilot[i] > 0 ? pilot[i] : Number.MIN_VALUE);
+ }
+ const g = Math.exp(logSum / n);
+ const hLocal = new Float64Array(n);
+ let hMax = h0;
+ for (let i = 0; i < n; i++) {
+ const ratio = g / (pilot[i] > 0 ? pilot[i] : Number.MIN_VALUE);
+ hLocal[i] = h0 * Math.pow(ratio, alpha);
+ if (hLocal[i] > hMax) hMax = hLocal[i];
+ }
+ // 4. Grid padded by the widest kernel's practical support so it isn't clipped.
+ const grid = autogrid1D(
+ data,
+ gaussianPracticalSupport(hMax),
+ numGridPoints,
+ 0.05,
+ );
+ // 5. Evaluate the adaptive density on the grid.
+ const y = new Array(numGridPoints).fill(0);
+ for (let gi = 0; gi < numGridPoints; gi++) {
+ const xg = grid[gi];
+ let s = 0;
+ for (let i = 0; i < n; i++) s += gaussianKernel1D(xg - data[i], hLocal[i]);
+ y[gi] = s / n;
+ }
+ return {
+ x: Array.from(grid),
+ y,
+ bandwidth: h0,
+ localBandwidths: Array.from(hLocal),
+ };
+}
+// ---------------------------------------------------------------------------
// argrelmax — port of scipy.signal.argrelmax(y, order=order)
//
// Returns indices i where y[i] is strictly greater than all y[i±j] for
@@ -621,23 +719,53 @@ export function areaFracs(x, y, boundaries) {
* Like fitKdeModes but takes (x, y) instead of raw data, so an interactive
* widget can re-fit modes on a slider change without recomputing the KDE.
*
- * Differences from fitKdeModes:
- * - Filters modes by KDE *area* (areaFracs) rather than the fraction of
- * raw data points falling in each bucket.
- * - Adds a minimum-separation guard: peaks closer than max(2, 5% of the x
- * range) collapse to the higher peak. This suppresses spurious modes
- * that KDE can produce on near-integer data (e.g. samples ∈ {0, 1}).
+ * How a mode is decided
+ * ---------------------
+ * Valley depth and integrated mass play distinct roles, mirroring how a human
+ * reads the curve:
+ * 1. Valley depth answers "are these two bumps *connected*?" — a shallow
+ * saddle means one mode with a wobble, a deep one means two basins.
+ * 2. Integrated mass (area under the bump, by the trapezoid rule) answers
+ * "is this bump *substantial enough to be a mode at all*?" This is the
+ * part the eye actually does — judging how much of the data lives under a
+ * hump — and it's what valley depth alone gets wrong. A single-sample
+ * outlier or a thin scatter of stragglers produces a tall, deep-valleyed
+ * spike that passes any depth test, yet carries almost no probability
+ * mass; a human dismisses it on sight. So a bump only survives as a mode
+ * when its area corresponds to at least `minSamples` real measurements
+ * (when the sample count n is known) and clears an absolute floor mdf.
+ *
+ * The height floor (mpf) is a companion to the mass test: it rejects wide,
+ * low scatter whose area would otherwise accumulate into a phantom mode even
+ * though it never rises into a real peak.
+ *
+ * Also keeps a minimum-separation guard: peaks closer than max(2, 5% of the x
+ * range) collapse to the higher peak, suppressing spurious modes that KDE can
+ * produce on near-integer data (e.g. samples ∈ {0, 1}).
*
* @param x - KDE x grid (uniform spacing)
* @param y - KDE density at each x
* @param vt - valley-depth threshold; the valley between two peaks must be
* shallower than vt × min(peak heights) for them to count as
* separate modes (0 = never split, 1 = always split)
- * @param mpf - minimum peak height as a fraction of the global max (default 0.05)
+ * @param n - number of underlying samples; enables the mass-to-sample-count
+ * floor (a mode must hold ≥ minSamples measurements). When omitted,
+ * only the absolute area floor mdf is applied.
+ * @param mpf - minimum peak height as a fraction of the global max (default 0.15)
+ * @param minSamples - minimum measurements a mode must contain, by integrated
+ * area × n (default 3); ignored when n is omitted
* @param mdf - minimum area fraction a mode must contain to be kept (default 0.05)
* @returns { peakLocs, boundaries } in the same units as x
*/
-export function fitModesFromKde(x, y, vt, mpf = 0.05, mdf = 0.05) {
+export function fitModesFromKde(
+ x,
+ y,
+ vt,
+ n,
+ mpf = 0.15,
+ minSamples = 3,
+ mdf = 0.05,
+) {
let yMax = 0;
for (let i = 0; i < y.length; i++) if (y[i] > yMax) yMax = y[i];
const peaks = argrelmax(y, 3).filter((i) => y[i] >= mpf * yMax);
@@ -646,8 +774,9 @@ export function fitModesFromKde(x, y, vt, mpf = 0.05, mdf = 0.05) {
for (let i = 1; i < y.length; i++) if (y[i] > y[gm]) gm = i;
return { peakLocs: [x[gm]], boundaries: [] };
}
- // Valley-depth filter — walk peaks left to right, keeping a peak only if
- // the valley between it and the previous kept peak is deep enough.
+ // Valley-depth filter — connectivity test. Walk peaks left to right, keeping
+ // a peak only if the valley between it and the previous kept peak is deep
+ // enough; otherwise merge it into whichever of the two is taller.
const good = [peaks[0]];
for (let k = 1; k < peaks.length; k++) {
const nxt = peaks[k];
@@ -669,28 +798,369 @@ export function fitModesFromKde(x, y, vt, mpf = 0.05, mdf = 0.05) {
}
return bs;
}
- // Area-fraction filter — drop modes whose KDE area is below mdf.
- const bs0 = computeBoundaries(good);
- const fr0 = areaFracs(x, y, bs0);
- const keep = good.map((_, i) => i).filter((i) => fr0[i] >= mdf);
- if (keep.length < 2) {
- const bp = good.reduce((a, b) => (y[a] > y[b] ? a : b));
+ // Mass filter — the integral test. Drop the least-massive mode one at a time,
+ // recomputing boundaries after each removal so the dropped bump's mass merges
+ // into the neighbour it actually belongs to (areaFracs re-buckets it), then
+ // re-test the remainder. A mode survives only when its area holds at least
+ // minSamples measurements (when n is known) and clears the absolute floor mdf.
+ const modes = good.slice();
+ while (modes.length > 1) {
+ const fr = areaFracs(x, y, computeBoundaries(modes));
+ let worst = 0;
+ for (let i = 1; i < fr.length; i++) if (fr[i] < fr[worst]) worst = i;
+ const enoughSamples = n === undefined || fr[worst] * n >= minSamples;
+ if (enoughSamples && fr[worst] >= mdf) break;
+ modes.splice(worst, 1);
+ }
+ if (modes.length < 2) {
+ const bp = modes.length
+ ? modes[0]
+ : good.reduce((a, b) => (y[a] > y[b] ? a : b));
return { peakLocs: [x[bp]], boundaries: [] };
}
- const fg = keep.map((i) => good[i]);
- const fb = computeBoundaries(fg);
- const locs = fg.map((i) => x[i]);
+ const locs = modes.map((i) => x[i]);
// Minimum-separation guard: KDE artefacts on near-integer data can put
// distinct peaks within a sample of each other. Collapse those to one.
const dataRange = x[x.length - 1] - x[0];
const minSep = Math.max(2, dataRange * 0.05);
for (let k = 1; k < locs.length; k++) {
if (locs[k] - locs[k - 1] < minSep) {
- const bestIdx = fg.reduce((a, b) => (y[a] > y[b] ? a : b));
+ const bestIdx = modes.reduce((a, b) => (y[a] > y[b] ? a : b));
return { peakLocs: [x[bestIdx]], boundaries: [] };
}
}
- return { peakLocs: locs, boundaries: fb };
+ return { peakLocs: locs, boundaries: computeBoundaries(modes) };
+}
+// ---------------------------------------------------------------------------
+// Gaussian-mixture mode detection (EM + BIC)
+//
+// Why a mixture model instead of carving the KDE at valley floors?
+// ----------------------------------------------------------------
+// Valley detection asks "where does the density dip", a one-parameter family
+// (the valley-depth slider) that often can't reproduce the grouping a human
+// would draw — in particular it gets the mode *count* wrong at every slider
+// setting, and a *diffuse* group (a slow path with high dispersion) has no
+// sharp KDE peak, so it shatters into ripples or is absorbed by the fast mode.
+//
+// A Gaussian mixture instead clusters the *samples*. Each mode is a component
+// with a weight (its probability mass — the integral the eye reads off), a
+// centre and a width, so:
+// - the number of modes is chosen by BIC (penalised likelihood), not a
+// threshold — fixing the "no slider value gives the right count" problem;
+// - a diffuse slow path is simply one wide-σ component, captured whole
+// rather than fragmented (and, because nothing smooths the density, a
+// valley between two real modes is never "filled in", unlike adaptive KDE);
+// - points are bucketed by maximum responsibility, so a boundary sits at the
+// Bayes-optimal crossing of two components, not at a KDE valley floor.
+//
+// EM is run once per K from a deterministic equal-count partition init (no
+// RNG, so results are reproducible), with a variance floor to avoid the
+// classic singular-component collapse on tightly-clustered samples.
+// ---------------------------------------------------------------------------
+function gaussPdf(x, mu, varr) {
+ return (
+ Math.exp(-((x - mu) * (x - mu)) / (2 * varr)) /
+ Math.sqrt(2 * Math.PI * varr)
+ );
+}
+// Init A — K contiguous equal-count chunks of the sorted data, one component
+// per chunk. Cheap and good when modes are comparable in size.
+function initChunks(sorted, K, varFloor) {
+ const n = sorted.length;
+ const means = new Array(K);
+ const vars = new Array(K);
+ const weights = new Array(K);
+ for (let k = 0; k < K; k++) {
+ const lo = Math.floor((k * n) / K);
+ const hi = Math.max(lo + 1, Math.floor(((k + 1) * n) / K));
+ let m = 0;
+ for (let i = lo; i < hi; i++) m += sorted[i];
+ m /= hi - lo;
+ let v = 0;
+ for (let i = lo; i < hi; i++) v += (sorted[i] - m) ** 2;
+ means[k] = m;
+ vars[k] = Math.max(v / (hi - lo), varFloor);
+ weights[k] = (hi - lo) / n;
+ }
+ return { means, vars, weights };
+}
+// Init B — deterministic k-means (farthest-point seeding + Lloyd), then derive
+// component params from the clusters. Robust when modes differ wildly in size
+// (the diffuse-slow case), where chunk init can straddle a cluster boundary.
+// Mirrors sklearn's default k-means init that our cross-check validated against.
+function initKmeans(data, sorted, K, varFloor) {
+ const n = data.length;
+ // Farthest-point seeding: start at the median, repeatedly add the point
+ // farthest (in value) from the nearest existing centre — deterministic, and
+ // it reliably plants one seed in each well-separated cluster.
+ const centres = [sorted[Math.floor(n / 2)]];
+ while (centres.length < K) {
+ let bestX = sorted[0];
+ let bestD = -1;
+ for (const x of sorted) {
+ let dmin = Infinity;
+ for (const c of centres) dmin = Math.min(dmin, Math.abs(x - c));
+ if (dmin > bestD) {
+ bestD = dmin;
+ bestX = x;
+ }
+ }
+ centres.push(bestX);
+ }
+ centres.sort((a, b) => a - b);
+ const assign = new Array(n).fill(0);
+ for (let iter = 0; iter < 50; iter++) {
+ let moved = false;
+ for (let i = 0; i < n; i++) {
+ let best = 0;
+ let bd = Infinity;
+ for (let k = 0; k < K; k++) {
+ const d = Math.abs(data[i] - centres[k]);
+ if (d < bd) {
+ bd = d;
+ best = k;
+ }
+ }
+ if (assign[i] !== best) {
+ assign[i] = best;
+ moved = true;
+ }
+ }
+ const sum = new Array(K).fill(0);
+ const cnt = new Array(K).fill(0);
+ for (let i = 0; i < n; i++) {
+ sum[assign[i]] += data[i];
+ cnt[assign[i]]++;
+ }
+ for (let k = 0; k < K; k++) if (cnt[k]) centres[k] = sum[k] / cnt[k];
+ if (!moved && iter > 0) break;
+ }
+ const means = centres.slice();
+ const vars = new Array(K).fill(0);
+ const weights = new Array(K).fill(0);
+ const cnt = new Array(K).fill(0);
+ for (let i = 0; i < n; i++) {
+ cnt[assign[i]]++;
+ vars[assign[i]] += (data[i] - means[assign[i]]) ** 2;
+ }
+ for (let k = 0; k < K; k++) {
+ weights[k] = cnt[k] / n;
+ vars[k] = Math.max(cnt[k] ? vars[k] / cnt[k] : varFloor, varFloor);
+ }
+ return { means, vars, weights };
+}
+// One EM run at fixed K from a supplied init {means, vars, weights}.
+function emGmm1D(data, init, varFloor, maxIter = 300) {
+ const n = data.length;
+ const K = init.means.length;
+ const means = init.means.slice();
+ const vars = init.vars.slice();
+ const weights = init.weights.slice();
+ const resp = Array.from({ length: n }, () => new Array(K).fill(0));
+ let prevLL = -Infinity;
+ let logL = -Infinity;
+ for (let iter = 0; iter < maxIter; iter++) {
+ // E-step
+ logL = 0;
+ for (let i = 0; i < n; i++) {
+ let denom = 0;
+ for (let k = 0; k < K; k++) {
+ const p = weights[k] * gaussPdf(data[i], means[k], vars[k]);
+ resp[i][k] = p;
+ denom += p;
+ }
+ if (denom > 0) {
+ for (let k = 0; k < K; k++) resp[i][k] /= denom;
+ logL += Math.log(denom);
+ } else {
+ for (let k = 0; k < K; k++) resp[i][k] = 1 / K;
+ logL -= 700;
+ }
+ }
+ // M-step
+ for (let k = 0; k < K; k++) {
+ let Nk = 0;
+ let mu = 0;
+ for (let i = 0; i < n; i++) {
+ Nk += resp[i][k];
+ mu += resp[i][k] * data[i];
+ }
+ if (Nk < 1e-12) continue;
+ mu /= Nk;
+ let v = 0;
+ for (let i = 0; i < n; i++) v += resp[i][k] * (data[i] - mu) ** 2;
+ means[k] = mu;
+ vars[k] = Math.max(v / Nk, varFloor);
+ weights[k] = Nk / n;
+ }
+ if (
+ iter > 0 &&
+ Math.abs(logL - prevLL) <= 1e-7 * (Math.abs(prevLL) + 1e-12)
+ )
+ break;
+ prevLL = logL;
+ }
+ return { means, vars, weights, logL };
+}
+// x in (a.mu, b.mu) where the two weighted components cross — the Bayes
+// decision boundary between adjacent modes. Scanned (robust) with linear
+// interpolation at the sign change.
+function componentCrossing(a, b) {
+ const lo = a.mu;
+ const hi = b.mu;
+ const va = Math.max(a.sigma * a.sigma, 1e-12);
+ const vb = Math.max(b.sigma * b.sigma, 1e-12);
+ const steps = 256;
+ let prev =
+ a.weight * gaussPdf(lo, a.mu, va) - b.weight * gaussPdf(lo, b.mu, vb);
+ for (let i = 1; i <= steps; i++) {
+ const x = lo + ((hi - lo) * i) / steps;
+ const diff =
+ a.weight * gaussPdf(x, a.mu, va) - b.weight * gaussPdf(x, b.mu, vb);
+ if (prev > 0 !== diff > 0) {
+ const xPrev = lo + ((hi - lo) * (i - 1)) / steps;
+ const t = prev / (prev - diff);
+ return xPrev + (x - xPrev) * t;
+ }
+ prev = diff;
+ }
+ return (lo + hi) / 2;
+}
+/**
+ * Detect modes by fitting a 1-D Gaussian mixture to the raw samples and
+ * selecting the number of components with BIC.
+ *
+ * @param data - raw sample values
+ * @param opts - { penaltyScale=1, maxComponents=5 }
+ * penaltyScale multiplies the BIC complexity penalty: >1 favours
+ * fewer modes, <1 favours more (the UI "mode sensitivity" knob).
+ * @returns { peakLocs, boundaries, fracs, components, nModes }
+ * peakLocs : component centres (sorted ascending)
+ * boundaries: Bayes crossings between adjacent components
+ * fracs : component weights (probability mass per mode)
+ * components: [{ mu, sigma, weight }] for drawing the mixture density
+ */
+export function fitGmmModes(data, opts = {}) {
+ const penaltyScale = opts.penaltyScale ?? 1;
+ const maxComponents = opts.maxComponents ?? 5;
+ const varFloorFrac = opts.varFloorFrac ?? 0.01;
+ const minSamples = opts.minSamples ?? 1.5;
+ const n = data.length;
+ const sorted = [...data].sort((a, b) => a - b);
+ if (n < 4 || sorted[0] === sorted[n - 1]) {
+ // Too few (or identical) samples to fit a mixture: report one mode at the
+ // mean, with the sample spread as its width (zero only when all equal).
+ const mean = n ? data.reduce((a, b) => a + b, 0) / n : 0;
+ const sigma = n
+ ? Math.sqrt(data.reduce((s, x) => s + (x - mean) ** 2, 0) / n)
+ : 0;
+ return {
+ peakLocs: [mean],
+ boundaries: [],
+ fracs: [1],
+ components: [{ mu: mean, sigma, weight: 1 }],
+ nModes: 1,
+ };
+ }
+ const span = sorted[n - 1] - sorted[0];
+ // Variance floor: the larger of a span fraction and the measurement
+ // resolution. The resolution term stops GMM from fitting one razor-thin
+ // component per rounding level on quantised data (e.g. ms-rounded timings) —
+ // a failure mode the scipy/sklearn reference shares, so it must be handled at
+ // the data level. resolution = smallest gap between distinct sample values.
+ let resolution = span;
+ for (let i = 1; i < n; i++) {
+ const gap = sorted[i] - sorted[i - 1];
+ if (gap > 0 && gap < resolution) resolution = gap;
+ }
+ const varFloor = Math.max(
+ (span * varFloorFrac) ** 2,
+ resolution * resolution,
+ 1e-12,
+ );
+ const Kmax = Math.max(1, Math.min(maxComponents, Math.ceil(n / 4)));
+ let best = null;
+ for (let K = 1; K <= Kmax; K++) {
+ // Try both inits and keep the higher-likelihood fit — this matches the
+ // robustness of sklearn's multi-restart k-means init (cross-checked) and
+ // avoids the bad local optima a single init falls into when modes differ
+ // greatly in size.
+ const inits =
+ K === 1
+ ? [initChunks(sorted, K, varFloor)]
+ : [
+ initChunks(sorted, K, varFloor),
+ initKmeans(data, sorted, K, varFloor),
+ ];
+ for (const init of inits) {
+ const fit = emGmm1D(data, init, varFloor, 300);
+ const params = 3 * K - 1; // K means + K variances + (K-1) free weights
+ const bic = -2 * fit.logL + penaltyScale * params * Math.log(n);
+ if (!best || bic < best.bic - 1e-9) best = { ...fit, K, bic };
+ }
+ }
+ // Components → modes: merge near-duplicate centres first, then drop
+ // near-empty ones. Merging before filtering lets two small close components
+ // combine their mass before being evaluated against minSamples.
+ let comps = best.means
+ .map((mu, k) => ({
+ mu,
+ sigma: Math.sqrt(best.vars[k]),
+ weight: best.weights[k],
+ }))
+ .filter((c) => Number.isFinite(c.mu))
+ .sort((a, b) => a.mu - b.mu);
+ if (!comps.length) {
+ const med = sorted[Math.floor(n / 2)];
+ comps = [{ mu: med, sigma: Math.sqrt(varFloor), weight: 1 }];
+ }
+ // minSep is relative to the data span so it is unit-independent.
+ const minSep = span * 0.02;
+ const merged = [comps[0]];
+ for (let i = 1; i < comps.length; i++) {
+ const p = merged[merged.length - 1];
+ const c = comps[i];
+ if (c.mu - p.mu < minSep) {
+ const w = p.weight + c.weight;
+ const mu = (p.mu * p.weight + c.mu * c.weight) / w;
+ const v =
+ (p.weight * (p.sigma ** 2 + (p.mu - mu) ** 2) +
+ c.weight * (c.sigma ** 2 + (c.mu - mu) ** 2)) /
+ w;
+ merged[merged.length - 1] = { mu, sigma: Math.sqrt(v), weight: w };
+ } else merged.push(c);
+ }
+ const filtered = merged.filter((c) => c.weight * n >= minSamples);
+ if (!filtered.length) {
+ const med = sorted[Math.floor(n / 2)];
+ filtered.push({ mu: med, sigma: Math.sqrt(varFloor), weight: 1 });
+ }
+ const wSum = filtered.reduce((s, c) => s + c.weight, 0) || 1;
+ for (const c of filtered) c.weight /= wSum;
+ const boundaries = [];
+ for (let i = 0; i < filtered.length - 1; i++)
+ boundaries.push(componentCrossing(filtered[i], filtered[i + 1]));
+ return {
+ peakLocs: filtered.map((c) => c.mu),
+ boundaries,
+ fracs: filtered.map((c) => c.weight),
+ components: filtered,
+ nModes: filtered.length,
+ };
+}
+/**
+ * Evaluate a Gaussian-mixture density on a grid, for drawing the fitted model
+ * over the KDE. Returns one density value per x.
+ */
+export function gmmDensity(components, x) {
+ const out = new Array(x.length);
+ for (let i = 0; i < x.length; i++) {
+ let s = 0;
+ for (const c of components)
+ s += c.weight * gaussPdf(x[i], c.mu, Math.max(c.sigma * c.sigma, 1e-12));
+ out[i] = s;
+ }
+ return out;
}
/**
* Assign single-letter labels to modes, with A = lowest peak location.
@@ -711,6 +1181,186 @@ export function assignLetters(locs) {
});
return out;
}
+
+/**
+ * Mode detection driven by the KDE curve rather than by GMM.
+ *
+ * Finds every local maximum in the KDE grid, then merges adjacent peaks whose
+ * separating valley is too shallow: if the valley height exceeds
+ * `valleyRatio * min(peakA.y, peakB.y)`, the two peaks are considered one
+ * visual mode and the shorter one is absorbed into the taller. The default
+ * ratio 0.5 means the valley must drop below half the lower peak's height to
+ * count as a genuine separation.
+ *
+ * This gives exactly one mode per visual bump in the KDE — no over-fitting,
+ * because the KDE bandwidth already defines the resolution at which bumps are
+ * distinguishable.
+ *
+ * @param {ArrayLike} kdeX - x coordinates of the KDE grid
+ * @param {ArrayLike} kdeY - corresponding density values
+ * @param {number[]} data - raw sample values (for fraction counts)
+ * @param {object} [opts]
+ * @param {number} [opts.valleyRatio=0.5]
+ * @returns {{ peakLocs: number[], boundaries: number[], fracs: number[], components: object[] }}
+ */
+export function fitKdePeakModes(kdeX, kdeY, data, opts = {}) {
+ const { valleyRatio = 0.5, minFrac = 0.05 } = opts;
+ const n = kdeX.length;
+ if (n < 3 || !data.length)
+ return { peakLocs: [], boundaries: [], fracs: [], components: [] };
+
+ // Collect all local maxima (strict on right so plateaux pick their left edge)
+ let peaks = [];
+ for (let i = 1; i < n - 1; i++) {
+ if (kdeY[i] >= kdeY[i - 1] && kdeY[i] > kdeY[i + 1]) {
+ peaks.push({ x: kdeX[i], y: kdeY[i], idx: i });
+ }
+ }
+ if (!peaks.length) {
+ // No local max (e.g. monotone or all-equal) → treat the global max as one peak
+ let mi = 0;
+ for (let i = 1; i < n; i++) if (kdeY[i] > kdeY[mi]) mi = i;
+ peaks = [{ x: kdeX[mi], y: kdeY[mi], idx: mi }];
+ }
+
+ // Merge adjacent peaks where the separating valley is too shallow.
+ let changed = true;
+ while (changed && peaks.length > 1) {
+ changed = false;
+ for (let i = 0; i < peaks.length - 1; i++) {
+ const pA = peaks[i];
+ const pB = peaks[i + 1];
+ let valleyY = Infinity;
+ for (let j = pA.idx + 1; j < pB.idx; j++) {
+ if (kdeY[j] < valleyY) valleyY = kdeY[j];
+ }
+ if (valleyY > valleyRatio * Math.min(pA.y, pB.y)) {
+ // Shallow valley — keep the taller peak, drop the shorter
+ peaks.splice(i, 2, pA.y >= pB.y ? pA : pB);
+ changed = true;
+ break;
+ }
+ }
+ }
+
+ // Boundaries: x-position of the KDE minimum between each adjacent peak pair
+ const boundaries = [];
+ for (let i = 0; i < peaks.length - 1; i++) {
+ let valleyY = Infinity;
+ let valleyX = (peaks[i].x + peaks[i + 1].x) / 2;
+ for (let j = peaks[i].idx + 1; j < peaks[i + 1].idx; j++) {
+ if (kdeY[j] < valleyY) {
+ valleyY = kdeY[j];
+ valleyX = kdeX[j];
+ }
+ }
+ boundaries.push(valleyX);
+ }
+
+ // Fracs: fraction of data points that fall in each mode's region
+ const initialRegionBounds = [-Infinity, ...boundaries, Infinity];
+ const fracs = peaks.map((_, i) => {
+ const lo = initialRegionBounds[i];
+ const hi = initialRegionBounds[i + 1];
+ return data.reduce((c, v) => c + (v > lo && v <= hi ? 1 : 0), 0) / data.length;
+ });
+
+ // Drop modes whose data fraction is below minFrac — they're outlier bumps,
+ // not real modes. Merge into the nearest neighbour by dropping the shared
+ // boundary; keep the taller peak as the representative.
+ let mergedPeaks = [...peaks];
+ let mergedBoundaries = [...boundaries];
+ let mergedFrags = [...fracs];
+ let mchanged = true;
+ while (mchanged && mergedPeaks.length > 1) {
+ mchanged = false;
+ let minIdx = -1;
+ let minFracVal = Infinity;
+ for (let i = 0; i < mergedFrags.length; i++) {
+ if (mergedFrags[i] < minFracVal) { minFracVal = mergedFrags[i]; minIdx = i; }
+ }
+ if (minFracVal >= minFrac) break;
+ // Merge minIdx with its nearest neighbour (left preferred, else right)
+ const mergeWith = minIdx > 0 ? minIdx - 1 : 1;
+ const lo = Math.min(minIdx, mergeWith);
+ const hi = Math.max(minIdx, mergeWith);
+ const keepPeak = mergedPeaks[lo].y >= mergedPeaks[hi].y ? mergedPeaks[lo] : mergedPeaks[hi];
+ mergedPeaks.splice(lo, 2, keepPeak);
+ mergedBoundaries.splice(lo, 1);
+ // Recompute fracs for merged region
+ const regionBounds2 = [-Infinity, ...mergedBoundaries, Infinity];
+ mergedFrags = mergedPeaks.map((_, i) => {
+ const lo2 = regionBounds2[i], hi2 = regionBounds2[i + 1];
+ return data.reduce((c, v) => c + (v > lo2 && v <= hi2 ? 1 : 0), 0) / data.length;
+ });
+ mchanged = true;
+ }
+ peaks = mergedPeaks;
+ boundaries.length = 0;
+ mergedBoundaries.forEach((b) => boundaries.push(b));
+ fracs.length = 0;
+ mergedFrags.forEach((f) => fracs.push(f));
+ const regionBounds = [-Infinity, ...boundaries, Infinity];
+
+ // Lightweight component descriptors (mean/sigma from the region's data)
+ // for downstream compatibility with matchModeLetters / gmmDensity callers.
+ const components = peaks.map((p, i) => {
+ const lo = regionBounds[i];
+ const hi = regionBounds[i + 1];
+ const region = data.filter((v) => v > lo && v <= hi);
+ const mean =
+ region.length ? region.reduce((s, v) => s + v, 0) / region.length : p.x;
+ const variance =
+ region.length > 1
+ ? region.reduce((s, v) => s + (v - mean) ** 2, 0) / region.length
+ : 1;
+ return { mean, sigma: Math.sqrt(variance), weight: fracs[i] };
+ });
+
+ return { peakLocs: peaks.map((p) => p.x), boundaries, fracs, components };
+}
+
+/**
+ * Cross-series mode matching. Assigns letters to base modes by position
+ * (A = leftmost), then uses `matchModes` to pair new modes with base modes.
+ * Matched new modes inherit the base letter; unmatched new modes get the next
+ * unused letter. This lets the caller see "Base A shifted to New A" vs.
+ * "New B is a novel mode".
+ *
+ * @param {number[]} baseLocs - peak positions for the base series
+ * @param {number[]} baseFracs - weight fractions for the base series
+ * @param {number[]} newLocs - peak positions for the new series
+ * @param {number[]} newFracs - weight fractions for the new series
+ * @returns {{ baseLetters: string[], newLetters: string[] }}
+ */
+export function matchModeLetters(baseLocs, baseFracs, newLocs, newFracs) {
+ const baseLetters = assignLetters(baseLocs);
+ if (!newLocs.length) return { baseLetters, newLetters: [] };
+ if (!baseLocs.length) return { baseLetters: [], newLetters: assignLetters(newLocs) };
+
+ const { pairs, un } = matchModes(baseLocs, baseFracs, newLocs, newFracs);
+
+ // Build a map from new-mode index → letter.
+ const newLetters = new Array(newLocs.length);
+ for (const [bi, ni] of pairs) {
+ newLetters[ni] = baseLetters[bi];
+ }
+
+ // Unmatched new modes get fresh letters not already used.
+ const usedLetters = new Set(baseLetters);
+ let nextCode = 65;
+ function freshLetter() {
+ while (usedLetters.has(String.fromCharCode(nextCode))) nextCode++;
+ const l = String.fromCharCode(nextCode++);
+ usedLetters.add(l);
+ return l;
+ }
+ for (const ni of un) {
+ newLetters[ni] = freshLetter();
+ }
+
+ return { baseLetters, newLetters };
+}
function popcount(x) {
let c = 0;
let v = x;