opentimstdf 1.2.6

//! mzML export for Bruker timsTOF `.d/` (TDF) bundles.
//!
//! Provides:
//!
//! * [`TdfSource`] - a [`openmassspec_core::SpectrumSource`] adapter over an
//!   open [`Reader`]. Use this when you want to feed timsTOF data into any
//!   `openmassspec-core`-shaped consumer (column store ingest, Arrow bridge,
//!   ...).
//! * [`write_mzml`] / [`write_indexed_mzml`] - convenience entry points
//!   that wrap the canonical writer in `openmassspec-core`.
//!
//! Frame -> spectrum projection:
//!
//! * **MS1 frames** (`msms_type == 0`): all mobility scans pooled into one
//!   MS1 spectrum, sorted by m/z. Native ID: `frame=F scan=1`.
//! * **PASEF DDA frames** (`msms_type == 8`): one MS2 spectrum per
//!   `PasefMsMsInfo` row, filtered to that row's mobility scan range, with
//!   precursor m/z / charge / intensity copied from the `Precursors` table.
//!   Native ID: `frame=F scan=A-B`.
//! * **diaPASEF frames** (`msms_type == 9`): one MS2 spectrum per
//!   `DiaWindow` in the frame, with the isolation window taken from the
//!   `DiaFrameMsMsWindows` table. Native ID: `frame=F scan=A-B`.
//!
//! `scan_number` is a monotonic 1-based counter across the whole bundle
//! (PASEF frames produce many spectra per frame, so frame ID alone is not
//! unique).
//!
//! Frames with `msms_type` other than 0/8/9 (e.g. PRM-PASEF = 10) are
//! skipped for now. `iter_spectra` decodes one frame at a time as the
//! iterator is driven; a frame that fails to decode is skipped rather than
//! aborting the whole run, per [`msc::SpectrumSource::iter_spectra`]'s
//! "skip silently" contract - the canonical writer trusts whatever the
//! iterator yields.

use std::io::Write;
use std::path::Path;

use openmassspec_core as msc;

use crate::error::Result;
use crate::{
    Calibration, DiaWindow, Frame, Metadata, PasefMsMsInfo, Peak, Precursor as TdfPrecursor, Reader,
};

const SOFTWARE_NAME: &str = "opentimstdf";
const SOFTWARE_VERSION: &str = env!("CARGO_PKG_VERSION");

fn source_file_format_cv() -> msc::CvTerm {
    // PSI-MS MS:1002817 = "Bruker TDF format".
    msc::CvTerm::new("MS:1002817", "Bruker TDF format")
}

fn native_id_format_cv() -> msc::CvTerm {
    // PSI-MS MS:1002818 = "Bruker TDF nativeID format".
    msc::CvTerm::new("MS:1002818", "Bruker TDF nativeID format")
}

/// Resolve a PSI-MS instrument CV term from the GlobalMetadata
/// `InstrumentName`. Falls back to the generic Bruker term when the model
/// string is not in the lookup table.
fn instrument_cv(meta: &Metadata) -> msc::CvTerm {
    let name = meta.instrument_name.as_str();
    let known: &[(&str, &str, &str)] = &[
        ("timsTOF SCP", "MS:1003229", "timsTOF SCP"),
        ("timsTOF HT", "MS:1003404", "timsTOF HT"),
        ("timsTOF Pro 2", "MS:1003230", "timsTOF Pro 2"),
        ("timsTOF Pro", "MS:1003005", "timsTOF Pro"),
        ("timsTOF fleX", "MS:1003124", "timsTOF fleX"),
        ("timsTOF", "MS:1003005", "timsTOF"),
        ("impact II", "MS:1002280", "impact II"),
        ("impact", "MS:1001581", "Bruker Daltonics impact series"),
        ("maXis II", "MS:1002281", "maXis II"),
        ("maXis", "MS:1001541", "Bruker Daltonics maXis series"),
    ];
    for (prefix, acc, term_name) in known {
        if name.starts_with(prefix) {
            return msc::CvTerm::new(acc, *term_name);
        }
    }
    msc::CvTerm::new("MS:1000122", "Bruker Daltonics instrument model")
}

/// Minimal structural check for `YYYY-MM-DDTHH:MM:SS[.fraction](Z|+HH:MM|-HH:MM)`.
///
/// Not a full RFC 3339 parser (no calendar/range validation) - just enough
/// to catch the one divergence that actually matters here: Bruker's
/// `AcquisitionDateTime` is documented as "ISO 8601", and ISO 8601 permits a
/// local time with no UTC offset at all, whereas RFC 3339 (what
/// [`msc::RunMetadata::start_timestamp`] promises) always requires one.
/// Every real `.d/` bundle inspected during this audit (schema versions
/// spanning multiple instruments/years) included a numeric offset, so this
/// is expected to pass in practice; the check exists so an unusual bundle
/// that omits the offset degrades to `None` instead of silently mislabeling
/// a local time as RFC 3339.
fn is_rfc3339(s: &str) -> bool {
    let b = s.as_bytes();
    if b.len() < 20 {
        return false;
    }
    let digits = |r: std::ops::Range<usize>| b[r].iter().all(u8::is_ascii_digit);
    if !(digits(0..4) && b[4] == b'-' && digits(5..7) && b[7] == b'-' && digits(8..10)) {
        return false;
    }
    if b[10] != b'T' && b[10] != b't' {
        return false;
    }
    if !(digits(11..13) && b[13] == b':' && digits(14..16) && b[16] == b':' && digits(17..19)) {
        return false;
    }
    let mut rest = &s[19..];
    if let Some(frac) = rest.strip_prefix('.') {
        let end = frac
            .find(|c: char| !c.is_ascii_digit())
            .unwrap_or(frac.len());
        if end == 0 {
            return false;
        }
        rest = &frac[end..];
    }
    if rest.eq_ignore_ascii_case("z") {
        return true;
    }
    let rb = rest.as_bytes();
    rb.len() == 6
        && (rb[0] == b'+' || rb[0] == b'-')
        && rb[1].is_ascii_digit()
        && rb[2].is_ascii_digit()
        && rb[3] == b':'
        && rb[4].is_ascii_digit()
        && rb[5].is_ascii_digit()
}

fn polarity_for(frame: &Frame) -> Option<msc::Polarity> {
    match frame.mz_calibration_id {
        1 => Some(msc::Polarity::Positive),
        2 => Some(msc::Polarity::Negative),
        _ => None,
    }
}

/// Project a slice of decoded peaks into the column arrays a
/// [`msc::SpectrumRecord`] expects, optionally filtered by mobility scan
/// range `[scan_lo, scan_hi)` (used for PASEF / diaPASEF sub-windows).
struct PeakArrays {
    mz: Vec<f64>,
    intensity: Vec<f32>,
    tic: f64,
    base_peak_mz: f64,
    base_peak_intensity: f64,
    low_mz: f64,
    high_mz: f64,
    inv_mobility: Option<f64>,
}

fn materialize_peaks(
    peaks: &[Peak],
    cal: &Calibration,
    scan_lo: Option<u32>,
    scan_hi: Option<u32>,
) -> Option<PeakArrays> {
    let mut filtered: Vec<(f64, f32, u32)> = Vec::new();
    for p in peaks {
        if let Some(lo) = scan_lo {
            if p.scan < lo {
                continue;
            }
        }
        if let Some(hi) = scan_hi {
            if p.scan >= hi {
                continue;
            }
        }
        let mz = cal.tof_to_mz(p.tof);
        filtered.push((mz, p.intensity as f32, p.scan));
    }
    if filtered.is_empty() {
        return None;
    }
    filtered.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));

    let mut mz = Vec::with_capacity(filtered.len());
    let mut intensity = Vec::with_capacity(filtered.len());
    let mut tic: f64 = 0.0;
    let mut bp_mz = filtered[0].0;
    let mut bp_int: f32 = 0.0;
    let mut scan_sum: u64 = 0;
    for (m, i, s) in &filtered {
        mz.push(*m);
        intensity.push(*i);
        tic += *i as f64;
        if *i > bp_int {
            bp_int = *i;
            bp_mz = *m;
        }
        scan_sum += *s as u64;
    }
    let mean_scan = scan_sum as f64 / filtered.len() as f64;
    let inv_mob = cal.scan_to_inv_mobility(mean_scan.round() as u32);

    let low_mz = filtered.first().map(|t| t.0).unwrap_or(0.0);
    let high_mz = filtered.last().map(|t| t.0).unwrap_or(0.0);

    Some(PeakArrays {
        mz,
        intensity,
        tic,
        base_peak_mz: bp_mz,
        base_peak_intensity: bp_int as f64,
        low_mz,
        high_mz,
        inv_mobility: Some(inv_mob),
    })
}

fn build_ms1(
    scan_number: u32,
    frame: &Frame,
    peaks: &[Peak],
    cal: &Calibration,
) -> msc::SpectrumRecord {
    let pa = materialize_peaks(peaks, cal, None, None).unwrap_or(PeakArrays {
        mz: Vec::new(),
        intensity: Vec::new(),
        tic: 0.0,
        base_peak_mz: 0.0,
        base_peak_intensity: 0.0,
        low_mz: 0.0,
        high_mz: 0.0,
        inv_mobility: None,
    });
    msc::SpectrumRecord {
        index: (scan_number as usize).saturating_sub(1),
        scan_number,
        native_id: format!("frame={} scan=1", frame.id),
        ms_level: 1,
        polarity: polarity_for(frame),
        scan_mode: Some(msc::ScanMode::Centroid),
        analyzer: Some(msc::Analyzer::TOFMS),
        filter: None,
        retention_time_sec: frame.time,
        total_ion_current: Some(pa.tic),
        base_peak_mz: Some(pa.base_peak_mz),
        base_peak_intensity: Some(pa.base_peak_intensity),
        low_mz: Some(pa.low_mz),
        high_mz: Some(pa.high_mz),
        ion_injection_time_ms: frame.accumulation_time,
        inv_mobility: pa.inv_mobility,
        precursor: None,
        mz: pa.mz,
        intensity: pa.intensity,
        inv_mobility_per_peak: None,
    }
}

fn build_pasef_ms2(
    scan_number: u32,
    frame: &Frame,
    info: &PasefMsMsInfo,
    tdf_prec: Option<&TdfPrecursor>,
    peaks: &[Peak],
    cal: &Calibration,
) -> msc::SpectrumRecord {
    let pa = materialize_peaks(
        peaks,
        cal,
        Some(info.scan_num_begin),
        Some(info.scan_num_end),
    )
    .unwrap_or(PeakArrays {
        mz: Vec::new(),
        intensity: Vec::new(),
        tic: 0.0,
        base_peak_mz: 0.0,
        base_peak_intensity: 0.0,
        low_mz: 0.0,
        high_mz: 0.0,
        inv_mobility: None,
    });
    let prec_mz = tdf_prec
        .and_then(|p| p.monoisotopic_mz)
        .or_else(|| tdf_prec.map(|p| p.largest_peak_mz));
    let precursor = Some(msc::PrecursorInfo {
        target_mz: Some(info.isolation_mz),
        selected_mz: prec_mz,
        isolation_width: Some(info.isolation_width),
        charge: tdf_prec.and_then(|p| p.charge).map(|c| c as i32),
        intensity: tdf_prec.map(|p| p.intensity),
        collision_energy: Some(info.collision_energy),
        ce_is_nce: false,
        precursor_native_id: None,
        activation: Some(msc::Activation::HCD),
        analyzer: Some(msc::Analyzer::TOFMS),
    });
    msc::SpectrumRecord {
        index: (scan_number as usize).saturating_sub(1),
        scan_number,
        native_id: format!(
            "frame={} scan={}-{}",
            frame.id, info.scan_num_begin, info.scan_num_end
        ),
        ms_level: 2,
        polarity: polarity_for(frame),
        scan_mode: Some(msc::ScanMode::Centroid),
        analyzer: Some(msc::Analyzer::TOFMS),
        filter: None,
        retention_time_sec: frame.time,
        total_ion_current: Some(pa.tic),
        base_peak_mz: Some(pa.base_peak_mz),
        base_peak_intensity: Some(pa.base_peak_intensity),
        low_mz: Some(pa.low_mz),
        high_mz: Some(pa.high_mz),
        ion_injection_time_ms: frame.accumulation_time,
        inv_mobility: pa.inv_mobility,
        precursor,
        mz: pa.mz,
        intensity: pa.intensity,
        inv_mobility_per_peak: None,
    }
}

fn build_dia_ms2(
    scan_number: u32,
    frame: &Frame,
    window: &DiaWindow,
    peaks: &[Peak],
    cal: &Calibration,
) -> msc::SpectrumRecord {
    let pa = materialize_peaks(
        peaks,
        cal,
        Some(window.scan_num_begin),
        Some(window.scan_num_end),
    )
    .unwrap_or(PeakArrays {
        mz: Vec::new(),
        intensity: Vec::new(),
        tic: 0.0,
        base_peak_mz: 0.0,
        base_peak_intensity: 0.0,
        low_mz: 0.0,
        high_mz: 0.0,
        inv_mobility: None,
    });
    let precursor = Some(msc::PrecursorInfo {
        target_mz: Some(window.isolation_mz),
        selected_mz: Some(window.isolation_mz),
        isolation_width: Some(window.isolation_width),
        charge: None,
        intensity: None,
        collision_energy: Some(window.collision_energy),
        ce_is_nce: false,
        precursor_native_id: None,
        activation: Some(msc::Activation::HCD),
        analyzer: Some(msc::Analyzer::TOFMS),
    });
    msc::SpectrumRecord {
        index: (scan_number as usize).saturating_sub(1),
        scan_number,
        native_id: format!(
            "frame={} scan={}-{}",
            frame.id, window.scan_num_begin, window.scan_num_end
        ),
        ms_level: 2,
        polarity: polarity_for(frame),
        scan_mode: Some(msc::ScanMode::Centroid),
        analyzer: Some(msc::Analyzer::TOFMS),
        filter: None,
        retention_time_sec: frame.time,
        total_ion_current: Some(pa.tic),
        base_peak_mz: Some(pa.base_peak_mz),
        base_peak_intensity: Some(pa.base_peak_intensity),
        low_mz: Some(pa.low_mz),
        high_mz: Some(pa.high_mz),
        ion_injection_time_ms: frame.accumulation_time,
        inv_mobility: pa.inv_mobility,
        precursor,
        mz: pa.mz,
        intensity: pa.intensity,
        inv_mobility_per_peak: None,
    }
}

/// `SpectrumSource` adapter over an open [`Reader`].
///
/// Construct with [`TdfSource::new`]; iterate via the trait's
/// [`iter_spectra`](openmassspec_core::SpectrumSource::iter_spectra) method.
///
/// `iter_spectra` decodes and projects one frame at a time as the returned
/// iterator is driven, so memory use is bounded by a single frame's worth
/// of spectra (PASEF/diaPASEF frames fan out into several spectra each,
/// which are queued and drained before the next frame is decoded) rather
/// than the whole run.
pub struct TdfSource<'a> {
    reader: &'a Reader,
    bundle_name: String,
    metadata: Metadata,
    calibration: Calibration,
    frames: Vec<Frame>,
}

impl<'a> TdfSource<'a> {
    /// Build a new source from an open [`Reader`] and the bundle's directory
    /// name (used as the `<sourceFile name="...">` in mzML output).
    pub fn new(reader: &'a Reader, bundle_name: impl Into<String>) -> Result<Self> {
        let metadata = reader.metadata()?;
        let calibration = reader.calibration()?;
        let frames = reader.frames()?;
        Ok(Self {
            reader,
            bundle_name: bundle_name.into(),
            metadata,
            calibration,
            frames,
        })
    }

    /// Open a bundle directory and build a source in one call.
    ///
    /// The reader is owned by the returned `OwnedTdfSource`.
    pub fn open<P: AsRef<Path>>(path: P) -> Result<OwnedTdfSource> {
        let path = path.as_ref();
        let reader = Reader::open(path)?;
        let metadata = reader.metadata()?;
        let calibration = reader.calibration()?;
        let frames = reader.frames()?;
        let bundle_name = path
            .file_name()
            .map(|n| n.to_string_lossy().into_owned())
            .unwrap_or_else(|| "bundle.d".into());
        Ok(OwnedTdfSource {
            reader,
            bundle_name,
            metadata,
            calibration,
            frames,
        })
    }
}

/// Self-contained variant of [`TdfSource`] that owns its [`Reader`]. Used
/// for the common case where the caller just wants `open() -> write_mzml`.
pub struct OwnedTdfSource {
    reader: Reader,
    bundle_name: String,
    metadata: Metadata,
    calibration: Calibration,
    frames: Vec<Frame>,
}

/// Project one frame into zero or more spectra, incrementing `scan_counter`
/// for each spectrum produced. Any decode failure - the frame's peaks, its
/// PASEF info rows, or its diaPASEF windows - causes that frame to be
/// skipped (returns an empty `Vec`) rather than aborting the caller's
/// iteration, matching [`msc::SpectrumSource::iter_spectra`]'s silent-skip
/// contract. A precursor lookup failure for a single PASEF MS2 is narrower:
/// it only drops that spectrum's precursor metadata, not the spectrum.
fn spectra_for_frame(
    reader: &Reader,
    frame: &Frame,
    calibration: &Calibration,
    scan_counter: &mut u32,
) -> Vec<msc::SpectrumRecord> {
    let Ok(peaks) = reader.decode_peaks(frame) else {
        return Vec::new();
    };
    match frame.msms_type {
        0 => {
            *scan_counter += 1;
            vec![build_ms1(*scan_counter, frame, &peaks, calibration)]
        }
        8 => {
            let Ok(infos) = reader.pasef_msms_info_for_frame(frame.id) else {
                return Vec::new();
            };
            let mut out = Vec::with_capacity(infos.len());
            for info in infos {
                let prec = reader.precursor(info.precursor_id).ok().flatten();
                *scan_counter += 1;
                out.push(build_pasef_ms2(
                    *scan_counter,
                    frame,
                    &info,
                    prec.as_ref(),
                    &peaks,
                    calibration,
                ));
            }
            out
        }
        9 => {
            let Ok(windows) = reader.dia_windows_for_frame(frame.id) else {
                return Vec::new();
            };
            let mut out = Vec::new();
            if let Some(group) = windows {
                for w in &group.windows {
                    *scan_counter += 1;
                    out.push(build_dia_ms2(*scan_counter, frame, w, &peaks, calibration));
                }
            }
            out
        }
        _ => Vec::new(),
    }
}

/// Build a lazy, frame-at-a-time spectrum iterator. Decodes and projects
/// one frame per call to `next()` that yields nothing from the pending
/// queue, so at most one frame's peaks (plus its fan-out of PASEF/diaPASEF
/// spectra) are held in memory at a time.
fn frame_iter<'s>(
    reader: &'s Reader,
    frames: &'s [Frame],
    calibration: &'s Calibration,
) -> impl Iterator<Item = msc::SpectrumRecord> + 's {
    let mut frame_idx = 0usize;
    let mut pending = std::collections::VecDeque::new();
    let mut scan_counter: u32 = 0;
    std::iter::from_fn(move || loop {
        if let Some(rec) = pending.pop_front() {
            return Some(rec);
        }
        let frame = frames.get(frame_idx)?;
        frame_idx += 1;
        pending.extend(spectra_for_frame(
            reader,
            frame,
            calibration,
            &mut scan_counter,
        ));
    })
}

fn run_metadata_for(meta: &Metadata, bundle_name: &str) -> msc::RunMetadata {
    msc::RunMetadata {
        source_file_name: bundle_name.to_string(),
        source_file_format: source_file_format_cv(),
        native_id_format: native_id_format_cv(),
        instrument: instrument_cv(meta),
        software_name: SOFTWARE_NAME.into(),
        software_version: SOFTWARE_VERSION.into(),
        // `Metadata::acquisition_date_time` is read straight from the
        // `AcquisitionDateTime` row in `analysis.tdf`'s GlobalMetadata table
        // (see `Reader::metadata` in reader.rs) and documented there as an
        // ISO 8601 string. `RunMetadata::start_timestamp` promises RFC 3339,
        // which additionally requires a UTC offset; `is_rfc3339` guards
        // against a bundle whose value omits one rather than passing it
        // through unchecked.
        start_timestamp: meta
            .acquisition_date_time
            .as_deref()
            .filter(|s| is_rfc3339(s))
            .map(str::to_string),
        mobility_array_kind: Some(msc::MobilityArrayKind::InverseReducedVsPerCm2),
    }
}

impl<'a> msc::SpectrumSource for TdfSource<'a> {
    fn run_metadata(&self) -> msc::RunMetadata {
        run_metadata_for(&self.metadata, &self.bundle_name)
    }
    fn iter_spectra<'s>(&'s mut self) -> Box<dyn Iterator<Item = msc::SpectrumRecord> + 's> {
        Box::new(frame_iter(self.reader, &self.frames, &self.calibration))
    }
}

impl msc::SpectrumSource for OwnedTdfSource {
    fn run_metadata(&self) -> msc::RunMetadata {
        run_metadata_for(&self.metadata, &self.bundle_name)
    }
    fn iter_spectra<'s>(&'s mut self) -> Box<dyn Iterator<Item = msc::SpectrumRecord> + 's> {
        Box::new(frame_iter(&self.reader, &self.frames, &self.calibration))
    }
}

/// Write a `.d/` bundle to mzML in one call.
///
/// Convenience wrapper that opens `bundle_dir`, projects every frame into
/// the appropriate MS1 / PASEF MS2 / diaPASEF MS2 spectra, and emits a
/// valid mzML 1.1.0 document via the canonical writer in `openmassspec-core`.
pub fn write_mzml<P: AsRef<Path>, W: Write>(bundle_dir: P, out: &mut W) -> Result<()> {
    let mut src = TdfSource::open(bundle_dir)?;
    msc::write_mzml(&mut src, out)?;
    Ok(())
}

/// Indexed-mzML equivalent of [`write_mzml`].
pub fn write_indexed_mzml<P: AsRef<Path>, W: Write>(bundle_dir: P, out: &mut W) -> Result<()> {
    let mut src = TdfSource::open(bundle_dir)?;
    msc::write_indexed_mzml(&mut src, out)?;
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;

    fn sample_metadata(acquisition_date_time: Option<&str>) -> Metadata {
        Metadata {
            schema_version_major: 3,
            schema_version_minor: 7,
            instrument_name: "timsTOF Pro".into(),
            acquisition_software: "timsControl".into(),
            acquisition_software_version: "2.0.18".into(),
            compression_type: 2,
            acquisition_date_time: acquisition_date_time.map(str::to_string),
        }
    }

    #[test]
    fn is_rfc3339_accepts_real_bruker_values() {
        // Pulled from GlobalMetadata.AcquisitionDateTime in real .d/ bundles
        // (multiple instruments/schema versions) - all include a numeric
        // UTC offset.
        for s in [
            "2018-08-21T20:40:14.356+02:00",
            "2022-02-22T23:47:02.147-08:00",
            "2026-03-30T22:54:48.394-07:00",
            "2019-08-24T18:27:02.345+02:00",
            "2026-01-01T00:00:00Z",
        ] {
            assert!(is_rfc3339(s), "expected {s} to be accepted");
        }
    }

    #[test]
    fn is_rfc3339_rejects_missing_offset_and_garbage() {
        for s in [
            "2019-01-17T09:14:39.730", // ISO 8601 local time, no offset
            "2019-01-17 09:14:39",     // space separator, no offset
            "not-a-timestamp",
            "",
        ] {
            assert!(!is_rfc3339(s), "expected {s:?} to be rejected");
        }
    }

    #[test]
    fn run_metadata_for_wires_up_start_timestamp() {
        let meta = sample_metadata(Some("2018-08-21T20:40:14.356+02:00"));
        let rm = run_metadata_for(&meta, "bundle.d");
        assert_eq!(
            rm.start_timestamp.as_deref(),
            Some("2018-08-21T20:40:14.356+02:00")
        );
    }

    #[test]
    fn run_metadata_for_none_when_absent() {
        let meta = sample_metadata(None);
        let rm = run_metadata_for(&meta, "bundle.d");
        assert_eq!(rm.start_timestamp, None);
    }

    #[test]
    fn run_metadata_for_none_when_not_rfc3339() {
        // Defensive path: don't claim RFC 3339 compliance for a value that
        // isn't, even though no real-world bundle observed so far hits this.
        let meta = sample_metadata(Some("2019-01-17T09:14:39.730"));
        let rm = run_metadata_for(&meta, "bundle.d");
        assert_eq!(rm.start_timestamp, None);
    }
}