feat: add asymmetric split gaussian shape

README.md

@@ -16,6 +16,7 @@ The current supported kinds of shapes:
 | Lorentzian Dispersive       |                                                                       <img src="https://tex.cheminfo.org/?tex=y%5Ccdot%5Cfrac%7B2%5Comega%5Cdelta%7D%7B4%5Cdelta%5E2%20%2B%20%5Comega%5E2%7D"/>                                                                        |
 | Generalized Lorentzian      | <img src="https://tex.cheminfo.org/?tex=y%5Ccdot%5Cleft%5B%281-%5Cgamma%29%5Ccdot%5Cfrac%7B1%7D%7B1%2Bu%7D%20%2B%20%5Cgamma%5Ccdot%5Cfrac%7B1%2B%5Cfrac%7Bu%7D%7B2%7D%7D%7B1%2Bu%2Bu%5E2%7D%5Cright%5D"/> |
 | Pseudo Voigt                | <img src="https://tex.cheminfo.org/v1/?tex=y%20%5Ccdot%5Cleft%5Bmu%5Ccdot%20exp%5Cleft%5B-%5Cfrac%7B1%7D%7B2%7D%5Cfrac%7B%5Cdelta%7D%7B%5Csigma%5E2%7D%5Cright%5D%20%2B%20%5Cleft(1%20-%20mu%5Cright)%5Ccdot%5Cfrac%7B%5Comega%5E2%7D%7B4%5Cdelta%20%2B%20%5Comega%5E2%7D%5Cright%5D"/> |
+| Split Gaussian (asymmetric) | Two gaussian halves sharing the apex: the lower-x half (`t ≤ x`) uses `fwhmLow`, the higher-x half (`t > x`) uses `fwhmHigh`. |
 
 where
 

node_modules

@@ -0,0 +1 @@
+/Users/lpatiny/git/mljs/peak-shape-generator/node_modules

src/index.ts

@@ -4,6 +4,7 @@ export * from './shapes/1d/lorentzianDispersive/LorentzianDispersive.ts';
 export * from './shapes/1d/pseudoVoigt/PseudoVoigt.ts';
 export * from './shapes/1d/pseudoVoigtTCH/PseudoVoigtTCH.ts';
 export * from './shapes/1d/generalizedLorentzian/GeneralizedLorentzian.ts';
+export * from './shapes/1d/splitGaussian/SplitGaussian.ts';
 export * from './shapes/2d/gaussian2D/Gaussian2D.ts';
 
 export * from './shapes/1d/getShape1D.ts';
@@ -17,6 +18,7 @@ export type {
   PseudoVoigtShape1D,
   PseudoVoigtTCHShape1D,
   Shape1D,
+  SplitGaussianShape1D,
 } from './shapes/1d/Shape1D.ts';
 export type { GaussianShape2D, Shape2D } from './shapes/2d/Shape2D.ts';
 export type {

src/shapes/1d/Shape1D.ts

@@ -3,6 +3,7 @@ import type { GeneralizedLorentzianClassOptions } from './generalizedLorentzian/
 import type { LorentzianClassOptions } from './lorentzian/Lorentzian.ts';
 import type { PseudoVoigtClassOptions } from './pseudoVoigt/PseudoVoigt.ts';
 import type { PseudoVoigtTCHClassOptions } from './pseudoVoigtTCH/PseudoVoigtTCH.ts';
+import type { SplitGaussianClassOptions } from './splitGaussian/SplitGaussian.ts';
 
 /**
  * kind of shape
@@ -31,10 +32,15 @@ export interface PseudoVoigtTCHShape1D extends PseudoVoigtTCHClassOptions {
   kind: 'pseudoVoigtTCH';
 }
 
+export interface SplitGaussianShape1D extends SplitGaussianClassOptions {
+  kind: 'splitGaussian';
+}
+
 export type Shape1D =
   | GaussianShape1D
   | LorentzianShape1D
   | PseudoVoigtShape1D
   | LorentzianDispersiveShape1D
   | GeneralizedLorentzianShape1D
-  | PseudoVoigtTCHShape1D;
+  | PseudoVoigtTCHShape1D
+  | SplitGaussianShape1D;

src/shapes/1d/Shape1DClass.ts

@@ -2,7 +2,8 @@ import type { DoubleArray } from 'cheminfo-types';
 
 import type { GetData1DOptions } from './GetData1DOptions.ts';
 
-export type Parameter = 'fwhm' | 'mu' | 'gamma' | 'fwhmG' | 'fwhmL';
+export type Parameter =
+  'fwhm' | 'mu' | 'gamma' | 'fwhmG' | 'fwhmL' | 'fwhmLow' | 'fwhmHigh';
 
 /**
  * The exact tuple of parameters a shape exposes. Constraining `T` to

src/shapes/1d/Shape1DInstance.ts

@@ -4,6 +4,7 @@ import type { Lorentzian } from './lorentzian/Lorentzian.ts';
 import type { LorentzianDispersive } from './lorentzianDispersive/LorentzianDispersive.ts';
 import type { PseudoVoigt } from './pseudoVoigt/PseudoVoigt.ts';
 import type { PseudoVoigtTCH } from './pseudoVoigtTCH/PseudoVoigtTCH.ts';
+import type { SplitGaussian } from './splitGaussian/SplitGaussian.ts';
 
 export interface Shape1DKindInstanceMap {
   gaussian: Gaussian;
@@ -12,6 +13,7 @@ export interface Shape1DKindInstanceMap {
   pseudoVoigtTCH: PseudoVoigtTCH;
   lorentzianDispersive: LorentzianDispersive;
   generalizedLorentzian: GeneralizedLorentzian;
+  splitGaussian: SplitGaussian;
 }
 
 export type Shape1DInstance<

src/shapes/1d/__tests__/getShape1D.test.ts

@@ -4,6 +4,7 @@ import type { Shape1D } from '../Shape1D.ts';
 import type { Shape1DInstance } from '../Shape1DInstance.ts';
 import { Gaussian } from '../gaussian/Gaussian.ts';
 import { getShape1D } from '../getShape1D.ts';
+import { SplitGaussian } from '../splitGaussian/SplitGaussian.ts';
 
 test('returns a Gaussian instance for gaussian input', () => {
   const shape = getShape1D({ kind: 'gaussian', fwhm: 10 });
@@ -13,6 +14,19 @@ test('returns a Gaussian instance for gaussian input', () => {
   expect(shape).toBeInstanceOf(Gaussian);
 });
 
+test('returns a SplitGaussian instance for splitGaussian input', () => {
+  const shape = getShape1D({
+    kind: 'splitGaussian',
+    fwhmLow: 10,
+    fwhmHigh: 30,
+  });
+
+  expectTypeOf(shape).toEqualTypeOf<SplitGaussian>();
+
+  expect(shape).toBeInstanceOf(SplitGaussian);
+  expect(shape.fwhm).toBe(20);
+});
+
 test('returns the broad union for widened Shape1D input', () => {
   const useGaussian = true as boolean;
   const shape: Shape1D = useGaussian

src/shapes/1d/getShape1D.ts

@@ -6,6 +6,7 @@ import { Lorentzian } from './lorentzian/Lorentzian.ts';
 import { LorentzianDispersive } from './lorentzianDispersive/LorentzianDispersive.ts';
 import { PseudoVoigt } from './pseudoVoigt/PseudoVoigt.ts';
 import { PseudoVoigtTCH } from './pseudoVoigtTCH/PseudoVoigtTCH.ts';
+import { SplitGaussian } from './splitGaussian/SplitGaussian.ts';
 
 /**
  * Generate an instance of a specific kind of shape.
@@ -30,6 +31,8 @@ export function getShape1D(shape: Shape1D): Shape1DInstance {
       return new LorentzianDispersive(shape);
     case 'generalizedLorentzian':
       return new GeneralizedLorentzian(shape);
+    case 'splitGaussian':
+      return new SplitGaussian(shape);
     default:
       throw new Error(`Unknown distribution ${kind as string}`);
   }

src/shapes/1d/splitGaussian/SplitGaussian.ts

@@ -0,0 +1,249 @@
+import { ROOT_PI_OVER_LN2 } from '../../../util/constants.ts';
+import type { GetData1DOptions } from '../GetData1DOptions.ts';
+import type {
+  ParameterTuple,
+  Shape1DClass,
+  Shape1DDerivative,
+} from '../Shape1DClass.ts';
+import {
+  gaussianDerivative,
+  gaussianFct,
+  gaussianFwhmToWidth,
+  gaussianWidthToFWHM,
+  getGaussianFactor,
+} from '../gaussian/Gaussian.ts';
+
+export interface SplitGaussianClassOptions {
+  /**
+   * Full width at half maximum of the lower-x half (x <= 0).
+   * @default 500
+   */
+  fwhmLow?: number;
+  /**
+   * Full width at half maximum of the higher-x half (x > 0).
+   * @default 500
+   */
+  fwhmHigh?: number;
+}
+
+interface CalculateSplitGaussianHeightOptions {
+  /**
+   * Full width at half maximum of the lower-x half.
+   * @default 500
+   */
+  fwhmLow?: number;
+  /**
+   * Full width at half maximum of the higher-x half.
+   * @default 500
+   */
+  fwhmHigh?: number;
+  /**
+   * @default 1
+   */
+  area?: number;
+}
+
+interface GetSplitGaussianAreaOptions {
+  /**
+   * The maximum intensity value of the shape.
+   * @default 1
+   */
+  height?: number;
+  /**
+   * Full width at half maximum of the lower-x half.
+   * @default 500
+   */
+  fwhmLow?: number;
+  /**
+   * Full width at half maximum of the higher-x half.
+   * @default 500
+   */
+  fwhmHigh?: number;
+}
+
+export class SplitGaussian implements Shape1DClass {
+  /**
+   * Full width at half maximum of the lower-x half (x <= 0).
+   * @default 500
+   */
+  public fwhmLow: number;
+  /**
+   * Full width at half maximum of the higher-x half (x > 0).
+   * @default 500
+   */
+  public fwhmHigh: number;
+
+  public constructor(options: SplitGaussianClassOptions = {}) {
+    const { fwhmLow = 500, fwhmHigh = 500 } = options;
+
+    this.fwhmLow = fwhmLow;
+    this.fwhmHigh = fwhmHigh;
+  }
+
+  /**
+   * Full width at half maximum of the peak. The half-maximum crossings are at
+   * `-fwhmLow / 2` and `fwhmHigh / 2`, so the width between them is the mean of
+   * both halves.
+   * @returns the full width at half maximum.
+   */
+  public get fwhm() {
+    return (this.fwhmLow + this.fwhmHigh) / 2;
+  }
+
+  /**
+   * Convert a full width at half maximum to the width between the inflection
+   * points. For this peak's own fwhm the result is exactly `σlow + σhigh`.
+   * @param fwhm - full width at half maximum. Defaults to the peak's fwhm.
+   * @returns the width between the inflection points.
+   */
+  public fwhmToWidth(fwhm = this.fwhm) {
+    return gaussianFwhmToWidth(fwhm);
+  }
+
+  /**
+   * Convert a width between the inflection points back to a full width at half
+   * maximum. A single width does not encode the asymmetry, so it cannot recover
+   * `fwhmLow` and `fwhmHigh` individually.
+   * @param width - width between the inflection points.
+   * @returns the corresponding full width at half maximum.
+   */
+  public widthToFWHM(width: number) {
+    return gaussianWidthToFWHM(width);
+  }
+
+  public fct(x: number) {
+    return splitGaussianFct(x, this.fwhmLow, this.fwhmHigh);
+  }
+
+  public getArea(
+    height = calculateSplitGaussianHeight({
+      fwhmLow: this.fwhmLow,
+      fwhmHigh: this.fwhmHigh,
+    }),
+  ) {
+    return getSplitGaussianArea({
+      fwhmLow: this.fwhmLow,
+      fwhmHigh: this.fwhmHigh,
+      height,
+    });
+  }
+
+  public getFactor(area?: number) {
+    return getGaussianFactor(area);
+  }
+
+  public getData(options: GetData1DOptions = {}) {
+    return getSplitGaussianData(this, options);
+  }
+
+  public calculateHeight(area = 1) {
+    return calculateSplitGaussianHeight({
+      fwhmLow: this.fwhmLow,
+      fwhmHigh: this.fwhmHigh,
+      area,
+    });
+  }
+
+  public getParameters(): ParameterTuple<['fwhmLow', 'fwhmHigh']> {
+    return ['fwhmLow', 'fwhmHigh'];
+  }
+
+  public derivative(x: number): Shape1DDerivative {
+    const { fct, dx, dFwhmLow, dFwhmHigh } = splitGaussianDerivative(
+      x,
+      this.fwhmLow,
+      this.fwhmHigh,
+    );
+    return { fct, dx, parameters: [dFwhmLow, dFwhmHigh] };
+  }
+}
+
+/**
+ * Calculate the peak height for a given area and both half-widths.
+ * @param options - fwhmLow, fwhmHigh and area.
+ * @returns the peak height.
+ */
+export function calculateSplitGaussianHeight(
+  options: CalculateSplitGaussianHeightOptions,
+) {
+  const { fwhmLow = 500, fwhmHigh = 500, area = 1 } = options;
+  return (4 * area) / ROOT_PI_OVER_LN2 / (fwhmLow + fwhmHigh);
+}
+
+/**
+ * Evaluate the split (asymmetric) gaussian function centered at x=0.
+ * The lower-x half (x <= 0) uses `fwhmLow`, the higher-x half (x > 0) uses `fwhmHigh`.
+ * @param x - position at which to evaluate.
+ * @param fwhmLow - full width at half maximum of the lower-x half.
+ * @param fwhmHigh - full width at half maximum of the higher-x half.
+ * @returns the intensity at x.
+ */
+export function splitGaussianFct(x: number, fwhmLow: number, fwhmHigh: number) {
+  return x <= 0 ? gaussianFct(x, fwhmLow) : gaussianFct(x, fwhmHigh);
+}
+
+/**
+ * Analytical value and partial derivatives of the split gaussian function centered at x=0.
+ * Each half's fwhm only affects its own side, so the off-side derivative is 0.
+ * @param x - position at which to evaluate.
+ * @param fwhmLow - full width at half maximum of the lower-x half.
+ * @param fwhmHigh - full width at half maximum of the higher-x half.
+ * @returns the value `fct` and its partial derivatives with respect to `x` (`dx`), `fwhmLow` (`dFwhmLow`) and `fwhmHigh` (`dFwhmHigh`).
+ */
+export function splitGaussianDerivative(
+  x: number,
+  fwhmLow: number,
+  fwhmHigh: number,
+) {
+  if (x <= 0) {
+    const { fct, dx, dFwhm } = gaussianDerivative(x, fwhmLow);
+    return { fct, dx, dFwhmLow: dFwhm, dFwhmHigh: 0 };
+  }
+  const { fct, dx, dFwhm } = gaussianDerivative(x, fwhmHigh);
+  return { fct, dx, dFwhmLow: 0, dFwhmHigh: dFwhm };
+}
+
+/**
+ * Calculate the area under a split gaussian peak.
+ * @param options - fwhmLow, fwhmHigh and height.
+ * @returns the area.
+ */
+export function getSplitGaussianArea(options: GetSplitGaussianAreaOptions) {
+  const { fwhmLow = 500, fwhmHigh = 500, height = 1 } = options;
+  return (height * ROOT_PI_OVER_LN2 * (fwhmLow + fwhmHigh)) / 4;
+}
+
+/**
+ * Generate an intensity array for a split gaussian shape.
+ * @param shape - split gaussian shape parameters (fwhmLow, fwhmHigh).
+ * @param options - sampling options (length, factor, height).
+ * @returns Float64Array of intensity values.
+ */
+export function getSplitGaussianData(
+  shape: SplitGaussianClassOptions = {},
+  options: GetData1DOptions = {},
+) {
+  const { fwhmLow = 500, fwhmHigh = 500 } = shape;
+
+  const {
+    factor = getGaussianFactor(),
+    height = calculateSplitGaussianHeight({ fwhmLow, fwhmHigh }),
+  } = options;
+  let { length } = options;
+
+  if (!length) {
+    length = Math.min(
+      Math.ceil(Math.max(fwhmLow, fwhmHigh) * factor),
+      2 ** 25 - 1,
+    );
+    if (length % 2 === 0) length++;
+  }
+
+  const center = (length - 1) / 2;
+  const data = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    data[i] = splitGaussianFct(i - center, fwhmLow, fwhmHigh) * height;
+  }
+
+  return data;
+}