fix: add generalizedLorentzian shape (#75)

src/index.ts

@@ -1,6 +1,7 @@
 export * from './shapes/1d/gaussian/Gaussian';
 export * from './shapes/1d/lorentzian/Lorentzian';
 export * from './shapes/1d/pseudoVoigt/PseudoVoigt';
+export * from './shapes/1d/generalizedLorentzian/GeneralizedLorentzian';
 export * from './shapes/2d/gaussian2D/Gaussian2D';
 
 export * from './shapes/1d/getShape1D';
@@ -11,6 +12,7 @@ export type {
   PseudoVoigtShape1D,
   GaussianShape1D,
   LorentzianShape1D,
+  GeneralizedLorentzianShape1D,
 } from './shapes/1d/Shape1D';
 export type { Shape2D, GaussianShape2D } from './shapes/2d/Shape2D';
 export type { Shape1DClass } from './shapes/1d/Shape1DClass';

src/shapes/1d/Shape1D.ts

@@ -1,4 +1,5 @@
 import type { GaussianClassOptions } from './gaussian/Gaussian';
+import { GeneralizedLorentzianClassOptions } from './generalizedLorentzian/GeneralizedLorentzian';
 import type { LorentzianClassOptions } from './lorentzian/Lorentzian';
 import type { PseudoVoigtClassOptions } from './pseudoVoigt/PseudoVoigt';
 
@@ -17,4 +18,13 @@ export interface PseudoVoigtShape1D extends PseudoVoigtClassOptions {
   kind: 'pseudoVoigt';
 }
 
-export type Shape1D = GaussianShape1D | LorentzianShape1D | PseudoVoigtShape1D;
+export interface GeneralizedLorentzianShape1D
+  extends GeneralizedLorentzianClassOptions {
+  kind: 'generalizedLorentzian';
+}
+
+export type Shape1D =
+  | GaussianShape1D
+  | LorentzianShape1D
+  | PseudoVoigtShape1D
+  | GeneralizedLorentzianShape1D;

src/shapes/1d/Shape1DClass.ts

@@ -2,7 +2,7 @@ import { DoubleArray } from 'cheminfo-types';
 
 import { GetData1DOptions } from './GetData1DOptions';
 
-export type Parameter = 'fwhm' | 'mu';
+export type Parameter = 'fwhm' | 'mu' | 'gamma';
 
 export interface Shape1DClass {
   fwhm: number;

src/shapes/1d/Shape1DInstance.ts

@@ -1,5 +1,10 @@
 import { Gaussian } from './gaussian/Gaussian';
+import { GeneralizedLorentzian } from './generalizedLorentzian/GeneralizedLorentzian';
 import { Lorentzian } from './lorentzian/Lorentzian';
 import { PseudoVoigt } from './pseudoVoigt/PseudoVoigt';
 
-export type Shape1DInstance = Gaussian | Lorentzian | PseudoVoigt;
+export type Shape1DInstance =
+  | Gaussian
+  | Lorentzian
+  | PseudoVoigt
+  | GeneralizedLorentzian;

src/shapes/1d/generalizedLorentzian/GeneralizedLorentzian.ts

@@ -0,0 +1,164 @@
+import { ROOT_THREE } from '../../../util/constants';
+import type { GetData1DOptions } from '../GetData1DOptions';
+import type { Parameter, Shape1DClass } from '../Shape1DClass';
+
+export interface GeneralizedLorentzianClassOptions {
+  /**
+   * Full width at half maximum.
+   * @default 500
+   */
+  fwhm?: number;
+  /**
+   * kurtosis parameter of the shape, between -1 to 2
+   * @default 1
+   */
+  gamma?: number;
+}
+
+interface GetGeneralizedLorentzianAreaOptions {
+  /**
+   * The maximum intensity value of the shape
+   * @default 1
+   */
+  height?: number;
+  /**
+   * Full width at half maximum.
+   * @default 500
+   */
+  fwhm?: number;
+  gamma?: number;
+  x?: number;
+}
+
+export class GeneralizedLorentzian implements Shape1DClass {
+  /**
+   * Full width at half maximum.
+   * @default 500
+   */
+  public fwhm: number;
+  /**
+   * kurtosis parameter of the shape, between -1 to 2
+   * @default 1
+   */
+  public gamma: number;
+
+  public constructor(options: GeneralizedLorentzianClassOptions = {}) {
+    const { fwhm = 500, gamma = 0.5 } = options;
+
+    this.fwhm = fwhm;
+    this.gamma = gamma;
+  }
+
+  public fwhmToWidth(fwhm = this.fwhm) {
+    return generalizedLorentzianFwhmToWidth(fwhm);
+  }
+
+  public widthToFWHM(width: number) {
+    return generalizedLorentzianWidthToFWHM(width);
+  }
+
+  public fct(x: number) {
+    return generalizedLorentzianFct(x, this.fwhm, this.gamma);
+  }
+
+  public getArea(height = 1) {
+    return getGeneralizedLorentzianArea({
+      fwhm: this.fwhm,
+      height,
+      gamma: this.gamma,
+    });
+  }
+
+  public getFactor(area?: number) {
+    return getGeneralizedLorentzianFactor(area);
+  }
+
+  public getData(options: GetData1DOptions = {}) {
+    return getGeneralizedLorentzianData(this, options);
+  }
+
+  public calculateHeight(area = 1) {
+    const { gamma, fwhm } = this;
+    return calculateGeneralizedLorentzianHeight({ fwhm, area, gamma });
+  }
+
+  public getParameters(): Parameter[] {
+    return ['fwhm', 'gamma'];
+  }
+}
+
+export const calculateGeneralizedLorentzianHeight = ({
+  fwhm = 1,
+  gamma = 1,
+  area = 1,
+}) => {
+  return area / fwhm / (3.14159 - 0.420894 * gamma) / 2;
+};
+
+/**
+ * expression of integral generated by Mathematica of the function
+ */
+export const getGeneralizedLorentzianArea = (
+  options: GetGeneralizedLorentzianAreaOptions,
+) => {
+  const { fwhm = 500, height = 1, gamma = 1 } = options;
+  return 2 * height * fwhm * (3.14159 - 0.420894 * gamma);
+};
+
+export const generalizedLorentzianFct = (
+  x: number,
+  fwhm: number,
+  gamma: number,
+) => {
+  const u = (x / 2 / fwhm) ** 2;
+  return (1 - gamma) / (1 + u) + (gamma * (1 + u / 2)) / (1 + u + u ** 2);
+};
+
+export const generalizedLorentzianWidthToFWHM = (width: number) => {
+  return width * ROOT_THREE;
+};
+
+export const generalizedLorentzianFwhmToWidth = (fwhm: number) => {
+  return fwhm / ROOT_THREE;
+};
+
+export const getGeneralizedLorentzianFactor = (area = 0.9999) => {
+  if (area >= 1) {
+    throw new Error('area should be (0 - 1)');
+  }
+  const halfResidual = (1 - area) * 0.5;
+  const quantileFunction = (p: number) => Math.tan(Math.PI * (p - 0.5));
+  return (
+    (quantileFunction(1 - halfResidual) - quantileFunction(halfResidual)) / 2
+  );
+};
+
+export type GetGeneralizedLorentzianData = GetData1DOptions & {
+  gamma?: number;
+};
+
+export const getGeneralizedLorentzianData = (
+  shape: GeneralizedLorentzianClassOptions = {},
+  options: GetGeneralizedLorentzianData = {},
+) => {
+  let { fwhm = 500, gamma = 1 } = shape;
+  let {
+    length,
+    factor = getGeneralizedLorentzianFactor(),
+    height = 1,
+  } = options;
+
+  if (!length) {
+    length = Math.min(Math.ceil(fwhm * factor), Math.pow(2, 25) - 1);
+    if (length % 2 === 0) length++;
+  }
+
+  const center = (length - 1) / 2;
+  const data = new Float64Array(length);
+  for (let i = 0; i <= center; i++) {
+    data[i] = generalizedLorentzianFct(i - center, fwhm, gamma) * height;
+    data[length - 1 - i] = data[i];
+  }
+
+  return data;
+};

src/shapes/1d/generalizedLorentzian/__tests__/GeneralizedLorentzian.shape.test.ts

@@ -0,0 +1,23 @@
+import { GeneralizedLorentzian } from '../GeneralizedLorentzian';
+
+describe('lorentzian', () => {
+  it('default factor area', () => {
+    const lorentzian = new GeneralizedLorentzian({ fwhm: 1.5, gamma: -1 });
+    const data = lorentzian.getData({
+      height: lorentzian.calculateHeight(),
+    });
+
+    const area = data.reduce((a, b) => a + b, 0);
+    expect(area).toBeCloseTo(1);
+
+    const data2 = lorentzian.getData({
+      height: lorentzian.calculateHeight(2),
+    });
+    const area2 = data2.reduce((a, b) => a + b, 0);
+    expect(area2).toBeCloseTo(2);
+    const expectedArea = 1;
+    const height = lorentzian.calculateHeight(expectedArea);
+    const computedArea = lorentzian.getArea(height);
+    expect(computedArea).toBeCloseTo(expectedArea, 2);
+  });
+});

src/shapes/1d/getShape1D.ts

@@ -1,6 +1,7 @@
 import type { Shape1D } from './Shape1D';
 import { Shape1DInstance } from './Shape1DInstance';
 import { Gaussian } from './gaussian/Gaussian';
+import { GeneralizedLorentzian } from './generalizedLorentzian/GeneralizedLorentzian';
 import { Lorentzian } from './lorentzian/Lorentzian';
 import { PseudoVoigt } from './pseudoVoigt/PseudoVoigt';
 
@@ -17,6 +18,8 @@ export function getShape1D(shape: Shape1D): Shape1DInstance {
       return new Lorentzian(shape);
     case 'pseudoVoigt':
       return new PseudoVoigt(shape);
+    case 'generalizedLorentzian':
+      return new GeneralizedLorentzian(shape);
     default: {
       throw Error(`Unknown distribution ${kind as string}`);
     }