feat: compute the jacobian analytically

.prettierignore

@@ -1,2 +1,4 @@
+CHANGELOG.md
 coverage
-CHANGELOG.md
+dist
+lib

README.md

@@ -7,7 +7,13 @@
 
 This is a spectra fitting package to optimize the position (x), max intensity (y),
 full width at half-maximum (FWHM = width) and the ratio of gaussian contribution (mu) if it's required.
-It supports three kinds of shapes:
+
+It supports the `gaussian`, `lorentzian`, `pseudoVoigt`, `pseudoVoigtTCH`,
+`lorentzianDispersive`, `generalizedLorentzian` and `splitGaussian` shapes of
+[ml-peak-shape-generator](https://github.com/mljs/peak-shape-generator).
+The `splitGaussian` shape is asymmetric: it is parameterized by `fwhmLow` and
+`fwhmHigh` instead of `fwhm`, and both are optimized independently.
+The three most common ones are:
 
 | Name         |                                                                                                                            Equation                                                                                                                             |
 | ------------ | :-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: |
@@ -42,72 +48,65 @@ const generator = new SpectrumGenerator({
   to: 1,
 });
 
-// by default the kind of shape is gaussian;
-generator.addPeak({ x: 0.5, y: 0.2 }, { fwhm: 0.2 });
+generator.addPeak(
+  { x: 0.5, y: 0.2 },
+  { shape: { kind: 'gaussian', fwhm: 0.2 } },
+);
 generator.addPeak(
   { x: -0.5, y: 0.2 },
-  {
-    shape: {
-      kind: 'lorentzian',
-      fwhm: 0.1,
-    },
-  },
+  { shape: { kind: 'lorentzian', fwhm: 0.1 } },
 );
 
 //points to fit {x, y};
-let data = generator.getSpectrum();
-console.log(JSON.stringify({ x: Array.from(data.x), y: Array.from(data.y) }));
-//the approximate values to be optimized, It could coming from a peak picking with ml-gsd
-let peaks = [
+const data = generator.getSpectrum();
+
+//the approximate values to be optimized, it could come from a peak picking with ml-gsd
+const peaks = [
   {
     x: -0.5,
     y: 0.22,
     shape: {
-      kind: 'gaussian',
+      kind: 'pseudoVoigt',
       fwhm: 0.25,
     },
   },
   {
     x: 0.52,
     y: 0.18,
     shape: {
-      kind: 'gaussian',
+      kind: 'pseudoVoigt',
       fwhm: 0.18,
     },
   },
 ];
 
-// the function receive an array of peak with {x, y, fwhm} as a guess
-// and return a list of objects
-let fittedParams = optimize(data, peaks, { shape: { kind: 'pseudoVoigt' } });
+// the function receives an array of peaks with {x, y, shape} as a guess
+// and returns the optimized peaks
+const fittedParams = optimize(data, peaks);
 
 console.log(fittedParams);
-const result = {
-  error: 0.12361588652854476,
-  iterations: 100,
-  peaks: [
-    {
-      x: -0.5000014532421942,
-      y: 0.19995307937326137,
-      shape: {
-        kind: 'pseudoVoigt',
-        fwhm: 0.10007670374735196,
-        mu: 0.004731136777288483,
-      },
-    },
-    {
-      x: 0.5001051783652894,
-      y: 0.19960010175400406,
-      shape: {
-        kind: 'pseudoVoigt',
-        fwhm: 0.19935932346969124,
-        mu: 1,
-      },
-    },
-  ],
-};
+// {
+//   error: 3.689192965926774e-8,
+//   iterations: 100,
+//   peaks: [
+//     {
+//       x: -0.49999998364851705,
+//       y: 0.2000001018716569,
+//       shape: { kind: 'pseudoVoigt', fwhm: 0.10000007214964078, mu: 0 },
+//     },
+//     {
+//       x: 0.500000000861883,
+//       y: 0.1999983363407256,
+//       shape: { kind: 'pseudoVoigt', fwhm: 0.20000322106167928, mu: 1 },
+//     },
+//   ],
+// }
 ```
 
+A `pseudoVoigt` fit recovers `mu = 0` for the lorentzian peak and `mu = 1` for
+the gaussian one. Each peak's own `shape.kind` wins over the `shape` passed in
+the options, which only applies to peaks that do not define one.
+
 ## Linked parameters
 
 You can link one parameter across multiple peaks so they share one optimization

eslint.config.js

@@ -1,4 +1,25 @@
 import { defineConfig, globalIgnores } from 'eslint/config';
-import cheminfo from 'eslint-config-cheminfo-typescript';
+import ts from 'eslint-config-cheminfo-typescript';
+import globals from 'globals';
 
-export default defineConfig([globalIgnores(['coverage', 'lib']), cheminfo]);
+export default defineConfig(
+  globalIgnores(['coverage', 'dist', 'lib']),
+  ts,
+  {
+    // Demo page, loaded by web/index.html as a module.
+    files: ['web/**'],
+    languageOptions: {
+      globals: globals.browser,
+    },
+  },
+  {
+    // Development scripts run with node; printing their results is the point.
+    files: ['script/**'],
+    languageOptions: {
+      globals: globals.nodeBuiltin,
+    },
+    rules: {
+      'no-console': 'off',
+    },
+  },
+);

package.json

@@ -14,11 +14,12 @@
     "build": "npm run tsc && cheminfo-build --entry lib/index.js --root SpectraFitting",
     "check-types": "tsc --noEmit",
     "clean": "rimraf coverage lib",
-    "eslint": "eslint src",
-    "eslint-fix": "eslint src --fix",
+    "dev": "vite web",
+    "eslint": "eslint .",
+    "eslint-fix": "eslint . --fix",
     "prepack": "npm run tsc",
-    "prettier": "prettier --check src",
-    "prettier-write": "prettier --write src",
+    "prettier": "prettier --check .",
+    "prettier-write": "prettier --write .",
     "test": "npm run test-only && npm run check-types && npm run eslint && npm run prettier",
     "test-only": "vitest run --coverage",
     "tsc": "npm run clean && npm run tsc-build",
@@ -46,22 +47,25 @@
   "dependencies": {
     "cheminfo-types": "^1.15.0",
     "ml-direct": "^1.0.0",
-    "ml-levenberg-marquardt": "^5.0.1",
-    "ml-peak-shape-generator": "^5.1.0",
-    "ml-spectra-processing": "^14.29.0"
+    "ml-levenberg-marquardt": "^5.1.0",
+    "ml-peak-shape-generator": "^5.5.0",
+    "ml-spectra-processing": "^14.34.0"
   },
   "devDependencies": {
-    "@types/node": "^25.9.3",
-    "@vitest/coverage-v8": "^4.1.8",
+    "@types/node": "^26.4.0",
+    "@vitest/coverage-v8": "^4.1.11",
     "@zakodium/tsconfig": "^1.0.5",
-    "cheminfo-build": "^1.3.2",
-    "eslint": "^9.28.0",
+    "chart.js": "^4.5.1",
+    "cheminfo-build": "^1.4.0",
+    "eslint": "^9.39.5",
     "eslint-config-cheminfo-typescript": "^22.1.0",
+    "globals": "^17.11.0",
     "jest-matcher-deep-close-to": "^3.0.2",
-    "prettier": "^3.8.4",
+    "prettier": "^3.9.6",
     "rimraf": "^6.1.3",
     "spectrum-generator": "^8.2.1",
     "typescript": "^6.0.3",
-    "vitest": "^4.1.8"
+    "vite": "^8.2.2",
+    "vitest": "^4.1.11"
   }
 }

script/benchmark.js

@@ -22,7 +22,9 @@ const data = generateSpectrum(peaks, {
 });
 
 let guess = structuredClone(peaks);
-guess.forEach((peak) => (peak.x += Math.random() / 10));
+for (const peak of guess) {
+  peak.x += Math.random() / 10;
+}
 
 let result = optimize(data, guess, {
   optimization: { options: { maxIterations: 10 } },

script/jacobianBenchmark.js

@@ -0,0 +1,102 @@
+// Compare the cost and the accuracy of the Levenberg-Marquardt fit.
+// Run with: node script/jacobianBenchmark.js
+
+import { generateSpectrum } from 'spectrum-generator';
+
+import { optimize } from '../src/index.ts';
+
+const scenarios = [
+  {
+    name: '1 gaussian, 101 points',
+    nbPoints: 101,
+    from: -1,
+    to: 1,
+    kind: 'gaussian',
+    peaks: [{ x: 0, y: 1, fwhm: 0.3 }],
+    repeats: 200,
+  },
+  {
+    name: '5 gaussians, 1024 points',
+    nbPoints: 1024,
+    from: -1,
+    to: 1,
+    kind: 'gaussian',
+    peaks: [
+      { x: -0.6, y: 1, fwhm: 0.08 },
+      { x: -0.3, y: 0.6, fwhm: 0.06 },
+      { x: 0, y: 0.9, fwhm: 0.1 },
+      { x: 0.35, y: 0.5, fwhm: 0.07 },
+      { x: 0.65, y: 0.8, fwhm: 0.09 },
+    ],
+    repeats: 20,
+  },
+  {
+    name: '10 pseudoVoigt, 3072 points',
+    nbPoints: 3072,
+    from: -1,
+    to: 1,
+    kind: 'pseudoVoigt',
+    peaks: Array.from({ length: 10 }, (_, i) => ({
+      x: -0.8 + i * 0.18,
+      y: 0.5 + (i % 3) * 0.2,
+      fwhm: 0.05 + (i % 4) * 0.01,
+    })),
+    repeats: 5,
+  },
+];
+
+function buildShape(kind, fwhm) {
+  return kind === 'pseudoVoigt' ? { kind, fwhm, mu: 0.5 } : { kind, fwhm };
+}
+
+for (const scenario of scenarios) {
+  const { name, nbPoints, from, to, kind, peaks, repeats } = scenario;
+
+  const truePeaks = peaks.map((peak) => ({
+    x: peak.x,
+    y: peak.y,
+    shape: buildShape(kind, peak.fwhm),
+  }));
+
+  const data = generateSpectrum(truePeaks, {
+    generator: { from, to, nbPoints, shape: { kind: 'gaussian' } },
+  });
+
+  // A deliberately offset starting guess, so the optimizer has real work to do.
+  const guess = peaks.map((peak) => ({
+    x: peak.x + 0.01,
+    y: peak.y * 0.9,
+    shape: buildShape(kind, peak.fwhm * 1.3),
+  }));
+
+  // warm up the JIT before timing
+  for (let i = 0; i < 3; i++) optimize(data, guess);
+
+  const start = performance.now();
+  let result;
+  for (let i = 0; i < repeats; i++) {
+    result = optimize(data, guess);
+  }
+  const elapsed = (performance.now() - start) / repeats;
+
+  let worstX = 0;
+  let worstFwhm = 0;
+  for (let i = 0; i < truePeaks.length; i++) {
+    worstX = Math.max(worstX, Math.abs(result.peaks[i].x - truePeaks[i].x));
+    worstFwhm = Math.max(
+      worstFwhm,
+      Math.abs(result.peaks[i].shape.fwhm - truePeaks[i].shape.fwhm),
+    );
+  }
+
+  console.log(
+    [
+      name.padEnd(28),
+      `${elapsed.toFixed(2).padStart(8)} ms`,
+      `error ${result.error.toExponential(3)}`,
+      `iter ${String(result.iterations).padStart(3)}`,
+      `dx ${worstX.toExponential(2)}`,
+      `dFwhm ${worstFwhm.toExponential(2)}`,
+    ].join('  '),
+  );
+}

src/__tests__/inputPeaksNotMutated.test.ts

@@ -0,0 +1,24 @@
+import type { DataXY } from 'cheminfo-types';
+import { generateSpectrum } from 'spectrum-generator';
+import { expect, test } from 'vitest';
+
+import { optimize } from '../index.ts';
+
+test('optimize does not modify the peaks it receives', () => {
+  const data: DataXY = generateSpectrum(
+    [{ x: -0.5, y: 0.001, shape: { kind: 'gaussian' as const, fwhm: 0.31 } }],
+    { generator: { from: -1, to: 1, nbPoints: 101 } },
+  );
+
+  const inputPeaks = [
+    { x: -0.52, y: 0.0009, shape: { kind: 'gaussian' as const, fwhm: 0.35 } },
+  ];
+
+  const result = optimize(data, inputPeaks);
+
+  expect(inputPeaks).toStrictEqual([
+    { x: -0.52, y: 0.0009, shape: { kind: 'gaussian', fwhm: 0.35 } },
+  ]);
+  expect(result.peaks[0].shape).not.toBe(inputPeaks[0].shape);
+  expect(result.peaks[0].x).toBeCloseTo(-0.5, 3);
+});