fix(xNoiseSanplot): use analytical expression for Raylegh inverse CDF (#289)

* fix(xNoiseSanplot): use analytical expression for Raylegh inverse cumulative distribution

* chore: reverse the use of callback for sort descending

* refactor: ensure typed array use

* refactor: avoid use let

* chore: add default values in tsdocs

* chore: test scaledFactor options

* refactor: enhance xNoiseSanPlot functionality and improve median calculation

Author: jobo322

package.json

@@ -108,7 +108,6 @@
     "fft.js": "^4.0.4",
     "is-any-array": "^2.0.1",
     "ml-matrix": "^6.12.0",
-    "ml-xsadd": "^3.0.1",
-    "spline-interpolator": "^1.0.0"
+    "ml-xsadd": "^3.0.1"
   }
 }

src/x/__tests__/xNoiseSanPlot.test.ts

@@ -1021,10 +1021,29 @@ test('get noise level', () => {
     ],
   };
 
-  const noise = xNoiseSanPlot(data.re, { magnitudeMode: false });
+  const noise = xNoiseSanPlot(data.re, { fixOffset: true });
   const noiseFromMagnitude = xNoiseSanPlot(reimAbsolute(data), {
     magnitudeMode: true,
   });
-  expect(noise.snr).toBeCloseTo(56.8, 1);
-  expect(noiseFromMagnitude.snr).toBeCloseTo(46.4, 1);
+  expect(noise.snr).toBeCloseTo(56.85, 1);
+  expect(noiseFromMagnitude.snr).toBeCloseTo(40.2, 1);
+
+  const noiseScaled = xNoiseSanPlot(data.re, {
+    scaleFactor: 2,
+  });
+  // the positive and negative noise level should be twice
+  expect(noiseScaled.positive).toBeCloseTo(noise.positive * 2, 0);
+  expect(noiseScaled.negative).toBeCloseTo(noise.negative * 2, 0);
+  //the signal to noise ratio should the same
+  expect(noiseScaled.snr).toBeCloseTo(noise.snr, 1);
+  //use a mask to ignore the bigger peak
+  const length = data.re.length;
+  const mask = Float64Array.from({ length }).fill(0);
+  mask.set(Float64Array.from({ length: 21 }).fill(1), length - 21);
+  const noiseWithoutBigPeaks = xNoiseSanPlot(data.re, {
+    mask,
+    fixOffset: true,
+  });
+  //the SNR should be less because the biggest peak is not present.
+  expect(noiseWithoutBigPeaks.snr).toBeLessThan(noise.snr);
 });

src/x/utils/__tests__/simpleNormInv.test.ts

@@ -0,0 +1,19 @@
+import { describe, it, expect } from 'vitest';
+
+import { simpleNormInvNumber } from '../simpleNormInv';
+
+describe('simpleNormInvNumber', () => {
+  it('should handle normal mode correctly', () => {
+    const result = simpleNormInvNumber(0.25);
+    expect(result).toBeCloseTo(-0.6744897495, 4);
+    const result2 = simpleNormInvNumber(0.025);
+    expect(result2).toBeCloseTo(-1.96, 4);
+  });
+
+  it('should handle magnitude mode correctly', () => {
+    const result = simpleNormInvNumber(0.9, { magnitudeMode: true });
+    // value take from Matlab raylinv function
+    // https://www.mathworks.com/help/stats/raylinv.html
+    expect(result).toBeCloseTo(-2.146, 4);
+  });
+});

src/x/utils/simpleNormInv.ts

@@ -0,0 +1,50 @@
+import type { NumberArray } from 'cheminfo-types';
+
+import erfcinv from './erfcinv';
+
+interface SimpleNormInvOptions {
+  /**
+   * If true the returns values will be calculated from the Rayleigh inverse cumulative distribution.
+   * it returns
+   * @default false
+   */
+  magnitudeMode?: boolean;
+}
+
+/**
+ * Applies a simple normalization inverse transformation to the input data.
+ * @param data - The input array of numbers to be transformed.
+ * @param options - Optional parameters for the transformation.
+ * @returns A new Float64Array containing the transformed data.
+ */
+export function simpleNormInv(
+  data: NumberArray,
+  options: SimpleNormInvOptions = {},
+): Float64Array {
+  const { magnitudeMode = false } = options;
+
+  const result = new Float64Array(data.length);
+  if (magnitudeMode) {
+    for (let i = 0; i < result.length; i++) {
+      result[i] = -Math.sqrt(-2 * Math.log(1 - data[i]));
+    }
+  } else {
+    for (let i = 0; i < result.length; i++) {
+      result[i] = -1 * Math.SQRT2 * erfcinv(2 * data[i]);
+    }
+  }
+  return result;
+}
+
+/**
+ * Convenience wrapper for single number processing by simpleNormInv function.
+ * @param data - The number to be normalized.
+ * @param options - The options for the normalization process.
+ * @returns The normalized number.
+ */
+export function simpleNormInvNumber(
+  data: number,
+  options: SimpleNormInvOptions = {},
+): number {
+  return simpleNormInv([data], options)[0];
+}