From f2aced986f6dc46798a5f015b2010572eca9c5ae Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Wed, 9 Sep 2026 16:39:13 +0200 Subject: [PATCH 1/9] feat: add xysFilter --- package-lock.json | 8 +- package.json | 3 +- .../__snapshots__/index.test.ts.snap | 3 + src/index.ts | 1 + src/xyArray/__tests__/xyArrayAlign.test.ts | 8 +- src/xyArray/xyArrayAlign.ts | 19 +- src/xys/__tests__/xysFilter.test.ts | 173 ++++++++++++++++++ src/xys/index.ts | 1 + src/xys/xysFilter.ts | 120 ++++++++++++ 9 files changed, 314 insertions(+), 22 deletions(-) create mode 100644 src/xys/__tests__/xysFilter.test.ts create mode 100644 src/xys/index.ts create mode 100644 src/xys/xysFilter.ts diff --git a/package-lock.json b/package-lock.json index 86fec393a..d7fd0314b 100644 --- a/package-lock.json +++ b/package-lock.json @@ -10,7 +10,7 @@ "license": "MIT", "dependencies": { "binary-search": "^1.3.6", - "cheminfo-types": "^1.15.0", + "cheminfo-types": "^1.16.0", "fft.js": "^4.0.4", "is-any-array": "^3.0.0", "ml-matrix": "^6.14.0", @@ -3848,9 +3848,9 @@ } }, "node_modules/cheminfo-types": { - "version": "1.15.0", - "resolved": "https://registry.npmjs.org/cheminfo-types/-/cheminfo-types-1.15.0.tgz", - "integrity": "sha512-shv45WN2u0yN9EHH1bisNrv+fy4Cw+eLM5lOoriP67mePrwbHZ1kJqg90C8GEU7K1A8gJsicEoVZHcuBbuul/w==", + "version": "1.16.0", + "resolved": "https://registry.npmjs.org/cheminfo-types/-/cheminfo-types-1.16.0.tgz", + "integrity": "sha512-j0MiXFnLRVEZg3vUr7SXm6F3HuDNVE5arfQyDMiy1w2fAcuOs3QRgVsw7aOIACvr0aqbvn85cucwt+yOu64IiQ==", "license": "MIT" }, "node_modules/ci-info": { diff --git a/package.json b/package.json index c292553a6..4b3398596 100644 --- a/package.json +++ b/package.json @@ -18,6 +18,7 @@ "./xy2": "./lib/xy2/index.js", "./xyArray": "./lib/xyArray/index.js", "./xyObject": "./lib/xyObject/index.js", + "./xys": "./lib/xys/index.js", "./zones": "./lib/zones/index.js" }, "files": [ @@ -40,7 +41,7 @@ }, "dependencies": { "binary-search": "^1.3.6", - "cheminfo-types": "^1.15.0", + "cheminfo-types": "^1.16.0", "fft.js": "^4.0.4", "is-any-array": "^3.0.0", "ml-matrix": "^6.14.0", diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index c44dcf1ec..fbeffd810 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -138,8 +138,10 @@ exports[`existence of exported functions 1`] = ` "xreimSortX", "xreimZeroFilling", "xyArrayAlign", + "xyArrayAlignByIntensity", "xyArrayAlignToFirst", "xyArrayMerge", + "xyArrayMergeByIntensity", "xyArrayWeightedMerge", "xyObjectBestPoints", "xyObjectCheck", @@ -154,6 +156,7 @@ exports[`existence of exported functions 1`] = ` "xyObjectSortX", "xyObjectSumY", "xyObjectToXY", + "xysFilter", "zonesNormalize", "zonesWithPoints", "zonesToXEquallySpaced", diff --git a/src/index.ts b/src/index.ts index 113946a47..d95d37ede 100644 --- a/src/index.ts +++ b/src/index.ts @@ -12,6 +12,7 @@ export * from './xreim/index.ts'; export * from './xyArray/index.ts'; export * from './xyObject/index.ts'; +export * from './xys/index.ts'; export * from './zones/index.ts'; diff --git a/src/xyArray/__tests__/xyArrayAlign.test.ts b/src/xyArray/__tests__/xyArrayAlign.test.ts index 379a91e38..241988fce 100644 --- a/src/xyArray/__tests__/xyArrayAlign.test.ts +++ b/src/xyArray/__tests__/xyArrayAlign.test.ts @@ -29,7 +29,11 @@ test('The y values must be present everywhere', () => { const result = xyArrayAlign(data, { delta: 0.15, requiredY: true }); expect(result).toStrictEqual({ - x: [3.025], - ys: [[1], [1], [2]], + x: Float64Array.from([3.025]), + ys: [ + Float64Array.from([1]), + Float64Array.from([1]), + Float64Array.from([2]), + ], }); }); diff --git a/src/xyArray/xyArrayAlign.ts b/src/xyArray/xyArrayAlign.ts index 64dbc1a74..e545c844a 100644 --- a/src/xyArray/xyArrayAlign.ts +++ b/src/xyArray/xyArrayAlign.ts @@ -1,6 +1,7 @@ import type { DataXY, NumberArray } from 'cheminfo-types'; import { xyJoinX } from '../xy/index.ts'; +import { xysFilter } from '../xys/index.ts'; import { getSlots } from './utils/getSlots.ts'; @@ -53,21 +54,9 @@ export function xyArrayAlign( } } - if (requiredY) return filterRequiredY(x, ys); + if (requiredY) { + return xysFilter({ x, ys }, { requiredY: true }); + } return { x, ys }; } - -function filterRequiredY(x: Float64Array, ys: Float64Array[]) { - const newX: number[] = []; - const newYs: number[][] = Array.from(ys, () => []); - for (let i = 0; i < x.length; i++) { - if (ys.every((y) => y[i] !== 0)) { - newX.push(x[i]); - for (let j = 0; j < ys.length; j++) { - newYs[j].push(ys[j][i]); - } - } - } - return { x: newX, ys: newYs }; -} diff --git a/src/xys/__tests__/xysFilter.test.ts b/src/xys/__tests__/xysFilter.test.ts new file mode 100644 index 000000000..8acd733d4 --- /dev/null +++ b/src/xys/__tests__/xysFilter.test.ts @@ -0,0 +1,173 @@ +import { expect, test } from 'vitest'; + +import { xyArrayAlign } from '../../xyArray/xyArrayAlign.ts'; +import { xysFilter } from '../xysFilter.ts'; + +/* + * Eight spectra, four slots. + * x=10 x=20 x=30 x=40 + * present 4 3 3 2 spectra + * longest 1 3 3 2 run of consecutive spectra + */ +function aligned() { + return { + x: [10, 20, 30, 40], + ys: [ + [1, 0, 0, 1], + [0, 0, 1, 1], + [1, 0, 1, 0], + [0, 0, 1, 0], + [1, 0, 0, 0], + [0, 1, 0, 0], + [1, 1, 0, 0], + [0, 1, 0, 0], + ], + }; +} + +test('without options nothing is removed', () => { + const data = aligned(); + const result = xysFilter(data); + + expect(Array.from(result.x)).toStrictEqual(data.x); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual(data.ys); +}); + +test('minY on its own removes nothing', () => { + const data = aligned(); + const result = xysFilter(data, { minY: 5 }); + + expect(Array.from(result.x)).toStrictEqual(data.x); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual(data.ys); +}); + +test('minNumberOfSpectra counts the spectra wherever they are', () => { + const result = xysFilter(aligned(), { minNumberOfSpectra: 3 }); + + expect(Array.from(result.x)).toStrictEqual([10, 20, 30]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [1, 0, 0], + [0, 0, 1], + [1, 0, 1], + [0, 0, 1], + [1, 0, 0], + [0, 1, 0], + [1, 1, 0], + [0, 1, 0], + ]); +}); + +test('minConsecutiveSpectra rejects the scattered slot the count kept', () => { + // x=10 is in four spectra but never in two neighbouring ones. + const result = xysFilter(aligned(), { minConsecutiveSpectra: 3 }); + + expect(Array.from(result.x)).toStrictEqual([20, 30]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [0, 0], + [0, 1], + [0, 1], + [0, 1], + [0, 0], + [1, 0], + [1, 0], + [1, 0], + ]); +}); + +test('two consecutive spectra are enough for the narrow slot', () => { + const result = xysFilter(aligned(), { minConsecutiveSpectra: 2 }); + + expect(Array.from(result.x)).toStrictEqual([20, 30, 40]); +}); + +test('the criteria are combined', () => { + const result = xysFilter(aligned(), { + minNumberOfSpectra: 3, + minConsecutiveSpectra: 2, + }); + + expect(Array.from(result.x)).toStrictEqual([20, 30]); +}); + +test('requiredY keeps only the slots every spectrum contains', () => { + const result = xysFilter( + { + x: [1, 2], + ys: [ + [1, 1], + [0, 1], + [3, 1], + ], + }, + { requiredY: true }, + ); + + expect(Array.from(result.x)).toStrictEqual([2]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([[1], [1], [1]]); +}); + +test('a negative peak counts as present until minY is set', () => { + const data = { + x: [1, 2], + ys: [ + [-5, 0], + [3, 4], + ], + }; + + expect(Array.from(xysFilter(data, { requiredY: true }).x)).toStrictEqual([1]); + // Setting minY compares the y itself, so the negative peak stops counting. + expect( + Array.from(xysFilter(data, { requiredY: true, minY: 0 }).x), + ).toStrictEqual([]); +}); + +test('minY sets how intense a peak must be to count', () => { + const data = { + x: [1, 2], + ys: [ + [5, 1], + [5, 1], + [5, 1], + ], + }; + + expect( + Array.from(xysFilter(data, { minNumberOfSpectra: 3 }).x), + ).toStrictEqual([1, 2]); + expect( + Array.from(xysFilter(data, { minNumberOfSpectra: 3, minY: 2 }).x), + ).toStrictEqual([1]); +}); + +test('agrees with the requiredY option of xyArrayAlign', () => { + const data = [ + { x: [1, 2, 3], y: [1, 1, 1] }, + { x: [0.1, 1.1, 2.1, 3.1, 4.1], y: [1, 1, 1, 1, 1] }, + { x: [2.9, 3.1, 3.9, 4.9], y: [1, 1, 1, 1] }, + ]; + const expected = xyArrayAlign(data, { delta: 0.15, requiredY: true }); + const result = xysFilter(xyArrayAlign(data, { delta: 0.15 }), { + requiredY: true, + }); + + expect(Array.from(result.x)).toStrictEqual(Array.from(expected.x)); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual( + expected.ys.map((y) => Array.from(y)), + ); +}); + +test('no slot survives', () => { + const result = xysFilter(aligned(), { minNumberOfSpectra: 8 }); + + expect(Array.from(result.x)).toStrictEqual([]); + expect(result.ys).toHaveLength(8); + expect(Array.from(result.ys[0])).toStrictEqual([]); +}); + +test('no spectra', () => { + const result = xysFilter({ x: [1, 2], ys: [] }, { requiredY: true }); + + expect(Array.from(result.x)).toStrictEqual([1, 2]); + expect(result.ys).toStrictEqual([]); +}); diff --git a/src/xys/index.ts b/src/xys/index.ts new file mode 100644 index 000000000..b52c468a7 --- /dev/null +++ b/src/xys/index.ts @@ -0,0 +1 @@ +export * from './xysFilter.ts'; diff --git a/src/xys/xysFilter.ts b/src/xys/xysFilter.ts new file mode 100644 index 000000000..0d614ef70 --- /dev/null +++ b/src/xys/xysFilter.ts @@ -0,0 +1,120 @@ +import type { DataXYs } from 'cheminfo-types'; + +export interface XYsFilterOptions { + /** + * The minimal y a spectrum must have in a slot to contain it, meant to be set + * above 0. When it is not set a spectrum contains a slot when its y there is + * not 0, so the negative peaks of a signed spectrum count as present. + * On its own this never removes a slot, it only changes what the criteria below + * count as contained. + * @default undefined + */ + minY?: number; + + /** + * The smallest number of spectra that must contain a slot. 0 does not filter. + * @default 0 + */ + minNumberOfSpectra?: number; + + /** + * The smallest number of neighbouring spectra that must all contain a slot, + * neighbouring meaning next to each other in `ys`. Made for spectra ordered by + * retention time, where a real compound is present in a run of consecutive + * scans while noise is scattered. 0 does not filter. + * @default 0 + */ + minConsecutiveSpectra?: number; + + /** + * If true, every spectrum must contain the slot. Same as setting + * `minNumberOfSpectra` to the number of spectra. + * @default false + */ + requiredY?: boolean; +} + +/** + * Removes the slots of an aligned set of spectra that too few spectra contain. + * + * `requiredY`, `minNumberOfSpectra` and `minConsecutiveSpectra` are the criteria, + * and `minY` only says what they count as contained. They may be combined, a slot + * is kept only when it passes all the ones that are set, so with no criterion set + * nothing is removed. They only ever remove slots, the y values are never changed. + * @param data - aligned spectra, as returned by the align methods. + * @param options - options. + * @returns the same spectra with the rejected slots removed. + */ +export function xysFilter( + data: DataXYs, + options: XYsFilterOptions = {}, +): DataXYs { + const { + minY, + minNumberOfSpectra = 0, + minConsecutiveSpectra = 0, + requiredY = false, + } = options; + + const { x, ys } = data; + const length = x.length; + const numberOfSpectra = ys.length; + const minNumber = requiredY ? numberOfSpectra : minNumberOfSpectra; + const hasMinY = minY !== undefined; + const threshold = minY ?? 0; + + const numbers = new Uint32Array(length); + let longestRuns: Uint32Array | undefined; + if (minConsecutiveSpectra > 0) { + longestRuns = new Uint32Array(length); + const runs = new Uint32Array(length); + for (let j = 0; j < numberOfSpectra; j++) { + const y = ys[j]; + for (let i = 0; i < length; i++) { + if (hasMinY ? y[i] > threshold : y[i] !== 0) { + numbers[i]++; + const run = runs[i] + 1; + runs[i] = run; + if (run > longestRuns[i]) longestRuns[i] = run; + } else { + runs[i] = 0; + } + } + } + } else if (minNumber > 0) { + for (let j = 0; j < numberOfSpectra; j++) { + const y = ys[j]; + for (let i = 0; i < length; i++) { + if (hasMinY ? y[i] > threshold : y[i] !== 0) numbers[i]++; + } + } + } + + let kept = 0; + for (let i = 0; i < length; i++) { + if (numbers[i] < minNumber) continue; + if (longestRuns && longestRuns[i] < minConsecutiveSpectra) continue; + kept++; + } + + const newX = new Float64Array(kept); + const positions = new Uint32Array(kept); + let at = 0; + for (let i = 0; i < length; i++) { + if (numbers[i] < minNumber) continue; + if (longestRuns && longestRuns[i] < minConsecutiveSpectra) continue; + newX[at] = x[i]; + positions[at] = i; + at++; + } + + const newYs = Array.from(ys, (y) => { + const newY = new Float64Array(kept); + for (let i = 0; i < kept; i++) { + newY[i] = y[positions[i]]; + } + return newY; + }); + + return { x: newX, ys: newYs }; +} From 69f1520da7d9fa30afcc0e6d2e171d275d109ffd Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Wed, 9 Sep 2026 16:59:28 +0200 Subject: [PATCH 2/9] feat: add xyArrayMaxY --- src/xyArray/__tests__/xyArrayMaxY.test.ts | 54 +++++++++++++++++++++++ src/xyArray/xyArrayMaxY.ts | 26 +++++++++++ 2 files changed, 80 insertions(+) create mode 100644 src/xyArray/__tests__/xyArrayMaxY.test.ts create mode 100644 src/xyArray/xyArrayMaxY.ts diff --git a/src/xyArray/__tests__/xyArrayMaxY.test.ts b/src/xyArray/__tests__/xyArrayMaxY.test.ts new file mode 100644 index 000000000..973d765c3 --- /dev/null +++ b/src/xyArray/__tests__/xyArrayMaxY.test.ts @@ -0,0 +1,54 @@ +import { expect, test } from 'vitest'; + +import { xyArrayMaxY } from '../xyArrayMaxY.ts'; + +test('the largest y of all the spectra', () => { + expect( + xyArrayMaxY([ + { x: [1, 2], y: [10, 50] }, + { x: [3], y: [30] }, + ]), + ).toBe(50); +}); + +test('empty spectra are skipped', () => { + expect( + xyArrayMaxY([ + { x: [], y: [] }, + { x: [1], y: [7] }, + { x: [], y: [] }, + ]), + ).toBe(7); +}); + +test('every y may be negative', () => { + expect( + xyArrayMaxY([ + { x: [1], y: [-5] }, + { x: [2], y: [-2] }, + ]), + ).toBe(-2); +}); + +test('typed arrays', () => { + expect( + xyArrayMaxY([ + { x: Float64Array.from([1]), y: Float64Array.from([3]) }, + { x: Float64Array.from([2]), y: Float64Array.from([9]) }, + ]), + ).toBe(9); +}); + +test('no point at all', () => { + expect(() => xyArrayMaxY([])).toThrow('can not process empty arrays'); + expect(() => xyArrayMaxY([{ x: [], y: [] }])).toThrow( + 'can not process empty arrays', + ); +}); + +test('the y must contain numbers', () => { + expect(() => + // @ts-expect-error the y are checked at runtime + xyArrayMaxY([{ x: [1], y: ['a'] }]), + ).toThrow('input must contain numbers'); +}); diff --git a/src/xyArray/xyArrayMaxY.ts b/src/xyArray/xyArrayMaxY.ts new file mode 100644 index 000000000..0ef1886c0 --- /dev/null +++ b/src/xyArray/xyArrayMaxY.ts @@ -0,0 +1,26 @@ +import type { DataXY } from 'cheminfo-types'; + +import { xMaxValue } from '../x/index.ts'; + +/** + * Returns the largest y of all the spectra. + * @param data - data. + * @returns the largest y. + */ +export function xyArrayMaxY(data: DataXY[]): number { + let max = Number.NEGATIVE_INFINITY; + let isEmpty = true; + + for (const spectrum of data) { + if (spectrum.y.length === 0) continue; + isEmpty = false; + const value = xMaxValue(spectrum.y); + if (value > max) max = value; + } + + if (isEmpty) { + throw new Error('can not process empty arrays'); + } + + return max; +} From 34da8e99fb9336d5270c1e7040d4923f5de4a065 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Wed, 9 Sep 2026 17:15:57 +0200 Subject: [PATCH 3/9] feat: add xyArrayAlignByIntensity --- .../__tests__/xyArrayFilterMinYValue.test.ts | 66 +++++++++++++++++++ src/xyArray/xyArrayFilterMinYValue.ts | 50 ++++++++++++++ 2 files changed, 116 insertions(+) create mode 100644 src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts create mode 100644 src/xyArray/xyArrayFilterMinYValue.ts diff --git a/src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts b/src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts new file mode 100644 index 000000000..920dc6d56 --- /dev/null +++ b/src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts @@ -0,0 +1,66 @@ +import { expect, test } from 'vitest'; + +import { xyArrayFilterMinYValue } from '../xyArrayFilterMinYValue.ts'; + +const data = [ + { x: [1, 2, 3], y: [100, 5, 1] }, + { x: [4, 5], y: [50, 2] }, +]; + +test('the threshold comes from the largest y of all the spectra', () => { + // 5% of 100, so the 2 of the second spectrum goes even though it is 4% of its own max + const result = xyArrayFilterMinYValue(data, 0.05); + + expect(result.map((spectrum) => Array.from(spectrum.x))).toStrictEqual([ + [1, 2], + [4], + ]); + expect(result.map((spectrum) => Array.from(spectrum.y))).toStrictEqual([ + [100, 5], + [50], + ]); +}); + +test('a spectrum may lose all of its points', () => { + const result = xyArrayFilterMinYValue(data, 0.6); + + expect(result.map((spectrum) => Array.from(spectrum.x))).toStrictEqual([ + [1], + [], + ]); +}); + +test('a spectrum may keep all of its points', () => { + const spectra = [ + { x: [1, 2], y: [100, 90] }, + { x: [3, 4], y: [50, 1] }, + ]; + const result = xyArrayFilterMinYValue(spectra, 0.05); + + expect(result[0].x).toStrictEqual([1, 2]); + expect(result[1].x).toStrictEqual([3]); +}); + +test('without a value nothing is filtered', () => { + expect(xyArrayFilterMinYValue(data)).toBe(data); +}); + +test('no data', () => { + expect(xyArrayFilterMinYValue([], 0.5)).toStrictEqual([]); + expect(xyArrayFilterMinYValue([{ x: [], y: [] }], 0.5)).toStrictEqual([ + { x: [], y: [] }, + ]); +}); + +test('the filtered spectra are plain arrays', () => { + const filtered = xyArrayFilterMinYValue( + [{ x: Float64Array.from([1, 2]), y: Float64Array.from([100, 1]) }], + 0.5, + ); + + // The overload gives number[], so the arrays can be pushed to. + filtered[0].x.push(3); + + expect(filtered[0].x).toStrictEqual([1, 3]); + expect(filtered[0].y).toStrictEqual([100]); +}); diff --git a/src/xyArray/xyArrayFilterMinYValue.ts b/src/xyArray/xyArrayFilterMinYValue.ts new file mode 100644 index 000000000..d15a7fe9e --- /dev/null +++ b/src/xyArray/xyArrayFilterMinYValue.ts @@ -0,0 +1,50 @@ +import type { DataXY } from 'cheminfo-types'; + +import { xyArrayMaxY } from './xyArrayMaxY.ts'; + +export function xyArrayFilterMinYValue( + data: DataXY[], + minRelativeYValue: number, +): Array>; +export function xyArrayFilterMinYValue( + data: DataXY[], + minRelativeYValue?: number, +): DataXY[]; + +/** + * Removes from every spectrum the points whose y is below a fraction of the + * largest y of the whole set. + * Unlike `xyFilterMinYValue`, which uses the largest y of the spectrum it filters, + * the threshold here is the same for all the spectra, so a weak spectrum may lose + * every one of its points. + * @param data - data. + * @param minRelativeYValue - the minimum relative value compare to the global Y max value. + * @returns the spectra without their smallest points, the data itself when there is nothing to filter on. + */ +export function xyArrayFilterMinYValue( + data: DataXY[], + minRelativeYValue?: number, +): DataXY[] { + if (minRelativeYValue === undefined) return data; + + let total = 0; + for (const spectrum of data) { + total += spectrum.y.length; + } + const threshold = total === 0 ? 0 : xyArrayMaxY(data) * minRelativeYValue; + + // The threshold is a number and not a callback, so the comparison is inlined + // rather than going through `xyFilter`, which costs 25% on 30 million points. + return data.map((spectrum) => { + const { x, y } = spectrum; + const newX: number[] = []; + const newY: number[] = []; + for (let i = 0; i < y.length; i++) { + if (y[i] >= threshold) { + newX.push(x[i]); + newY.push(y[i]); + } + } + return { x: newX, y: newY }; + }); +} From 2061f6401f3ddc18a627ff104701291fc2617544 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Wed, 9 Sep 2026 17:19:23 +0200 Subject: [PATCH 4/9] feat: add isLittleEndian --- src/utils/__tests__/isLittleEndian.test.ts | 14 ++++++++++++++ src/utils/index.ts | 1 + src/utils/isLittleEndian.ts | 11 +++++++++++ 3 files changed, 26 insertions(+) create mode 100644 src/utils/__tests__/isLittleEndian.test.ts create mode 100644 src/utils/isLittleEndian.ts diff --git a/src/utils/__tests__/isLittleEndian.test.ts b/src/utils/__tests__/isLittleEndian.test.ts new file mode 100644 index 000000000..34c59cf4e --- /dev/null +++ b/src/utils/__tests__/isLittleEndian.test.ts @@ -0,0 +1,14 @@ +import { expect, test } from 'vitest'; + +import { isLittleEndian } from '../isLittleEndian.ts'; + +test('matches the byte order seen through a DataView', () => { + const view = new DataView(new ArrayBuffer(2)); + new Uint16Array(view.buffer)[0] = 1; + + expect(isLittleEndian()).toBe(view.getUint16(0, true) === 1); +}); + +test('is stable across calls', () => { + expect(isLittleEndian()).toBe(isLittleEndian()); +}); diff --git a/src/utils/index.ts b/src/utils/index.ts index d7b069db1..e1c100ebb 100644 --- a/src/utils/index.ts +++ b/src/utils/index.ts @@ -6,6 +6,7 @@ export { clearFFTCache, setFFTCacheMaxSize } from './fftCache.ts'; export * from './getCombinations.ts'; export * from './getCombinationsIterator.ts'; export * from './getRescaler.ts'; +export * from './isLittleEndian.ts'; export * from './isPowerOfTwo.ts'; export * from './nextPowerOfTwo.ts'; export * from './recursiveResolve.ts'; diff --git a/src/utils/isLittleEndian.ts b/src/utils/isLittleEndian.ts new file mode 100644 index 000000000..7807f9aa3 --- /dev/null +++ b/src/utils/isLittleEndian.ts @@ -0,0 +1,11 @@ +const LITTLE_ENDIAN = new Uint8Array(new Uint16Array([1]).buffer)[0] === 1; + +/** + * Check whether the host stores multi-byte values least significant byte first. + * Needed whenever a typed array is reinterpreted through another view, because + * the position of the words of a value then depends on the platform. + * @returns true on a little-endian host. + */ +export function isLittleEndian(): boolean { + return LITTLE_ENDIAN; +} From d412ecbbf8932620f65c7ba903a31d106b882575 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Thu, 10 Sep 2026 07:25:38 +0200 Subject: [PATCH 5/9] feat: add xGetSortOrder and xGetApproximateSortOrder Radix sort over the bits of each double, so the cost is a fixed number of linear passes instead of n log n comparisons. Measured against a comparator sort of an index array (benchmark.js, minSamples 30, random doubles, ns per element, node 26): n comparator xGetSortOrder xGetApproximateSortOrder 50 69.8 78.8 (0.9x) 46.3 (1.5x) 1000 122.5 25.4 (4.8x) 13.9 (8.8x) 100000 230.9 24.0 (9.6x) 9.3 (24.8x) 1000000 290.4 23.0 (12.6x) 9.6 (30.3x) xGetSortOrder keys on all 64 bits and orders exactly, NaN, infinities, subnormals and the two zeros included. xGetApproximateSortOrder keys on the high word only, halving the passes, and orders to 20 bits of relative precision. The digit is 8 bits below 65536 values and 16 above, which keeps the O(radix) counting pass from dominating short arrays. Co-Authored-By: Claude Opus 5 (1M context) --- .../__snapshots__/index.test.ts.snap | 5 +- .../xGetApproximateSortOrder.test.ts | 119 +++++++++++ src/x/__tests__/xGetSortOrder.test.ts | 139 +++++++++++++ src/x/index.ts | 2 + src/x/utils/radixSortOrder.ts | 184 ++++++++++++++++++ src/x/xGetApproximateSortOrder.ts | 35 ++++ src/x/xGetSortOrder.ts | 30 +++ 7 files changed, 512 insertions(+), 2 deletions(-) create mode 100644 src/x/__tests__/xGetApproximateSortOrder.test.ts create mode 100644 src/x/__tests__/xGetSortOrder.test.ts create mode 100644 src/x/utils/radixSortOrder.ts create mode 100644 src/x/xGetApproximateSortOrder.ts create mode 100644 src/x/xGetSortOrder.ts diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index fbeffd810..72ae0ecf3 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -33,7 +33,9 @@ exports[`existence of exported functions 1`] = ` "xEnsureFloat64", "xEqualIntegrationVectorSimilarity", "xFindClosestIndex", + "xGetApproximateSortOrder", "xGetFromToIndex", + "xGetSortOrder", "xGetTargetIndex", "xHilbertTransform", "xHistogram", @@ -138,10 +140,8 @@ exports[`existence of exported functions 1`] = ` "xreimSortX", "xreimZeroFilling", "xyArrayAlign", - "xyArrayAlignByIntensity", "xyArrayAlignToFirst", "xyArrayMerge", - "xyArrayMergeByIntensity", "xyArrayWeightedMerge", "xyObjectBestPoints", "xyObjectCheck", @@ -198,6 +198,7 @@ exports[`existence of exported functions 1`] = ` "getCombinations", "getCombinationsIterator", "getRescaler", + "isLittleEndian", "isPowerOfTwo", "nextPowerOfTwo", "recursiveResolve", diff --git a/src/x/__tests__/xGetApproximateSortOrder.test.ts b/src/x/__tests__/xGetApproximateSortOrder.test.ts new file mode 100644 index 000000000..1075c0d00 --- /dev/null +++ b/src/x/__tests__/xGetApproximateSortOrder.test.ts @@ -0,0 +1,119 @@ +import { XSadd } from 'ml-xsadd'; +import { expect, test } from 'vitest'; + +import { xGetApproximateSortOrder } from '../xGetApproximateSortOrder.ts'; + +/** + * Applies an order to the values. + * @param values - the values. + * @param order - the indices, in order. + * @returns the reordered values. + */ +function apply(values: number[], order: Uint32Array): number[] { + return Array.from(order, (index) => values[index]); +} + +test('ascending order of well separated values', () => { + const values = [3, 1, 4, 1, 5, 9, 2, 6]; + + expect(apply(values, xGetApproximateSortOrder(values))).toStrictEqual([ + 1, 1, 2, 3, 4, 5, 6, 9, + ]); +}); + +test('descending order of well separated values', () => { + const values = [3, 1, 4, 1, 5]; + + expect( + apply(values, xGetApproximateSortOrder(values, { descending: true })), + ).toStrictEqual([5, 4, 3, 1, 1]); +}); + +test('negative values order below positive ones', () => { + const values = [1, -1, 0, -3.5, 2.5, -0.5]; + + expect(apply(values, xGetApproximateSortOrder(values))).toStrictEqual([ + -3.5, -1, -0.5, 0, 1, 2.5, + ]); +}); + +test('equal values keep their input order', () => { + const values = [2, 1, 2, 1, 2, 1]; + + expect(Array.from(xGetApproximateSortOrder(values))).toStrictEqual([ + 1, 3, 5, 0, 2, 4, + ]); +}); + +test('empty and single element arrays', () => { + expect(xGetApproximateSortOrder([])).toStrictEqual(new Uint32Array(0)); + expect(Array.from(xGetApproximateSortOrder([42]))).toStrictEqual([0]); +}); + +test('values closer than the key resolution keep their input order', () => { + // The high word ignores the low 32 bits of the mantissa, so these three share a key. + const values = [1, 1 + 2 ** -40, 1 + 2 ** -30]; + + expect( + Array.from(xGetApproximateSortOrder(values, { descending: true })), + ).toStrictEqual([0, 1, 2]); +}); + +test('orders to 20 bits of relative precision on the 16 bit digit path', () => { + const length = 70000; + const values = new Float64Array(length); + for (let i = 0; i < length; i++) { + values[i] = Math.abs(Math.sin(i)) * 1e6 + 1; + } + const ordered = Array.from( + xGetApproximateSortOrder(values), + (i) => values[i], + ); + + expect(ordered).toHaveLength(length); + + let worst = 0; + for (let i = 1; i < length; i++) { + if (ordered[i] < ordered[i - 1]) { + worst = Math.max(worst, (ordered[i - 1] - ordered[i]) / ordered[i - 1]); + } + } + + expect(worst).toBeLessThan(2 ** -20); +}); + +test('stays within the key resolution of a native descending sort', () => { + const length = 5000; + const values = new Float64Array(length); + for (let i = 0; i < length; i++) { + values[i] = Math.abs(Math.cos(i)) * 1e4 + 1; + } + const ordered = Float64Array.from( + xGetApproximateSortOrder(values, { descending: true }), + (i) => values[i], + ); + const native = values.toSorted().toReversed(); + let worst = 0; + for (let i = 0; i < length; i++) { + worst = Math.max(worst, Math.abs(ordered[i] - native[i]) / native[i]); + } + + expect(worst).toBeLessThan(2 ** -20); +}); + +test('matches a comparator sort on seeded values spanning many magnitudes', () => { + const random = new XSadd(7); + const values = new Float64Array(10000); + for (let i = 0; i < values.length; i++) { + values[i] = (random.random() - 0.5) * 10 ** ((random.random() - 0.5) * 40); + } + + const expected = Array.from( + { length: values.length }, + (_, index) => index, + ).toSorted((first, second) => values[second] - values[first]); + + expect( + Array.from(xGetApproximateSortOrder(values, { descending: true })), + ).toStrictEqual(expected); +}); diff --git a/src/x/__tests__/xGetSortOrder.test.ts b/src/x/__tests__/xGetSortOrder.test.ts new file mode 100644 index 000000000..2bcddf112 --- /dev/null +++ b/src/x/__tests__/xGetSortOrder.test.ts @@ -0,0 +1,139 @@ +import { expect, test } from 'vitest'; + +import { xGetSortOrder } from '../xGetSortOrder.ts'; + +/** + * Applies an order to the values. + * @param values - the values. + * @param order - the indices, in order. + * @returns the reordered values. + */ +function apply(values: number[], order: Uint32Array): number[] { + return Array.from(order, (index) => values[index]); +} + +test('ascending order of a plain array', () => { + const values = [3, 1, 4, 1, 5, 9, 2, 6]; + + expect(Array.from(xGetSortOrder(values))).toStrictEqual([ + 1, 3, 6, 0, 2, 4, 7, 5, + ]); + expect(apply(values, xGetSortOrder(values))).toStrictEqual([ + 1, 1, 2, 3, 4, 5, 6, 9, + ]); +}); + +test('descending order', () => { + const values = [3, 1, 4, 1, 5]; + + expect( + apply(values, xGetSortOrder(values, { descending: true })), + ).toStrictEqual([5, 4, 3, 1, 1]); +}); + +test('negative values order below positive ones', () => { + const values = [1, -1, 0, -3.5, 2.5, -0.5]; + + expect(apply(values, xGetSortOrder(values))).toStrictEqual([ + -3.5, -1, -0.5, 0, 1, 2.5, + ]); +}); + +test('equal values keep their input order', () => { + const values = [2, 1, 2, 1, 2, 1]; + + expect(Array.from(xGetSortOrder(values))).toStrictEqual([1, 3, 5, 0, 2, 4]); +}); + +test('equal values keep their input order when descending', () => { + const values = [2, 1, 2, 1, 2, 1]; + + expect(Array.from(xGetSortOrder(values, { descending: true }))).toStrictEqual( + [0, 2, 4, 1, 3, 5], + ); +}); + +test('empty and single element arrays', () => { + expect(xGetSortOrder([])).toStrictEqual(new Uint32Array(0)); + expect(Array.from(xGetSortOrder([42]))).toStrictEqual([0]); +}); + +test('accepts a Float64Array without modifying it', () => { + const values = Float64Array.from([3, 1, 2]); + + expect(Array.from(xGetSortOrder(values))).toStrictEqual([1, 2, 0]); + expect(values).toStrictEqual(Float64Array.from([3, 1, 2])); +}); + +test('orders infinities, both zeros and subnormals like the native sort', () => { + const values = [ + 3, + -Infinity, + 0, + -0, + Infinity, + -3, + Number.MIN_VALUE, + -Number.MIN_VALUE, + 5e-324, + ]; + const native = Float64Array.from(values).toSorted(); + + expect(Float64Array.from(apply(values, xGetSortOrder(values)))).toStrictEqual( + native, + ); +}); + +test('NaN sorts last, like the native sort', () => { + const values = [3, Number.NaN, 1, 2]; + + expect(apply(values, xGetSortOrder(values))).toStrictEqual([ + 1, + 2, + 3, + Number.NaN, + ]); +}); + +test('matches the native sort on the 8 bit digit path', () => { + const values = new Float64Array(5000); + for (let i = 0; i < values.length; i++) { + values[i] = Math.sin(i) * 1e6; + } + const sorted = values.toSorted(); + const order = xGetSortOrder(values); + + expect(order).toHaveLength(5000); + expect(Float64Array.from(order, (index) => values[index])).toStrictEqual( + sorted, + ); +}); + +test('matches the native sort on the 16 bit digit path', () => { + const length = 70000; + const values = new Float64Array(length); + for (let i = 0; i < length; i++) { + values[i] = Math.sin(i) * 1e6 - 5e5; + } + const sorted = values.toSorted(); + const order = xGetSortOrder(values); + + expect(order).toHaveLength(length); + expect(Float64Array.from(order, (index) => values[index])).toStrictEqual( + sorted, + ); +}); + +test('descending is the exact reverse of ascending when values are distinct', () => { + const values = new Float64Array(3000); + for (let i = 0; i < values.length; i++) { + values[i] = Math.cos(i) * 1234.5; + } + const ascending = Float64Array.from(xGetSortOrder(values), (i) => values[i]); + const descending = Float64Array.from( + xGetSortOrder(values, { descending: true }), + (i) => values[i], + ); + + expect(descending).toStrictEqual(ascending.toReversed()); +}); diff --git a/src/x/index.ts b/src/x/index.ts index d185ae2f8..2a0d24c84 100644 --- a/src/x/index.ts +++ b/src/x/index.ts @@ -23,7 +23,9 @@ export * from './xDotProduct.ts'; export * from './xEnsureFloat64.ts'; export * from './xEqualIntegrationVectorSimilarity.ts'; export * from './xFindClosestIndex.ts'; +export * from './xGetApproximateSortOrder.ts'; export * from './xGetFromToIndex.ts'; +export * from './xGetSortOrder.ts'; export * from './xGetTargetIndex.ts'; export * from './xHilbertTransform.ts'; export * from './xHistogram.ts'; diff --git a/src/x/utils/radixSortOrder.ts b/src/x/utils/radixSortOrder.ts new file mode 100644 index 000000000..65f05427c --- /dev/null +++ b/src/x/utils/radixSortOrder.ts @@ -0,0 +1,184 @@ +import type { NumberArray } from 'cheminfo-types'; + +import { isLittleEndian } from '../../utils/isLittleEndian.ts'; + +// Internal helpers: not re-exported from `x/index.ts`, so they stay out of the public API. + +// The counting pass costs O(radix) whatever the length, so the digit stays narrow until the +// array is long enough to amortize a wide one. +const SMALL_DIGIT_BITS = 8; +const LARGE_DIGIT_BITS = 16; +const LARGE_ARRAY_LENGTH = 1 << LARGE_DIGIT_BITS; + +/** + * Orders the values by their full 64 bits, as an LSD radix sort. + * @param array - the values to order. + * @param descending - order from the largest to the smallest. + * @returns the indices of the values, in order. + */ +export function radixSortOrderExact( + array: NumberArray, + descending: boolean, +): Uint32Array { + const values = ensureFloat64(array); + const length = values.length; + let order = new Uint32Array(length); + if (length < 2) return order; + + const digitBits = getDigitBits(length); + const radix = 1 << digitBits; + const digitMask = radix - 1; + const words = new Uint32Array(values.buffer, values.byteOffset, length * 2); + const flip = descending ? 0xffffffff : 0; + const highWord = isLittleEndian() ? 1 : 0; + const lowWord = 1 - highWord; + + let keysLow = new Uint32Array(length); + let keysHigh = new Uint32Array(length); + for (let i = 0; i < length; i++) { + const high = words[i * 2 + highWord]; + const low = words[i * 2 + lowWord]; + // Monotonic unsigned key of the double: a negative flips every bit, a positive only the + // sign bit. Complementing that key reverses the order. + if ((high & 0x80000000) === 0) { + keysHigh[i] = high ^ 0x80000000 ^ flip; + keysLow[i] = low ^ flip; + } else { + keysHigh[i] = ~high ^ flip; + keysLow[i] = ~low ^ flip; + } + order[i] = i; + } + + let nextOrder = new Uint32Array(length); + let nextLow = new Uint32Array(length); + let nextHigh = new Uint32Array(length); + const counts = new Uint32Array(radix); + + for (let shift = 0; shift < 64; shift += digitBits) { + const onHigh = shift >= 32; + const wordShift = onHigh ? shift - 32 : shift; + const keys = onHigh ? keysHigh : keysLow; + + counts.fill(0); + for (let i = 0; i < length; i++) { + counts[(keys[i] >>> wordShift) & digitMask]++; + } + // A digit shared by every value would make this pass the identity. + if (counts[(keys[0] >>> wordShift) & digitMask] === length) continue; + cumulate(counts); + + for (let i = 0; i < length; i++) { + const low = keysLow[i]; + const high = keysHigh[i]; + const at = counts[((onHigh ? high : low) >>> wordShift) & digitMask]++; + nextOrder[at] = order[i]; + nextLow[at] = low; + nextHigh[at] = high; + } + // A destructuring swap makes JavaScriptCore give up the typed array specialization of + // these variables, which costs 8x on the whole sort. + let swap = order; + order = nextOrder; + nextOrder = swap; + swap = keysLow; + keysLow = nextLow; + nextLow = swap; + swap = keysHigh; + keysHigh = nextHigh; + nextHigh = swap; + } + + return order; +} + +/** + * Orders the values by the high word of each double, as an LSD radix sort. Values differing by + * less than their 2^-20 relative resolution share a key and keep their input order. + * @param array - the values to order. + * @param descending - order from the largest to the smallest. + * @returns the indices of the values, in order. + */ +export function radixSortOrderApproximate( + array: NumberArray, + descending: boolean, +): Uint32Array { + const values = ensureFloat64(array); + const length = values.length; + let order = new Uint32Array(length); + if (length < 2) return order; + + const digitBits = getDigitBits(length); + const radix = 1 << digitBits; + const digitMask = radix - 1; + const words = new Uint32Array(values.buffer, values.byteOffset, length * 2); + const flip = descending ? 0xffffffff : 0; + const highWord = isLittleEndian() ? 1 : 0; + + let keys = new Uint32Array(length); + for (let i = 0; i < length; i++) { + const high = words[i * 2 + highWord]; + keys[i] = ((high & 0x80000000) === 0 ? high ^ 0x80000000 : ~high) ^ flip; + order[i] = i; + } + + let nextOrder = new Uint32Array(length); + let nextKeys = new Uint32Array(length); + const counts = new Uint32Array(radix); + + for (let shift = 0; shift < 32; shift += digitBits) { + counts.fill(0); + for (let i = 0; i < length; i++) { + counts[(keys[i] >>> shift) & digitMask]++; + } + if (counts[(keys[0] >>> shift) & digitMask] === length) continue; + cumulate(counts); + + for (let i = 0; i < length; i++) { + const key = keys[i]; + const at = counts[(key >>> shift) & digitMask]++; + nextOrder[at] = order[i]; + nextKeys[at] = key; + } + // See radixSortOrderExact: a destructuring swap is 8x slower on JavaScriptCore. + let swap = order; + order = nextOrder; + nextOrder = swap; + swap = keys; + keys = nextKeys; + nextKeys = swap; + } + + return order; +} + +/** + * Turns the counts of each digit into the position its first entry goes to. + * @param counts - counts per digit, replaced in place by the start offsets. + */ +function cumulate(counts: Uint32Array): void { + let start = 0; + for (let i = 0; i < counts.length; i++) { + const count = counts[i]; + counts[i] = start; + start += count; + } +} + +/** + * Picks the digit width, keeping the radix no wider than the array it scans. + * @param length - number of values to order. + * @returns the number of bits of one digit. + */ +function getDigitBits(length: number): number { + return length >= LARGE_ARRAY_LENGTH ? LARGE_DIGIT_BITS : SMALL_DIGIT_BITS; +} + +/** + * Returns the values as a Float64Array, without copying one that already is. + * @param array - values. + * @returns the values, readable as doubles. + */ +function ensureFloat64(array: NumberArray): Float64Array { + return array instanceof Float64Array ? array : Float64Array.from(array); +} diff --git a/src/x/xGetApproximateSortOrder.ts b/src/x/xGetApproximateSortOrder.ts new file mode 100644 index 000000000..49aa665fa --- /dev/null +++ b/src/x/xGetApproximateSortOrder.ts @@ -0,0 +1,35 @@ +import type { NumberArray } from 'cheminfo-types'; + +import { radixSortOrderApproximate } from './utils/radixSortOrder.ts'; + +export interface XGetApproximateSortOrderOptions { + /** + * Order the values from the largest to the smallest. + * @default false + */ + descending?: boolean; +} + +/** + * Returns the indices that approximately order the values, without moving them. + * Only the high word of each double is used as the sort key, which halves the work of + * `xGetSortOrder` but orders to 20 bits of relative precision instead of exactly: two values + * closer than 2^-20 relative share a key and keep their input order. Reading the result back as + * values therefore does not give a monotonic sequence — on random data about one adjacent pair + * in fourteen is transposed, by at most 2^-20 relative. + * Never feed this order to anything that assumes sorted values, a binary search in particular; + * use `xGetSortOrder` there. It is meant for ranking by magnitude, such as taking the n most + * intense points of a spectrum, where that resolution is far below the noise. + * The sort is stable, so equal values keep their input order. + * @param array - values to order. + * @param options - options. + * @returns the indices of the values, in approximate order. + */ +export function xGetApproximateSortOrder( + array: NumberArray, + options: XGetApproximateSortOrderOptions = {}, +): Uint32Array { + const { descending = false } = options; + + return radixSortOrderApproximate(array, descending); +} diff --git a/src/x/xGetSortOrder.ts b/src/x/xGetSortOrder.ts new file mode 100644 index 000000000..a2a2dcc59 --- /dev/null +++ b/src/x/xGetSortOrder.ts @@ -0,0 +1,30 @@ +import type { NumberArray } from 'cheminfo-types'; + +import { radixSortOrderExact } from './utils/radixSortOrder.ts'; + +export interface XGetSortOrderOptions { + /** + * Order the values from the largest to the smallest. + * @default false + */ + descending?: boolean; +} + +/** + * Returns the indices that order the values, without moving them. + * The sort is a radix sort over the full 64 bits of each double: it costs a fixed number of + * linear passes instead of the n log n comparisons of `Array.prototype.sort`, and orders + * exactly, NaN, infinities, subnormals and the two zeros included. + * The sort is stable, so equal values keep their input order. + * @param array - values to order. + * @param options - options. + * @returns the indices of the values, in order. + */ +export function xGetSortOrder( + array: NumberArray, + options: XGetSortOrderOptions = {}, +): Uint32Array { + const { descending = false } = options; + + return radixSortOrderExact(array, descending); +} From 2fd32f105af9aa208e2f2b262ca5bf5e9131530b Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Thu, 10 Sep 2026 08:02:02 +0200 Subject: [PATCH 6/9] feat: add xDoubleTypedArrayLength A typed array cannot be resized, so a buffer of unknown size has to allocate a longer one and copy into it. Doubling keeps the total cost of filling it linear, and an empty array becomes one element long so that doubling makes progress. Co-Authored-By: Claude Opus 5 (1M context) --- .../__snapshots__/index.test.ts.snap | 1 + .../__tests__/xDoubleTypedArrayLength.test.ts | 33 +++++++++++++++++++ src/x/index.ts | 1 + src/x/xDoubleTypedArrayLength.ts | 22 +++++++++++++ 4 files changed, 57 insertions(+) create mode 100644 src/x/__tests__/xDoubleTypedArrayLength.test.ts create mode 100644 src/x/xDoubleTypedArrayLength.ts diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index 72ae0ecf3..357628eb1 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -30,6 +30,7 @@ exports[`existence of exported functions 1`] = ` "xDistributionStats", "xDivide", "xDotProduct", + "xDoubleTypedArrayLength", "xEnsureFloat64", "xEqualIntegrationVectorSimilarity", "xFindClosestIndex", diff --git a/src/x/__tests__/xDoubleTypedArrayLength.test.ts b/src/x/__tests__/xDoubleTypedArrayLength.test.ts new file mode 100644 index 000000000..6e4981505 --- /dev/null +++ b/src/x/__tests__/xDoubleTypedArrayLength.test.ts @@ -0,0 +1,33 @@ +import { expect, test } from 'vitest'; + +import { xDoubleTypedArrayLength } from '../xDoubleTypedArrayLength.ts'; + +test('keeps the values and zeroes the added positions', () => { + expect( + xDoubleTypedArrayLength(Float64Array.from([1.5, -2, 3])), + ).toStrictEqual(Float64Array.from([1.5, -2, 3, 0, 0, 0])); +}); + +test('returns an array of the same kind', () => { + expect(xDoubleTypedArrayLength(Uint32Array.from([7, 8]))).toStrictEqual( + Uint32Array.from([7, 8, 0, 0]), + ); + expect(xDoubleTypedArrayLength(Int16Array.from([-4]))).toStrictEqual( + Int16Array.from([-4, 0]), + ); + expect(xDoubleTypedArrayLength(Uint8Array.from([255]))).toStrictEqual( + Uint8Array.from([255, 0]), + ); +}); + +test('an empty array becomes one element long, so doubling progresses', () => { + expect(xDoubleTypedArrayLength(new Float64Array(0))).toStrictEqual( + new Float64Array(1), + ); +}); + +test('does not return the array it was given', () => { + const array = Float64Array.from([1, 2]); + + expect(xDoubleTypedArrayLength(array)).not.toBe(array); +}); diff --git a/src/x/index.ts b/src/x/index.ts index 2a0d24c84..0e3cdcd34 100644 --- a/src/x/index.ts +++ b/src/x/index.ts @@ -20,6 +20,7 @@ export * from './xCumulative.ts'; export * from './xDistributionStats.ts'; export * from './xDivide.ts'; export * from './xDotProduct.ts'; +export * from './xDoubleTypedArrayLength.ts'; export * from './xEnsureFloat64.ts'; export * from './xEqualIntegrationVectorSimilarity.ts'; export * from './xFindClosestIndex.ts'; diff --git a/src/x/xDoubleTypedArrayLength.ts b/src/x/xDoubleTypedArrayLength.ts new file mode 100644 index 000000000..3576ff9c6 --- /dev/null +++ b/src/x/xDoubleTypedArrayLength.ts @@ -0,0 +1,22 @@ +import type { NumberArray } from 'cheminfo-types'; + +type TypedNumberArray = Exclude; + +/** + * Copies the values into an array of the same kind of twice the length, the added positions + * left at zero. A typed array cannot be resized, so a buffer that fills up has to allocate a + * longer one and copy into it; doubling is what keeps the total cost of filling a buffer of + * unknown size linear. An empty array becomes one element long, so that doubling it repeatedly + * makes progress. + * @param array - values to copy. + * @returns an array of the same kind, twice as long, holding the same values. + */ +export function xDoubleTypedArrayLength( + array: ArrayType, +): ArrayType { + const Constructor = array.constructor as new (length: number) => ArrayType; + const grown = new Constructor(Math.max(1, array.length * 2)); + grown.set(array); + + return grown; +} From 4a799f2f55ae0b777b7a08882e1f821a859f2656 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Thu, 10 Sep 2026 08:43:39 +0200 Subject: [PATCH 7/9] feat: add xyArrayAlignByIntensity and xyArrayMergeByIntensity Align and merge centroided spectra around the most intense peaks: the strongest peak of the whole set opens a slot and takes every peak within delta of it, the strongest of the rest opens the next, and so on. A slot therefore never grows wider than 2 * delta, where a chain of close neighbours can. Merge keeps only the common peaks, so it scales where the dense alignment matrix does not. mergeSortedXY picks its strategy from the number of spectra: a binary heap costs O(log k) per point, ordering the concatenated points costs the same whatever k is, and the two meet at thirty spectra (benchmark/mergeSortedXY.ts, 2M points, ns per point): spectra heap sorting 5 13.6 23.2 15 21.1 23.8 30 25.9 23.9 100 33.1 21.7 1000 55.3 22.1 5000 67.8 23.4 30000 spectra of 1000 peaks merge into 7141 slots in 2.1s. Co-Authored-By: Claude Opus 5 (1M context) --- benchmark/mergeSortedXY.ts | 195 ++++++++++++++++++ .../__snapshots__/index.test.ts.snap | 4 + .../__tests__/xyArrayAlignByIntensity.test.ts | 146 +++++++++++++ .../__tests__/xyArrayMergeByIntensity.test.ts | 116 +++++++++++ src/xyArray/index.ts | 4 + .../utils/__tests__/mergeSortedXY.test.ts | 126 +++++++++++ src/xyArray/utils/getIntensitySlots.ts | 183 ++++++++++++++++ src/xyArray/utils/mergeSortedXY.ts | 177 ++++++++++++++++ src/xyArray/utils/sortIntensitySlots.ts | 54 +++++ src/xyArray/xyArrayAlignByIntensity.ts | 56 +++++ src/xyArray/xyArrayMergeByIntensity.ts | 49 +++++ 11 files changed, 1110 insertions(+) create mode 100644 benchmark/mergeSortedXY.ts create mode 100644 src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts create mode 100644 src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts create mode 100644 src/xyArray/utils/__tests__/mergeSortedXY.test.ts create mode 100644 src/xyArray/utils/getIntensitySlots.ts create mode 100644 src/xyArray/utils/mergeSortedXY.ts create mode 100644 src/xyArray/utils/sortIntensitySlots.ts create mode 100644 src/xyArray/xyArrayAlignByIntensity.ts create mode 100644 src/xyArray/xyArrayMergeByIntensity.ts diff --git a/benchmark/mergeSortedXY.ts b/benchmark/mergeSortedXY.ts new file mode 100644 index 000000000..c69a532c4 --- /dev/null +++ b/benchmark/mergeSortedXY.ts @@ -0,0 +1,195 @@ +/* eslint-disable no-console */ +import Benchmark from 'benchmark'; +import type { DataXY } from 'cheminfo-types'; +import { XSadd } from 'ml-xsadd'; + +import { xGetSortOrder } from '../src/x/xGetSortOrder.ts'; + +// The two strategies of mergeSortedXY, copied here so one process compares them on +// the same data. The heap costs O(log k) per point, ordering the concatenated +// points costs the same whatever the number of spectra. + +const TOTAL_POINTS = 2_000_000; + +/** + * Counts the points of every spectrum. + * @param data - spectra. + * @returns the number of points of all of them. + */ +function getTotalLength(data: DataXY[]): number { + let total = 0; + for (const spectrum of data) { + total += spectrum.x.length; + } + return total; +} + +/** + * Merges by repeatedly taking the smallest x of a binary heap of the spectra. + * @param data - spectra, each one with ascending x values. + * @returns every point of every spectrum, ascending by x. + */ +function mergeByHeap(data: DataXY[]) { + const total = getTotalLength(data); + const x = new Float64Array(total); + const y = new Float64Array(total); + const spectrumIndex = new Uint32Array(total); + + const ids = new Uint32Array(data.length); + const keys = new Float64Array(data.length); + const positions = new Uint32Array(data.length); + let size = 0; + for (let i = 0; i < data.length; i++) { + const spectrum = data[i]; + if (spectrum.x.length === 0) continue; + const key = spectrum.x[0]; + let child = size++; + while (child > 0) { + const parent = (child - 1) >> 1; + if (keys[parent] <= key) break; + keys[child] = keys[parent]; + ids[child] = ids[parent]; + child = parent; + } + keys[child] = key; + ids[child] = i; + } + + let at = 0; + while (size > 0) { + const id = ids[0]; + const spectrum = data[id]; + const position = positions[id]; + const currentX = keys[0]; + x[at] = currentX; + y[at] = spectrum.y[position]; + spectrumIndex[at] = id; + at++; + + const next = position + 1; + positions[id] = next; + let key: number; + let newId: number; + if (next < spectrum.x.length) { + key = spectrum.x[next]; + newId = id; + } else { + size--; + if (size === 0) break; + key = keys[size]; + newId = ids[size]; + } + + let parent = 0; + for (;;) { + const left = 2 * parent + 1; + if (left >= size) break; + const right = left + 1; + const child = right < size && keys[right] < keys[left] ? right : left; + if (keys[child] >= key) break; + keys[parent] = keys[child]; + ids[parent] = ids[child]; + parent = child; + } + keys[parent] = key; + ids[parent] = newId; + } + + return { x, y, spectrumIndex }; +} + +/** + * Merges by concatenating the spectra and ordering the result by x. + * @param data - spectra, each one with ascending x values. + * @returns every point of every spectrum, ascending by x. + */ +function mergeBySorting(data: DataXY[]) { + const total = getTotalLength(data); + const concatenatedX = new Float64Array(total); + const concatenatedY = new Float64Array(total); + const concatenatedIndex = new Uint32Array(total); + let at = 0; + for (let i = 0; i < data.length; i++) { + const spectrum = data[i]; + concatenatedX.set(spectrum.x, at); + concatenatedY.set(spectrum.y, at); + concatenatedIndex.fill(i, at, at + spectrum.x.length); + at += spectrum.x.length; + } + + const order = xGetSortOrder(concatenatedX); + const x = new Float64Array(total); + const y = new Float64Array(total); + const spectrumIndex = new Uint32Array(total); + for (let i = 0; i < total; i++) { + const from = order[i]; + x[i] = concatenatedX[from]; + y[i] = concatenatedY[from]; + spectrumIndex[i] = concatenatedIndex[from]; + } + + return { x, y, spectrumIndex }; +} + +/** + * Builds centroided spectra sharing a peak list, each one drifting slightly. + * @param numberOfSpectra - how many spectra. + * @param peaksPerSpectrum - how many peaks in each. + * @returns the spectra, each with ascending x values. + */ +function makeSpectra( + numberOfSpectra: number, + peaksPerSpectrum: number, +): DataXY[] { + const { random } = new XSadd(42); + const spacing = 1900 / peaksPerSpectrum; + const data: DataXY[] = []; + for (let spectrum = 0; spectrum < numberOfSpectra; spectrum++) { + const x = new Float64Array(peaksPerSpectrum); + const y = new Float64Array(peaksPerSpectrum); + for (let i = 0; i < peaksPerSpectrum; i++) { + x[i] = 100 + i * spacing + (random() - 0.5) * spacing * 0.4; + y[i] = 10 ** (random() * 5); + } + data.push({ x, y }); + } + return data; +} + +console.log(`${TOTAL_POINTS} points, split over a growing number of spectra`); +for (const numberOfSpectra of [5, 15, 30, 100, 1000, 5000]) { + const peaksPerSpectrum = Math.round(TOTAL_POINTS / numberOfSpectra); + const data = makeSpectra(numberOfSpectra, peaksPerSpectrum); + const total = numberOfSpectra * peaksPerSpectrum; + + const fromHeap = mergeByHeap(data); + const fromSorting = mergeBySorting(data); + let identical = true; + for (let i = 0; i < total; i++) { + if (fromHeap.x[i] !== fromSorting.x[i]) { + identical = false; + break; + } + } + + const nanoseconds: Record = {}; + new Benchmark.Suite() + .add('heap', () => mergeByHeap(data), { minSamples: 30 }) + .add('sorting', () => mergeBySorting(data), { minSamples: 30 }) + .on('cycle', (event: Benchmark.Event) => { + const { target } = event; + const { stats, name } = target; + if (!stats || !name) return; + nanoseconds[name] = (stats.mean * 1e9) / total; + }) + .run(); + + const { heap, sorting } = nanoseconds; + console.log( + `${String(numberOfSpectra).padStart(5)} spectra x ${String(peaksPerSpectrum).padStart(6)} peaks` + + ` same x: ${identical}` + + ` heap ${heap.toFixed(1).padStart(5)} ns/point` + + ` sorting ${sorting.toFixed(1).padStart(5)} ns/point` + + ` ${heap > sorting ? `sorting ${(heap / sorting).toFixed(2)}x` : `heap ${(sorting / heap).toFixed(2)}x`}`, + ); +} diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index 357628eb1..8bdabba27 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -141,8 +141,12 @@ exports[`existence of exported functions 1`] = ` "xreimSortX", "xreimZeroFilling", "xyArrayAlign", + "xyArrayAlignByIntensity", "xyArrayAlignToFirst", + "xyArrayFilterMinYValue", + "xyArrayMaxY", "xyArrayMerge", + "xyArrayMergeByIntensity", "xyArrayWeightedMerge", "xyObjectBestPoints", "xyObjectCheck", diff --git a/src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts b/src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts new file mode 100644 index 000000000..ea4ce7a29 --- /dev/null +++ b/src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts @@ -0,0 +1,146 @@ +import { expect, test } from 'vitest'; + +import { xyArrayAlign } from '../xyArrayAlign.ts'; +import { xyArrayAlignByIntensity } from '../xyArrayAlignByIntensity.ts'; + +test('the most intense peak takes the peaks of the other spectra', () => { + const data = [ + { x: [10, 10.4, 10.8], y: [1, 10, 1] }, + { x: [10.2, 10.6], y: [2, 2] }, + ]; + const result = xyArrayAlignByIntensity(data, { delta: 0.25 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([10, 10.4, 10.8]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [1, 10, 1], + [0, 4, 0], + ]); +}); + +test('same length spectra, simple integers', () => { + const data = [ + { x: [1, 2, 3], y: [1, 1, 1] }, + { x: [0.1, 1.1, 2.1, 3.1, 4.1], y: [1, 1, 1, 1, 1] }, + { x: [2.9, 3.1, 3.9, 4.9], y: [1, 1, 1, 1] }, + ]; + const result = xyArrayAlignByIntensity(data, { delta: 0.15 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([ + 0.1, 1.05, 2.05, 2.95, 3.1, 3.9, 4.1, 4.9, + ]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [0, 1, 1, 1, 0, 0, 0, 0], + [1, 1, 1, 0, 1, 0, 1, 0], + [0, 0, 0, 1, 1, 1, 0, 1], + ]); +}); + +test('a slot never grows wider than 2 * delta', () => { + // Each x is within delta of the previous one, so a chain of close neighbours + // ends up in a single slot spanning the whole range. + const x = Array.from({ length: 21 }, (_, index) => 100 + index * 0.004); + const y = x.map(() => 1); + const chained = xyArrayAlign([{ x, y }], { delta: 0.005 }); + + expect(chained.x).toHaveLength(1); + expect(chained.x[0]).toBeCloseTo(100.04, 10); + + const result = xyArrayAlignByIntensity([{ x, y }], { delta: 0.005 }); + + // Each slot holds the peak that opened it and the one 0.004 further, never more. + expect(Array.from(result.x)).toBeDeepCloseTo([ + 100.002, 100.01, 100.018, 100.026, 100.034, 100.042, 100.05, 100.058, + 100.066, 100.074, 100.08, + ]); + expect(Array.from(result.ys[0])).toStrictEqual([ + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, + ]); +}); + +test('delta may be a function of x', () => { + const data = [ + { x: [100, 100.5, 1000], y: [1, 1, 1] }, + { x: [100.6, 1005], y: [1, 1] }, + ]; + const result = xyArrayAlignByIntensity(data, { delta: (x) => x * 0.006 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([100.36666666666666, 1002.5]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [2, 1], + [1, 1], + ]); +}); + +test('no data', () => { + expect(xyArrayAlignByIntensity([])).toStrictEqual({ + x: new Float64Array(0), + ys: [], + }); + expect(xyArrayAlignByIntensity([{ x: [], y: [] }])).toStrictEqual({ + x: new Float64Array(0), + ys: [new Float64Array(0)], + }); +}); + +test('x values must be ascending', () => { + expect(() => + xyArrayAlignByIntensity([ + { x: [1, 2], y: [1, 1] }, + { x: [3, 2], y: [1, 1] }, + ]), + ).toThrow('x values of the spectrum 1 must be ascending'); +}); + +// The four spectra share the same x, but the big peak sits at 1 for the first two +// and at 1.5 for the last two. +const twoCompetingPeaks = [ + { x: [1, 1.5, 10], y: [100, 1, 1] }, + { x: [1, 1.5, 10], y: [10, 1, 1] }, + { x: [1, 1.5, 10], y: [1, 10, 1] }, + { x: [1, 1.5, 10], y: [1, 100, 1] }, +]; + +test('competing peaks closer than delta end up in one slot', () => { + // The 100 at x = 1 opens the first slot, and delta = 1 reaches 1.5, so every + // point of both peaks joins it: 224 of intensity at (112 * 1 + 112 * 1.5) / 224. + const result = xyArrayAlignByIntensity(twoCompetingPeaks); + + expect(Array.from(result.x)).toBeDeepCloseTo([1.25, 10]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [101, 1], + [11, 1], + [11, 1], + [101, 1], + ]); +}); + +test('competing peaks further apart than delta keep their own slot', () => { + // 1.5 is now out of reach of the slot opened at 1, so each peak keeps its own + // and the spectra come back exactly as they went in. + const result = xyArrayAlignByIntensity(twoCompetingPeaks, { delta: 0.3 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([1, 1.5, 10]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [100, 1, 1], + [10, 1, 1], + [1, 10, 1], + [1, 100, 1], + ]); +}); + +test('the weighted centre is pulled towards the more intense spectrum', () => { + // Every peak of the second spectrum sits 0.5 above its counterpart in the first and + // within delta, so each pair shares a slot. The centre is the intensity weighted + // average, so it lands near the strong spectrum rather than halfway: a peak of 10 at + // 1 and one of 1 at 1.5 give (10 * 1 + 1 * 1.5) / 11. The peaks at 10 coincide. + const result = xyArrayAlignByIntensity([ + { x: [1, 5, 10], y: [10, 10, 1] }, + { x: [1.5, 5.5, 10], y: [1, 1, 1] }, + ]); + + expect(Array.from(result.x)).toBeDeepCloseTo([11.5 / 11, 55.5 / 11, 10]); + expect(result.ys.map((y) => Array.from(y))).toStrictEqual([ + [10, 10, 1], + [1, 1, 1], + ]); +}); diff --git a/src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts b/src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts new file mode 100644 index 000000000..6c798582e --- /dev/null +++ b/src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts @@ -0,0 +1,116 @@ +import type { DataXY } from 'cheminfo-types'; +import { XSadd } from 'ml-xsadd'; +import { expect, test } from 'vitest'; + +import { xyArrayMergeByIntensity } from '../xyArrayMergeByIntensity.ts'; + +test('the x of a slot is the intensity weighted average of its peaks', () => { + const data = [ + { x: [100.002, 200.01], y: [10, 20] }, + { x: [100.001, 200.02], y: [11, 21] }, + ]; + const result = xyArrayMergeByIntensity(data, { delta: 0.02 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([ + (100.002 * 10 + 100.001 * 11) / 21, + (200.01 * 20 + 200.02 * 21) / 41, + ]); + expect(Array.from(result.y)).toStrictEqual([21, 41]); + expect(Array.from(result.number)).toStrictEqual([2, 2]); + expect(Array.from(result.from)).toStrictEqual([100.001, 200.01]); + expect(Array.from(result.to)).toStrictEqual([100.002, 200.02]); +}); + +test('the strongest peak takes its neighbours before the weaker ones do', () => { + // 5 is closer to 4 than to 10, but 10 is the most intense peak of the set and + // takes it first, which leaves 4 alone. + const data = [{ x: [4, 5, 10], y: [3, 1, 100] }]; + const result = xyArrayMergeByIntensity(data, { delta: 5.5 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([4, (5 * 1 + 10 * 100) / 101]); + expect(Array.from(result.y)).toStrictEqual([3, 101]); + expect(Array.from(result.number)).toStrictEqual([1, 2]); +}); + +test('a slot may hold several peaks of the same spectrum', () => { + const data = [{ x: [1, 1.1, 1.2], y: [1, 5, 1] }]; + const result = xyArrayMergeByIntensity(data, { delta: 0.2 }); + + expect(Array.from(result.x)).toBeDeepCloseTo([1.1]); + expect(Array.from(result.y)).toStrictEqual([7]); + expect(Array.from(result.number)).toStrictEqual([3]); +}); + +test('no data', () => { + const result = xyArrayMergeByIntensity([]); + + expect(Array.from(result.x)).toStrictEqual([]); + expect(Array.from(result.y)).toStrictEqual([]); +}); + +test('many drifting spectra keep the total intensity and stay within delta', () => { + const delta = 0.005; + const random = new XSadd(42); + const trueMasses = Array.from( + { length: 200 }, + (_, index) => 100 + index * 0.5, + ); + const data: DataXY[] = []; + let totalIntensity = 0; + for (let spectrum = 0; spectrum < 500; spectrum++) { + const x: number[] = []; + const y: number[] = []; + for (const mass of trueMasses) { + if (random.random() < 0.3) continue; + x.push(mass + (random.random() - 0.5) * 0.004); + const intensity = 1 + random.random() * 1000; + y.push(intensity); + totalIntensity += intensity; + } + data.push({ x, y }); + } + + const result = xyArrayMergeByIntensity(data, { delta }); + + expect(result.x).toHaveLength(200); + + let sum = 0; + let numberOfPeaks = 0; + for (let i = 0; i < result.x.length; i++) { + sum += result.y[i]; + numberOfPeaks += result.number[i]; + + expect(result.to[i] - result.from[i]).toBeLessThanOrEqual(2 * delta); + expect(result.x[i]).toBeCloseTo(trueMasses[i], 2); + } + + for (let i = 1; i < result.x.length; i++) { + expect(result.x[i]).toBeGreaterThan(result.x[i - 1]); + } + + expect(sum).toBeCloseTo(totalIntensity, 6); + expect(numberOfPeaks).toBe(data.reduce((total, s) => total + s.x.length, 0)); +}); + +test('slots with no intensity keep the x of the peak that opened them', () => { + const result = xyArrayMergeByIntensity([{ x: [1, 1.05, 5], y: [0, 0, 0] }], { + delta: 0.1, + }); + + expect(Array.from(result.x)).toStrictEqual([1, 5]); + expect(Array.from(result.y)).toStrictEqual([0, 0]); + expect(Array.from(result.number)).toStrictEqual([2, 1]); +}); + +test('more slots than the initially allocated ones', () => { + const length = 5000; + const x = Array.from({ length }, (_, index) => index * 10); + const y = Array.from({ length }, (_, index) => (index % 7) + 1); + + const result = xyArrayMergeByIntensity([{ x, y }], { delta: 1 }); + + expect(result.x).toHaveLength(length); + expect(Array.from(result.x)).toStrictEqual(x); + expect(Array.from(result.y)).toStrictEqual(y); + expect(Array.from(result.number)).toStrictEqual(x.map(() => 1)); +}); diff --git a/src/xyArray/index.ts b/src/xyArray/index.ts index a084fb3dd..aea413682 100644 --- a/src/xyArray/index.ts +++ b/src/xyArray/index.ts @@ -1,4 +1,8 @@ export * from './xyArrayAlign.ts'; +export * from './xyArrayAlignByIntensity.ts'; export * from './xyArrayAlignToFirst.ts'; +export * from './xyArrayFilterMinYValue.ts'; +export * from './xyArrayMaxY.ts'; export * from './xyArrayMerge.ts'; +export * from './xyArrayMergeByIntensity.ts'; export * from './xyArrayWeightedMerge.ts'; diff --git a/src/xyArray/utils/__tests__/mergeSortedXY.test.ts b/src/xyArray/utils/__tests__/mergeSortedXY.test.ts new file mode 100644 index 000000000..2685a9a73 --- /dev/null +++ b/src/xyArray/utils/__tests__/mergeSortedXY.test.ts @@ -0,0 +1,126 @@ +import { expect, test } from 'vitest'; + +import { mergeSortedXY } from '../mergeSortedXY.ts'; + +test('merges the points of every spectrum ascending by x', () => { + const merged = mergeSortedXY([ + { x: [1, 4, 7], y: [10, 40, 70] }, + { x: [], y: [] }, + { x: [2, 4.5, 9], y: [20, 45, 90] }, + { x: [0], y: [0] }, + ]); + + expect(Array.from(merged.x)).toStrictEqual([0, 1, 2, 4, 4.5, 7, 9]); + expect(Array.from(merged.y)).toStrictEqual([0, 10, 20, 40, 45, 70, 90]); + expect(Array.from(merged.spectrumIndex)).toStrictEqual([3, 0, 2, 0, 2, 0, 2]); +}); + +test('points sharing an x are all kept', () => { + const merged = mergeSortedXY([ + { x: [4], y: [40] }, + { x: [4], y: [41] }, + ]); + + expect(Array.from(merged.x)).toStrictEqual([4, 4]); + expect(Array.from(merged.y).toSorted((a, b) => a - b)).toStrictEqual([ + 40, 41, + ]); + expect( + Array.from(merged.spectrumIndex).toSorted((a, b) => a - b), + ).toStrictEqual([0, 1]); +}); + +test('no data', () => { + const merged = mergeSortedXY([]); + + expect(Array.from(merged.x)).toStrictEqual([]); + expect(Array.from(merged.spectrumIndex)).toStrictEqual([]); +}); + +test('x values must be ascending', () => { + expect(() => mergeSortedXY([{ x: [3, 1], y: [1, 1] }])).toThrow( + 'x values of the spectrum 0 must be ascending', + ); +}); + +// Thirty spectra or more are merged by ordering the concatenated points instead of +// by the heap, so the cases below have to reach that many to exercise it. +const MANY_SPECTRA = 30; + +/** + * Builds spectra whose x values interleave, one point of each per unit. + * @param numberOfSpectra - how many spectra. + * @returns the spectra. + */ +function interleavedSpectra(numberOfSpectra: number) { + return Array.from({ length: numberOfSpectra }, (_, spectrum) => ({ + x: [spectrum / numberOfSpectra, 1 + spectrum / numberOfSpectra], + y: [spectrum, 100 + spectrum], + })); +} + +test('many spectra are merged ascending by x', () => { + const merged = mergeSortedXY(interleavedSpectra(MANY_SPECTRA)); + + expect(merged.x).toHaveLength(2 * MANY_SPECTRA); + expect(Array.from(merged.x)).toStrictEqual( + Array.from(merged.x).toSorted((a, b) => a - b), + ); + // Each spectrum contributes its low point, in order, then its high point. + expect(Array.from(merged.spectrumIndex)).toStrictEqual([ + ...Array.from({ length: MANY_SPECTRA }, (_, i) => i), + ...Array.from({ length: MANY_SPECTRA }, (_, i) => i), + ]); + expect(Array.from(merged.y)).toStrictEqual([ + ...Array.from({ length: MANY_SPECTRA }, (_, i) => i), + ...Array.from({ length: MANY_SPECTRA }, (_, i) => 100 + i), + ]); +}); + +test('both strategies agree on the same points', () => { + // Empty spectra contribute no point and keep the indices of the real ones, so + // padding to the threshold sends the very same data through the other strategy. + const spectra = interleavedSpectra(10); + const padded = [ + ...spectra, + ...Array.from({ length: MANY_SPECTRA - spectra.length }, () => ({ + x: [] as number[], + y: [] as number[], + })), + ]; + + expect(mergeSortedXY(padded)).toStrictEqual(mergeSortedXY(spectra)); +}); + +test('many spectra still reject descending x values', () => { + const spectra = interleavedSpectra(MANY_SPECTRA); + spectra[7].x = [5, 2]; + + expect(() => mergeSortedXY(spectra)).toThrow( + 'x values of the spectrum 7 must be ascending', + ); +}); + +test('many spectra keep every point sharing an x', () => { + const spectra = Array.from({ length: MANY_SPECTRA }, (_, spectrum) => ({ + x: [4], + y: [spectrum], + })); + const merged = mergeSortedXY(spectra); + + expect(Array.from(merged.x)).toStrictEqual( + Array.from({ length: MANY_SPECTRA }, () => 4), + ); + expect(Array.from(merged.y).toSorted((a, b) => a - b)).toStrictEqual( + Array.from({ length: MANY_SPECTRA }, (_, i) => i), + ); +}); + +test('many spectra including empty ones', () => { + const spectra = interleavedSpectra(MANY_SPECTRA); + spectra[3] = { x: [], y: [] }; + const merged = mergeSortedXY(spectra); + + expect(merged.x).toHaveLength(2 * MANY_SPECTRA - 2); + expect(Array.from(merged.spectrumIndex)).not.toContain(3); +}); diff --git a/src/xyArray/utils/getIntensitySlots.ts b/src/xyArray/utils/getIntensitySlots.ts new file mode 100644 index 000000000..166adbadc --- /dev/null +++ b/src/xyArray/utils/getIntensitySlots.ts @@ -0,0 +1,183 @@ +import { xDoubleTypedArrayLength } from '../../x/xDoubleTypedArrayLength.ts'; +import { xGetApproximateSortOrder } from '../../x/xGetApproximateSortOrder.ts'; + +import type { MergedXY } from './mergeSortedXY.ts'; +import { sortIntensitySlots } from './sortIntensitySlots.ts'; + +const FREE = 0xffffffff; + +export interface GetIntensitySlotsOptions { + /** + * A point joins a slot when its x is within `delta` of the x of the most intense + * point of that slot. It may also be a function that allows to change `delta` + * depending on the x value. + * @default 1 + */ + delta?: ((arg: number) => number) | number; +} + +export interface IntensitySlots { + /** + * The intensity weighted average of the x of the points of each slot, ascending. + */ + x: Float64Array; + /** + * The sum of the y of the points of each slot. + */ + y: Float64Array; + /** + * The smallest x of each slot. + */ + from: Float64Array; + /** + * The largest x of each slot. + */ + to: Float64Array; + /** + * The number of points of each slot. + */ + number: Uint32Array; + /** + * For each merged point, the index of the slot it belongs to. + */ + slotIndex: Uint32Array; +} + +/** + * Groups points around the most intense ones. + * The most intense point of the whole set opens the first slot and takes every + * point within `delta` of it, the most intense of the remaining points opens the + * next one, and so on until every point belongs to exactly one slot. Unlike a + * chain of close neighbours, a slot can therefore never grow wider than `2 * delta`. + * @param merged - all the points of all the spectra, ascending by x. + * @param options - options. + * @returns the slots, ascending by x. + */ +export function getIntensitySlots( + merged: MergedXY, + options: GetIntensitySlotsOptions = {}, +): IntensitySlots { + const { delta = 1 } = options; + const getDelta = typeof delta === 'function' ? delta : () => delta; + const { x, y } = merged; + const length = x.length; + + const slotIndex = new Uint32Array(length).fill(FREE); + // Union-find over the still free positions, so that claiming a point is done + // once and the scan of a window never walks over the points already taken. + const nextFree = new Uint32Array(length + 1); + const previousFree = new Uint32Array(length + 1); + for (let i = 0; i <= length; i++) { + nextFree[i] = i; + previousFree[i] = i; + } + + let capacity = 1024; + let sumY = new Float64Array(capacity); + let sumXY = new Float64Array(capacity); + let seeds = new Float64Array(capacity); + let from = new Float64Array(capacity); + let to = new Float64Array(capacity); + let number = new Uint32Array(capacity); + let numberOfSlots = 0; + + const order = xGetApproximateSortOrder(y, { descending: true }); + for (let i = 0; i < length; i++) { + const seed = order[i]; + if (slotIndex[seed] !== FREE) continue; + + if (numberOfSlots === capacity) { + capacity *= 2; + sumY = xDoubleTypedArrayLength(sumY); + sumXY = xDoubleTypedArrayLength(sumXY); + seeds = xDoubleTypedArrayLength(seeds); + from = xDoubleTypedArrayLength(from); + to = xDoubleTypedArrayLength(to); + number = xDoubleTypedArrayLength(number); + } + + const slot = numberOfSlots++; + const center = x[seed]; + const currentDelta = getDelta(center); + const lowerLimit = center - currentDelta; + const upperLimit = center + currentDelta; + + let slotSumY = 0; + let slotSumXY = 0; + let slotNumber = 0; + let slotFrom = center; + let slotTo = center; + + let position = seed; + while (position < length && x[position] <= upperLimit) { + slotIndex[position] = slot; + nextFree[position] = position + 1; + previousFree[position + 1] = position; + slotSumY += y[position]; + slotSumXY += x[position] * y[position]; + slotNumber++; + slotTo = x[position]; + position = findNext(nextFree, position + 1); + } + + position = findPrevious(previousFree, seed); + while (position > 0 && x[position - 1] >= lowerLimit) { + const claimed = position - 1; + slotIndex[claimed] = slot; + nextFree[claimed] = claimed + 1; + previousFree[claimed + 1] = claimed; + slotSumY += y[claimed]; + slotSumXY += x[claimed] * y[claimed]; + slotNumber++; + slotFrom = x[claimed]; + position = findPrevious(previousFree, claimed); + } + + sumY[slot] = slotSumY; + sumXY[slot] = slotSumXY; + seeds[slot] = center; + from[slot] = slotFrom; + to[slot] = slotTo; + number[slot] = slotNumber; + } + + return sortIntensitySlots({ + numberOfSlots, + sumY, + sumXY, + seeds, + from, + to, + number, + slotIndex, + }); +} + +/** + * Follows the union-find to the first free position at or after `position`. + * @param nextFree - union-find of the free positions. + * @param position - where to start looking. + * @returns the first free position, at most the length of the data. + */ +function findNext(nextFree: Uint32Array, position: number): number { + while (nextFree[position] !== position) { + nextFree[position] = nextFree[nextFree[position]]; + position = nextFree[position]; + } + return position; +} + +/** + * Follows the union-find to the first boundary at or before `position` whose + * preceding point is still free. + * @param previousFree - union-find of the free positions. + * @param position - where to start looking. + * @returns the boundary, 0 when no free point remains on the left. + */ +function findPrevious(previousFree: Uint32Array, position: number): number { + while (previousFree[position] !== position) { + previousFree[position] = previousFree[previousFree[position]]; + position = previousFree[position]; + } + return position; +} diff --git a/src/xyArray/utils/mergeSortedXY.ts b/src/xyArray/utils/mergeSortedXY.ts new file mode 100644 index 000000000..890897eeb --- /dev/null +++ b/src/xyArray/utils/mergeSortedXY.ts @@ -0,0 +1,177 @@ +import type { DataXY } from 'cheminfo-types'; + +import { xGetSortOrder } from '../../x/xGetSortOrder.ts'; + +// Below this many spectra the heap wins, above it the sort does. The heap costs +// O(log k) per point while ordering the concatenated points is independent of k, +// and the two meet around thirty spectra. +const MANY_SPECTRA = 30; + +export interface MergedXY { + /** + * All the x values of all the spectra, ascending. + */ + x: Float64Array; + /** + * The y value matching each x. + */ + y: Float64Array; + /** + * The index, in `data`, of the spectrum each point comes from. + */ + spectrumIndex: Uint32Array; +} + +/** + * Merges the points of many spectra into one array of points ascending by x. + * Each spectrum must already have ascending x values. + * @param data - spectra, each one with ascending x values. + * @returns every point of every spectrum, ascending by x. + */ +export function mergeSortedXY(data: DataXY[]): MergedXY { + return data.length < MANY_SPECTRA ? mergeByHeap(data) : mergeBySorting(data); +} + +/** + * Merges the spectra by repeatedly taking the smallest x of a binary heap of the + * spectra, which reads each spectrum in order and costs O(log k) per point. + * @param data - spectra, each one with ascending x values. + * @returns every point of every spectrum, ascending by x. + */ +function mergeByHeap(data: DataXY[]): MergedXY { + const total = getTotalLength(data); + + const x = new Float64Array(total); + const y = new Float64Array(total); + const spectrumIndex = new Uint32Array(total); + + // Binary heap of the spectra, keyed by the x value each one is currently on. + const ids = new Uint32Array(data.length); + const keys = new Float64Array(data.length); + const positions = new Uint32Array(data.length); + let size = 0; + for (let i = 0; i < data.length; i++) { + const spectrum = data[i]; + if (spectrum.x.length === 0) continue; + const key = spectrum.x[0]; + let child = size++; + while (child > 0) { + const parent = (child - 1) >> 1; + if (keys[parent] <= key) break; + keys[child] = keys[parent]; + ids[child] = ids[parent]; + child = parent; + } + keys[child] = key; + ids[child] = i; + } + + let at = 0; + while (size > 0) { + const id = ids[0]; + const spectrum = data[id]; + const position = positions[id]; + const currentX = keys[0]; + x[at] = currentX; + y[at] = spectrum.y[position]; + spectrumIndex[at] = id; + at++; + + const next = position + 1; + positions[id] = next; + + let key: number; + let newId: number; + if (next < spectrum.x.length) { + key = spectrum.x[next]; + if (key < currentX) { + throw new Error(`x values of the spectrum ${id} must be ascending`); + } + newId = id; + } else { + size--; + if (size === 0) break; + key = keys[size]; + newId = ids[size]; + } + + let parent = 0; + for (;;) { + const left = 2 * parent + 1; + if (left >= size) break; + const right = left + 1; + const child = right < size && keys[right] < keys[left] ? right : left; + if (keys[child] >= key) break; + keys[parent] = keys[child]; + ids[parent] = ids[child]; + parent = child; + } + keys[parent] = key; + ids[parent] = newId; + } + + return { x, y, spectrumIndex }; +} + +/** + * Merges the spectra by concatenating them and ordering the result by x, which + * costs the same whatever the number of spectra. + * @param data - spectra, each one with ascending x values. + * @returns every point of every spectrum, ascending by x. + */ +function mergeBySorting(data: DataXY[]): MergedXY { + const total = getTotalLength(data); + + const concatenatedX = new Float64Array(total); + const concatenatedY = new Float64Array(total); + const concatenatedIndex = new Uint32Array(total); + let at = 0; + for (let i = 0; i < data.length; i++) { + const spectrum = data[i]; + checkAscending(spectrum.x, i); + concatenatedX.set(spectrum.x, at); + concatenatedY.set(spectrum.y, at); + concatenatedIndex.fill(i, at, at + spectrum.x.length); + at += spectrum.x.length; + } + + const order = xGetSortOrder(concatenatedX); + const x = new Float64Array(total); + const y = new Float64Array(total); + const spectrumIndex = new Uint32Array(total); + for (let i = 0; i < total; i++) { + const from = order[i]; + x[i] = concatenatedX[from]; + y[i] = concatenatedY[from]; + spectrumIndex[i] = concatenatedIndex[from]; + } + + return { x, y, spectrumIndex }; +} + +/** + * Counts the points of every spectrum. + * @param data - spectra. + * @returns the number of points of all of them. + */ +function getTotalLength(data: DataXY[]): number { + let total = 0; + for (const spectrum of data) { + total += spectrum.x.length; + } + return total; +} + +/** + * Throws when the x values of a spectrum are not ascending, which the heap + * notices on its own but ordering the points would silently accept. + * @param x - the x values of one spectrum. + * @param id - the index of the spectrum, for the message. + */ +function checkAscending(x: DataXY['x'], id: number): void { + for (let i = 1; i < x.length; i++) { + if (x[i] < x[i - 1]) { + throw new Error(`x values of the spectrum ${id} must be ascending`); + } + } +} diff --git a/src/xyArray/utils/sortIntensitySlots.ts b/src/xyArray/utils/sortIntensitySlots.ts new file mode 100644 index 000000000..84aaec432 --- /dev/null +++ b/src/xyArray/utils/sortIntensitySlots.ts @@ -0,0 +1,54 @@ +import { xGetSortOrder } from '../../x/xGetSortOrder.ts'; + +import type { IntensitySlots } from './getIntensitySlots.ts'; + +export interface RawIntensitySlots { + numberOfSlots: number; + sumY: Float64Array; + sumXY: Float64Array; + seeds: Float64Array; + from: Float64Array; + to: Float64Array; + number: Uint32Array; + slotIndex: Uint32Array; +} + +/** + * Turns the slots built by decreasing intensity into slots ascending by x. + * @param raw - the slots, by decreasing intensity of the point that opened them. + * @returns the slots, ascending by x, and the slot of each merged point. + */ +export function sortIntensitySlots(raw: RawIntensitySlots): IntensitySlots { + const { numberOfSlots, sumY, sumXY, seeds, slotIndex } = raw; + + const centers = new Float64Array(numberOfSlots); + for (let i = 0; i < numberOfSlots; i++) { + // A slot with no intensity at all has no weighted average, it stays where + // the point that opened it was. + centers[i] = sumY[i] === 0 ? seeds[i] : sumXY[i] / sumY[i]; + } + + const bySlot = xGetSortOrder(centers); + + const x = new Float64Array(numberOfSlots); + const y = new Float64Array(numberOfSlots); + const from = new Float64Array(numberOfSlots); + const to = new Float64Array(numberOfSlots); + const number = new Uint32Array(numberOfSlots); + const rankOfSlot = new Uint32Array(numberOfSlots); + for (let rank = 0; rank < numberOfSlots; rank++) { + const slot = bySlot[rank]; + rankOfSlot[slot] = rank; + x[rank] = centers[slot]; + y[rank] = sumY[slot]; + from[rank] = raw.from[slot]; + to[rank] = raw.to[slot]; + number[rank] = raw.number[slot]; + } + + for (let i = 0; i < slotIndex.length; i++) { + slotIndex[i] = rankOfSlot[slotIndex[i]]; + } + + return { x, y, from, to, number, slotIndex }; +} diff --git a/src/xyArray/xyArrayAlignByIntensity.ts b/src/xyArray/xyArrayAlignByIntensity.ts new file mode 100644 index 000000000..beb1b2d04 --- /dev/null +++ b/src/xyArray/xyArrayAlignByIntensity.ts @@ -0,0 +1,56 @@ +import type { DataXY, DataXYs } from 'cheminfo-types'; + +import { getIntensitySlots } from './utils/getIntensitySlots.ts'; +import { mergeSortedXY } from './utils/mergeSortedXY.ts'; + +export interface XYArrayAlignByIntensityOptions { + /** + * A peak joins a slot when its x is within `delta` of the x of the most intense + * peak of that slot, so a slot is never wider than `2 * delta`. It may also be a + * function that allows to change `delta` depending on the x value. + * @default 1 + */ + delta?: ((arg: number) => number) | number; +} + +/** + * Aligns spectra on the most intense peaks, which is meant for centroided data + * whose x values drift slightly from one spectrum to the next. + * + * The most intense peak of the whole set defines the first common x, and takes + * every peak within `delta` of it. The most intense of the remaining peaks defines + * the next one, and so on until every peak has been placed. The reported x of a + * slot is the intensity weighted average of the peaks it holds, and its y is their + * sum, so no intensity is created or lost. Because a slot only ever reaches + * `delta` away from the peak that opened it, dense data cannot chain into a single + * wide slot the way `xyArrayAlign` does. + * + * The x values of each spectrum must be ascending. + * + * The result is a dense matrix: it holds `data.length * x.length` numbers, which + * is out of reach for tens of thousands of spectra. Use + * `xyArrayMergeByIntensity` to get the common peaks alone, and `xysFilter` to + * drop the slots too few spectra contain. + * @param data - data + * @param options - options. + * @returns the common `x` axis and one aligned `y` array per input spectrum. + */ +export function xyArrayAlignByIntensity( + data: DataXY[], + options: XYArrayAlignByIntensityOptions = {}, +): DataXYs { + const { delta = 1 } = options; + + const merged = mergeSortedXY(data); + const slots = getIntensitySlots(merged, { delta }); + + const x = slots.x; + const ys = Array.from(data, () => new Float64Array(x.length)); + const { slotIndex } = slots; + const { spectrumIndex, y: mergedY } = merged; + for (let i = 0; i < slotIndex.length; i++) { + ys[spectrumIndex[i]][slotIndex[i]] += mergedY[i]; + } + + return { x, ys }; +} diff --git a/src/xyArray/xyArrayMergeByIntensity.ts b/src/xyArray/xyArrayMergeByIntensity.ts new file mode 100644 index 000000000..a66d4e0e7 --- /dev/null +++ b/src/xyArray/xyArrayMergeByIntensity.ts @@ -0,0 +1,49 @@ +import type { DataXY } from 'cheminfo-types'; + +import type { IntensitySlots } from './utils/getIntensitySlots.ts'; +import { getIntensitySlots } from './utils/getIntensitySlots.ts'; +import { mergeSortedXY } from './utils/mergeSortedXY.ts'; + +export interface XYArrayMergeByIntensityOptions { + /** + * A peak joins a slot when its x is within `delta` of the x of the most intense + * peak of that slot, so a slot is never wider than `2 * delta`. It may also be a + * function that allows to change `delta` depending on the x value. + * @default 1 + */ + delta?: ((arg: number) => number) | number; +} + +/** + * The merged peaks: every field of a slot except which slot each input point landed in. + */ +export type XYArrayMergeByIntensityResult = Omit; + +/** + * Merges the peaks of many spectra around the most intense ones. + * + * This is `xyArrayAlignByIntensity` without the per spectrum y arrays: the + * most intense peak of the whole set takes every peak within `delta` of it, the + * most intense of the remaining peaks opens the next slot, and so on. Every peak + * ends up in exactly one slot, so the sum of `y` is the sum of the intensities of + * all the spectra. + * + * The x values of each spectrum must be ascending. Only the merged peaks are held + * in memory, so it scales to tens of millions of peaks where a dense alignment + * matrix does not. + * @param data - data + * @param options - options. + * @returns the common peaks, ascending by x. + */ +export function xyArrayMergeByIntensity( + data: DataXY[], + options: XYArrayMergeByIntensityOptions = {}, +): XYArrayMergeByIntensityResult { + const { delta = 1 } = options; + + const { x, y, from, to, number } = getIntensitySlots(mergeSortedXY(data), { + delta, + }); + + return { x, y, from, to, number }; +} From de88c4450f7e6aaa365229a5fa0899d8201ef599 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Thu, 10 Sep 2026 15:49:19 +0200 Subject: [PATCH 8/9] chore: add a benchmark for the shapes xGetSortOrder can order node 26 / bun 1.3, ns per point at n = 1M: permuting x and y 555 -> 29, ordering points 445 -> 28, the nth largest y 70 -> 24, sorted values 70 -> 20. --- benchmark/xGetSortOrder.ts | 302 +++++++++++++++++++++++++++++++++++++ 1 file changed, 302 insertions(+) create mode 100644 benchmark/xGetSortOrder.ts diff --git a/benchmark/xGetSortOrder.ts b/benchmark/xGetSortOrder.ts new file mode 100644 index 000000000..519cf2c2a --- /dev/null +++ b/benchmark/xGetSortOrder.ts @@ -0,0 +1,302 @@ +/* eslint-disable no-console */ +import Benchmark from 'benchmark'; +import type { DataXY } from 'cheminfo-types'; +import { XSadd } from 'ml-xsadd'; + +import type { Point } from '../src/types/index.ts'; +import { xGetApproximateSortOrder } from '../src/x/xGetApproximateSortOrder.ts'; +import { xGetSortOrder } from '../src/x/xGetSortOrder.ts'; + +// The four shapes of sorting the library does, each one with the implementation currently +// in src copied here, so one process compares them on the same data: +// 1. permute two parallel arrays (xySortX, xyUniqueX) +// 2. read the nth largest value (xyGetNMaxY, xyFilterTopYValues) +// 3. materialize the sorted values (getSortedFloat64, xNoiseSanPlot, getSlots) +// 4. order an array of points (xyObjectSortX, xyObjectBestPoints) + +// The library keeps the comparator below RADIX_ORDER_MIN_LENGTH: the crossover measured on +// V8 and JavaScriptCore, for all three input shapes, is around 150 values. +const SIZES = [1_000, 10_000, 100_000, 1_000_000]; + +/** + * A spectrum whose x values are in random order, as xySortX gets them when the points come + * from several sources. + * @param size - number of points. + * @returns unordered x and the y that goes with them. + */ +function makeShuffledSpectrum(size: number): DataXY { + 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 * 0.01 + 100; + y[i] = 10 ** (random() * 5); + } + for (let i = size - 1; i > 0; i--) { + const j = Math.floor(random() * (i + 1)); + let swap = x[i]; + x[i] = x[j]; + x[j] = swap; + swap = y[i]; + y[i] = y[j]; + y[j] = swap; + } + return { x, y }; +} + +// 1. Permute two parallel arrays ------------------------------------------------------------- + +/** + * xySortX as it is today: one object per point, sorted by a comparator, then scattered back. + * @param data - points, x in any order. + * @returns the same points, ascending by x. + */ +function sortXYByObjects(data: DataXY): DataXY { + const { x, y } = data; + const xyObject = Array.from(x, (val, index) => ({ + x: val, + y: y[index], + })); + xyObject.sort((a, b) => a.x - b.x); + + const response = { + x: new Float64Array(x.length), + y: new Float64Array(y.length), + }; + for (let i = 0; i < x.length; i++) { + response.x[i] = xyObject[i].x; + response.y[i] = xyObject[i].y; + } + return response; +} + +/** + * The same permutation read off xGetSortOrder, with no object allocated. + * @param data - points, x in any order. + * @returns the same points, ascending by x. + */ +function sortXYByOrder(data: DataXY): DataXY { + const { x, y } = data; + const length = x.length; + const order = xGetSortOrder(x); + const sortedX = new Float64Array(length); + const sortedY = new Float64Array(length); + for (let i = 0; i < length; i++) { + const at = order[i]; + sortedX[i] = x[at]; + sortedY[i] = y[at]; + } + return { x: sortedX, y: sortedY }; +} + +// 2. Read the nth largest value -------------------------------------------------------------- + +/** + * The threshold xyGetNMaxY uses today: a sorted copy of y, reversed, read at n - 1. + * @param y - intensities. + * @param n - number of points to keep. + * @returns the nth largest intensity. + */ +function thresholdBySortedCopy(y: Float64Array, n: number): number { + const floatY = Float64Array.from(y); + floatY.sort(); + floatY.reverse(); + return floatY[n - 1]; +} + +/** + * The same threshold read through the order, leaving y where it is. + * @param y - intensities. + * @param n - number of points to keep. + * @returns the nth largest intensity. + */ +function thresholdByOrder(y: Float64Array, n: number): number { + const order = xGetSortOrder(y, { descending: true }); + return y[order[n - 1]]; +} + +/** + * The same threshold, ranking only by the high word of each intensity. + * @param y - intensities. + * @param n - number of points to keep. + * @returns the nth largest intensity, to 2^-20 relative. + */ +function thresholdByApproximateOrder(y: Float64Array, n: number): number { + const order = xGetApproximateSortOrder(y, { descending: true }); + return y[order[n - 1]]; +} + +// 3. Materialize the sorted values ----------------------------------------------------------- + +/** + * getSortedFloat64 as it is today. + * @param array - values. + * @returns the values, ascending. + */ +function sortedValuesByTypedSort(array: Float64Array): Float64Array { + const sorted = Float64Array.from(array); + sorted.sort(); + return sorted; +} + +/** + * The same values gathered through the order. + * @param array - values. + * @returns the values, ascending. + */ +function sortedValuesByOrder(array: Float64Array): Float64Array { + const length = array.length; + const order = xGetSortOrder(array); + const sorted = new Float64Array(length); + for (let i = 0; i < length; i++) { + sorted[i] = array[order[i]]; + } + return sorted; +} + +// 4. Order an array of points ---------------------------------------------------------------- + +/** + * xyObjectSortX as it is today, on a copy so the benchmark never re-sorts sorted points. + * @param points - points, x in any order. + * @returns the points, ascending by x. + */ +function sortPointsByComparator(points: Point[]): Point[] { + const sorted = points.slice(); + sorted.sort((a, b) => a.x - b.x); + return sorted; +} + +/** + * The same points gathered through the order of their x values. + * @param points - points, x in any order. + * @returns the points, ascending by x. + */ +function sortPointsByOrder(points: Point[]): Point[] { + const length = points.length; + const x = new Float64Array(length); + for (let i = 0; i < length; i++) { + x[i] = points[i].x; + } + const order = xGetSortOrder(x); + const sorted = new Array(length); + for (let i = 0; i < length; i++) { + sorted[i] = points[order[i]]; + } + // xyObjectSortX orders in place, so the gathered points are written back over the input. + const inPlace = points.slice(); + for (let i = 0; i < length; i++) { + inPlace[i] = sorted[i]; + } + return inPlace; +} + +/** + * Runs one A/B suite and prints ns per element for each case. + * @param title - what is being compared. + * @param size - number of elements one call goes through. + * @param cases - the implementations to compare, the current one first. + */ +function compare( + title: string, + size: number, + cases: Array<[string, () => unknown]>, +): void { + const nanoseconds: Record = {}; + const rmes: Record = {}; + const suite = new Benchmark.Suite(); + for (const [name, run] of cases) { + suite.add(name, run, { minSamples: 30 }); + } + suite + .on('cycle', (event: Benchmark.Event) => { + const { stats, name } = event.target; + if (!stats || !name) return; + nanoseconds[name] = (stats.mean * 1e9) / size; + rmes[name] = stats.rme; + }) + .run(); + + const reference = nanoseconds[cases[0][0]]; + const line = cases + .map(([name]) => { + const ns = nanoseconds[name]; + const speedup = + name === cases[0][0] ? '' : ` ${(reference / ns).toFixed(2)}x`; + return `${name} ${ns.toFixed(2)} ns/pt (±${rmes[name].toFixed(1)}%)${speedup}`; + }) + .join(' | '); + console.log(` ${title.padEnd(22)} ${line}`); +} + +for (const size of SIZES) { + const data = makeShuffledSpectrum(size); + const { x, y } = data; + const n = Math.max(1, Math.round(size / 100)); + + const byObjects = sortXYByObjects(data); + const byOrder = sortXYByOrder(data); + let samePermutation = true; + for (let i = 0; i < size; i++) { + if (byObjects.x[i] !== byOrder.x[i] || byObjects.y[i] !== byOrder.y[i]) { + samePermutation = false; + break; + } + } + + const exactThreshold = thresholdBySortedCopy(y, n); + const orderThreshold = thresholdByOrder(y, n); + const approximateThreshold = thresholdByApproximateOrder(y, n); + + const typedSorted = sortedValuesByTypedSort(x); + const orderSorted = sortedValuesByOrder(x); + let sameValues = true; + for (let i = 0; i < size; i++) { + if (typedSorted[i] !== orderSorted[i]) { + sameValues = false; + break; + } + } + + const points: Point[] = new Array(size); + for (let i = 0; i < size; i++) { + points[i] = { x: x[i], y: y[i] }; + } + const byComparator = sortPointsByComparator(points); + const pointsByOrder = sortPointsByOrder(points); + let samePoints = true; + for (let i = 0; i < size; i++) { + if (byComparator[i] !== pointsByOrder[i]) { + samePoints = false; + break; + } + } + + console.log(`\nn = ${size}`); + compare('permute x and y', size, [ + ['objects', () => sortXYByObjects(data)], + ['order', () => sortXYByOrder(data)], + ]); + console.log(` same points: ${samePermutation}`); + + compare(`${n}th largest y`, size, [ + ['sorted copy', () => thresholdBySortedCopy(y, n)], + ['order', () => thresholdByOrder(y, n)], + ['approximate order', () => thresholdByApproximateOrder(y, n)], + ]); + console.log( + ` threshold: sorted copy ${exactThreshold} | order ${orderThreshold} | approximate ${approximateThreshold}`, + ); + + compare('sorted values', size, [ + ['typed sort', () => sortedValuesByTypedSort(x)], + ['order + gather', () => sortedValuesByOrder(x)], + ]); + console.log(` same values: ${sameValues}`); + + compare('order points', size, [ + ['comparator', () => sortPointsByComparator(points)], + ['order', () => sortPointsByOrder(points)], + ]); + console.log(` same order: ${samePoints}`); +} From 502f73ea51e90c4bfc6e6092620669da020f8f60 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Thu, 10 Sep 2026 15:49:23 +0200 Subject: [PATCH 9/9] perf: order xySortX and xyObjectSortX through xGetSortOrder ns per point on node, from n = 1000 to n = 1M: xySortX 118 -> 41 and 555 -> 29, xyObjectSortX 102 -> 42 and 445 -> 28. Below 256 values the comparator still wins, so it stays. --- src/x/utils/radixSortOrder.ts | 8 +++++++ src/xy/__tests__/xySortX.test.ts | 19 ++++++++++++++++ src/xy/xySortX.ts | 23 +++++++++++++++---- src/xyObject/__tests__/xyObjectSortX.test.ts | 19 ++++++++++++++++ src/xyObject/xyObjectSortX.ts | 24 +++++++++++++++++++- 5 files changed, 87 insertions(+), 6 deletions(-) diff --git a/src/x/utils/radixSortOrder.ts b/src/x/utils/radixSortOrder.ts index 65f05427c..8ca8a0d73 100644 --- a/src/x/utils/radixSortOrder.ts +++ b/src/x/utils/radixSortOrder.ts @@ -10,6 +10,14 @@ const SMALL_DIGIT_BITS = 8; const LARGE_DIGIT_BITS = 16; const LARGE_ARRAY_LENGTH = 1 << LARGE_DIGIT_BITS; +/** + * Length from which ordering through the radix beats a comparator, for a caller that has to + * build the keys and read the values back through the order. The measured crossover is around + * 150 values on both V8 and JavaScriptCore, whether the values come from a typed array, a + * plain array or the x of an array of points. + */ +export const RADIX_ORDER_MIN_LENGTH = 1 << SMALL_DIGIT_BITS; + /** * Orders the values by their full 64 bits, as an LSD radix sort. * @param array - the values to order. diff --git a/src/xy/__tests__/xySortX.test.ts b/src/xy/__tests__/xySortX.test.ts index 6bb690460..618c25313 100644 --- a/src/xy/__tests__/xySortX.test.ts +++ b/src/xy/__tests__/xySortX.test.ts @@ -1,3 +1,4 @@ +import { XSadd } from 'ml-xsadd'; import { expect, test } from 'vitest'; import { xySortX } from '../xySortX.ts'; @@ -60,3 +61,21 @@ test('typed XY arrays', () => { y: Float64Array.from([8, 9, 1, 7, 2, 0, 5, 10, 6]), }); }); + +// The short arrays above are ordered by a comparator, the long one by xGetSortOrder: both +// have to give the same stable order, ties included. +test('same stable order whatever the length', () => { + for (const length of [100, 1000]) { + const { random } = new XSadd(42); + const x = Array.from({ length }, () => Math.round(random() * 20)); + const y = Array.from({ length }, (value, index) => index); + const expected = Array.from({ length }, (value, index) => index).toSorted( + (a, b) => x[a] - x[b], + ); + + const result = xySortX({ x, y }); + + expect(result.x).toStrictEqual(Float64Array.from(expected, (at) => x[at])); + expect(result.y).toStrictEqual(Float64Array.from(expected, (at) => y[at])); + } +}); diff --git a/src/xy/xySortX.ts b/src/xy/xySortX.ts index baf0c9e29..04eb6948c 100644 --- a/src/xy/xySortX.ts +++ b/src/xy/xySortX.ts @@ -1,10 +1,12 @@ import type { DataXY } from 'cheminfo-types'; -import { xIsMonotonic } from '../x/index.ts'; +import { xGetSortOrder, xIsMonotonic } from '../x/index.ts'; +import { RADIX_ORDER_MIN_LENGTH } from '../x/utils/radixSortOrder.ts'; /** * This function performs a quick sort of the x array while transforming the y array to preserve the coordinates. * @param data - object that contains property x (Array) and y (Array). + * @returns the points, ascending by x. */ export function xySortX(data: DataXY): DataXY { const { x, y } = data; @@ -22,16 +24,27 @@ export function xySortX(data: DataXY): DataXY { }; } + const response = { + x: new Float64Array(x.length), + y: new Float64Array(y.length), + }; + + if (x.length >= RADIX_ORDER_MIN_LENGTH) { + const order = xGetSortOrder(x); + for (let i = 0; i < order.length; i++) { + const at = order[i]; + response.x[i] = x[at]; + response.y[i] = y[at]; + } + return response; + } + const xyObject = Array.from(x, (val, index) => ({ x: val, y: y[index], })); xyObject.sort((a, b) => a.x - b.x); - const response = { - x: new Float64Array(x.length), - y: new Float64Array(y.length), - }; for (let i = 0; i < x.length; i++) { response.x[i] = xyObject[i].x; response.y[i] = xyObject[i].y; diff --git a/src/xyObject/__tests__/xyObjectSortX.test.ts b/src/xyObject/__tests__/xyObjectSortX.test.ts index f2780f27e..e3c0882f6 100644 --- a/src/xyObject/__tests__/xyObjectSortX.test.ts +++ b/src/xyObject/__tests__/xyObjectSortX.test.ts @@ -1,3 +1,4 @@ +import { XSadd } from 'ml-xsadd'; import { expect, test } from 'vitest'; import { xyObjectSortX } from '../xyObjectSortX.ts'; @@ -17,3 +18,21 @@ test('xyObjectSortX', () => { { x: 3, y: 3 }, ]); }); + +// The short array above is ordered by a comparator, the long one by xGetSortOrder: both have +// to give the same stable order, and both in place. +test('same stable order whatever the length', () => { + for (const length of [100, 1000]) { + const { random } = new XSadd(42); + const points = Array.from({ length }, (value, index) => ({ + x: Math.round(random() * 20), + y: index, + })); + const expected = points.toSorted((a, b) => a.x - b.x); + + const result = xyObjectSortX(points); + + expect(result).toBe(points); + expect(result).toStrictEqual(expected); + } +}); diff --git a/src/xyObject/xyObjectSortX.ts b/src/xyObject/xyObjectSortX.ts index 95787951c..07d235934 100644 --- a/src/xyObject/xyObjectSortX.ts +++ b/src/xyObject/xyObjectSortX.ts @@ -1,4 +1,6 @@ import type { Point } from '../types/index.ts'; +import { xGetSortOrder } from '../x/index.ts'; +import { RADIX_ORDER_MIN_LENGTH } from '../x/utils/radixSortOrder.ts'; /** * Sorts an array of points in-place. @@ -6,6 +8,26 @@ import type { Point } from '../types/index.ts'; * @returns sorted array of points {x,y}. */ export function xyObjectSortX(points: Point[]): Point[] { - points.sort((a, b) => a.x - b.x); + const { length } = points; + + if (length < RADIX_ORDER_MIN_LENGTH) { + points.sort((a, b) => a.x - b.x); + return points; + } + + const x = new Float64Array(length); + for (let i = 0; i < length; i++) { + x[i] = points[i].x; + } + + const order = xGetSortOrder(x); + const sorted = new Array(length); + for (let i = 0; i < length; i++) { + sorted[i] = points[order[i]]; + } + for (let i = 0; i < length; i++) { + points[i] = sorted[i]; + } + return points; }