-
Notifications
You must be signed in to change notification settings - Fork 8
/
optimizePeaksWithLogs.ts
125 lines (110 loc) · 3.68 KB
/
optimizePeaksWithLogs.ts
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
import type { DataXY, PeakXYWidth } from 'cheminfo-types';
import { getShape1D } from 'ml-peak-shape-generator';
import { optimize } from 'ml-spectra-fitting';
import { xGetFromToIndex } from 'ml-spectra-processing';
import { GSDPeakOptimized } from '../GSDPeakOptimized';
import { MakeMandatory } from '../utils/MakeMandatory';
import { addMissingShape } from '../utils/addMissingShape';
import { groupPeaks } from '../utils/groupPeaks';
import { OptimizePeaksOptions } from './optimizePeaks';
export interface Peak extends PeakXYWidth {
id?: string;
}
export type GSDPeakOptimizedID = MakeMandatory<GSDPeakOptimized, 'id'>;
type GSDPeakOptimizedIDOrNot<T extends Peak> = T extends {
id: string;
}
? GSDPeakOptimizedID
: GSDPeakOptimized;
/**
* Optimize the position (x), max intensity (y), full width at half maximum (fwhm)
* and the ratio of gaussian contribution (mu) if it's required. It currently supports three kind of shapes: gaussian, lorentzian and pseudovoigt
* @param data - An object containing the x and y data to be fitted.
* @param peakList - A list of initial parameters to be optimized. e.g. coming from a peak picking [{x, y, width}].
*/
export function optimizePeaksWithLogs<T extends Peak>(
data: DataXY,
peakList: T[],
options: OptimizePeaksOptions = {},
): { logs: any[]; optimizedPeaks: Array<GSDPeakOptimizedIDOrNot<T>> } {
const {
fromTo = {},
baseline,
shape = { kind: 'gaussian' },
groupingFactor = 1,
factorLimits = 2,
optimization = {
kind: 'lm',
options: {
timeout: 10,
},
},
}: OptimizePeaksOptions = options;
/*
The optimization algorithm will take some group of peaks.
We can not simply optimize everything because there would be too many variables to optimize
and it would be too time consuming.
*/
const groups = groupPeaks(peakList, { factor: groupingFactor });
const logs: any[] = [];
const results: Array<GSDPeakOptimizedIDOrNot<T>> = [];
groups.forEach((peakGroup) => {
const start = Date.now();
// In order to make optimization we will add fwhm and shape on all the peaks
const peaks = addMissingShape(peakGroup, { shape });
const firstPeak = peaks[0];
const lastPeak = peaks[peaks.length - 1];
const {
from = firstPeak.x - firstPeak.width * factorLimits,
to = lastPeak.x + lastPeak.width * factorLimits,
} = fromTo;
const { fromIndex, toIndex } = xGetFromToIndex(data.x, { from, to });
const x =
data.x instanceof Float64Array
? data.x.subarray(fromIndex, toIndex)
: data.x.slice(fromIndex, toIndex);
const y =
data.y instanceof Float64Array
? data.y.subarray(fromIndex, toIndex)
: data.y.slice(fromIndex, toIndex);
const log = {
range: { from, to },
parameters: optimization,
groupSize: peakGroup.length,
time: Date.now() - start,
};
if (x.length > 5) {
const {
iterations,
error,
peaks: optimizedPeaks,
} = optimize({ x, y }, peaks, {
shape,
baseline,
optimization,
});
for (let i = 0; i < peaks.length; i++) {
results.push({
...optimizedPeaks[i],
width: getShape1D(peaks[i].shape).fwhmToWidth(
optimizedPeaks[i].shape.fwhm,
),
} as GSDPeakOptimizedIDOrNot<T>);
}
logs.push({
...log,
iterations,
error,
message: 'optimization successful',
});
} else {
results.push(...(peaks as Array<GSDPeakOptimizedIDOrNot<T>>));
logs.push({
...log,
iterations: 0,
message: 'x length too small for optimization',
});
}
});
return { logs, optimizedPeaks: results };
}