feat: reintroduce an improved noiseLevel option

Author: lpatiny

README.md

@@ -12,6 +12,7 @@
 Preprocessing of the data involves the following parameters
 
 - `maxCriteria`: search either for maxima or minima. We will invert the data and the results if searching for a minima
+- `noiseLevel`: specifies the noise level. All the peaks bellow this value (or above in case of maxCriteria=false) are ignored. By default the noiseLevel will be set to the median + 3 x sd. This is a good value when not too many peaks are present in the spectrum.
 - `sgOptions`: Savitzky-Golay filter that is used to smooth the data for the calculation of the derivatives
 - `smoothY`: If this value is true the SG filter is not only applied during the calculation of the derivatives but also on the original data
 

package.json

@@ -63,16 +63,16 @@
     "cheminfo-build": "^1.1.11",
     "eslint": "^8.9.0",
     "eslint-config-cheminfo-typescript": "^10.3.0",
-    "eslint-plugin-jest": "^26.1.0",
+    "eslint-plugin-jest": "^26.1.1",
     "esm": "^3.2.25",
     "jest": "^27.5.1",
     "jest-matcher-deep-close-to": "^3.0.2",
-    "mf-global": "^1.4.20",
+    "mf-global": "^1.4.22",
     "ml-stat": "^1.3.3",
     "nodemon": "^2.0.15",
     "prettier": "^2.5.1",
     "rimraf": "^3.0.2",
-    "spectrum-generator": "^6.0.4",
+    "spectrum-generator": "^7.0.0",
     "ts-jest": "^27.1.3",
     "typescript": "^4.5.5",
     "xy-parser": "^4.0.1"
@@ -82,6 +82,6 @@
     "ml-peak-shape-generator": "^4.0.3",
     "ml-savitzky-golay-generalized": "4.0.0",
     "ml-spectra-fitting": "^3.0.4",
-    "ml-spectra-processing": "^10.0.0"
+    "ml-spectra-processing": "^10.1.2"
   }
 }

src/__tests__/gaussian-noise.test.ts

@@ -0,0 +1,70 @@
+import type { DataXY } from 'cheminfo-types';
+import { toBeDeepCloseTo, toMatchCloseTo } from 'jest-matcher-deep-close-to';
+
+import { gsd } from '../gsd';
+
+expect.extend({ toBeDeepCloseTo, toMatchCloseTo });
+
+// eslint-disable-next-line @typescript-eslint/no-var-requires
+const { generateSpectrum } = require('spectrum-generator');
+
+describe('smooth:false option', () => {
+  const peaks = [
+    { x: -0.5, y: 1, width: 0.05 },
+    { x: 0.5, y: 1, width: 0.05 },
+  ];
+
+  const data: DataXY = generateSpectrum(peaks, {
+    generator: {
+      from: -1,
+      to: 1,
+      nbPoints: 101,
+    },
+    peaks: {
+      factor: 6,
+    },
+    noise: {
+      percent: 5,
+    },
+  });
+
+  it('positive maxima peaks', () => {
+    let peakList = gsd(data);
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: 1.131 },
+      { x: 0.5, y: 1.05 },
+    ]);
+  });
+
+  it('negative maxima peaks', () => {
+    let peakList = gsd({ x: data.x, y: data.y.map((value) => value - 2) }, {});
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: -0.868 },
+      { x: 0.5, y: -0.95 },
+    ]);
+  });
+
+  it('Negative peaks', () => {
+    // we check negative peaks
+    let peakList = gsd(
+      { x: data.x, y: data.y.map((value) => -value) },
+      { maxCriteria: false },
+    );
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: -1.131 },
+      { x: 0.5, y: -1.05 },
+    ]);
+  });
+
+  it('minima peaks', () => {
+    // we check negative peaks
+    let peakList = gsd(
+      { x: data.x, y: data.y.map((value) => 1 - value) },
+      { maxCriteria: false },
+    );
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: -0.131 },
+      { x: 0.5, y: -0.05 },
+    ]);
+  });
+});

src/__tests__/gaussian-noiseParam.test.ts

@@ -0,0 +1,87 @@
+import type { DataXY } from 'cheminfo-types';
+import { toBeDeepCloseTo, toMatchCloseTo } from 'jest-matcher-deep-close-to';
+
+import { gsd } from '../gsd';
+
+expect.extend({ toBeDeepCloseTo, toMatchCloseTo });
+
+// eslint-disable-next-line @typescript-eslint/no-var-requires
+const { generateSpectrum } = require('spectrum-generator');
+
+describe('smooth:false option', () => {
+  const peaks = [
+    { x: -0.5, y: 1, width: 0.05 },
+    { x: 0.5, y: 1, width: 0.05 },
+  ];
+
+  const data: DataXY = generateSpectrum(peaks, {
+    generator: {
+      from: -1,
+      to: 1,
+      nbPoints: 101,
+    },
+    peaks: {
+      factor: 6,
+    },
+    noise: {
+      percent: 5,
+    },
+  });
+
+  it('positive maxima peaks but noiseLevel over the peaks', () => {
+    let peakList = gsd(data, { noiseLevel: 1.2 });
+    expect(peakList).toStrictEqual([]);
+  });
+
+  it('positive maxima peaks', () => {
+    let peakList = gsd(data, { noiseLevel: 0.5 });
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: 1.131 },
+      { x: 0.5, y: 1.05 },
+    ]);
+  });
+
+  it('negative maxima peaks', () => {
+    let peakList = gsd(
+      { x: data.x, y: data.y.map((value) => value - 2) },
+      { noiseLevel: -1.5 },
+    );
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: -0.868 },
+      { x: 0.5, y: -0.95 },
+    ]);
+  });
+
+  it('Negative peaks', () => {
+    // we check negative peaks
+    let peakList = gsd(
+      { x: data.x, y: data.y.map((value) => -value) },
+      { maxCriteria: false, noiseLevel: -0.5 },
+    );
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: -1.131 },
+      { x: 0.5, y: -1.05 },
+    ]);
+  });
+
+  it('Negative peaks with noiseLevel too low', () => {
+    // we check negative peaks
+    let peakList = gsd(
+      { x: data.x, y: data.y.map((value) => -value) },
+      { maxCriteria: false, noiseLevel: -1.2 },
+    );
+    expect(peakList).toStrictEqual([]);
+  });
+
+  it('minima peaks', () => {
+    // we check negative peaks
+    let peakList = gsd(
+      { x: data.x, y: data.y.map((value) => 1 - value) },
+      { maxCriteria: false, noiseLevel: 0.5 },
+    );
+    expect(peakList).toMatchCloseTo([
+      { x: -0.5, y: -0.131 },
+      { x: 0.5, y: -0.05 },
+    ]);
+  });
+});

src/__tests__/zero.test.ts

@@ -6,7 +6,7 @@ describe('Simple test cases', () => {
     let y = [];
     expect(() => {
       gsd({ x, y });
-    }).toThrow('Window size is higher than the data length 9>0');
+    }).toThrow('input must not be empty');
   });
 
   it('length = 2', () => {

src/gsd.ts

@@ -5,6 +5,7 @@ import {
   xIsMonotoneIncreasing,
   xMinValue,
   xMaxValue,
+  xNoiseStandardDeviation,
 } from 'ml-spectra-processing';
 
 import { GSDPeak } from './GSDPeak';
@@ -30,11 +31,14 @@ export interface GSDOptions {
    */
   maxCriteria?: boolean;
   /**
-   * Threshold to determine if a given peak should be considered as a noise
+   * Peak under the noiseLevel (or over in case of maxCriteria=false) are ignored.
+   */
+  noiseLevel?: number;
+  /**
+   * Minimal height of small peaks based on the ratio between min and max
    * @default 0.00025
    */
   minMaxRatio?: number;
-
   /**
    * Use a quadratic optimizations with the peak and its 3 closest neighbors
    * @default false
@@ -56,6 +60,7 @@ export function gsd(data: DataXY, options: GSDOptions = {}): GSDPeak[] {
       windowSize: 9,
       polynomial: 3,
     },
+    noiseLevel,
     smoothY = false,
     maxCriteria = true,
     minMaxRatio = 0.00025,
@@ -72,11 +77,35 @@ export function gsd(data: DataXY, options: GSDOptions = {}): GSDPeak[] {
   // we can assume it to be equally spaced variable
   let equallySpaced = xIsEquallySpaced(x);
 
+  if (noiseLevel === undefined) {
+    if (equallySpaced) {
+      const noiseInfo = xNoiseStandardDeviation(y);
+      if (maxCriteria) {
+        noiseLevel = noiseInfo.median + 1.5 * noiseInfo.sd;
+      } else {
+        noiseLevel = -noiseInfo.median + 1.5 * noiseInfo.sd;
+      }
+    } else {
+      noiseLevel = 0;
+    }
+  } else {
+    if (maxCriteria === false) {
+      noiseLevel *= -1;
+    }
+  }
+
   if (maxCriteria === false) {
     for (let i = 0; i < y.length; i++) {
       y[i] *= -1;
     }
   }
+  if (noiseLevel) {
+    for (let i = 0; i < y.length; i++) {
+      if (y[i] < noiseLevel) {
+        y[i] = noiseLevel;
+      }
+    }
+  }
 
   let yData = y;
   let dY, ddY;