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,
+);
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",
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/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<Array<string | number>>,
   dictionary: Record<string, number>,
 ): 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[][];
 }
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;
     }
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<ArrayBuffer> {
+  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.