1use itertools::multizip;
2use mzdeisotope::{scorer::ScoreType, solution::DeconvolvedSolutionPeak};
3use num_traits::Zero;
4use std::{
5 borrow::Cow,
6 boxed::Box,
7 cmp::Ordering,
8 ops::{Bound, RangeBounds},
9};
10use tracing::warn;
11
12use mzpeaks::{
13 feature::{ChargedFeature, FeatureView, TimeArray},
14 peak::MZPoint,
15 prelude::*,
16 IonMobility, Mass, MZ,
17};
18
19use mzsignal::{
20 feature_mapping::graph::ChargeAwareFeatureMerger,
21 feature_statistics::{FitPeaksOn, PeakFitArgs},
22};
23
24use mzdata::{
25 prelude::*,
26 spectrum::{BinaryArrayMap, BinaryDataArrayType},
27 utils::mass_charge_ratio,
28};
29use mzdata::{
30 spectrum::bindata::{
31 ArrayRetrievalError, ArrayType, BinaryArrayMap3D, BuildArrayMap3DFrom, BuildFromArrayMap3D,
32 DataArray,
33 },
34 utils::neutral_mass,
35};
36
37#[derive(Default, Debug, Clone)]
38#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
39pub struct MZPointSeries {
40 mz: Vec<f64>,
41 intensity: Vec<f32>,
42}
43
44impl<'a> MZPointSeries {
45 pub fn new(mz: Vec<f64>, intensity: Vec<f32>) -> Self {
46 Self { mz, intensity }
47 }
48
49 pub fn push<T: CoordinateLike<MZ> + IntensityMeasurement>(&mut self, pt: T) {
50 self.push_raw(pt.mz(), pt.intensity());
51 }
52
53 pub fn push_raw(&mut self, mz: f64, intensity: f32) {
54 self.mz.push(mz);
55 self.intensity.push(intensity);
56 }
57
58 pub fn split_at(&self, i: usize) -> (Self, Self) {
59 let mz_a = self.mz[..i].to_vec();
60 let mz_b = self.mz[i..].to_vec();
61
62 let inten_a = self.intensity[..i].to_vec();
63 let inten_b = self.intensity[i..].to_vec();
64
65 (Self::new(mz_a, inten_a), Self::new(mz_b, inten_b))
66 }
67
68 pub fn slice<I: RangeBounds<usize> + Clone>(&self, bounds: I) -> Self {
69 let start = match bounds.start_bound() {
70 Bound::Included(i) | Bound::Excluded(i) => *i,
71 Bound::Unbounded => 0,
72 };
73
74 let end = match bounds.end_bound() {
75 Bound::Included(i) => *i + 1,
76 Bound::Excluded(i) => *i,
77 Bound::Unbounded => self.mz.len(),
78 };
79
80 Self::new(
81 self.mz[start..end].to_vec(),
82 self.intensity[start..end].to_vec(),
83 )
84 }
85
86 pub fn len(&self) -> usize {
87 self.mz.len()
88 }
89
90 pub fn is_empty(&self) -> bool {
91 self.mz.is_empty()
92 }
93
94 pub fn at(&self, index: usize) -> Option<MZPoint> {
95 if index < self.len() {
96 Some(MZPoint::new(self.mz[index], self.intensity[index]))
97 } else {
98 None
99 }
100 }
101
102 pub fn iter(&self) -> std::iter::Zip<std::slice::Iter<'_, f64>, std::slice::Iter<'_, f32>> {
103 self.mz.iter().zip(self.intensity.iter())
104 }
105
106 pub fn as_feature_view<Y>(&'a self, time: &'a [f64]) -> FeatureView<'a, MZ, Y> {
107 let start = 0;
108 let end = time.len().min(self.len());
109
110 FeatureView::new(&self.mz[start..end], time, &self.intensity[start..end])
111 }
112}
113
114#[derive(Debug, Clone)]
115#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
116pub struct DeconvolvedSolutionFeature<Y: Clone> {
117 inner: ChargedFeature<Mass, Y>,
118 pub score: ScoreType,
119 pub scores: Vec<ScoreType>,
120 envelope: Box<[MZPointSeries]>,
121}
122
123impl<'a, Y: Clone + 'a> FitPeaksOn<'a> for DeconvolvedSolutionFeature<Y> {}
124
125impl<'a, Y: Clone + 'a> From<&'a DeconvolvedSolutionFeature<Y>> for PeakFitArgs<'a, 'a> {
126 fn from(value: &'a DeconvolvedSolutionFeature<Y>) -> Self {
127 let view = value.as_inner();
128 let view = view.as_view();
129 let view = view.into_inner().0;
130 let (_, y, z) = view.into_inner();
131 PeakFitArgs::new(Cow::Borrowed(y), Cow::Borrowed(z))
132 }
133}
134
135impl<Y: Clone> Default for DeconvolvedSolutionFeature<Y> {
136 fn default() -> Self {
137 Self {
138 inner: ChargedFeature::empty(0),
139 score: Default::default(),
140 scores: Default::default(),
141 envelope: Default::default(),
142 }
143 }
144}
145
146impl<Y: Clone> AsPeakIter for DeconvolvedSolutionFeature<Y> {
147 type Peak = DeconvolvedSolutionPeak;
148
149 type Iter<'a>
150 = PeakIter<'a, Y>
151 where
152 Self: 'a;
153
154 fn iter_peaks(&self) -> Self::Iter<'_> {
155 self.iter_peaks()
156 }
157}
158
159impl<Y: Clone> BuildFromPeak<DeconvolvedSolutionPeak> for DeconvolvedSolutionFeature<Y> {
160 fn push_peak(&mut self, value: &DeconvolvedSolutionPeak, time: f64) {
161 self.push_peak(value, time);
162 }
163}
164
165impl<Y: Clone> KnownChargeMut for DeconvolvedSolutionFeature<Y> {
166 fn charge_mut(&mut self) -> &mut i32 {
167 &mut self.inner.charge
168 }
169}
170
171impl<Y: Clone> DeconvolvedSolutionFeature<Y> {
172 pub fn new(
173 inner: ChargedFeature<Mass, Y>,
174 score: ScoreType,
175 scores: Vec<ScoreType>,
176 envelope: Box<[MZPointSeries]>,
177 ) -> Self {
178 Self {
179 inner,
180 score,
181 scores,
182 envelope,
183 }
184 }
185
186 pub fn as_inner(&self) -> &ChargedFeature<Mass, Y> {
187 &self.inner
188 }
189
190 pub fn into_inner(
191 self,
192 ) -> (
193 ChargedFeature<Mass, Y>,
194 Vec<ScoreType>,
195 Box<[MZPointSeries]>,
196 ) {
197 (self.inner, self.scores, self.envelope)
198 }
199
200 pub fn len(&self) -> usize {
201 self.inner.len()
202 }
203
204 pub fn is_empty(&self) -> bool {
205 self.inner.is_empty()
206 }
207
208 pub fn iter(&self) -> mzpeaks::feature::Iter<'_, Mass, Y> {
209 self.inner.iter()
210 }
211
212 pub fn iter_mut(&mut self) -> mzpeaks::feature::IterMut<'_, Mass, Y> {
213 self.inner.iter_mut()
214 }
215
216 pub fn iter_peaks(&self) -> PeakIter<'_, Y> {
217 PeakIter::new(self)
218 }
219
220 pub fn iter_envelope(&self) -> EnvelopeIter<'_, Y> {
221 EnvelopeIter::new(self)
222 }
223
224 pub fn push<T: CoordinateLike<Mass> + IntensityMeasurement>(&mut self, pt: &T, time: f64) {
225 self.inner.push(pt, time);
226 self.scores.push(0.0);
227 }
228
229 pub fn push_peak(&mut self, peak: &DeconvolvedSolutionPeak, time: f64) {
230 let n_before = self.inner.len();
231 self.inner.push_raw(peak.neutral_mass, time, peak.intensity);
232 let n_after = self.len();
233 let did_resize = n_after != n_before;
234 if did_resize {
235 self.scores.push(peak.score);
236 if self.envelope.is_empty() {
237 let mut env_set = Vec::new();
238 for pt in peak.envelope.iter() {
239 let mut series = MZPointSeries::default();
240 series.push(pt);
241 env_set.push(series);
242 }
243 self.envelope = env_set.into_boxed_slice();
244 } else {
245 self.envelope
246 .iter_mut()
247 .zip(peak.envelope.iter())
248 .for_each(|(ev, pt)| ev.push(pt));
249 }
250 } else {
251 let q = self.scores.last_mut().unwrap();
252 *q = peak.score.max(*q);
253 self.envelope
254 .iter_mut()
255 .zip(peak.envelope.iter())
256 .for_each(|(ev, pt)| {
257 let last_int = *ev.intensity.last().unwrap();
258 let int = last_int + pt.intensity;
259 let last_mz = *ev.mz.last().unwrap();
260 let mz =
261 ((pt.mz * pt.intensity as f64) + (last_mz * last_int as f64)) / int as f64;
262 *ev.mz.last_mut().unwrap() = mz;
263 *ev.intensity.last_mut().unwrap() = int
264 });
265 }
266 }
267
268 pub fn envelope(&self) -> Vec<FeatureView<'_, MZ, Y>> {
269 let times = self.inner.time_view();
270 self.envelope
271 .iter()
272 .map(|s| s.as_feature_view(times))
273 .collect()
274 }
275}
276
277impl<Y0: Clone> AsRef<ChargedFeature<Mass, Y0>> for DeconvolvedSolutionFeature<Y0> {
278 fn as_ref(&self) -> &ChargedFeature<Mass, Y0> {
279 &self.inner
280 }
281}
282
283impl<Y0: Clone> AsMut<ChargedFeature<Mass, Y0>> for DeconvolvedSolutionFeature<Y0> {
284 fn as_mut(&mut self) -> &mut ChargedFeature<Mass, Y0> {
285 &mut self.inner
286 }
287}
288
289impl<Y0: Clone> FeatureLikeMut<Mass, Y0> for DeconvolvedSolutionFeature<Y0> {
290 fn clear(&mut self) {
291 self.inner.clear();
292 for e in self.envelope.iter_mut() {
293 e.intensity.clear();
294 e.mz.clear();
295 }
296 }
297
298 fn iter_mut(&mut self) -> impl Iterator<Item = (&mut f64, &mut f64, &mut f32)> {
299 <ChargedFeature<Mass, Y0> as FeatureLikeMut<Mass, Y0>>::iter_mut(&mut self.inner)
300 }
301
302 fn push<T: CoordinateLike<Mass> + IntensityMeasurement>(&mut self, pt: &T, time: f64) {
303 <ChargedFeature<Mass, Y0> as FeatureLikeMut<Mass, Y0>>::push(&mut self.inner, pt, time)
304 }
305
306 fn push_raw(&mut self, x: f64, y: f64, z: f32) {
307 <ChargedFeature<Mass, Y0> as FeatureLikeMut<Mass, Y0>>::push_raw(&mut self.inner, x, y, z)
308 }
309}
310
311impl<Y0: Clone> TimeInterval<Y0> for DeconvolvedSolutionFeature<Y0> {
312 fn apex_time(&self) -> Option<f64> {
313 <ChargedFeature<Mass, Y0> as TimeInterval<Y0>>::apex_time(&self.inner)
314 }
315
316 fn area(&self) -> f32 {
317 <ChargedFeature<Mass, Y0> as TimeInterval<Y0>>::area(&self.inner)
318 }
319
320 fn end_time(&self) -> Option<f64> {
321 <ChargedFeature<Mass, Y0> as TimeInterval<Y0>>::end_time(&self.inner)
322 }
323
324 fn start_time(&self) -> Option<f64> {
325 <ChargedFeature<Mass, Y0> as TimeInterval<Y0>>::start_time(&self.inner)
326 }
327
328 fn iter_time(&self) -> impl Iterator<Item = f64> {
329 <ChargedFeature<Mass, Y0> as TimeInterval<Y0>>::iter_time(&self.inner)
330 }
331
332 fn find_time(&self, time: f64) -> (Option<usize>, f64) {
333 <ChargedFeature<Mass, Y0> as TimeInterval<Y0>>::find_time(&self.inner, time)
334 }
335}
336
337impl<Y0: Clone> TimeArray<Y0> for DeconvolvedSolutionFeature<Y0> {
338 fn time_view(&self) -> &[f64] {
339 self.inner.time_view()
340 }
341
342 fn intensity_view(&self) -> &[f32] {
343 self.inner.intensity_view()
344 }
345}
346
347impl<Y0: Clone> FeatureLike<Mass, Y0> for DeconvolvedSolutionFeature<Y0> {
348 fn len(&self) -> usize {
349 <ChargedFeature<Mass, Y0> as FeatureLike<Mass, Y0>>::len(&self.inner)
350 }
351
352 fn iter(&self) -> impl Iterator<Item = (f64, f64, f32)> {
353 <ChargedFeature<Mass, Y0> as FeatureLike<Mass, Y0>>::iter(&self.inner)
354 }
355}
356
357impl<Y: Clone> PartialOrd for DeconvolvedSolutionFeature<Y> {
358 fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
359 if self == other {
360 return Some(Ordering::Equal);
361 }
362 match self.neutral_mass().total_cmp(&other.neutral_mass()) {
363 Ordering::Equal => {}
364 x => return Some(x),
365 };
366 self.start_time()
367 .partial_cmp(&other.start_time())
368 .map(|c| c.then(self.score.total_cmp(&other.score)))
369 }
370}
371
372impl<Y: Clone> PartialEq for DeconvolvedSolutionFeature<Y> {
373 fn eq(&self, other: &Self) -> bool {
374 self.inner == other.inner && self.score == other.score
375 }
376}
377
378impl<Y0: Clone> CoordinateLike<Mass> for DeconvolvedSolutionFeature<Y0> {
379 fn coordinate(&self) -> f64 {
380 <ChargedFeature<Mass, Y0> as CoordinateLike<Mass>>::coordinate(&self.inner)
381 }
382}
383
384impl<Y0: Clone> CoordinateLike<MZ> for DeconvolvedSolutionFeature<Y0> {
385 fn coordinate(&self) -> f64 {
386 <ChargedFeature<Mass, Y0> as CoordinateLike<MZ>>::coordinate(&self.inner)
387 }
388}
389
390impl<Y0: Clone> KnownCharge for DeconvolvedSolutionFeature<Y0> {
391 fn charge(&self) -> i32 {
392 <ChargedFeature<Mass, Y0> as KnownCharge>::charge(&self.inner)
393 }
394}
395
396impl<Y0: Clone> IntensityMeasurement for DeconvolvedSolutionFeature<Y0> {
397 fn intensity(&self) -> f32 {
398 <ChargedFeature<Mass, Y0> as IntensityMeasurement>::intensity(&self.inner)
399 }
400}
401
402impl<Y: Clone> SplittableFeatureLike<'_, Mass, Y> for DeconvolvedSolutionFeature<Y> {
403 type ViewType = DeconvolvedSolutionFeature<Y>;
404
405 fn split_at_time(&self, point: f64) -> (Self::ViewType, Self::ViewType) {
406 if let (Some(idx), _) = self.find_time(point) {
407 let (before, after) = self.inner.split_at_time(point);
408 let mut envelope_before = Vec::new();
409 let mut envelope_after = Vec::new();
410 for (env_before_i, env_after_i) in self.envelope.iter().map(|e| {
411 let (a, b) = e.split_at(idx);
412 (a.to_owned(), b.to_owned())
413 }) {
414 envelope_before.push(env_before_i);
415 envelope_after.push(env_after_i);
416 }
417 (
418 Self::new(
419 before.to_owned(),
420 self.score,
421 self.scores[..idx].to_vec(),
422 envelope_before.into_boxed_slice(),
423 ),
424 Self::new(
425 after.to_owned(),
426 self.score,
427 self.scores[idx..].to_vec(),
428 envelope_after.into_boxed_slice(),
429 ),
430 )
431 } else {
432 let mut envelope_before = Vec::new();
433 let mut envelope_after = Vec::new();
434 for (env_before_i, env_after_i) in self
435 .envelope
436 .iter()
437 .map(|_| (MZPointSeries::default(), MZPointSeries::default()))
438 {
439 envelope_before.push(env_before_i);
440 envelope_after.push(env_after_i);
441 }
442 (
443 Self::new(
444 ChargedFeature::empty(self.charge()),
445 self.score,
446 Vec::new(),
447 envelope_before.into_boxed_slice(),
448 ),
449 Self::new(
450 ChargedFeature::empty(self.charge()),
451 self.score,
452 Vec::new(),
453 envelope_after.into_boxed_slice(),
454 ),
455 )
456 }
457 }
458
459 fn split_at(&self, point: usize) -> (Self::ViewType, Self::ViewType) {
460 let (before, after) = self.inner.split_at(point);
461 let mut envelope_before = Vec::new();
462 let mut envelope_after = Vec::new();
463 for (env_before_i, env_after_i) in self.envelope.iter().map(|e| {
464 let (a, b) = e.split_at(point);
465 (a.to_owned(), b.to_owned())
466 }) {
467 envelope_before.push(env_before_i);
468 envelope_after.push(env_after_i);
469 }
470 (
471 Self::new(
472 before.to_owned(),
473 self.score,
474 self.scores[..point].to_vec(),
475 envelope_before.into_boxed_slice(),
476 ),
477 Self::new(
478 after.to_owned(),
479 self.score,
480 self.scores[point..].to_vec(),
481 envelope_after.into_boxed_slice(),
482 ),
483 )
484 }
485
486 fn slice<I: std::ops::RangeBounds<usize> + Clone>(&self, bounds: I) -> Self::ViewType {
487 let inner = self.inner.slice(bounds.clone()).to_owned();
488 let envelope: Vec<_> = self
489 .envelope
490 .iter()
491 .map(|e| e.slice(bounds.clone()).to_owned())
492 .collect();
493
494 let start = match bounds.start_bound() {
495 Bound::Included(i) => *i,
496 Bound::Excluded(i) => *i,
497 Bound::Unbounded => 0,
498 };
499 let end = match bounds.end_bound() {
500 Bound::Included(i) => *i + 1,
501 Bound::Excluded(i) => *i,
502 Bound::Unbounded => self.scores.len(),
503 };
504
505 let scores = self.scores[start..end].to_vec();
506
507 Self::new(inner, self.score, scores, envelope.into_boxed_slice())
508 }
509}
510
511pub struct PeakIter<'a, Y: Clone> {
512 feature: &'a DeconvolvedSolutionFeature<Y>,
513 i: usize,
514}
515
516impl<'a, Y: Clone> PeakIter<'a, Y> {
517 pub fn new(feature: &'a DeconvolvedSolutionFeature<Y>) -> Self {
518 Self { feature, i: 0 }
519 }
520}
521
522impl<Y: Clone> Iterator for PeakIter<'_, Y> {
523 type Item = (DeconvolvedSolutionPeak, f64);
524
525 fn next(&mut self) -> Option<Self::Item> {
526 let i = self.i;
527 if i < self.feature.len() {
528 let (mass, time, inten) = self.feature.at(i).unwrap();
529 let score = self.feature.scores[i];
530 let mut env = Vec::with_capacity(self.feature.envelope.len());
531 env.extend(self.feature.envelope.iter().enumerate().map(|(j, e)| {
532 e.at(i).unwrap_or_else(|| {
533 let sizes: Vec<_> = self.feature.envelope.iter().map(|e| e.len()).collect();
534 warn!(
535 "Envelope {j} of length {} does not have an {i}th element, on a feature of length {} [{sizes:?}]",
536 e.len(),
537 self.feature.len(),
538 );
539 MZPoint::new(e.mz.first().copied().unwrap_or_default(), 1.0)
540 })
541 }));
542 let peak = DeconvolvedSolutionPeak::new(
543 mass,
544 inten,
545 self.feature.charge(),
546 0,
547 score,
548 Box::new(env),
549 );
550 self.i += 1;
551 Some((peak, time))
552 } else {
553 None
554 }
555 }
556}
557
558pub struct EnvelopeIter<'a, Y: Clone> {
559 feature: &'a DeconvolvedSolutionFeature<Y>,
560 i: usize,
561}
562
563impl<'a, Y: Clone> EnvelopeIter<'a, Y> {
564 pub fn new(feature: &'a DeconvolvedSolutionFeature<Y>) -> Self {
565 Self { feature, i: 0 }
566 }
567}
568
569impl<Y: Clone> Iterator for EnvelopeIter<'_, Y> {
570 type Item = (f64, Box<[(f64, f32)]>);
571
572 fn next(&mut self) -> Option<Self::Item> {
573 let i = self.i;
574 if i < self.feature.len() {
575 let (_, time, _) = self.feature.at(i).unwrap();
576 let env = self
577 .feature
578 .envelope
579 .iter()
580 .map(|e| {
581 let pt = e.at(i).unwrap();
582 (pt.mz, pt.intensity)
583 })
584 .collect();
585 self.i += 1;
586 Some((time, env))
587 } else {
588 None
589 }
590 }
591}
592
593pub(crate) fn reflow_feature<Y: Clone + Default>(
595 feature: DeconvolvedSolutionFeature<Y>,
596) -> DeconvolvedSolutionFeature<Y> {
597 let mut sink = DeconvolvedSolutionFeature::default();
598 *sink.charge_mut() = feature.charge();
599 sink.score = feature.score;
600
601 for (peak, time) in feature.iter_peaks() {
602 sink.push_peak(&peak, time);
603 }
604 sink
605}
606
607pub type FeatureMerger<Y> = ChargeAwareFeatureMerger<Mass, Y, DeconvolvedSolutionFeature<Y>>;
608
609const DECONVOLUTION_SCORE_ARRAY_NAME: &str = "deconvolution score array";
610const SUMMARY_SCORE_ARRAY_NAME: &str = "summary deconvolution score array";
611const ISOTOPIC_ENVELOPE_ARRAY_NAME: &str = "isotopic envelopes array";
612const FEATURE_IDENTIFIER_ARRAY_NAME: &str = "feature identifier array";
613
614impl BuildArrayMapFrom for DeconvolvedSolutionFeature<IonMobility> {
615 fn arrays_included(&self) -> Option<Vec<ArrayType>> {
616 Some(vec![
617 ArrayType::MZArray,
618 ArrayType::IntensityArray,
619 ArrayType::ChargeArray,
620 ArrayType::DeconvolutedIonMobilityArray,
621 ArrayType::nonstandard(SUMMARY_SCORE_ARRAY_NAME),
622 ArrayType::nonstandard(DECONVOLUTION_SCORE_ARRAY_NAME),
623 ArrayType::nonstandard(FEATURE_IDENTIFIER_ARRAY_NAME),
624 ])
625 }
626
627 fn as_arrays(source: &[Self]) -> BinaryArrayMap {
628 let m = source.len();
629 let n: usize = source.iter().map(|f| f.len()).sum();
630 let n_envelope = source
631 .iter()
632 .map(|f| f.envelope.iter().map(|e| e.len() + 2).sum::<usize>())
633 .sum::<usize>()
634 * 2;
635
636 let mut mz_array: Vec<u8> = Vec::with_capacity(n * BinaryDataArrayType::Float64.size_of());
637
638 let mut intensity_array: Vec<u8> =
639 Vec::with_capacity(n * BinaryDataArrayType::Float32.size_of());
640
641 let mut ion_mobility_array: Vec<u8> =
642 Vec::with_capacity(n * BinaryDataArrayType::Float64.size_of());
643
644 let mut charge_array: Vec<u8> =
645 Vec::with_capacity(n * BinaryDataArrayType::Int32.size_of());
646
647 let mut score_array: Vec<u8> =
648 Vec::with_capacity(n * BinaryDataArrayType::Float32.size_of());
649
650 let mut summary_score_array: Vec<u8> =
651 Vec::with_capacity(m * BinaryDataArrayType::Float32.size_of());
652
653 let mut marker_array: Vec<u8> =
654 Vec::with_capacity(n * BinaryDataArrayType::Int32.size_of());
655
656 let mut isotopic_envelope_array: Vec<u8> =
657 Vec::with_capacity(n_envelope * BinaryDataArrayType::Float32.size_of());
658
659 let mut acc = Vec::with_capacity(n);
660 source.iter().enumerate().for_each(|(i, f)| {
661 let envelope = &f.envelope;
662 let n_env = envelope.len() as f32;
663 let m_env = envelope.first().map(|f| f.len()).unwrap_or_default() as f32;
664
665 isotopic_envelope_array.extend_from_slice(&n_env.to_le_bytes());
666 isotopic_envelope_array.extend_from_slice(&m_env.to_le_bytes());
667 for env in envelope.iter() {
668 for (x, y) in env.iter() {
669 isotopic_envelope_array.extend_from_slice(&((*x) as f32).to_le_bytes());
670 isotopic_envelope_array.extend_from_slice(&(*y).to_le_bytes());
671 }
672 isotopic_envelope_array.extend_from_slice(&(0.0f32).to_le_bytes());
673 isotopic_envelope_array.extend_from_slice(&(0.0f32).to_le_bytes());
674 }
675
676 summary_score_array.extend(f.score.to_le_bytes());
677 f.iter().enumerate().for_each(|(j, (mass, im, inten))| {
678 acc.push((
679 mass_charge_ratio(mass, f.charge()),
680 im,
681 inten,
682 f.charge(),
683 f.scores[j],
684 i,
685 ))
686 })
687 });
688
689 acc.sort_by(
690 |(mz_a, im_a, _, _, _, key_a), (mz_b, im_b, _, _, _, key_b)| {
691 mz_a.total_cmp(mz_b)
692 .then(im_a.total_cmp(im_b))
693 .then(key_a.cmp(key_b))
694 },
695 );
696
697 for (mz, im, inten, charge, score, key) in acc.iter() {
698 mz_array.extend(mz.to_le_bytes());
699 intensity_array.extend(inten.to_le_bytes());
700 ion_mobility_array.extend(im.to_le_bytes());
701 charge_array.extend(charge.to_le_bytes());
702 score_array.extend(score.to_le_bytes());
703 marker_array.extend((*key as i32).to_le_bytes());
704 }
705
706 let mut map = BinaryArrayMap::default();
707 map.add(DataArray::wrap(
708 &ArrayType::MZArray,
709 BinaryDataArrayType::Float64,
710 mz_array,
711 ));
712 map.add(DataArray::wrap(
713 &ArrayType::IntensityArray,
714 BinaryDataArrayType::Float32,
715 intensity_array,
716 ));
717 map.add(DataArray::wrap(
718 &ArrayType::ChargeArray,
719 BinaryDataArrayType::Int32,
720 charge_array,
721 ));
722 map.add(DataArray::wrap(
723 &ArrayType::DeconvolutedIonMobilityArray,
724 BinaryDataArrayType::Float64,
725 ion_mobility_array,
726 ));
727 map.add(DataArray::wrap(
728 &ArrayType::nonstandard(SUMMARY_SCORE_ARRAY_NAME),
729 BinaryDataArrayType::Float32,
730 summary_score_array,
731 ));
732 map.add(DataArray::wrap(
733 &ArrayType::nonstandard(DECONVOLUTION_SCORE_ARRAY_NAME),
734 BinaryDataArrayType::Float32,
735 score_array,
736 ));
737 map.add(DataArray::wrap(
738 &ArrayType::nonstandard(FEATURE_IDENTIFIER_ARRAY_NAME),
739 BinaryDataArrayType::Int32,
740 marker_array,
741 ));
742 map.add(DataArray::wrap(
743 &ArrayType::nonstandard(ISOTOPIC_ENVELOPE_ARRAY_NAME),
744 BinaryDataArrayType::Float32,
745 isotopic_envelope_array,
746 ));
747
748 map
749 }
750}
751
752impl BuildArrayMap3DFrom for DeconvolvedSolutionFeature<IonMobility> {}
753
754impl BuildFromArrayMap for DeconvolvedSolutionFeature<IonMobility> {
755 fn arrays_required() -> Option<Vec<ArrayType>> {
756 Some(vec![
757 ArrayType::MZArray,
758 ArrayType::IntensityArray,
759 ArrayType::ChargeArray,
760 ArrayType::DeconvolutedIonMobilityArray,
761 ArrayType::nonstandard(SUMMARY_SCORE_ARRAY_NAME),
762 ArrayType::nonstandard(DECONVOLUTION_SCORE_ARRAY_NAME),
763 ArrayType::nonstandard(FEATURE_IDENTIFIER_ARRAY_NAME),
764 ])
765 }
766
767 fn try_from_arrays(arrays: &BinaryArrayMap) -> Result<Vec<Self>, ArrayRetrievalError> {
768 let arrays_3d = arrays.try_into()?;
769 Self::try_from_arrays_3d(&arrays_3d)
770 }
771}
772
773impl BuildFromArrayMap3D for DeconvolvedSolutionFeature<IonMobility> {
774 fn try_from_arrays_3d(arrays: &BinaryArrayMap3D) -> Result<Vec<Self>, ArrayRetrievalError> {
775 let key = ArrayType::nonstandard(FEATURE_IDENTIFIER_ARRAY_NAME);
776 let mut n: usize = 0;
777 for (_, arr) in arrays.iter() {
778 if arr.is_empty() {
779 continue;
780 }
781 if let Some(arr) = arr.get(&key) {
782 if let Some(i) = arr.iter_i32()?.map(|i| i as usize).max() {
783 n = n.max(i);
784 }
785 }
786 }
787
788 let isotopic_envelope_array_key = ArrayType::nonstandard(ISOTOPIC_ENVELOPE_ARRAY_NAME);
789 let score_array_key = ArrayType::nonstandard(DECONVOLUTION_SCORE_ARRAY_NAME);
790 let summary_score_array_key = ArrayType::nonstandard(SUMMARY_SCORE_ARRAY_NAME);
791
792 if n == 0 {
793 return Ok(Vec::new());
794 }
795
796 let mut index = Vec::with_capacity(n + 1);
797 index.resize(n + 1, Self::default());
798
799 for (im, arr) in arrays.iter() {
800 if arr.is_empty() {
801 continue;
802 }
803
804 let mz_array = arr.mzs()?;
805 let intensity_array = arr.intensities()?;
806 let charge_array = arr.charges()?;
807 let scores_array = arr
808 .get(&score_array_key)
809 .ok_or_else(|| ArrayRetrievalError::NotFound(score_array_key.clone()))?
810 .to_f32()?;
811 let marker_array = arr
812 .get(&key)
813 .ok_or_else(|| ArrayRetrievalError::NotFound(key.clone()))?
814 .to_i32()?;
815
816 for (mz, inten, charge, score, key_i) in multizip((
817 mz_array.iter(),
818 intensity_array.iter(),
819 charge_array.iter(),
820 scores_array.iter(),
821 marker_array.iter(),
822 )) {
823 let f = &mut index[(*key_i) as usize];
824 if f.is_empty() {
825 f.inner.charge = *charge;
826 }
827 f.score += *score;
828 f.push_raw(neutral_mass(*mz, *charge), im, *inten);
829 f.scores.push(*score);
830 }
831 }
832
833 if let Some(isotopic_envelopes_array) =
834 arrays.additional_arrays.get(&isotopic_envelope_array_key)
835 {
836 let mut isotopic_envelopes = Vec::new();
837 let isotopic_envelopes_array = isotopic_envelopes_array.to_f32()?;
838 let mut chunks = isotopic_envelopes_array.chunks_exact(2);
839 while let Some((n_traces, _trace_size)) = chunks.next().map(|header| {
840 let n_traces = header[0];
841 let trace_size = header[1];
842 (n_traces as usize, trace_size as usize)
843 }) {
844 let mut traces: Vec<MZPointSeries> = Vec::with_capacity(n_traces);
845 let mut current_trace = MZPointSeries::default();
846 while traces.len() < n_traces {
847 for block in chunks.by_ref() {
848 let mz = block[0] as f64;
849 let intensity = block[1];
850 if mz.is_zero() && intensity.is_zero() {
851 break;
852 }
853 current_trace.push(MZPoint::new(mz, intensity));
854 }
855 if current_trace.is_empty() {
856 warn!("Empty trace detected");
857 }
858 traces.push(current_trace);
859 current_trace = MZPointSeries::default();
860 }
861 isotopic_envelopes.push(Some(traces));
862 }
863 for (i, f) in index.iter_mut().enumerate() {
864 f.envelope = isotopic_envelopes[i].take().unwrap().into_boxed_slice();
865 }
866 }
867
868 if let Some(summary_score_array) = arrays.additional_arrays.get(&summary_score_array_key) {
869 let summary_scores = summary_score_array.to_f32()?;
870 for (score, f) in summary_scores.iter().zip(index.iter_mut()) {
871 f.score = *score;
872 }
873 }
874 Ok(index)
875 }
876}