feat: add xyRollingCircle that provide the rolling circle transform +…

.npmrc

@@ -1 +1 @@
-package-lock=false
+ignore-scripts=true

benchmark/xyRollingCircle.ts

@@ -0,0 +1,45 @@
+/* eslint-disable no-console */
+import Benchmark from 'benchmark';
+import { XSadd } from 'ml-xsadd';
+
+import { xyRollingCircle } from '../src/xy/xyRollingCircle.ts';
+import { xyRollingCircleTransform } from '../src/xy/xyRollingCircleTransform.ts';
+
+const size = 100000;
+
+// Deterministic, reproducible input: a slowly varying baseline plus noise.
+const { random } = new XSadd(42);
+const x = new Float64Array(size);
+const y = new Float64Array(size);
+for (let i = 0; i < size; i++) {
+  x[i] = i;
+  y[i] = Math.sin(i / 500) * 10 + random();
+}
+const data = { x, y };
+
+// What the touched points cost: xyRollingCircle does the same scan as the
+// transform, plus one extra pass over each window to mark the supporting points.
+for (const radius of [1.5, 50]) {
+  const { points } = xyRollingCircle(data, { radius });
+  console.log(
+    `\nradius ${radius}, n = ${size}, ${points.length} touched points`,
+  );
+
+  new Benchmark.Suite()
+    .add('transform', () => xyRollingCircleTransform(data, { radius }), {
+      minSamples: 30,
+    })
+    .add('with points', () => xyRollingCircle(data, { radius }), {
+      minSamples: 30,
+    })
+    .on('cycle', (event: Benchmark.Event) => {
+      const { target } = event;
+      const stats = target.stats;
+      if (!stats) return;
+      const nsPerPoint = (stats.mean * 1e9) / size;
+      console.log(
+        `${target.name?.padEnd(10)} ${(1 / stats.mean).toFixed(1).padStart(7)} ops/sec | ${nsPerPoint.toFixed(1)} ns/point | rme ${stats.rme.toFixed(1)}% | ${stats.sample.length} samples`,
+      );
+    })
+    .run();
+}

benchmark/xyRollingCircleMarking.ts

@@ -0,0 +1,124 @@
+/* eslint-disable no-console */
+/*
+ * Marking the touched points needs `y[j] + yScale * sqrt(radius2 - deltaX²) >= limit`.
+ * That comparison can be squared into `radius2 - deltaX² >= needed²` to get rid
+ * of the square root, which looks like an obvious win. It is not:
+ *
+ *              node 24              bun 1.3
+ * radius  1.5  15.5 -> 16.7 ns/pt   11.0 -> 11.3 ns/pt
+ * radius 50   247.6 -> 256.8 ns/pt  177.2 -> 158.2 ns/pt
+ *
+ * so it is a wash: slower on V8 in both cases, slower on JSC for a narrow window
+ * and faster only for a wide one. The square root stays, it is simpler.
+ *
+ * Measure it here and not through `xyRollingCircle`: that one allocates one
+ * object per touched point (~70000 of them at radius 1.5), and the resulting GC
+ * pressure makes its timings bimodal, 26 or 37 ns/point for the very same code.
+ */
+import Benchmark from 'benchmark';
+import { XSadd } from 'ml-xsadd';
+
+const size = 100000;
+const RELATIVE_TOLERANCE = 1e-10;
+
+// Deterministic, reproducible input: a slowly varying baseline plus noise.
+const { random } = new XSadd(42);
+const x = new Float64Array(size);
+const y = new Float64Array(size);
+for (let i = 0; i < size; i++) {
+  x[i] = i;
+  y[i] = Math.sin(i / 500) * 10 + random();
+}
+
+/**
+ * Marks the touched points with one square root per point.
+ * @param radius - radius of the circle.
+ * @returns the mask of the touched points.
+ */
+function markWithSqrt(radius: number): Uint8Array {
+  const radius2 = radius * radius;
+  const touching = new Uint8Array(size);
+  let fromX = 0;
+  let toX = 0;
+  for (let i = 0; i < size; i++) {
+    const x0 = x[i];
+    while (x[fromX] - x0 < -radius) fromX++;
+    while (toX + 1 < size && x[toX + 1] - x0 <= radius) toX++;
+    let yShift = y[i] + radius;
+    for (let j = fromX; j <= toX; j++) {
+      const deltaX = x[j] - x0;
+      const current = y[j] + Math.sqrt(radius2 - deltaX * deltaX);
+      if (current > yShift) yShift = current;
+    }
+    const limit = yShift - RELATIVE_TOLERANCE * (Math.abs(yShift) + radius);
+    for (let j = fromX; j <= toX; j++) {
+      const deltaX = x[j] - x0;
+      if (y[j] + Math.sqrt(radius2 - deltaX * deltaX) >= limit) {
+        touching[j] = 1;
+      }
+    }
+  }
+  return touching;
+}
+
+/**
+ * Marks the touched points with the squared comparison, no square root.
+ * @param radius - radius of the circle.
+ * @returns the mask of the touched points.
+ */
+function markWithSquare(radius: number): Uint8Array {
+  const radius2 = radius * radius;
+  const touching = new Uint8Array(size);
+  let fromX = 0;
+  let toX = 0;
+  for (let i = 0; i < size; i++) {
+    const x0 = x[i];
+    while (x[fromX] - x0 < -radius) fromX++;
+    while (toX + 1 < size && x[toX + 1] - x0 <= radius) toX++;
+    let yShift = y[i] + radius;
+    for (let j = fromX; j <= toX; j++) {
+      const deltaX = x[j] - x0;
+      const current = y[j] + Math.sqrt(radius2 - deltaX * deltaX);
+      if (current > yShift) yShift = current;
+    }
+    const limit = yShift - RELATIVE_TOLERANCE * (Math.abs(yShift) + radius);
+    for (let j = fromX; j <= toX; j++) {
+      const needed = limit - y[j];
+      if (needed <= 0) {
+        touching[j] = 1;
+        continue;
+      }
+      const deltaX = x[j] - x0;
+      if (radius2 - deltaX * deltaX >= needed * needed) {
+        touching[j] = 1;
+      }
+    }
+  }
+  return touching;
+}
+
+for (const radius of [1.5, 50]) {
+  const withSqrt = markWithSqrt(radius);
+  const withSquare = markWithSquare(radius);
+  let touched = 0;
+  for (let i = 0; i < size; i++) {
+    if (withSqrt[i] !== withSquare[i]) {
+      throw new Error(`masks differ at ${i}`);
+    }
+    touched += withSqrt[i];
+  }
+  console.log(`\nradius ${radius}, n = ${size}, ${touched} touched points`);
+
+  new Benchmark.Suite()
+    .add('sqrt', () => markWithSqrt(radius), { minSamples: 30 })
+    .add('squared', () => markWithSquare(radius), { minSamples: 30 })
+    .on('cycle', (event: Benchmark.Event) => {
+      const { target } = event;
+      const stats = target.stats;
+      if (!stats) return;
+      console.log(
+        `${target.name?.padEnd(8)} ${(1 / stats.mean).toFixed(1).padStart(7)} ops/sec | ${((stats.mean * 1e9) / size).toFixed(1)} ns/point | rme ${stats.rme.toFixed(1)}% | ${stats.sample.length} samples`,
+      );
+    })
+    .run();
+}

benchmark/xyRollingCircleTransform.ts

@@ -0,0 +1,117 @@
+/* eslint-disable no-console */
+import Benchmark from 'benchmark';
+import type { NumberArray } from 'cheminfo-types';
+import { XSadd } from 'ml-xsadd';
+
+import { xFindClosestIndex } from '../src/x/xFindClosestIndex.ts';
+
+const size = 100000;
+
+// Deterministic, reproducible input: a slowly varying baseline plus noise.
+const { random } = new XSadd(42);
+const x = new Float64Array(size);
+const y = new Float64Array(size);
+for (let i = 0; i < size; i++) {
+  x[i] = i;
+  y[i] = Math.sin(i / 500) * 10 + random();
+}
+
+/**
+ * The window bounds as they were located *before*: two binary searches per
+ * center, plus a guard because the closest index may fall outside the window.
+ * @param radius - radius of the circle.
+ * @returns the y center of the circle for each x value.
+ */
+function withFindClosestIndex(radius: number): Float64Array {
+  const radius2 = radius * radius;
+  const yCenters = new Float64Array(size);
+  for (let i = 0; i < size; i++) {
+    const x0 = x[i];
+    const fromX = xFindClosestIndex(x, x0 - radius);
+    const toX = xFindClosestIndex(x, x0 + radius);
+    let yShift = y[i] + radius;
+    for (let j = fromX; j <= toX; j++) {
+      const deltaX = x[j] - x0;
+      if (Math.abs(deltaX) > radius) {
+        continue;
+      }
+      const current = y[j] + Math.sqrt(radius2 - deltaX * deltaX);
+      if (current > yShift) {
+        yShift = current;
+      }
+    }
+    yCenters[i] = yShift;
+  }
+  return yCenters;
+}
+
+/**
+ * The window bounds as they are located *now*: both ends only move forward, so
+ * the bounds are exact and no per-element guard is needed.
+ * @param radius - radius of the circle.
+ * @returns the y center of the circle for each x value.
+ */
+function withSlidingWindow(radius: number): Float64Array {
+  const radius2 = radius * radius;
+  const yCenters = new Float64Array(size);
+  let fromX = 0;
+  let toX = 0;
+  for (let i = 0; i < size; i++) {
+    const x0 = x[i];
+    while (x[fromX] - x0 < -radius) {
+      fromX++;
+    }
+    while (toX + 1 < size && x[toX + 1] - x0 <= radius) {
+      toX++;
+    }
+    let yShift = y[i] + radius;
+    for (let j = fromX; j <= toX; j++) {
+      const deltaX = x[j] - x0;
+      const current = y[j] + Math.sqrt(radius2 - deltaX * deltaX);
+      if (current > yShift) {
+        yShift = current;
+      }
+    }
+    yCenters[i] = yShift;
+  }
+  return yCenters;
+}
+
+function checksum(array: NumberArray): number {
+  let sum = 0;
+  for (const value of array) {
+    sum += value;
+  }
+  return sum;
+}
+
+for (const radius of [1.5, 50]) {
+  const before = withFindClosestIndex(radius);
+  const after = withSlidingWindow(radius);
+  for (let i = 0; i < before.length; i++) {
+    if (before[i] !== after[i]) {
+      throw new Error(`values differ at ${i}: ${before[i]} vs ${after[i]}`);
+    }
+  }
+
+  const points = 2 * Math.floor(radius) + 1;
+  console.log(
+    `\nradius ${radius} (~${points} points per window), n = ${size}, checksum ${checksum(after).toFixed(6)}`,
+  );
+
+  new Benchmark.Suite()
+    .add('xFindClosestIndex', () => withFindClosestIndex(radius), {
+      minSamples: 30,
+    })
+    .add('sliding window', () => withSlidingWindow(radius), { minSamples: 30 })
+    .on('cycle', (event: Benchmark.Event) => {
+      const { target } = event;
+      const stats = target.stats;
+      if (!stats) return;
+      const nsPerPoint = (stats.mean * 1e9) / size;
+      console.log(
+        `${target.name?.padEnd(18)} ${(1 / stats.mean).toFixed(1).padStart(7)} ops/sec | ${nsPerPoint.toFixed(1)} ns/point | rme ${stats.rme.toFixed(1)}% | ${stats.sample.length} samples`,
+      );
+    })
+    .run();
+}