implement Whitaker smoother and all the utilities to be used outside

src/__tests__/__snapshots__/index.test.ts.snap

@@ -70,6 +70,7 @@ exports[`test existence of exported functions 1`] = `
   "xSum",
   "xUniqueSorted",
   "xVariance",
+  "xWhitakerSmoother",
   "xyAlign",
   "xyCheck",
   "xyCovariance",
@@ -140,9 +141,11 @@ exports[`test existence of exported functions 1`] = `
   "matrixBoxPlot",
   "matrixCenterZMean",
   "matrixCheck",
+  "matrixCholeskySolver",
   "matrixClone",
   "matrixColumnsCorrelation",
   "matrixCreateEmpty",
+  "matrixCuthillMckee",
   "matrixHistogram",
   "matrixMaxAbsoluteZ",
   "matrixMedian",
@@ -167,5 +170,6 @@ exports[`test existence of exported functions 1`] = `
   "nextPowerOfTwo",
   "recursiveResolve",
   "stringify",
+  "calculateAdaptiveWeights",
 ]
 `;

src/matrix/__tests__/matrixCholeskySolver.test.ts

@@ -0,0 +1,51 @@
+import type { NumberArray } from 'cheminfo-types';
+import { expect, test } from 'vitest';
+
+import { addWeights } from '../../utils/addWeights';
+import { createSystemMatrix } from '../../utils/createSystemMatrix';
+import { xSequentialFillFromTo } from '../../x';
+import { matrixCholeskySolver } from '../matrixCholeskySolver';
+import { matrixCuthillMckee } from '../matrixCuthillMckee';
+
+type Func = (b: NumberArray) => NumberArray;
+test('solve a least square system', () => {
+  const x = xSequentialFillFromTo({ from: 0, to: Math.PI, size: 101 });
+  const noise = x.map(() => Math.random() * 0.1 - 0.05);
+  const y = x.map((value, index) => Math.cos(value) + noise[index]);
+  const weights = new Float64Array(y.length).fill(1);
+  y[50] = 0.9; // add outliers
+
+  const lambda = 20;
+  const dimension = x.length;
+  const upperTriangularNonZeros = createSystemMatrix(dimension, lambda);
+
+  const weighted = addWeights(upperTriangularNonZeros, y, weights);
+
+  const cho = matrixCholeskySolver(weighted.leftHandSide, dimension) as Func;
+  expect(cho).not.toBeNull();
+  const smoothed = cho(weighted.rightHandSide);
+  expect(smoothed[50]).toBeLessThan(0.2);
+  expect(smoothed[50]).toBeGreaterThan(-0.2);
+
+  //ignore the outlier, it implicates the smooth should pass closer to zero.
+  weights[50] = 0;
+  const weighted2 = addWeights(upperTriangularNonZeros, y, weights);
+  const cho2 = matrixCholeskySolver(weighted.leftHandSide, dimension) as Func;
+  expect(cho2).not.toBeNull();
+  const smoothed2 = cho2(weighted2.rightHandSide);
+  expect(smoothed2[50]).toBeLessThan(smoothed[50]);
+
+  const permutationEncodedArray = matrixCuthillMckee(
+    upperTriangularNonZeros,
+    dimension,
+  );
+  const cho3 = matrixCholeskySolver(
+    weighted2.leftHandSide,
+    dimension,
+    permutationEncodedArray,
+  ) as Func;
+
+  expect(cho3).not.toBeNull();
+  const smoothed3 = cho2(weighted2.rightHandSide);
+  expect(smoothed3[50]).toEqual(smoothed3[50]);
+});

src/matrix/__tests__/matrixCuthillMckee.test.ts

@@ -0,0 +1,67 @@
+import type { NumberArray } from 'cheminfo-types';
+import { expect, test } from 'vitest';
+
+import { matrixCuthillMckee } from '../matrixCuthillMckee';
+
+test('permutation to reduce bandwidth', () => {
+  const dimension = 10;
+  const matrix = [
+    [2, 2, 1.5],
+    [1, 1, 1],
+    [1, 4, 0.02],
+    [5, 5, 1.2],
+    [7, 7, 1.6],
+    [4, 4, 2.6],
+    [3, 3, 1.1],
+    [4, 7, 0.09],
+    [4, 6, 0.16],
+    [0, 0, 1.7],
+    [4, 8, 0.52],
+    [0, 8, 0.13],
+    [6, 6, 1.3],
+    [7, 8, 0.11],
+    [4, 9, 0.53],
+    [8, 8, 1.4],
+    [9, 9, 3.1],
+    [1, 9, 0.01],
+    [6, 9, 0.56],
+  ];
+
+  const permutationEncodedArray = matrixCuthillMckee(matrix, dimension);
+
+  const reducedBandwidthMatrix = permut(
+    matrix,
+    permutationEncodedArray,
+    dimension,
+  );
+  expect(bandwidth(reducedBandwidthMatrix)).toBeLessThan(bandwidth(matrix));
+});
+
+function permut(
+  nonZerosArray: NumberArray[],
+  permutationEncoded: NumberArray,
+  dimension: number,
+) {
+  const pinv = new Array(dimension);
+  for (let k = 0; k < dimension; k++) {
+    pinv[permutationEncoded[k]] = k;
+  }
+  const mt: NumberArray[] = new Array(nonZerosArray.length);
+  for (let a = 0; a < nonZerosArray.length; ++a) {
+    const [r, c, value] = nonZerosArray[a];
+    const [ar, ac] = [pinv[r], pinv[c]];
+    mt[a] = ac < ar ? [ac, ar, value] : [ar, ac, value];
+  }
+  return mt;
+}
+
+function bandwidth(matrix: NumberArray[]) {
+  let bandwidth = 0;
+  for (const [i, j] of matrix) {
+    const distance = Math.abs(i - j);
+    if (distance > bandwidth) {
+      bandwidth = distance;
+    }
+  }
+  return bandwidth;
+}

src/matrix/index.ts

@@ -4,9 +4,11 @@ export * from './matrixAutoCorrelation';
 export * from './matrixBoxPlot';
 export * from './matrixCenterZMean';
 export * from './matrixCheck';
+export * from './matrixCholeskySolver';
 export * from './matrixClone';
 export * from './matrixColumnsCorrelation';
 export * from './matrixCreateEmpty';
+export * from './matrixCuthillMckee';
 export * from './matrixHistogram';
 export * from './matrixMaxAbsoluteZ';
 export * from './matrixMedian';