From 3f3c382308df69cd7b9465bfce6e9a820d2463f4 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Mon, 20 Jul 2026 18:47:27 +0200 Subject: [PATCH 1/5] feat: add xyEqualIntegrationVector and xEqualIntegrationVectorSimilarity Split a spectrum in parts of equal integration rather than on centers of mass. A center of mass averages positions, so its splits land between the peaks instead of on them (3.32 and 6.68 for peaks at 2.5/5/7.5), and every point is weighted by its own x value, which lets the noise of a long baseline drag the result. The new splits follow the integration, so they point at the signal and hardly move with a noise as high as the peaks. xyMassCenterVector and xMassCenterVectorSimilarity are deprecated in favor of them and keep their exact behaviour, the latter now delegating so that the similarity algorithm exists only once. --- .../__snapshots__/index.test.ts.snap | 2 + .../xEqualIntegrationVectorSimilarity.test.ts | 95 +++++++++++++ src/x/index.ts | 1 + src/x/xEqualIntegrationVectorSimilarity.ts | 94 +++++++++++++ src/x/xMassCenterVectorSimilarity.ts | 65 ++------- .../xyEqualIntegrationVector.test.ts | 130 ++++++++++++++++++ src/xy/index.ts | 1 + src/xy/utils/cumulativeDensity.ts | 108 +++++++++++++++ src/xy/xyEqualIntegrationVector.ts | 90 ++++++++++++ src/xy/xyMassCenterVector.ts | 3 + 10 files changed, 532 insertions(+), 57 deletions(-) create mode 100644 src/x/__tests__/xEqualIntegrationVectorSimilarity.test.ts create mode 100644 src/x/xEqualIntegrationVectorSimilarity.ts create mode 100644 src/xy/__tests__/xyEqualIntegrationVector.test.ts create mode 100644 src/xy/utils/cumulativeDensity.ts create mode 100644 src/xy/xyEqualIntegrationVector.ts diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index ea521371d..7dd000e18 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -31,6 +31,7 @@ exports[`existence of exported functions 1`] = ` "xDivide", "xDotProduct", "xEnsureFloat64", + "xEqualIntegrationVectorSimilarity", "xFindClosestIndex", "xGetFromToIndex", "xGetTargetIndex", @@ -88,6 +89,7 @@ exports[`existence of exported functions 1`] = ` "xyCumulativeDistributionStatistics", "xyEnsureFloat64", "xyEnsureGrowingX", + "xyEqualIntegrationVector", "xyEquallySpaced", "xyExtract", "xyFilter", diff --git a/src/x/__tests__/xEqualIntegrationVectorSimilarity.test.ts b/src/x/__tests__/xEqualIntegrationVectorSimilarity.test.ts new file mode 100644 index 000000000..82375485b --- /dev/null +++ b/src/x/__tests__/xEqualIntegrationVectorSimilarity.test.ts @@ -0,0 +1,95 @@ +import { expect, test } from 'vitest'; + +import { xyEqualIntegrationVector } from '../../xy/index.ts'; +import { xEqualIntegrationVectorSimilarity } from '../xEqualIntegrationVectorSimilarity.ts'; + +// a single peak, the spectrum starts and ends on the baseline +const data = { x: [1, 2, 3, 4, 5, 6, 7], y: [0, 1, 2, 1, 0, 0, 0] }; +// exactly the same spectrum, shifted by 2 and padded with zeros +const shifted = { + x: [1, 2, 3, 4, 5, 6, 7, 8, 9], + y: [0, 0, 0, 1, 2, 1, 0, 0, 0], +}; + +test('basic', () => { + const vector = xyEqualIntegrationVector(data); + + expect(xEqualIntegrationVectorSimilarity(vector, vector)).toBeCloseTo(1); +}); + +test('a global shift is only penalized at the first level', () => { + const vector1 = xyEqualIntegrationVector(data, { depth: 2 }); + const vector2 = xyEqualIntegrationVector(shifted, { depth: 2 }); + + // zero padding does not change the integration, so the vector is exactly translated + expect(vector2).toBeDeepCloseTo( + [...vector1].map((value) => value + 2), + 10, + ); + // only the first of the 2 levels differs once the tree is recentered + expect(xEqualIntegrationVectorSimilarity(vector1, vector2)).toBeCloseTo(0.5); +}); + +test('a global shift, deeper tree', () => { + const similarity = xEqualIntegrationVectorSimilarity( + xyEqualIntegrationVector(data, { depth: 3 }), + xyEqualIntegrationVector(shifted, { depth: 3 }), + ); + + expect(similarity).toBeCloseTo(2 / 3); +}); + +test('without recentering a global shift matches nothing', () => { + const vector1 = xyEqualIntegrationVector(data, { depth: 3 }); + const vector2 = xyEqualIntegrationVector(shifted, { depth: 3 }); + const clone1 = vector1.slice(); + const clone2 = vector2.slice(); + + expect( + xEqualIntegrationVectorSimilarity(vector1, vector2, { recenter: false }), + ).toBe(0); + // the inputs must never be modified + expect(vector1).toStrictEqual(clone1); + expect(vector2).toStrictEqual(clone2); +}); + +test('a partial shift, only the second half moves', () => { + const data1 = { + x: [1, 2, 3, 4, 5, 6, 7, 8, 9], + y: [1, 2, 2, 1, 0, 1, 2, 2, 1], + }; + const data2 = { + x: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], + y: [1, 2, 2, 1, 0, 0, 0, 1, 2, 2, 1], + }; + + expect( + xEqualIntegrationVectorSimilarity( + xyEqualIntegrationVector(data1, { depth: 2 }), + xyEqualIntegrationVector(data2, { depth: 2 }), + ), + ).toBeCloseTo(0.75); + expect( + xEqualIntegrationVectorSimilarity( + xyEqualIntegrationVector(data1, { depth: 3 }), + xyEqualIntegrationVector(data2, { depth: 3 }), + ), + ).toBeCloseTo(5 / 6); +}); + +test('a kind function for similarity', () => { + const similarity = xEqualIntegrationVectorSimilarity( + xyEqualIntegrationVector(data, { depth: 3 }), + xyEqualIntegrationVector(shifted, { depth: 3 }), + { similarityFct: (a, b) => (b - a < 0.1 ? 1 : (b - a) / (b + a)) }, + ); + + // the shift of 2 scores (b - a) / (b + a) = 0.25 on the first level, the two others match + expect(similarity).toBeCloseTo(0.25 / 3 + 1 / 3 + 1 / 3); +}); + +test('the array length must be a power of 2 minus 1', () => { + expect(() => + xEqualIntegrationVectorSimilarity([1, 2, 3], [1, 2, 3, 4]), + ).toThrow('the array length is not a power of 2 minus 1'); +}); diff --git a/src/x/index.ts b/src/x/index.ts index 2c2b3cb3c..85949fa41 100644 --- a/src/x/index.ts +++ b/src/x/index.ts @@ -18,6 +18,7 @@ export * from './xDistributionStats.ts'; export * from './xDivide.ts'; export * from './xDotProduct.ts'; export * from './xEnsureFloat64.ts'; +export * from './xEqualIntegrationVectorSimilarity.ts'; export * from './xFindClosestIndex.ts'; export * from './xGetFromToIndex.ts'; export * from './xGetTargetIndex.ts'; diff --git a/src/x/xEqualIntegrationVectorSimilarity.ts b/src/x/xEqualIntegrationVectorSimilarity.ts new file mode 100644 index 000000000..4cd8ad1eb --- /dev/null +++ b/src/x/xEqualIntegrationVectorSimilarity.ts @@ -0,0 +1,94 @@ +import type { NumberArray } from 'cheminfo-types'; + +export interface XEqualIntegrationVectorSimilarityOptions { + /** + * Function that based on the difference between the two values, return a similarity score between 0 and 1 + * The default one matches values that are identical up to the precision of a float64. The splits are the + * solution of a dichotomic search, so two spectra that are the same up to a shift do not give bit to bit + * identical values. + * @default (a, b) => (Math.abs(a - b) <= 1e-9 * Math.max(1, Math.abs(a), Math.abs(b)) ? 1 : 0) + */ + similarityFct?: (a: number, b: number) => number; + + /** + * Should we recenter the tree ? + * @default: true + */ + recenter?: boolean; +} + +/** + * Check the similarity between array created by xyEqualIntegrationVector + * @param array1 - first equal integration vector. + * @param array2 - second equal integration vector. + * @param options - options. + * @returns similarity score between 0 and 1. + */ +export function xEqualIntegrationVectorSimilarity( + array1: NumberArray, + array2: NumberArray, + options: XEqualIntegrationVectorSimilarityOptions = {}, +): number { + const { + recenter = true, + similarityFct = (a: number, b: number) => + Math.abs(a - b) <= 1e-9 * Math.max(1, Math.abs(a), Math.abs(b)) ? 1 : 0, + } = options; + const depth1 = getDepth(array1); + const depth2 = getDepth(array2); + const depth = Math.min(depth1, depth2); + + // need to copy the array because we shift the data in place if recenter is true + if (recenter) { + array1 = array1.slice(); + } + + let similarity = 0; + // we will compare level by level + // and recenter the array at each level + + for (let level = 0; level < depth; level++) { + const maxSimilarity = 1 / depth / (1 << level); + + for (let slot = 0; slot < 1 << level; slot++) { + const index = (1 << level) - 1 + slot; + const value1 = array1[index]; + const value2 = array2[index]; + similarity += similarityFct(value1, value2) * maxSimilarity; + if (recenter) { + shiftSubtree(array1, depth, level, slot, value2 - value1); + } + } + } + return similarity; +} + +function shiftSubtree( + array: NumberArray, + depth: number, + level: number, + slot: number, + shift: number, +) { + for (let currentLevel = level; currentLevel < depth; currentLevel++) { + const levelSlotShift = slot * (1 << (currentLevel - level)); + const levelShift = (1 << currentLevel) - 1; + const levelSlotSize = 1 << (currentLevel - level); + for ( + let slotIndex = levelSlotShift; + slotIndex < levelSlotShift + levelSlotSize; + slotIndex++ + ) { + const index = levelShift + slotIndex; + array[index] += shift; + } + } +} + +function getDepth(array: NumberArray) { + const depth = Math.log2(array.length + 1); + if (depth % 1 !== 0) { + throw new Error('the array length is not a power of 2 minus 1'); + } + return depth; +} diff --git a/src/x/xMassCenterVectorSimilarity.ts b/src/x/xMassCenterVectorSimilarity.ts index 4c4d04c9c..69291729a 100644 --- a/src/x/xMassCenterVectorSimilarity.ts +++ b/src/x/xMassCenterVectorSimilarity.ts @@ -1,5 +1,7 @@ import type { NumberArray } from 'cheminfo-types'; +import { xEqualIntegrationVectorSimilarity } from './xEqualIntegrationVectorSimilarity.ts'; + interface XMassCenterVectorSimilarityOptions { /** * Function that based on the difference between the two values, return a similarity score between 0 and 1 @@ -16,6 +18,8 @@ interface XMassCenterVectorSimilarityOptions { /** * Check the similarity between array created by xyMassCenterVector + * @deprecated The center of mass trees it compares are misleading, see `xyMassCenterVector`. + * Use `xEqualIntegrationVectorSimilarity` together with `xyEqualIntegrationVector`. * @param array1 - first mass center vector. * @param array2 - second mass center vector. * @param options - options. @@ -30,61 +34,8 @@ export function xMassCenterVectorSimilarity( recenter = true, similarityFct = (a: number, b: number) => (a === b ? 1 : 0), } = options; - const depth1 = getDepth(array1); - const depth2 = getDepth(array2); - const depth = Math.min(depth1, depth2); - - // need to copy the array because we shift the data in place if recenter is true - if (recenter) { - array1 = array1.slice(); - } - - let similarity = 0; - // we will compare level by level - // and recenter the array at each level - - for (let level = 0; level < depth; level++) { - const maxSimilarity = 1 / depth / (1 << level); - - for (let slot = 0; slot < 1 << level; slot++) { - const index = (1 << level) - 1 + slot; - const value1 = array1[index]; - const value2 = array2[index]; - similarity += similarityFct(value1, value2) * maxSimilarity; - if (recenter) { - shiftSubtree(array1, depth, level, slot, value2 - value1); - } - } - } - return similarity; -} - -function shiftSubtree( - array: NumberArray, - depth: number, - level: number, - slot: number, - shift: number, -) { - for (let currentLevel = level; currentLevel < depth; currentLevel++) { - const levelSlotShift = slot * (1 << (currentLevel - level)); - const levelShift = (1 << currentLevel) - 1; - const levelSlotSize = 1 << (currentLevel - level); - for ( - let slotIndex = levelSlotShift; - slotIndex < levelSlotShift + levelSlotSize; - slotIndex++ - ) { - const index = levelShift + slotIndex; - array[index] += shift; - } - } -} - -function getDepth(array: NumberArray) { - const depth = Math.log2(array.length + 1); - if (depth % 1 !== 0) { - throw new Error('the array length is not a power of 2 minus 1'); - } - return depth; + return xEqualIntegrationVectorSimilarity(array1, array2, { + recenter, + similarityFct, + }); } diff --git a/src/xy/__tests__/xyEqualIntegrationVector.test.ts b/src/xy/__tests__/xyEqualIntegrationVector.test.ts new file mode 100644 index 000000000..3c1a60dcb --- /dev/null +++ b/src/xy/__tests__/xyEqualIntegrationVector.test.ts @@ -0,0 +1,130 @@ +import { SpectrumGenerator } from 'spectrum-generator'; +import { expect, test } from 'vitest'; + +import { xyEqualIntegrationVector } from '../xyEqualIntegrationVector.ts'; +import { xyIntegration } from '../xyIntegration.ts'; + +/** + * 3 peaks of the same height at 2.5, 5 and 7.5, on a 0 to 10 window. + * @param percent - percentage of uniform noise added to the peaks. + * @param seed - seed of the pseudo random generator. + * @returns the spectrum. + */ +function threePeaks(percent = 0, seed = 1) { + const generator = new SpectrumGenerator({ from: 0, to: 10, nbPoints: 10001 }); + generator.addPeak({ x: 2.5, y: 1, width: 0.1 }); + generator.addPeak({ x: 5, y: 1, width: 0.1 }); + generator.addPeak({ x: 7.5, y: 1, width: 0.1 }); + if (percent > 0) { + generator.addNoise({ distribution: 'uniform', seed, percent }); + } + return generator.getSpectrum(); +} + +test('a constant density is split in equal parts', () => { + const data = { x: [1, 2, 3, 4, 5], y: [1, 1, 1, 1, 1] }; + + // 1/2, then 1/4 and 3/4, then 1/8, 3/8, 5/8 and 7/8 of the width + expect(xyEqualIntegrationVector(data, { depth: 3 })).toBeDeepCloseTo([ + 3, 2, 4, 1.5, 2.5, 3.5, 4.5, + ]); +}); + +test('the integration is really the same on both sides', () => { + const generator = new SpectrumGenerator({ + from: 0, + to: 10, + nbPoints: 10001, + }); + generator.addPeak({ x: 2, y: 3, width: 1 }); + generator.addPeak({ x: 7, y: 1, width: 0.6 }); + const data = generator.getSpectrum(); + const [middle, firstQuarter, thirdQuarter] = xyEqualIntegrationVector(data, { + depth: 2, + }); + + // xyIntegration snaps its bounds to the closest point, hence the tolerance + expect(xyIntegration(data, { from: 0, to: middle })).toBeCloseTo( + xyIntegration(data, { from: middle, to: 10 }), + 2, + ); + // each half is split in two quarters that also have the same integration + expect(xyIntegration(data, { from: 0, to: firstQuarter })).toBeCloseTo( + xyIntegration(data, { from: firstQuarter, to: middle }), + 2, + ); + expect(xyIntegration(data, { from: middle, to: thirdQuarter })).toBeCloseTo( + xyIntegration(data, { from: thirdQuarter, to: 10 }), + 2, + ); +}); + +test('the splits land on the peaks, not between them', () => { + // a center of mass reports 3.32 and 6.68 on the very same data, which is bare baseline + expect(xyEqualIntegrationVector(threePeaks(), { depth: 2 })).toBeDeepCloseTo( + [5, 2.5, 7.5], + 1, + ); +}); + +test('3 peaks buried in noise, the splits hardly move', () => { + for (const seed of [1, 2, 3]) { + expect( + xyEqualIntegrationVector(threePeaks(20, seed), { depth: 2 }), + ).toBeDeepCloseTo([5, 2.5, 7.5], 1); + + // now with a noise as high as the peaks themselves + const result = xyEqualIntegrationVector(threePeaks(100, seed), { + depth: 2, + }); + const expected = [5, 2.5, 7.5]; + + for (let i = 0; i < expected.length; i++) { + expect(Math.abs(result[i] - expected[i])).toBeLessThan(0.1); + } + } +}); + +test('a long noisy baseline does not drag the result', () => { + const generator = new SpectrumGenerator({ from: 0, to: 40, nbPoints: 4001 }); + generator.addPeak({ x: 7, y: 1, width: 0.3 }); + const clean = generator.getSpectrum(); + generator.addNoise({ distribution: 'uniform', seed: 42, percent: 10 }); + const noisy = generator.getSpectrum(); + + // a single peak at 7 in a window that is 0..40, so the noise has a long lever arm + expect(xyEqualIntegrationVector(clean, { depth: 1 })[0]).toBeCloseTo(7, 3); + expect(xyEqualIntegrationVector(noisy, { depth: 1 })[0]).toBeCloseTo(7, 2); +}); + +test('range selection with from / to', () => { + const data = { + x: [1, 2, 3, 4, 5, 6, 7, 8, 9], + y: [1, 1, 1, 1, 1, 0, 0, 0, 0], + }; + + expect( + xyEqualIntegrationVector(data, { depth: 2, from: 1, to: 5 }), + ).toStrictEqual( + xyEqualIntegrationVector( + { x: [1, 2, 3, 4, 5], y: [1, 1, 1, 1, 1] }, + { depth: 2 }, + ), + ); +}); + +test('the window may be wider than the data', () => { + const data = { x: [4, 5, 6], y: [0, 1, 0] }; + const result = xyEqualIntegrationVector(data, { depth: 1, from: 0, to: 10 }); + + // the peak is symmetric around 5, the padding carries no integration + expect(result[0]).toBeCloseTo(5, 6); +}); + +test('range of less than 2 points', () => { + const data = { x: [1, 2, 3], y: [1, 1, 1] }; + + expect(() => + xyEqualIntegrationVector(data, { fromIndex: 1, toIndex: 1 }), + ).toThrow('the selected range must contain at least 2 points'); +}); diff --git a/src/xy/index.ts b/src/xy/index.ts index 472429d3f..51141849a 100644 --- a/src/xy/index.ts +++ b/src/xy/index.ts @@ -4,6 +4,7 @@ export * from './xyCovariance.ts'; export * from './xyCumulativeDistributionStatistics.ts'; export * from './xyEnsureFloat64.ts'; export * from './xyEnsureGrowingX.ts'; +export * from './xyEqualIntegrationVector.ts'; export * from './xyEquallySpaced.ts'; export * from './xyExtract.ts'; export * from './xyFilter.ts'; diff --git a/src/xy/utils/cumulativeDensity.ts b/src/xy/utils/cumulativeDensity.ts new file mode 100644 index 000000000..93801ce4f --- /dev/null +++ b/src/xy/utils/cumulativeDensity.ts @@ -0,0 +1,108 @@ +import type { DataXY, NumberArray } from 'cheminfo-types'; + +import { xGetFromToIndex } from '../../x/index.ts'; +import type { XGetFromToIndexOptions } from '../../x/xGetFromToIndex.ts'; + +export interface CumulativeDensity { + x: NumberArray; + y: NumberArray; + fromIndex: number; + toIndex: number; + /** integration from the first point of the range, on every point of the range */ + integral: Float64Array; + /** beginning of the window, that may be outside of the data */ + beginX: number; + /** end of the window, that may be outside of the data */ + endX: number; +} + +/** + * The data seen as a piecewise linear density, with the integration precomputed on every point. + * @param data - object that contains property x (an ordered increasing array) and y (an array). + * @param options - options. + * @returns the cumulated integration over the selected window. + */ +export function getCumulativeDensity( + data: DataXY, + options: XGetFromToIndexOptions = {}, +): CumulativeDensity { + const { x, y } = data; + const { fromIndex, toIndex } = xGetFromToIndex(x, options); + if (toIndex - fromIndex < 1) { + throw new Error('the selected range must contain at least 2 points'); + } + + const { beginX, endX } = getBoundaries(x, options, fromIndex, toIndex); + const length = toIndex - fromIndex + 1; + const integral = new Float64Array(length); + for (let i = 1; i < length; i++) { + const previousIndex = fromIndex + i - 1; + const currentIndex = fromIndex + i; + integral[i] = + integral[i - 1] + + ((y[previousIndex] + y[currentIndex]) / 2) * + (x[currentIndex] - x[previousIndex]); + } + return { x, y, fromIndex, toIndex, integral, beginX, endX }; +} + +/** + * Integration from the beginning of the range up to a specific x value, interpolating inside + * the interval that contains it. + * @param cumulative - cumulated integration of the density. + * @param target - x value. + * @returns the integration. + */ +export function evaluateCumulative( + cumulative: CumulativeDensity, + target: number, +): number { + const { x, y, fromIndex, toIndex, integral } = cumulative; + if (target <= x[fromIndex]) return 0; + const lastIndex = toIndex - fromIndex; + if (target >= x[toIndex]) return integral[lastIndex]; + + let low = 0; + let high = lastIndex; + while (high - low > 1) { + const middle = low + ((high - low) >> 1); + if (x[fromIndex + middle] <= target) { + low = middle; + } else { + high = middle; + } + } + + const startX = x[fromIndex + low]; + const startY = y[fromIndex + low]; + const width = target - startX; + const targetY = + startY + + ((y[fromIndex + low + 1] - startY) * width) / + (x[fromIndex + low + 1] - startX); + return integral[low] + ((startY + targetY) / 2) * width; +} + +/** + * Boundaries of the window, that may be wider than the data itself. + * @param x - array of x values. + * @param options - options. + * @param fromIndex - first index of the selected range. + * @param toIndex - last index of the selected range. + * @returns the x values delimiting the window. + */ +function getBoundaries( + x: NumberArray, + options: XGetFromToIndexOptions, + fromIndex: number, + toIndex: number, +) { + const { from, to } = options; + const firstX = x[fromIndex]; + const lastX = x[toIndex]; + let beginX = + options.fromIndex === undefined && from !== undefined ? from : firstX; + let endX = options.toIndex === undefined && to !== undefined ? to : lastX; + if (beginX > endX) [beginX, endX] = [endX, beginX]; + return { beginX: Math.min(beginX, firstX), endX: Math.max(endX, lastX) }; +} diff --git a/src/xy/xyEqualIntegrationVector.ts b/src/xy/xyEqualIntegrationVector.ts new file mode 100644 index 000000000..75f496a21 --- /dev/null +++ b/src/xy/xyEqualIntegrationVector.ts @@ -0,0 +1,90 @@ +import type { DataXY } from 'cheminfo-types'; + +import type { XGetFromToIndexOptions } from '../x/xGetFromToIndex.ts'; + +import type { CumulativeDensity } from './utils/cumulativeDensity.ts'; +import { + evaluateCumulative, + getCumulativeDensity, +} from './utils/cumulativeDensity.ts'; +import { xyCheck } from './xyCheck.ts'; + +export interface XYEqualIntegrationVectorOptions extends XGetFromToIndexOptions { + /** + * Number of levels of the binary tree. The returned vector has `2 ** depth - 1` entries. + * @default 5 + */ + depth?: number; +} + +/** + * Split a spectrum in parts that all have the same integration and return the x values of the splits. + * The first entry splits the spectrum in two halves of equal integration, the two next ones split each half + * again, etc. This is the same as looking for the 1/2, then 1/4 and 3/4, then 1/8, 3/8, 5/8 and 7/8 fractions + * of the total integration, which is how it is computed. + * This describes a spectrum by a small vector and should provide an efficient way to store and search + * similar spectra, comparing the vectors with `xEqualIntegrationVectorSimilarity`. + * A split never lands between the peaks: it follows the integration, so it points at the signal itself. The + * noise of a long baseline does not drag it either, because a point is never weighted by its own x value. + * The data is considered as a piecewise linear density, so the result varies continuously with the y values + * and hardly depends on the sampling density. + * The window may extend beyond the data, the missing points being considered to have a y value of zero. + * @param data - object that contains property x (an ordered increasing array) and y (an array). + * @param options - options. + * @returns array of x values splitting the spectrum in parts of equal integration. + */ +export function xyEqualIntegrationVector( + data: DataXY, + options: XYEqualIntegrationVectorOptions = {}, +): Float64Array { + xyCheck(data, { minLength: 2 }); + const { depth = 5 } = options; + const cumulative = getCumulativeDensity(data, options); + const { beginX, endX } = cumulative; + const total = evaluateCumulative(cumulative, endX); + + const results = new Float64Array((1 << depth) - 1); + let index = 0; + for (let level = 0; level < depth; level++) { + const nbSlots = 1 << level; + for (let slot = 0; slot < nbSlots; slot++) { + const fraction = (2 * slot + 1) / (2 * nbSlots); + results[index++] = findFraction( + cumulative, + beginX, + endX, + total * fraction, + ); + } + } + + return results; +} + +/** + * Dichotomic search of the x value for which the integration reaches a target. + * @param cumulative - cumulated integration of the density. + * @param beginX - beginning of the window. + * @param endX - end of the window. + * @param target - integration to reach. + * @returns the x value. + */ +function findFraction( + cumulative: CumulativeDensity, + beginX: number, + endX: number, + target: number, +): number { + let low = beginX; + let high = endX; + // the interval is halved 60 times, far below the precision of a float64 + for (let i = 0; i < 60; i++) { + const middle = (low + high) / 2; + if (evaluateCumulative(cumulative, middle) < target) { + low = middle; + } else { + high = middle; + } + } + return (low + high) / 2; +} diff --git a/src/xy/xyMassCenterVector.ts b/src/xy/xyMassCenterVector.ts index f40cbbe02..9cf84a541 100644 --- a/src/xy/xyMassCenterVector.ts +++ b/src/xy/xyMassCenterVector.ts @@ -12,6 +12,9 @@ export interface XYMassCenterVectorOptions { * We will calculate a vector containing center of mass of DataXY as well as center of mass of both parts, etc. * This approach allows to efficiently represent spectra like XPS, NMR, etc. It should provide an extremely efficient * way to store and search similar spectra. + * @deprecated A center of mass averages positions, so a split lands between the peaks instead of on them, and + * every point is weighted by its own x value, which lets the noise of a long baseline drag the result. + * Use `xyEqualIntegrationVector`, that splits the spectrum in parts of equal integration. * @param data - object that contains property x (an ordered increasing array) and y (an array). * @param options - options. * @returns array of centers of mass. From 67467f18a3c55c3e3c69d04650b5826711c97b9a Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Mon, 20 Jul 2026 18:47:30 +0200 Subject: [PATCH 2/5] chore: do not lint .claude eslint . was descending into the agent worktrees under .claude and reporting errors on stale copies of the sources. The directory is already gitignored. --- eslint.config.js | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/eslint.config.js b/eslint.config.js index 00cc43d24..ec43b0709 100644 --- a/eslint.config.js +++ b/eslint.config.js @@ -1,4 +1,7 @@ import { defineConfig, globalIgnores } from 'eslint/config'; import cheminfo from 'eslint-config-cheminfo-typescript'; -export default defineConfig(globalIgnores(['coverage', 'lib']), cheminfo); +export default defineConfig( + globalIgnores(['.claude', 'coverage', 'lib']), + cheminfo, +); From 85096bcff6d24c4b9572cf7d7111056ad062c6a5 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Mon, 20 Jul 2026 18:49:50 +0200 Subject: [PATCH 3/5] fix: matrixApplyNumericalEncoding returned a matrix of zeros It started from matrixCreateEmpty instead of the initial matrix, so no value was ever a string, nothing was encoded and the numeric columns were lost. Clone the initial matrix and encode it in place. The test only checked that every value was a number, which a matrix of zeros satisfies, so it now asserts the expected values. --- .../matrixApplyNumericalDecoding.test.ts | 13 ++++++++++++ src/matrix/matrixApplyNumericalEncoding.ts | 21 ++++++++----------- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/matrix/__tests__/matrixApplyNumericalDecoding.test.ts b/src/matrix/__tests__/matrixApplyNumericalDecoding.test.ts index 07757703c..402f80e62 100644 --- a/src/matrix/__tests__/matrixApplyNumericalDecoding.test.ts +++ b/src/matrix/__tests__/matrixApplyNumericalDecoding.test.ts @@ -28,4 +28,17 @@ test('should return an array of numbers', () => { const nonNumbers = matrix.flat().filter((value) => typeof value !== 'number'); expect(nonNumbers).toHaveLength(0); + + // numeric columns are preserved and strings are mapped through the dictionary + expect(matrix).toStrictEqual([ + [78, 198, 220, 147], + [81, 202, 200, 183], + [88, 221, 222, 177], + [78, 223, 210, 159], + [82, 222, 202, 177], + [86, 224, 218, 175], + [78, 223, 225, 175], + [76, 223, 226, 149], + [96, 200, 216, 192], + ]); }); diff --git a/src/matrix/matrixApplyNumericalEncoding.ts b/src/matrix/matrixApplyNumericalEncoding.ts index 7991a7ae4..c5fd91cf6 100644 --- a/src/matrix/matrixApplyNumericalEncoding.ts +++ b/src/matrix/matrixApplyNumericalEncoding.ts @@ -1,6 +1,6 @@ import { xMaxValue } from '../x/index.ts'; -import { matrixCreateEmpty } from './matrixCreateEmpty.ts'; +import { matrixClone } from './matrixClone.ts'; /** * Numerically encodes the strings in the matrix with an encoding dictionary. @@ -12,11 +12,7 @@ export function matrixApplyNumericalEncoding( matrixInitial: Array>, dictionary: Record, ): number[][] { - const matrix = matrixCreateEmpty({ - nbRows: matrixInitial.length, - nbColumns: matrixInitial[0].length, - ArrayConstructor: Array, - }); + const matrix = matrixClone(matrixInitial); const arrayOfValues: number[] = []; for (const key in dictionary) { @@ -25,18 +21,19 @@ export function matrixApplyNumericalEncoding( let k = xMaxValue(arrayOfValues); for (const row of matrix) { - for (let j = 0; j < matrix[0].length; j++) { - if (typeof row[j] === 'string') { - if (row[j] in dictionary) { - row[j] = dictionary[row[j]]; + for (let j = 0; j < row.length; j++) { + const value = row[j]; + if (typeof value === 'string') { + if (value in dictionary) { + row[j] = dictionary[value]; } else { k++; - dictionary[row[j]] = k; + dictionary[value] = k; row[j] = k; } } } } - return matrix; + return matrix as number[][]; } From c061747605b5826f0cfde5d72a76f92bbc03a0a9 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Mon, 20 Jul 2026 18:49:53 +0200 Subject: [PATCH 4/5] refactor: compute the rescaling factor once in matrixZPivotRescale The factor does not depend on the column, it was recomputed on every iteration. --- src/matrix/matrixZPivotRescale.ts | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/matrix/matrixZPivotRescale.ts b/src/matrix/matrixZPivotRescale.ts index 2d11ca818..9f6447e33 100644 --- a/src/matrix/matrixZPivotRescale.ts +++ b/src/matrix/matrixZPivotRescale.ts @@ -41,10 +41,9 @@ export function matrixZPivotRescale< const newMatrix = matrixCreateEmpty({ nbColumns, nbRows, ArrayConstructor }); const currentMax = matrixMaxAbsoluteZ(matrix); + const factor = max / currentMax; for (let column = 0; column < nbColumns; column++) { - const factor = max / currentMax; - for (let row = 0; row < nbRows; row++) { newMatrix[row][column] = matrix[row][column] * factor; } From fc6405cad7cd243baf208a37406fc19c444410e8 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Mon, 20 Jul 2026 18:50:49 +0200 Subject: [PATCH 5/5] chore: update dependencies --- package.json | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/package.json b/package.json index 47d110328..0fe93cdfc 100644 --- a/package.json +++ b/package.json @@ -43,23 +43,23 @@ "cheminfo-types": "^1.15.0", "fft.js": "^4.0.4", "is-any-array": "^3.0.0", - "ml-matrix": "^6.13.0", + "ml-matrix": "^6.14.0", "ml-xsadd": "^3.0.1" }, "devDependencies": { - "@types/node": "^25.9.3", - "@vitest/coverage-v8": "^4.1.9", + "@types/node": "^26.1.1", + "@vitest/coverage-v8": "^4.1.10", "@zakodium/tsconfig": "^1.0.5", "cheminfo-build": "^1.3.2", "eslint": "^9.39.2", "eslint-config-cheminfo-typescript": "^22.1.0", "jest-matcher-deep-close-to": "^3.0.2", - "ml-spectra-fitting": "^6.1.0", - "prettier": "^3.8.4", + "ml-spectra-fitting": "^6.2.1", + "prettier": "^3.9.5", "rimraf": "^6.1.3", "spectrum-generator": "^8.2.1", "typescript": "^6.0.3", - "vitest": "^4.1.9" + "vitest": "^4.1.10" }, "repository": { "type": "git",