diff --git a/examples/xy-peaks-plotly.html b/examples/xy-peaks-plotly.html new file mode 100644 index 0000000..0a37f32 --- /dev/null +++ b/examples/xy-peaks-plotly.html @@ -0,0 +1,102 @@ + + + + + + XY Data Visualization with Peaks + + + + +

XY Data Visualization with Peaks

+ +
+ + + diff --git a/package.json b/package.json index 8c47f2b..aa0bf3c 100644 --- a/package.json +++ b/package.json @@ -3,7 +3,9 @@ "version": "13.0.1", "description": "Global Spectral Deconvolution", "type": "module", - "exports": "./lib/index.js", + "exports": { + ".": "./lib/index.js" + }, "files": [ "lib", "src" diff --git a/src/XIndex.ts b/src/XIndex.ts new file mode 100644 index 0000000..0397385 --- /dev/null +++ b/src/XIndex.ts @@ -0,0 +1,4 @@ +export interface XIndex { + x: number; + index: number; +} diff --git a/src/__tests__/data/broadNotLorentzianShape.json b/src/__tests__/data/broadNotLorentzianShape.json new file mode 100644 index 0000000..a2146b1 --- /dev/null +++ b/src/__tests__/data/broadNotLorentzianShape.json @@ -0,0 +1,528 @@ +{ + "x": [ + 7.4723387985080425, + 7.472644400345208, + 7.472950002182374, + 7.47325560401954, + 7.473561205856706, + 7.473866807693872, + 7.47417240953104, + 7.474478011368205, + 7.474783613205371, + 7.475089215042537, + 7.475394816879703, + 7.475700418716871, + 7.476006020554037, + 7.4763116223912025, + 7.476617224228368, + 7.476922826065534, + 7.477228427902702, + 7.477534029739868, + 7.477839631577034, + 7.4781452334142, + 7.4784508352513654, + 7.478756437088531, + 7.479062038925699, + 7.479367640762865, + 7.479673242600031, + 7.479978844437197, + 7.4802844462743625, + 7.48059004811153, + 7.480895649948696, + 7.481201251785862, + 7.481506853623028, + 7.481812455460194, + 7.482118057297361, + 7.482423659134527, + 7.482729260971693, + 7.483034862808859, + 7.483340464646025, + 7.483646066483191, + 7.483951668320358, + 7.484257270157524, + 7.48456287199469, + 7.484868473831856, + 7.485174075669022, + 7.48547967750619, + 7.4857852793433555, + 7.486090881180521, + 7.486396483017687, + 7.486702084854853, + 7.487007686692021, + 7.487313288529187, + 7.4876188903663525, + 7.487924492203518, + 7.488230094040684, + 7.488535695877852, + 7.488841297715018, + 7.489146899552184, + 7.48945250138935, + 7.4897581032265155, + 7.490063705063681, + 7.490369306900849, + 7.490674908738015, + 7.490980510575181, + 7.491286112412347, + 7.4915917142495125, + 7.49189731608668, + 7.492202917923846, + 7.492508519761012, + 7.492814121598178, + 7.493119723435344, + 7.493425325272511, + 7.493730927109677, + 7.494036528946843, + 7.494342130784009, + 7.494647732621175, + 7.494953334458341, + 7.495258936295508, + 7.495564538132674, + 7.49587013996984, + 7.496175741807006, + 7.496481343644172, + 7.49678694548134, + 7.4970925473185055, + 7.497398149155671, + 7.497703750992837, + 7.498009352830003, + 7.498314954667171, + 7.498620556504337, + 7.498926158341503, + 7.499231760178668, + 7.499537362015834, + 7.499842963853, + 7.500148565690168, + 7.500454167527334, + 7.5007597693645, + 7.5010653712016655, + 7.501370973038831, + 7.501676574875999, + 7.501982176713165, + 7.502287778550331, + 7.502593380387497, + 7.502898982224663, + 7.50320458406183, + 7.503510185898996, + 7.503815787736162, + 7.504121389573328, + 7.504426991410494, + 7.50473259324766, + 7.505038195084827, + 7.505343796921993, + 7.505649398759159, + 7.505955000596325, + 7.506260602433491, + 7.5065662042706585, + 7.506871806107824, + 7.50717740794499, + 7.507483009782156, + 7.507788611619322, + 7.50809421345649, + 7.5083998152936555, + 7.508705417130821, + 7.509011018967987, + 7.509316620805153, + 7.509622222642319, + 7.509927824479487, + 7.510233426316653, + 7.5105390281538185, + 7.510844629990984, + 7.51115023182815, + 7.511455833665318, + 7.511761435502484, + 7.51206703733965, + 7.5123726391768155, + 7.512678241013981, + 7.512983842851149, + 7.513289444688315, + 7.513595046525481, + 7.513900648362647, + 7.514206250199813, + 7.51451185203698, + 7.514817453874146, + 7.515123055711312, + 7.515428657548478, + 7.515734259385644, + 7.51603986122281, + 7.516345463059977, + 7.516651064897143, + 7.516956666734309, + 7.517262268571475, + 7.517567870408641, + 7.5178734722458085, + 7.518179074082974, + 7.51848467592014, + 7.518790277757306, + 7.519095879594472, + 7.51940148143164, + 7.5197070832688055, + 7.520012685105971, + 7.520318286943137, + 7.520623888780303, + 7.520929490617469, + 7.521235092454637, + 7.521540694291803, + 7.5218462961289685, + 7.522151897966134, + 7.5224574998033, + 7.522763101640468, + 7.523068703477634, + 7.5233743053148, + 7.5236799071519656, + 7.523985508989131, + 7.524291110826299, + 7.524596712663465, + 7.524902314500631, + 7.525207916337797, + 7.525513518174963, + 7.5258191200121285, + 7.526124721849296, + 7.526430323686462, + 7.526735925523628, + 7.527041527360794, + 7.52734712919796, + 7.527652731035127, + 7.527958332872293, + 7.528263934709459, + 7.528569536546625, + 7.528875138383791, + 7.5291807402209585, + 7.529486342058124, + 7.52979194389529, + 7.530097545732456, + 7.530403147569622, + 7.530708749406788, + 7.531014351243956, + 7.5313199530811215, + 7.531625554918287, + 7.531931156755453, + 7.532236758592619, + 7.532542360429787, + 7.532847962266953, + 7.5331535641041185, + 7.533459165941284, + 7.53376476777845, + 7.534070369615618, + 7.534375971452784, + 7.53468157328995, + 7.534987175127116, + 7.5352927769642815, + 7.535598378801447, + 7.535903980638615, + 7.536209582475781, + 7.536515184312947, + 7.536820786150113, + 7.5371263879872785, + 7.537431989824446, + 7.537737591661612, + 7.538043193498778, + 7.538348795335944, + 7.53865439717311, + 7.538959999010277, + 7.539265600847443, + 7.539571202684609, + 7.539876804521775, + 7.540182406358941, + 7.5404880081961085, + 7.540793610033274, + 7.54109921187044, + 7.541404813707606, + 7.541710415544772, + 7.542016017381938, + 7.542321619219106, + 7.5426272210562715, + 7.542932822893437, + 7.543238424730603, + 7.543544026567769, + 7.543849628404937, + 7.544155230242103, + 7.5444608320792685, + 7.544766433916434, + 7.5450720357536, + 7.545377637590768, + 7.545683239427934, + 7.5459888412651, + 7.546294443102266, + 7.5466000449394315, + 7.546905646776597, + 7.547211248613765, + 7.547516850450931, + 7.547822452288097, + 7.548128054125263, + 7.5484336559624285, + 7.548739257799596, + 7.549044859636762, + 7.549350461473928, + 7.549656063311094, + 7.54996166514826, + 7.550267266985427, + 7.550572868822593, + 7.550878470659759, + 7.551184072496925, + 7.551489674334091, + 7.551795276171257 + ], + "y": [ + 1306.60205078125, + 1316.4267578125, + 1298.068359375, + 1279.7294921875, + 1269.0380859375, + 1316.6865234375, + 1381.2861328125, + 1466.94091796875, + 1491.8486328125, + 1472.84521484375, + 1495.01513671875, + 1521.20263671875, + 1587.6630859375, + 1615.69384765625, + 1686.85498046875, + 1802.880859375, + 1813.015625, + 1831.43310546875, + 1928.32861328125, + 2002.90185546875, + 2135.2509765625, + 2199.71728515625, + 2180.5654296875, + 2228.65185546875, + 2365.5224609375, + 2547.42919921875, + 2754.6455078125, + 2957.9697265625, + 3092.22314453125, + 3159.67578125, + 3237.07666015625, + 3390.29931640625, + 3604.94140625, + 3874.33984375, + 4123.71923828125, + 4323.8447265625, + 4552.0673828125, + 4852.697265625, + 5175.55322265625, + 5480.50830078125, + 5913.18017578125, + 6489.40234375, + 7052.8466796875, + 7674.63525390625, + 8481.373046875, + 9317.92138671875, + 10245.50927734375, + 11324.927734375, + 12558.79052734375, + 13925.37255859375, + 15488.56689453125, + 17312.24609375, + 19441.28125, + 21826.310546875, + 24467.40869140625, + 27407.7412109375, + 30665.35595703125, + 34308.72607421875, + 38178.13330078125, + 42196.86279296875, + 46236.62060546875, + 50480.478515625, + 54883.96044921875, + 59254.93603515625, + 63547.91064453125, + 67728.46142578125, + 71609.025390625, + 75196.11279296875, + 78395.541015625, + 81394.177734375, + 84191.90283203125, + 86623.78466796875, + 88665.4931640625, + 90193.92724609375, + 91034.48046875, + 91189.37744140625, + 90775.009765625, + 89560.7890625, + 87490.92822265625, + 84628.22900390625, + 81128.63330078125, + 77149.76904296875, + 72808.2236328125, + 68147.27392578125, + 63241.89111328125, + 58421.8359375, + 53892.9072265625, + 49788.8193359375, + 46081.00439453125, + 42807.052734375, + 40016.07568359375, + 37686.0634765625, + 35667.06298828125, + 33707.78955078125, + 31797.22412109375, + 29907.927734375, + 28139.4150390625, + 26466.216796875, + 24941.85107421875, + 23679.0244140625, + 22733.46630859375, + 22261.59619140625, + 22171.8505859375, + 22421.84619140625, + 23085.67724609375, + 24106.60107421875, + 25475.21875, + 27213.99560546875, + 29341.9248046875, + 31801.29931640625, + 34442.64990234375, + 37406.04638671875, + 40851.873046875, + 44775.18896484375, + 49205.67138671875, + 54332.30078125, + 60151.326171875, + 66667.83056640625, + 73801.1640625, + 81454.25732421875, + 89616.77587890625, + 98286.294921875, + 107285.6669921875, + 116457.1611328125, + 125637.12939453125, + 134819.685546875, + 143939.326171875, + 153032.85205078125, + 161814.50244140625, + 170015.36181640625, + 177581.64990234375, + 184322.77197265625, + 190284.6640625, + 195430.74951171875, + 199672.4599609375, + 203138.01220703125, + 205604.5126953125, + 206656.5712890625, + 206060.24658203125, + 203634.84912109375, + 199552.1015625, + 193720.3203125, + 186250.6494140625, + 177471.28662109375, + 167524.4423828125, + 156932.49462890625, + 146119.2763671875, + 135190.9990234375, + 124356.01708984375, + 113892.16796875, + 103942.537109375, + 94506.7900390625, + 85638.806640625, + 77227.6640625, + 69358.2138671875, + 62056.341796875, + 55381.75390625, + 49310.04443359375, + 43903.56640625, + 39209.615234375, + 35184.61083984375, + 31799.74755859375, + 29135.54248046875, + 27089.86962890625, + 25590.6484375, + 24496.0576171875, + 23812.060546875, + 23432.0693359375, + 23420.697265625, + 23833.1416015625, + 24631.55126953125, + 25758.27880859375, + 27138.60791015625, + 28771.58544921875, + 30646.6318359375, + 32727.5146484375, + 34940.54443359375, + 37492.97412109375, + 40401.79052734375, + 43614.154296875, + 47106.962890625, + 51068.90576171875, + 55614.79638671875, + 60719.81494140625, + 66111.685546875, + 71577.94482421875, + 77063.2587890625, + 82714.38720703125, + 88476.27880859375, + 94079.37158203125, + 99398.3291015625, + 104297.75390625, + 108684.2509765625, + 112438.51416015625, + 115592.18896484375, + 118275.28466796875, + 120441.60400390625, + 122054.8115234375, + 123088.75927734375, + 123563.40673828125, + 123476.16064453125, + 122837.2861328125, + 121380.98486328125, + 118908.73388671875, + 115344.431640625, + 110720.8330078125, + 105186.533203125, + 98792.193359375, + 91580.96435546875, + 83727.458984375, + 75512.0498046875, + 67268.73876953125, + 59480.33544921875, + 52338.44775390625, + 45879.71630859375, + 40337.5048828125, + 35710.47265625, + 31768.3076171875, + 28442.650390625, + 25571.68212890625, + 22897.7041015625, + 20434.2900390625, + 18228.9013671875, + 16149.2568359375, + 14312.8369140625, + 12744.9580078125, + 11439.06103515625, + 10448.3583984375, + 9676.95947265625, + 8966.53125, + 8342.06005859375, + 7766.158203125, + 7124.79150390625, + 6662.2666015625, + 6360.26513671875, + 6044.79296875, + 5686.46533203125, + 5385.580078125, + 5077.259765625, + 4776.38232421875, + 4497.7119140625, + 4256.52734375, + 4087.5634765625, + 3841.275390625, + 3657.32470703125, + 3507.29638671875, + 3243.6328125, + 3045.40185546875, + 2979.54638671875, + 2833.37744140625, + 2674.27587890625, + 2560.095703125, + 2401.0625, + 2326.259765625, + 2215.7490234375, + 2123.59814453125, + 2030.19580078125, + 1949.86474609375, + 1890.271484375, + 1864.92919921875, + 1926.38427734375 + ] +} diff --git a/src/__tests__/gaussian-with-shoulders.test.ts b/src/__tests__/gaussian-with-shoulders.test.ts new file mode 100644 index 0000000..c1b2a50 --- /dev/null +++ b/src/__tests__/gaussian-with-shoulders.test.ts @@ -0,0 +1,37 @@ +import type { DataXY } from 'cheminfo-types'; +import { generateSpectrum } from 'spectrum-generator'; +import { describe, expect, it } from 'vitest'; + +import { gsd } from '../gsd.ts'; + +describe('gaussian overlapping', () => { + const expectedPeaks = [ + { x: 0.1, y: 0.3, width: 0.05 }, + { x: -0.1, y: 0.3, width: 0.05 }, + { x: 0, y: 1, width: 0.1 }, + ]; + + const data: DataXY = generateSpectrum(expectedPeaks, { + generator: { + from: -1, + to: 1, + nbPoints: 1001, + }, + peakOptions: { + factor: 6, // need a high factor so that we don't detect the end of the simulated peak + }, + }); + + it('1 peak should detected by first derivative algorithm', () => { + const peaks = gsd(data, { peakDetectionAlgorithm: 'first' }); + expect(peaks).toHaveLength(1); + expect(peaks[0].x).toBeCloseTo(0); + expect(peaks[0].y).toBeCloseTo(1); + }); + it('3 peak should detected by auto algorithm', () => { + const peaks = gsd(data, { peakDetectionAlgorithm: 'auto' }); + expect(peaks).toHaveLength(3); + expect(peaks[1].x).toBeCloseTo(0); + expect(peaks[1].y).toBeCloseTo(1); + }); +}); diff --git a/src/__tests__/noPerfectShape.test.ts b/src/__tests__/noPerfectShape.test.ts new file mode 100644 index 0000000..ca6743d --- /dev/null +++ b/src/__tests__/noPerfectShape.test.ts @@ -0,0 +1,45 @@ +import { readFileSync } from 'fs'; + +import { describe, expect, it } from 'vitest'; + +import { gsd } from '../gsd.ts'; + +describe('Global spectra deconvolution NMR spectra', () => { + // Test case obtained from Pag 443, Chap 8. + it('Should have 5 peaks', () => { + const spectrum: { x: number[]; y: number[] } = JSON.parse( + readFileSync(`${__dirname}/data/broadNotLorentzianShape.json`, 'utf-8'), + ); + const sgOptions = { windowSize: 7, polynomial: 3 }; + const peaks = gsd( + { x: spectrum.x, y: spectrum.y }, + { + minMaxRatio: 0.05, + noiseLevel: 1500, + smoothY: false, + realTopDetection: false, + maxCriteria: true, + sgOptions, + }, + ); + + expect(peaks).toHaveLength(3); + expect(peaks).toMatchObject([ + { + x: 7.49587013996984, + y: 89560.7890625, + width: 0.006723240417656484, + }, + { + x: 7.514817453874146, + y: 203634.84912109375, + width: 0.007028842254822365, + }, + { + x: 7.534375971452784, + y: 118908.73388671875, + width: 0.007028842254822365, + }, + ]); + }); +}); diff --git a/src/algorithms/PeakData.ts b/src/algorithms/PeakData.ts new file mode 100644 index 0000000..63c1ff3 --- /dev/null +++ b/src/algorithms/PeakData.ts @@ -0,0 +1,11 @@ +import type { NumberArray } from 'cheminfo-types'; + +export interface PeakData { + x: NumberArray; + y: NumberArray; + yData: NumberArray; + dY: NumberArray; + ddY: NumberArray; + yThreshold: number; + dX: number; +} diff --git a/src/algorithms/autoAlgorithm.ts b/src/algorithms/autoAlgorithm.ts new file mode 100644 index 0000000..a01a04c --- /dev/null +++ b/src/algorithms/autoAlgorithm.ts @@ -0,0 +1,96 @@ +import type { NumberArray } from 'cheminfo-types'; + +import type { GSDPeakID } from '../gsd.ts'; + +import { getMinMaxIntervalsDy } from './getMinMaxIntervals.ts'; +import { tryMatchOneIntervalWithMinData } from './tryMatchOneIntervalWithMinData.ts'; + +export function autoAlgorithm(input: { + x: NumberArray; + y: NumberArray; + yData: NumberArray; + dY: NumberArray; + ddY: NumberArray; + yThreshold: number; + dX: number; +}) { + const { x, y, yData, dY, ddY, dX, yThreshold } = input; + + const minddY: number[] = []; + const crossDy: number[] = []; + const { intervalL, intervalR } = getMinMaxIntervalsDy(y, x, dY, dX); + + for (let i = 1; i < y.length - 1; ++i) { + if ((dY[i] < 0 && dY[i + 1] > 0) || (dY[i] > 0 && dY[i + 1] < 0)) { + // push the index of the element closer to zero + crossDy.push(Math.abs(dY[i]) < Math.abs(dY[i + 1]) ? i : i + 1); + } + // Handle exact zero + if ( + dY[i] === 0 && + dY[i] < Math.abs(dY[i + 1]) && + dY[i] < Math.abs(dY[i - 1]) + ) { + crossDy.push(i); + } + + // Minimum in second derivative + if (ddY[i] < ddY[i - 1] && ddY[i] < ddY[i + 1]) { + minddY.push(i); + } + } + + const peaks: GSDPeakID[] = []; + let [lastK, lastJ] = [-1, -1]; + for (let i = 0; i < intervalL.length; i++) { + const intervalWidth = (intervalR[i].x - intervalL[i].x) / 2; + const intervalCenter = (intervalR[i].x + intervalL[i].x) / 2; + + let yIndex = -1; + let match = tryMatchOneIntervalWithMinData({ + x, + yData, + lastK, + yThreshold, + intervalWidth, + intervalCenter, + minData: crossDy, + }); + lastK = match.lastIndex; + if (match.possible !== -1) { + yIndex = crossDy[match.possible]; + } else { + match = tryMatchOneIntervalWithMinData({ + x, + yData, + yThreshold, + lastK: lastJ, + intervalWidth, + intervalCenter, + minData: minddY, + }); + if (match.possible !== -1) { + yIndex = minddY[match.possible]; + } + lastJ = match.lastIndex; + } + + if (yIndex !== -1) { + const width = Math.abs(intervalR[i].x - intervalL[i].x); + peaks.push({ + id: crypto.randomUUID(), + x: x[yIndex], + y: y[yIndex], + width, + index: yIndex, + ddY: ddY[yIndex], + inflectionPoints: { + from: intervalL[i], + to: intervalR[i], + }, + }); + } + } + + return peaks; +} diff --git a/src/algorithms/firstDerivative.ts b/src/algorithms/firstDerivative.ts new file mode 100644 index 0000000..5af6564 --- /dev/null +++ b/src/algorithms/firstDerivative.ts @@ -0,0 +1,35 @@ +import type { PeakData } from './PeakData.ts'; +import { getMinMaxIntervalsDy } from './getMinMaxIntervals.ts'; +import { getPeakFromIntervals } from './getPeaksFromIntervals.ts'; + +export function firstDerivative(input: PeakData) { + const { x, y, yData, dY, ddY, dX, yThreshold } = input; + + const crossDy: number[] = []; + const { intervalL, intervalR } = getMinMaxIntervalsDy(y, x, dY, dX); + + for (let i = 1; i < y.length - 1; ++i) { + if ((dY[i] < 0 && dY[i + 1] > 0) || (dY[i] > 0 && dY[i + 1] < 0)) { + // push the index of the element closer to zero + crossDy.push(Math.abs(dY[i]) < Math.abs(dY[i + 1]) ? i : i + 1); + } + // Handle exact zero + if ( + dY[i] === 0 && + dY[i] < Math.abs(dY[i + 1]) && + dY[i] < Math.abs(dY[i - 1]) + ) { + crossDy.push(i); + } + } + + return getPeakFromIntervals({ + minData: crossDy, + intervalL, + intervalR, + x, + yData, + yThreshold, + ddY, + }); +} diff --git a/src/algorithms/getMinMaxIntervals.ts b/src/algorithms/getMinMaxIntervals.ts new file mode 100644 index 0000000..f6fe382 --- /dev/null +++ b/src/algorithms/getMinMaxIntervals.ts @@ -0,0 +1,47 @@ +import type { NumberArray } from 'cheminfo-types'; + +import type { XIndex } from '../XIndex.ts'; + +export function getMinMaxIntervalsDy( + y: NumberArray, + x: NumberArray, + dY: NumberArray, + dX: number, +) { + let lastMax: XIndex | null = null; + let lastMin: XIndex | null = null; + const intervalL: XIndex[] = []; + const intervalR: XIndex[] = []; + for (let i = 1; i < y.length - 1; ++i) { + if ( + (dY[i] < dY[i - 1] && dY[i] <= dY[i + 1]) || + (dY[i] <= dY[i - 1] && dY[i] < dY[i + 1]) + ) { + lastMin = { + x: x[i], + index: i, + }; + if (dX > 0 && lastMax !== null) { + intervalL.push(lastMax); + intervalR.push(lastMin); + } + } + + // Maximum in first derivative + if ( + (dY[i] >= dY[i - 1] && dY[i] > dY[i + 1]) || + (dY[i] > dY[i - 1] && dY[i] >= dY[i + 1]) + ) { + lastMax = { + x: x[i], + index: i, + }; + if (dX < 0 && lastMin !== null) { + intervalL.push(lastMax); + intervalR.push(lastMin); + } + } + } + + return { intervalL, intervalR }; +} diff --git a/src/algorithms/getPeaksFromIntervals.ts b/src/algorithms/getPeaksFromIntervals.ts new file mode 100644 index 0000000..b1310fc --- /dev/null +++ b/src/algorithms/getPeaksFromIntervals.ts @@ -0,0 +1,51 @@ +import type { XIndex } from '../XIndex.ts'; +import type { GSDPeakID } from '../gsd.ts'; + +import type { PeakData } from './PeakData.ts'; +import { tryMatchOneIntervalWithMinData } from './tryMatchOneIntervalWithMinData.ts'; + +export function getPeakFromIntervals( + options: Pick & { + intervalR: XIndex[]; + intervalL: XIndex[]; + minData: number[]; + }, +) { + let lastK = -1; + const peaks: GSDPeakID[] = []; + const { x, ddY, yData, yThreshold, intervalR, intervalL, minData } = options; + + for (let i = 0; i < intervalL.length; i++) { + const intervalWidth = (intervalR[i].x - intervalL[i].x) / 2; + const intervalCenter = (intervalR[i].x + intervalL[i].x) / 2; + const { possible = -1, lastIndex } = tryMatchOneIntervalWithMinData({ + x, + lastK, + minData, + yThreshold, + intervalWidth, + intervalCenter, + yData, + }); + + if (possible !== -1) { + const centerIndex = minData[possible]; + const width = Math.abs(intervalR[i].x - intervalL[i].x); + peaks.push({ + id: crypto.randomUUID(), + x: x[centerIndex], + y: yData[centerIndex], + width, + index: centerIndex, + ddY: ddY[centerIndex], + inflectionPoints: { + from: intervalL[i], + to: intervalR[i], + }, + }); + } + lastK = lastIndex; + } + + return peaks; +} diff --git a/src/algorithms/secondDerivative.ts b/src/algorithms/secondDerivative.ts new file mode 100644 index 0000000..9e18ef3 --- /dev/null +++ b/src/algorithms/secondDerivative.ts @@ -0,0 +1,37 @@ +import type { NumberArray } from 'cheminfo-types'; + +import { getMinMaxIntervalsDy } from './getMinMaxIntervals.ts'; +import { getPeakFromIntervals } from './getPeaksFromIntervals.ts'; + +export function secondDerivative(input: { + x: NumberArray; + y: NumberArray; + yData: NumberArray; + dY: NumberArray; + ddY: NumberArray; + yThreshold: number; + dX: number; +}) { + const { x, y, yData, dY, ddY, dX, yThreshold } = input; + + const minddY: number[] = []; + const { intervalL, intervalR } = getMinMaxIntervalsDy(y, x, dY, dX); + + // By the intermediate value theorem We cannot find 2 consecutive maximum or minimum + for (let i = 1; i < y.length - 1; ++i) { + // Minimum in second derivative + if (ddY[i] < ddY[i - 1] && ddY[i] < ddY[i + 1]) { + minddY.push(i); + } + } + + return getPeakFromIntervals({ + minData: minddY, + intervalL, + intervalR, + x, + yData, + yThreshold, + ddY, + }); +} diff --git a/src/algorithms/tryMatchOneIntervalWithMinData.ts b/src/algorithms/tryMatchOneIntervalWithMinData.ts new file mode 100644 index 0000000..a904987 --- /dev/null +++ b/src/algorithms/tryMatchOneIntervalWithMinData.ts @@ -0,0 +1,50 @@ +import type { NumberArray } from 'cheminfo-types'; + +interface TryMatchOneIntervalWithMinDataOptions { + x: NumberArray; + lastK: number; + minData: number[]; + yThreshold: number; + intervalWidth: number; + intervalCenter: number; + yData: NumberArray; +} + +export function tryMatchOneIntervalWithMinData( + options: TryMatchOneIntervalWithMinDataOptions, +) { + const { + x, + lastK, + minData, + yThreshold, + intervalWidth, + intervalCenter, + yData, + } = options; + + let minDistance = Number.POSITIVE_INFINITY; + let possible = -1; + let newLastIndex = lastK; + for (let k = newLastIndex + 1; k < minData.length; k++) { + const centerIndex = minData[k]; + if (yData[centerIndex] <= yThreshold) { + continue; + } + + const deltaX = x[centerIndex]; + const currentDistance = Math.abs(deltaX - intervalCenter); + + if (currentDistance < intervalWidth) { + if (currentDistance < minDistance) { + possible = k; + } + newLastIndex = k; + } + + if (currentDistance >= minDistance) break; + minDistance = currentDistance; + } + + return { lastIndex: newLastIndex, possible }; +} diff --git a/src/gsd.ts b/src/gsd.ts index 4a14672..8bfc592 100644 --- a/src/gsd.ts +++ b/src/gsd.ts @@ -9,6 +9,9 @@ import { } from 'ml-spectra-processing'; import type { GSDPeak } from './GSDPeak.ts'; +import { autoAlgorithm } from './algorithms/autoAlgorithm.ts'; +import { firstDerivative } from './algorithms/firstDerivative.ts'; +import { secondDerivative } from './algorithms/secondDerivative.ts'; import type { MakeMandatory } from './utils/MakeMandatory.ts'; import { optimizeTop } from './utils/optimizeTop.ts'; @@ -45,7 +48,18 @@ export interface GSDOptions { * @default false */ realTopDetection?: boolean; + /** + * Algorithm used for peak detection: + * - 'first': Uses the first derivative to detect peaks, detecting only + * the peaks which first derivative cross the zero. + * - 'second': Uses the second derivative to detect peaks (inflection points). + * - 'auto': Automatically selects the peaks by checking the zero crossing of the first derivative + * or the local minima in the second derivative. + * @default 'second' + */ + peakDetectionAlgorithm?: 'first' | 'second' | 'auto'; } + export type GSDPeakID = MakeMandatory; /** @@ -67,6 +81,7 @@ export function gsd(data: DataXY, options: GSDOptions = {}): GSDPeakID[] { maxCriteria = true, minMaxRatio = 0.00025, realTopDetection = false, + peakDetectionAlgorithm = 'second', } = options; const { x } = data; let { y } = data; @@ -78,10 +93,10 @@ export function gsd(data: DataXY, options: GSDOptions = {}): GSDPeakID[] { // If the max difference between delta x is less than 5%, then, // we can assume it to be equally spaced variable - const equallySpaced = xIsEquallySpaced(x); + const isEquallySpaced = xIsEquallySpaced(x); if (noiseLevel === undefined) { - if (equallySpaced) { + if (isEquallySpaced) { const noiseInfo = xNoiseStandardDeviation(y); if (maxCriteria) { noiseLevel = noiseInfo.median + 1.5 * noiseInfo.sd; @@ -108,7 +123,7 @@ export function gsd(data: DataXY, options: GSDOptions = {}): GSDPeakID[] { } } - const xValue = equallySpaced ? x[1] - x[0] : x; + const xValue = isEquallySpaced ? x[1] - x[0] : x; const yData = smoothY ? sgg(y, xValue, { @@ -117,114 +132,31 @@ export function gsd(data: DataXY, options: GSDOptions = {}): GSDPeakID[] { }) : y; + const { min: minY, max: maxY } = xMinMaxValues(yData); + if (minY > maxY || minY === maxY) return []; + const dY = sgg(y, xValue, { ...sgOptions, derivative: 1, }); + const ddY = sgg(y, xValue, { ...sgOptions, derivative: 2, }); - const { min: minY, max: maxY } = xMinMaxValues(yData); - - if (minY > maxY || minY === maxY) return []; - const yThreshold = minY + (maxY - minY) * minMaxRatio; const dX = x[1] - x[0]; - interface XIndex { - x: number; - index: number; - } - - let lastMax: XIndex | null = null; - let lastMin: XIndex | null = null; - const minddY: number[] = []; - const intervalL: XIndex[] = []; - const intervalR: XIndex[] = []; - - // By the intermediate value theorem We cannot find 2 consecutive maximum or minimum - for (let i = 1; i < yData.length - 1; ++i) { - if ( - (dY[i] < dY[i - 1] && dY[i] <= dY[i + 1]) || - (dY[i] <= dY[i - 1] && dY[i] < dY[i + 1]) - ) { - lastMin = { - x: x[i], - index: i, - }; - if (dX > 0 && lastMax !== null) { - intervalL.push(lastMax); - intervalR.push(lastMin); - } - } - - // Maximum in first derivative - if ( - (dY[i] >= dY[i - 1] && dY[i] > dY[i + 1]) || - (dY[i] > dY[i - 1] && dY[i] >= dY[i + 1]) - ) { - lastMax = { - x: x[i], - index: i, - }; - if (dX < 0 && lastMin !== null) { - intervalL.push(lastMax); - intervalR.push(lastMin); - } - } - - // Minimum in second derivative - if (ddY[i] < ddY[i - 1] && ddY[i] < ddY[i + 1]) { - minddY.push(i); - } - } - - let lastK = -1; - - const peaks: GSDPeakID[] = []; - for (const minddYIndex of minddY) { - const deltaX = x[minddYIndex]; - let possible = -1; - let k = lastK + 1; - let minDistance = Number.POSITIVE_INFINITY; - let currentDistance = 0; - while (possible === -1 && k < intervalL.length) { - currentDistance = Math.abs( - deltaX - (intervalL[k].x + intervalR[k].x) / 2, - ); - if (currentDistance < (intervalR[k].x - intervalL[k].x) / 2) { - possible = k; - lastK = k; - } - ++k; - - // Not getting closer? - if (currentDistance >= minDistance) { - break; - } - minDistance = currentDistance; - } - - if (possible !== -1) { - if (yData[minddYIndex] > yThreshold) { - const width = Math.abs(intervalR[possible].x - intervalL[possible].x); - peaks.push({ - id: crypto.randomUUID(), - x: deltaX, - y: yData[minddYIndex], - width, - index: minddYIndex, - ddY: ddY[minddYIndex], - inflectionPoints: { - from: intervalL[possible], - to: intervalR[possible], - }, - }); - } - } + const peakData = { x, y, yData, dY, ddY, dX, yThreshold }; + let peaks: GSDPeakID[] = []; + if (peakDetectionAlgorithm === 'first') { + peaks = firstDerivative(peakData); + } else if (peakDetectionAlgorithm === 'second') { + peaks = secondDerivative(peakData); + } else { + peaks = autoAlgorithm(peakData); } if (realTopDetection) {