feat: equal integration vector

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,
+);

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",

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",

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],
+  ]);
 });

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[][];
 }

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;
     }

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');
+});

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';

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;
+}