feat: add support for selecting peak shape algorithm (lorentzian or gaussian)

Author: hamed-musallam

app/scripts/nmr-cli/src/index.ts

@@ -25,20 +25,21 @@ Options for 'parse-spectra' command:
   -p, --path          Directory path
   -s, --capture-snapshot   Capture snapshot
 
-
+ 
 Options for 'predict' command:
-  -n, --nucleus       Predicted nucleus  "1H" or "13C"  (required)
+  -ps,--peakShape     Peak shape algorithm (default: "lorentzian") choices: ["gaussian", "lorentzian"]
+  -n, --nucleus       Predicted nucleus, choices: ["1H","13C"]    (required)
   -i, --id            Input ID (default: 1)
   -t, --type          NMR type (default: "nmr;1H;1d")
   -s, --shifts        Chemical shifts (default: "1")
       --solvent       NMR solvent (default: "Dimethylsulphoxide-D6 (DMSO-D6, C2D6SO)")
   -m, --molText       MOL text (required)
       --from          From in (ppm)
       --to            To in (ppm)
-      --nbPoints      Number of points (default: 128)
+      --nbPoints      Number of points (default: 1024)
       --lineWidth     Line width (default: 1)
       --frequency     NMR frequency (MHz) (default: 400)
-      --tolerance     Tolerance (default: 0.001)
+      --tolerance     Tolerance to group peaks with close shift (default: 0.001)
 
 
 

app/scripts/nmr-cli/src/prediction/generatePredictedSpectrumData.ts

@@ -8,23 +8,21 @@ export interface ShiftsItem {
   spheres: number
 }
 
-interface LorentzianOptions {
+type PeakShape = 'gaussian' | 'lorentzian'
+
+interface PeakShapeOptions {
   x: number
   fwhm: number
 }
 
-function lorentzian2(options: LorentzianOptions) {
-  const { x, fwhm } = options
-  return fwhm ** 2 / (4 * x ** 2 + fwhm ** 2)
-}
-
-interface GenerateSpectrumOptions {
+export interface GenerateSpectrumOptions {
   from?: number
   to?: number
   nbPoints?: number
   lineWidth?: number
   frequency?: number
   tolerance?: number
+  peakShape?: PeakShape
 }
 
 interface GroupItem {
@@ -38,7 +36,37 @@ interface Data1D {
   re: number[]
 }
 
-const getLorentzianFactor = (area = 0.9999) => {
+const GAUSSIAN_EXP_FACTOR = -4 * Math.LN2
+
+function lorentzian(options: PeakShapeOptions) {
+  const { x, fwhm } = options
+  return fwhm ** 2 / (4 * x ** 2 + fwhm ** 2)
+}
+
+function gaussian(options: PeakShapeOptions) {
+  const { x, fwhm } = options
+  return Math.exp(GAUSSIAN_EXP_FACTOR * Math.pow(x / fwhm, 2))
+}
+
+function erfinv(x: number): number {
+  let a = 0.147
+  if (x === 0) return 0
+  let ln1MinusXSqrd = Math.log(1 - x * x)
+  let lnEtcBy2Plus2 = ln1MinusXSqrd / 2 + 2 / (Math.PI * a)
+  let firstSqrt = Math.sqrt(lnEtcBy2Plus2 ** 2 - ln1MinusXSqrd / a)
+  let secondSqrt = Math.sqrt(firstSqrt - lnEtcBy2Plus2)
+  return secondSqrt * (x > 0 ? 1 : -1)
+}
+
+function getGaussianFactor(area = 0.9999) {
+  if (area >= 1) {
+    throw new Error('area should be (0 - 1)')
+  }
+
+  return Math.sqrt(2) * erfinv(area)
+}
+
+function getLorentzianFactor(area = 0.9999) {
   if (area >= 1) {
     throw new Error('area should be (0 - 1)')
   }
@@ -49,6 +77,22 @@ const getLorentzianFactor = (area = 0.9999) => {
   )
 }
 
+function peakShapeFunction(options: PeakShapeOptions, shape: PeakShape) {
+  if (shape === 'lorentzian') {
+    return lorentzian(options)
+  }
+
+  return gaussian(options)
+}
+
+function getPeakShapeFactor(area: number, shape: PeakShape) {
+  if (shape === 'lorentzian') {
+    return getLorentzianFactor(area)
+  }
+
+  return getGaussianFactor(area)
+}
+
 function groupEquivalentShifts(shifts: ShiftsItem[], tolerance = 0.001) {
   const groups: GroupItem[] = []
 
@@ -81,6 +125,7 @@ export function generatePredictedSpectrumData(
     frequency = 400,
     lineWidth = 1,
     tolerance = 0.001,
+    peakShape = 'lorentzian',
   } = options
 
   if (!shifts || shifts.length === 0) return []
@@ -107,24 +152,22 @@ export function generatePredictedSpectrumData(
   const stepSize = (to - from) / (nbPoints - 1)
   const groupedShifts = groupEquivalentShifts(acceptedShifts, tolerance)
 
-  console.log(groupedShifts)
-
-  const limit = (lineWidth * getLorentzianFactor(0.99)) / frequency
-  console.log(limit)
+  const limit = (lineWidth * getPeakShapeFactor(0.99, peakShape)) / frequency
   for (let i = 0; i < nbPoints; i++) {
     const x = from + i * stepSize
     let intensity = 0
     for (const { prediction, count } of groupedShifts) {
       if (Math.abs(x - prediction) <= limit) {
         intensity +=
-          lorentzian2({ x: x - prediction, fwhm: lineWidth / frequency }) *
-          count
+          peakShapeFunction(
+            { x: x - prediction, fwhm: lineWidth / frequency },
+            peakShape
+          ) * count
       }
     }
-    if (intensity > 0) {
-      data.x.push(x)
-      data.re.push(intensity)
-    }
+
+    data.x.push(x)
+    data.re.push(intensity)
   }
 
   return data