diff --git a/README.md b/README.md index 0592977..b4fc094 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,7 @@ The current supported kinds of shapes: | Lorentzian Dispersive | | | Generalized Lorentzian | | | Pseudo Voigt | | +| 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 diff --git a/node_modules b/node_modules new file mode 120000 index 0000000..42d64e5 --- /dev/null +++ b/node_modules @@ -0,0 +1 @@ +/Users/lpatiny/git/mljs/peak-shape-generator/node_modules \ No newline at end of file diff --git a/src/index.ts b/src/index.ts index edaa9f8..584055d 100644 --- a/src/index.ts +++ b/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 { diff --git a/src/shapes/1d/Shape1D.ts b/src/shapes/1d/Shape1D.ts index 79135a0..bc20efd 100644 --- a/src/shapes/1d/Shape1D.ts +++ b/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; diff --git a/src/shapes/1d/Shape1DClass.ts b/src/shapes/1d/Shape1DClass.ts index 71a74a5..bf410de 100644 --- a/src/shapes/1d/Shape1DClass.ts +++ b/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 diff --git a/src/shapes/1d/Shape1DInstance.ts b/src/shapes/1d/Shape1DInstance.ts index bddaebb..b1a396d 100644 --- a/src/shapes/1d/Shape1DInstance.ts +++ b/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< diff --git a/src/shapes/1d/__tests__/getShape1D.test.ts b/src/shapes/1d/__tests__/getShape1D.test.ts index 11001f6..fd10bbb 100644 --- a/src/shapes/1d/__tests__/getShape1D.test.ts +++ b/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(); + + 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 diff --git a/src/shapes/1d/getShape1D.ts b/src/shapes/1d/getShape1D.ts index d1bd013..2100348 100644 --- a/src/shapes/1d/getShape1D.ts +++ b/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}`); } diff --git a/src/shapes/1d/splitGaussian/SplitGaussian.ts b/src/shapes/1d/splitGaussian/SplitGaussian.ts new file mode 100644 index 0000000..dd5221f --- /dev/null +++ b/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; +} diff --git a/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.derivative.test.ts b/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.derivative.test.ts new file mode 100644 index 0000000..29982f4 --- /dev/null +++ b/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.derivative.test.ts @@ -0,0 +1,81 @@ +import { expect, test } from 'vitest'; + +import { gaussianDerivative } from '../../gaussian/Gaussian.ts'; +import { + SplitGaussian, + splitGaussianDerivative, + splitGaussianFct, +} from '../SplitGaussian.ts'; + +const fwhmLow = 0.2; +const fwhmHigh = 0.4; +const h = 1e-6; + +test('splitGaussianDerivative matches numerical derivatives on each side', () => { + for (const x of [-0.4, -0.1, 0.05, 0.25]) { + const { fct, dx, dFwhmLow, dFwhmHigh } = splitGaussianDerivative( + x, + fwhmLow, + fwhmHigh, + ); + + expect(fct).toBeCloseTo(splitGaussianFct(x, fwhmLow, fwhmHigh), 12); + + const numericalDx = + (splitGaussianFct(x + h, fwhmLow, fwhmHigh) - + splitGaussianFct(x - h, fwhmLow, fwhmHigh)) / + (2 * h); + + expect(dx).toBeCloseTo(numericalDx, 6); + + const numericalDFwhmLow = + (splitGaussianFct(x, fwhmLow + h, fwhmHigh) - + splitGaussianFct(x, fwhmLow - h, fwhmHigh)) / + (2 * h); + + expect(dFwhmLow).toBeCloseTo(numericalDFwhmLow, 6); + + const numericalDFwhmHigh = + (splitGaussianFct(x, fwhmLow, fwhmHigh + h) - + splitGaussianFct(x, fwhmLow, fwhmHigh - h)) / + (2 * h); + + expect(dFwhmHigh).toBeCloseTo(numericalDFwhmHigh, 6); + } +}); + +test('each side is the gaussian of its own half', () => { + const lower = splitGaussianDerivative(-0.1, fwhmLow, fwhmHigh); + const higher = splitGaussianDerivative(0.1, fwhmLow, fwhmHigh); + + expect(lower.dFwhmLow).toBe(gaussianDerivative(-0.1, fwhmLow).dFwhm); + expect(lower.dFwhmHigh).toBe(0); + expect(higher.dFwhmHigh).toBe(gaussianDerivative(0.1, fwhmHigh).dFwhm); + expect(higher.dFwhmLow).toBe(0); +}); + +test('the apex is a smooth maximum', () => { + const { fct, dx, dFwhmLow, dFwhmHigh } = splitGaussianDerivative( + 0, + fwhmLow, + fwhmHigh, + ); + + expect(fct).toBe(1); + expect(dx).toBeCloseTo(0, 12); + expect(dFwhmLow).toBeCloseTo(0, 12); + expect(dFwhmHigh).toBeCloseTo(0, 12); +}); + +test('SplitGaussian.derivative returns parameters in getParameters() order', () => { + const shape = new SplitGaussian({ fwhmLow, fwhmHigh }); + const { fct, dx, parameters } = shape.derivative(0.05); + const expected = splitGaussianDerivative(0.05, fwhmLow, fwhmHigh); + + expect(shape.getParameters()).toStrictEqual(['fwhmLow', 'fwhmHigh']); + expect(fct).toBeCloseTo(expected.fct, 12); + expect(dx).toBeCloseTo(expected.dx, 12); + expect(parameters).toHaveLength(2); + expect(parameters[0]).toBeCloseTo(expected.dFwhmLow, 12); + expect(parameters[1]).toBeCloseTo(expected.dFwhmHigh, 12); +}); diff --git a/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.fct.test.ts b/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.fct.test.ts new file mode 100644 index 0000000..d947db8 --- /dev/null +++ b/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.fct.test.ts @@ -0,0 +1,36 @@ +import { expect, test } from 'vitest'; + +import { SplitGaussian, splitGaussianFct } from '../SplitGaussian.ts'; + +test('SplitGaussian.fct is asymmetric around the apex', () => { + const shape = new SplitGaussian({ fwhmLow: 0.2, fwhmHigh: 0.4 }); + + expect(shape.fct(0)).toBeCloseTo(1); + // lower-x half reaches half-max at -fwhmLow/2 + expect(shape.fct(-0.1)).toBeCloseTo(0.5); + // higher-x half reaches half-max at +fwhmHigh/2 + expect(shape.fct(0.2)).toBeCloseTo(0.5); + // same distance, different intensity -> asymmetry + expect(shape.fct(-0.1)).toBeLessThan(shape.fct(0.1)); +}); + +test('splitGaussianFct with equal halves matches a symmetric gaussian', () => { + expect(splitGaussianFct(-0.1, 0.2, 0.2)).toBeCloseTo(0.5); + expect(splitGaussianFct(0.1, 0.2, 0.2)).toBeCloseTo(0.5); +}); + +test('fwhm is the distance between both half-maximum crossings', () => { + const shape = new SplitGaussian({ fwhmLow: 200, fwhmHigh: 600 }); + + expect(shape.fct(-100)).toBeCloseTo(0.5, 12); + expect(shape.fct(300)).toBeCloseTo(0.5, 12); + expect(shape.fwhm).toBe(400); +}); + +test('both halves default to 500', () => { + const shape = new SplitGaussian(); + + expect(shape.fwhmLow).toBe(500); + expect(shape.fwhmHigh).toBe(500); + expect(shape.fwhm).toBe(500); +}); diff --git a/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.shape.test.ts b/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.shape.test.ts new file mode 100644 index 0000000..c5b1d3d --- /dev/null +++ b/src/shapes/1d/splitGaussian/__tests__/SplitGaussian.shape.test.ts @@ -0,0 +1,84 @@ +import { expect, test } from 'vitest'; + +import { ROOT_PI_OVER_LN2 } from '../../../../util/constants.ts'; +import { getGaussianData } from '../../gaussian/Gaussian.ts'; +import { + SplitGaussian, + calculateSplitGaussianHeight, + getSplitGaussianArea, +} from '../SplitGaussian.ts'; + +test('height 1, apex at the center', () => { + const shape = new SplitGaussian({ fwhmLow: 10, fwhmHigh: 30 }); + const data = shape.getData({ height: 1 }); + + const center = (data.length - 1) / 2; + + expect(data[center]).toBe(1); + // higher-x half is wider, so it decays slower than the lower-x half + expect(data[center + 5]).toBeGreaterThan(data[center - 5]); + + const area = data.reduce((a, b) => a + b, 0); + + expect(area).toBeCloseTo((ROOT_PI_OVER_LN2 * (10 + 30)) / 4, 2); + expect(shape.getParameters()).toStrictEqual(['fwhmLow', 'fwhmHigh']); +}); + +test('normalized area is close to 1', () => { + const shape = new SplitGaussian({ fwhmLow: 40, fwhmHigh: 60 }); + const data = shape.getData(); + + const area = data.reduce((a, b) => a + b, 0); + + // the window is sized on the wider half, so it covers more than the 0.9999 factor + expect(area).toBeCloseTo(1, 5); + expect(shape.getArea()).toBeCloseTo(1, 12); +}); + +test('equal halves reduce to a symmetric gaussian area', () => { + const shape = new SplitGaussian({ fwhmLow: 50, fwhmHigh: 50 }); + + expect(shape.fwhm).toBe(50); + expect(getSplitGaussianArea({ fwhmLow: 50, fwhmHigh: 50 })).toBeCloseTo( + (ROOT_PI_OVER_LN2 * 50) / 2, + 6, + ); +}); + +test('equal halves produce the same data as a gaussian', () => { + const data = new SplitGaussian({ fwhmLow: 50, fwhmHigh: 50 }).getData(); + + expect(data).toStrictEqual(getGaussianData({ fwhm: 50 })); +}); + +test('length option is honored, apex stays at the center', () => { + const shape = new SplitGaussian({ fwhmLow: 10, fwhmHigh: 30 }); + const data = shape.getData({ length: 101, height: 1 }); + + expect(data).toHaveLength(101); + expect(data[50]).toBe(1); + expect(data[40]).toBe(shape.fct(-10)); + expect(data[60]).toBe(shape.fct(10)); +}); + +test('height calculation is consistent with the area', () => { + const shape = new SplitGaussian({ fwhmLow: 100, fwhmHigh: 300 }); + const height = shape.calculateHeight(); + const expected = calculateSplitGaussianHeight({ + fwhmLow: 100, + fwhmHigh: 300, + area: 1, + }); + + expect(height).toBeCloseTo(expected, 6); + expect( + getSplitGaussianArea({ fwhmLow: 100, fwhmHigh: 300, height }), + ).toBeCloseTo(1, 6); +}); + +test('change height should change area proportionally', () => { + const shape = new SplitGaussian({ fwhmLow: 100, fwhmHigh: 200 }); + const area = shape.getArea(1); + + expect(shape.getArea(2)).toBeCloseTo(2 * area, 4); +});