corems.mass_spectra.calc.lc_calc

   1import math
   2import time
   3from itertools import product
   4
   5import numpy as np
   6import pandas as pd
   7import warnings
   8import scipy
   9
  10if not hasattr(np, "trapezoid"):  # numpy < 2.0
  11    np.trapezoid = np.trapz
  12from ripser import ripser
  13from scipy import sparse
  14from scipy.spatial import KDTree
  15from sklearn.svm import SVR
  16from sklearn.cluster import AgglomerativeClustering
  17import matplotlib.pyplot as plt
  18from tqdm import tqdm
  19
  20from corems.chroma_peak.factory.chroma_peak_classes import LCMSMassFeature
  21from corems.encapsulation.plot_utils import _finalize_plot
  22from corems.encapsulation.spawn_pool import spawn_pool
  23from corems.mass_spectra.calc import SignalProcessing as sp
  24from corems.mass_spectra.factory.chromat_data import EIC_Data
  25from corems.mass_spectrum.input.numpyArray import ms_from_array_profile
  26
  27warnings.filterwarnings("ignore", category=RuntimeWarning)
  28
  29
  30def find_closest(A, target):
  31    """Find the index of closest value in A to each value in target.
  32
  33    Parameters
  34    ----------
  35    A : :obj:`~numpy.array`
  36        The array to search (blueprint). A must be sorted.
  37    target : :obj:`~numpy.array`
  38        The array of values to search for. target must be sorted.
  39
  40    Returns
  41    -------
  42    :obj:`~numpy.array`
  43        The indices of the closest values in A to each value in target.
  44    """
  45    idx = A.searchsorted(target)
  46    idx = np.clip(idx, 1, len(A) - 1)
  47    left = A[idx - 1]
  48    right = A[idx]
  49    idx -= target - left < right - target
  50    return idx
  51
  52
  53def _mz_grid_centers(mz_desc, radius):
  54    """m/z bin centers, highest intensity first.
  55
  56    ``mz_desc`` is unique and ordered from highest intensity to lowest. A value
  57    is kept when no kept higher-intensity value lies within ``radius``. That is
  58    the parent rule of the distance-pair grid, without building the pair matrix.
  59    """
  60    mz = np.asarray(mz_desc, dtype=np.float64)
  61    n = mz.size
  62    if n == 0:
  63        return mz
  64    rank_order = np.argsort(mz, kind="mergesort")
  65    mz_by_rank = mz[rank_order]
  66    rank_of = np.empty(n, dtype=np.intp)
  67    rank_of[rank_order] = np.arange(n)
  68    bit = np.zeros(n + 1, dtype=np.int32)
  69    keep = np.zeros(n, dtype=bool)
  70
  71    def prefix(end):
  72        total = 0
  73        while end:
  74            total += int(bit[end])
  75            end -= end & -end
  76        return total
  77
  78    def mark(rank):
  79        node = rank + 1
  80        while node <= n:
  81            bit[node] += 1
  82            node += node & -node
  83
  84    for i, value in enumerate(mz):
  85        lo = int(np.searchsorted(mz_by_rank, value - radius, side="left"))
  86        hi = int(np.searchsorted(mz_by_rank, value + radius, side="right"))
  87        if prefix(hi) - prefix(lo) > 0:
  88            continue
  89        mark(int(rank_of[i]))
  90        keep[i] = True
  91    return mz[keep]
  92
  93
  94class LCCalculations:
  95    """Methods for performing LC calculations on mass spectra data.
  96
  97    Notes
  98    -----
  99    This class is intended to be used as a mixin for the LCMSBase class.
 100
 101    Methods
 102    -------
 103    * get_max_eic(eic_data).
 104        Returns the maximum EIC value from the given EIC data. A static method.
 105    * smooth_tic(tic).
 106        Smooths the TIC data using the specified smoothing method and settings.
 107    * eic_centroid_detector(rt, eic, max_eic).
 108        Performs EIC centroid detection on the given EIC data.
 109    * find_nearest_scan(rt).
 110        Finds the nearest scan to the given retention time.
 111    * get_average_mass_spectrum(scan_list, apex_scan, spectrum_mode="profile", ms_level=1, auto_process=True, use_parser=False, perform_checks=True, polarity=None).
 112        Returns an averaged mass spectrum object.
 113    * find_mass_features(ms_level=1).
 114        Find regions of interest for a given MS level (default is MS1).
 115    * integrate_mass_features(drop_if_fail=False, ms_level=1).
 116        Integrate mass features of interest and extracts EICs.
 117    * find_c13_mass_features().
 118        Evaluate mass features and mark likely C13 isotopes.
 119    * deconvolute_ms1_mass_features().
 120        Deconvolute mass features' ms1 mass spectra.
 121    """
 122
 123    @staticmethod
 124    def get_max_eic(eic_data: dict):
 125        """Returns the maximum EIC value from the given EIC data.
 126
 127        Notes
 128        -----
 129        This is a static method.
 130
 131        Parameters
 132        ----------
 133        eic_data : dict
 134            A dictionary containing EIC data.
 135
 136        Returns
 137        -------
 138        float
 139            The maximum EIC value.
 140        """
 141        max_eic = 0
 142        for eic_data in eic_data.values():
 143            ind_max_eic = max(eic_data.get("EIC"))
 144            max_eic = ind_max_eic if ind_max_eic > max_eic else max_eic
 145
 146        return max_eic
 147
 148    def smooth_tic(self, tic):
 149        """Smooths the TIC or EIC data using the specified smoothing method and settings.
 150
 151        Parameters
 152        ----------
 153        tic : numpy.ndarray
 154            The TIC (or EIC) data to be smoothed.
 155
 156        Returns
 157        -------
 158        numpy.ndarray
 159            The smoothed TIC data.
 160        """
 161        implemented_smooth_method = self.parameters.lc_ms.implemented_smooth_method
 162
 163        pol_order = self.parameters.lc_ms.savgol_pol_order
 164
 165        window_len = self.parameters.lc_ms.smooth_window
 166
 167        window = self.parameters.lc_ms.smooth_method
 168
 169        return sp.smooth_signal(
 170            tic, window_len, window, pol_order, implemented_smooth_method
 171        )
 172
 173    def eic_centroid_detector(self, rt, eic, max_eic, apex_indexes=[]):
 174        """Performs EIC centroid detection on the given EIC data.
 175
 176        Parameters
 177        ----------
 178        rt : numpy.ndarray
 179            The retention time data.
 180        eic : numpy.ndarray
 181            The EIC data.
 182        max_eic : float
 183            The maximum EIC value.
 184        apex_indexes : list, optional
 185            The apexes of the EIC peaks. Defaults to [], which means that the apexes will be calculated by the function.
 186
 187        Returns
 188        -------
 189        numpy.ndarray
 190            The indexes of left, apex, and right limits as a generator.
 191        """
 192
 193        max_prominence = self.parameters.lc_ms.peak_max_prominence_percent
 194
 195        max_height = self.parameters.lc_ms.peak_height_max_percent
 196
 197        signal_threshold = self.parameters.lc_ms.eic_signal_threshold
 198
 199        min_peak_datapoints = self.parameters.lc_ms.min_peak_datapoints
 200
 201        peak_derivative_threshold = self.parameters.lc_ms.peak_derivative_threshold
 202
 203        include_indexes = sp.peak_picking_first_derivative(
 204            domain=rt,
 205            signal=eic,
 206            max_height=max_height,
 207            max_prominence=max_prominence,
 208            max_signal=max_eic,
 209            min_peak_datapoints=min_peak_datapoints,
 210            peak_derivative_threshold=peak_derivative_threshold,
 211            signal_threshold=signal_threshold,
 212            correct_baseline=False,
 213            plot_res=False,
 214            apex_indexes=apex_indexes,
 215        )
 216        # include_indexes is a generator of tuples (left_index, apex_index, right_index)
 217        include_indexes = list(include_indexes)
 218        # Add check to make sure that there are at least 1/2 of min_peak_datapoints on either side of the apex
 219        indicies = [x for x in include_indexes]
 220        for idx in indicies:
 221            if (idx[1] - idx[0] < min_peak_datapoints / 2) or (
 222                idx[2] - idx[1] < min_peak_datapoints / 2
 223            ):
 224                include_indexes.remove(idx)
 225        return include_indexes
 226
 227    def find_nearest_scan(self, rt):
 228        """Finds the nearest scan to the given retention time.
 229
 230        Parameters
 231        ----------
 232        rt : float
 233            The retention time (in minutes) to find the nearest scan for.
 234
 235        Returns
 236        -------
 237        int
 238            The scan number of the nearest scan.
 239        """
 240        array_rt = np.array(self.retention_time)
 241
 242        scan_index = (np.abs(array_rt - rt)).argmin()
 243
 244        real_scan = self.scans_number[scan_index]
 245
 246        return real_scan
 247
 248    def add_peak_metrics(self, remove_by_metrics=True, induced_features=False):
 249        """Add peak metrics to the mass features.
 250
 251        This function calculates the peak metrics for each mass feature and adds them to the mass feature objects.
 252
 253        Parameters
 254        ----------
 255        remove_by_metrics : bool, optional
 256            If True, remove mass features based on their peak metrics such as S/N, Gaussian similarity,
 257            dispersity index, and noise score. Default is True, which checks the setting in the processing parameters.
 258            If False, peak metrics are calculated but no mass features are removed, regardless of the setting in the processing parameters.
 259        induced_features : bool, optional
 260            Whether the mass features to be integrated were induced. Default is False.
 261        """
 262        # Check that at least some mass features have eic data
 263        if induced_features:
 264            mf_dict_values = self.induced_mass_features.values()
 265        else:
 266            mf_dict_values = self.mass_features.values()
 267
 268        if not any([mf._eic_data is not None for mf in mf_dict_values]):
 269            raise ValueError(
 270                "No mass features have EIC data. Run integrate_mass_features first."
 271            )
 272
 273        for mass_feature in mf_dict_values:
 274            # Check if the mass feature has been integrated
 275            if mass_feature._eic_data is not None and mass_feature.area is not None:
 276                # Calculate peak metrics
 277                mass_feature.calc_half_height_width()
 278                mass_feature.calc_tailing_factor()
 279                mass_feature.calc_dispersity_index()
 280                mass_feature.calc_gaussian_similarity()
 281                mass_feature.calc_noise_score()
 282
 283        # Remove mass features by peak metrics if designated in parameters
 284        if (
 285            self.parameters.lc_ms.remove_mass_features_by_peak_metrics
 286            and remove_by_metrics
 287        ):
 288            self._remove_mass_features_by_peak_metrics(
 289                induced_features=induced_features
 290            )
 291
 292    def get_average_mass_spectrum(
 293        self,
 294        scan_list,
 295        apex_scan,
 296        spectrum_mode="profile",
 297        ms_level=1,
 298        auto_process=True,
 299        use_parser=False,
 300        perform_checks=True,
 301        polarity=None,
 302        ms_params=None,
 303    ):
 304        """Returns an averaged mass spectrum object
 305
 306        Parameters
 307        ----------
 308        scan_list : list
 309            List of scan numbers to average.
 310        apex_scan : int
 311            Number of the apex scan
 312        spectrum_mode : str, optional
 313            The spectrum mode to use. Defaults to "profile". Not that only "profile" mode is supported for averaging.
 314        ms_level : int, optional
 315            The MS level to use. Defaults to 1.
 316        auto_process : bool, optional
 317            If True, the averaged mass spectrum will be auto-processed. Defaults to True.
 318        use_parser : bool, optional
 319            If True, the mass spectra will be obtained from the parser. Defaults to False.
 320        perform_checks : bool, optional
 321            If True, the function will check if the data are within the ms_unprocessed dictionary and are the correct mode. Defaults to True. Only set to False if you are sure the data are profile, and (if not using the parser) are in the ms_unprocessed dictionary!  ms_unprocessed dictionary also must be indexed on scan
 322        polarity : int, optional
 323            The polarity of the mass spectra (1 or -1). If not set, the polarity will be determined from the dataset. Defaults to None. (fastest if set to -1 or 1)
 324        ms_params : MSParameters, optional
 325            The mass spectrum parameters to use. If not set (None), the globally set parameters will be used. Defaults to None.
 326
 327        Returns
 328        -------
 329        MassSpectrumProfile
 330            The averaged mass spectrum object.
 331
 332        Raises
 333        ------
 334        ValueError
 335            If the spectrum mode is not "profile".
 336            If the MS level is not found in the unprocessed mass spectra dictionary.
 337            If not all scan numbers are found in the unprocessed mass spectra dictionary.
 338        """
 339        if perform_checks:
 340            if spectrum_mode != "profile":
 341                raise ValueError("Averaging only supported for profile mode")
 342
 343        if polarity is None:
 344            # set polarity to -1 if negative mode, 1 if positive mode (for mass spectrum creation)
 345            if self.polarity == "negative":
 346                polarity = -1
 347            elif self.polarity == "positive":
 348                polarity = 1
 349            else:
 350                raise ValueError(
 351                    "Polarity not set for dataset, must be a set containing either 'positive' or 'negative'"
 352                )
 353
 354        # if not using_parser, check that scan numbers are in _ms_unprocessed
 355        if not use_parser:
 356            if perform_checks:
 357                # Set index to scan for faster lookup
 358                ms_df = (
 359                    self._ms_unprocessed[ms_level]
 360                    .copy()
 361                    .set_index("scan", drop=False)
 362                    .sort_index()
 363                )
 364                my_ms_df = ms_df.loc[scan_list]
 365                # Check that all scan numbers are in the ms_df
 366                if not all(np.isin(scan_list, ms_df.index)):
 367                    raise ValueError(
 368                        "Not all scan numbers found in the unprocessed mass spectra dictionary"
 369                    )
 370            else:
 371                my_ms_df = (
 372                    pd.DataFrame({"scan": scan_list})
 373                    .set_index("scan")
 374                    .join(self._ms_unprocessed[ms_level], how="left")
 375                )
 376
 377        if use_parser:
 378            parser = self.spectra_parser
 379            ms_list = [
 380                parser.get_mass_spectrum_from_scan(
 381                    x, spectrum_mode=spectrum_mode, auto_process=False
 382                )
 383                for x in scan_list
 384            ]
 385            ms_mz = [x._mz_exp for x in ms_list]
 386            ms_int = [x._abundance for x in ms_list]
 387            my_ms_df = []
 388            for i in np.arange(len(ms_mz)):
 389                my_ms_df.append(
 390                    pd.DataFrame(
 391                        {"mz": ms_mz[i], "intensity": ms_int[i], "scan": scan_list[i]}
 392                    )
 393                )
 394            my_ms_df = pd.concat(my_ms_df)
 395
 396        if not self.check_if_grid(my_ms_df):
 397            my_ms_df = self.grid_data(my_ms_df)
 398
 399        my_ms_ave = my_ms_df.groupby("mz")["intensity"].sum().reset_index()
 400
 401        ms = ms_from_array_profile(
 402            my_ms_ave.mz,
 403            my_ms_ave.intensity,
 404            self.file_location,
 405            polarity=polarity,
 406            auto_process=False,
 407        )
 408
 409        # Set the mass spectrum parameters, auto-process if auto_process is True, and add to the dataset
 410        if ms is not None:
 411            if ms_params is not None:
 412                ms.parameters = ms_params
 413            ms.scan_number = apex_scan
 414            if auto_process:
 415                ms.process_mass_spec()
 416        return ms
 417
 418    def _resolve_lcms_peak_picking_method(self, ms_level=1):
 419        """Resolve and validate the 2D peak picking method against scan formats.
 420
 421        Parameters
 422        ----------
 423        ms_level : int, optional
 424            MS level whose ``scan_df`` formats are checked. Default is 1.
 425
 426        Returns
 427        -------
 428        str
 429            Concrete method: ``persistent homology`` or
 430            ``centroided_persistent_homology``. If the stored method is
 431            ``by_format``, it is chosen from scan formats (all profile or all
 432            centroid).
 433
 434        Raises
 435        ------
 436        ValueError
 437            If the method is not implemented, ``by_format`` cannot be resolved
 438            (empty, mixed, or unknown formats), or an explicit method does
 439            not match the scan formats at ``ms_level``.
 440        """
 441        method = self.parameters.lc_ms.peak_picking_method
 442        msx_scan_df = self.scan_df[self.scan_df["ms_level"] == ms_level]
 443        values = list(msx_scan_df["ms_format"])
 444
 445        if method == "by_format":
 446            if len(values) == 0:
 447                raise ValueError(
 448                    "Cannot resolve by_format peak picking: no scans found "
 449                    "for the requested MS level."
 450                )
 451            if all(v == "profile" for v in values):
 452                return "persistent homology"
 453            if all(v == "centroid" for v in values):
 454                return "centroided_persistent_homology"
 455            found = sorted({str(v) for v in values})
 456            raise ValueError(
 457                "Cannot resolve by_format peak picking: MS scans are mixed "
 458                f"or unknown format (found {found}). Set peak_picking_method "
 459                "explicitly to 'persistent homology' or "
 460                "'centroided_persistent_homology'."
 461            )
 462
 463        if method == "persistent homology":
 464            if len(values) == 0 or not all(v == "profile" for v in values):
 465                raise ValueError(
 466                    "MS{} scans are not profile mode, which is required for persistent homology peak picking.".format(
 467                        ms_level
 468                    )
 469                )
 470            return method
 471
 472        if method == "centroided_persistent_homology":
 473            if len(values) == 0 or not all(v == "centroid" for v in values):
 474                raise ValueError(
 475                    "MS{} scans are not centroid mode, which is required for persistent homology centroided peak picking.".format(
 476                        ms_level
 477                    )
 478                )
 479            return method
 480
 481        raise ValueError("Peak picking method not implemented")
 482
 483    def find_mass_features(
 484        self,
 485        ms_level=1,
 486        grid=True,
 487        assign_ms2_scans=False,
 488        ms2_scan_filter=None,
 489        targeted_search=False,
 490        target_search_dict=None,
 491        accumulate_features=False,
 492    ):
 493        """Find mass features within an LCMSBase object
 494
 495        Wrapper that dispatches to a 2D peak picking implementation.
 496
 497        If ``parameters.lc_ms.peak_picking_method`` is ``by_format``, the method is
 498        chosen from ``scan_df.ms_format`` at the requested MS level (all
 499        profile → persistent homology; all centroid → centroided persistent
 500        homology). A successful resolve overwrites
 501        ``peak_picking_method`` with the concrete method so saved parameters
 502        record what was used.
 503
 504        Parameters
 505        ----------
 506        ms_level : int, optional
 507            The MS level to use for peak picking Default is 1.
 508        grid : bool, optional
 509            If True, will regrid the data before running the persistent homology calculations (after checking if the data is gridded),
 510            used for persistent homology peak picking for profile data only. Default is True.
 511        assign_ms2_scans : bool, optional
 512            If True, assign MS2 scan numbers to mass features after peak picking.
 513            This populates the ms2_scan_numbers attribute on each mass feature, which enables
 514            choosing representative features based on MS2 availability. Default is False.
 515        ms2_scan_filter : str or None, optional
 516            Filter string for MS2 scans when assign_ms2_scans is True (e.g., 'hcd').
 517            If None, all MS2 scans are considered. Default is None.
 518        targeted_search : bool, optional
 519            If True, perform targeted mass feature search using the target_search_dict.
 520            This mode filters data to only m/z and RT windows of interest and bypasses
 521            intensity and persistence thresholds. Default is False.
 522        target_search_dict : dict or None, optional
 523            Dictionary containing target search parameters. Required if targeted_search is True.
 524            Must contain:
 525                - 'target_mz_list': list of target m/z values
 526                - 'target_rt_list': list of target retention times (in minutes)
 527                - 'mz_tolerance_ppm': m/z tolerance in ppm
 528                - 'rt_tolerance': retention time tolerance (in minutes)
 529            Optionally can contain:
 530                - 'type': type label for mass features (e.g., "internal standard")
 531                  If not provided, defaults to "targeted"
 532            Default is None.
 533        accumulate_features : bool, optional
 534            If True, new mass features will be added to existing features rather than replacing them.
 535            This allows multiple sequential calls to find_mass_features to build up a combined set.
 536            Default is False (replace existing features for backwards compatibility).
 537
 538        Raises
 539        ------
 540        ValueError
 541            If no MS level data is found on the object.
 542            If ``by_format`` cannot resolve a method (empty, mixed, or unknown
 543            ``ms_format`` values at the requested MS level).
 544            If persistent homology peak picking is attempted on non-profile mode data.
 545            If centroided persistent homology is attempted on non-centroid mode data.
 546            If data is not gridded and grid is False.
 547            If peak picking method is not implemented.
 548            If targeted_search is True but target_search_dict is None or invalid.
 549
 550        Returns
 551        -------
 552        None, but assigns the mass_features and eics attributes to the object.
 553
 554        """
 555        # Validate targeted search parameters
 556        if targeted_search:
 557            if target_search_dict is None:
 558                raise ValueError(
 559                    "target_search_dict must be provided when targeted_search is True"
 560                )
 561            required_keys = [
 562                "target_mz_list",
 563                "target_rt_list",
 564                "mz_tolerance_ppm",
 565                "rt_tolerance",
 566            ]
 567            for key in required_keys:
 568                if key not in target_search_dict:
 569                    raise ValueError(f"target_search_dict must contain '{key}'")
 570            if len(target_search_dict["target_mz_list"]) != len(
 571                target_search_dict["target_rt_list"]
 572            ):
 573                raise ValueError(
 574                    "target_mz_list and target_rt_list must have the same length"
 575                )
 576
 577        requested_method = self.parameters.lc_ms.peak_picking_method
 578        pp_method = self._resolve_lcms_peak_picking_method(ms_level=ms_level)
 579        if requested_method == "by_format":
 580            self.parameters.lc_ms.peak_picking_method = pp_method
 581
 582        if targeted_search:
 583            mf_type = target_search_dict.get("type", "targeted")
 584        else:
 585            mf_type = "untargeted"
 586
 587        if pp_method == "persistent homology":
 588            self.find_mass_features_ph(
 589                ms_level=ms_level,
 590                grid=grid,
 591                targeted_search=targeted_search,
 592                target_search_dict=target_search_dict,
 593                mf_type=mf_type,
 594                accumulate_features=accumulate_features,
 595            )
 596        elif pp_method == "centroided_persistent_homology":
 597            self.find_mass_features_ph_centroid(
 598                ms_level=ms_level,
 599                targeted_search=targeted_search,
 600                target_search_dict=target_search_dict,
 601                mf_type=mf_type,
 602                accumulate_features=accumulate_features,
 603            )
 604        else:
 605            raise ValueError("Peak picking method not implemented")
 606
 607        # Cluster mass features to remove redundant features
 608        self.cluster_mass_features(drop_children=True)
 609
 610        # Optionally assign MS2 scan numbers to mass features during peak picking
 611        # This helps with choosing representative features that have MS2 data
 612        if assign_ms2_scans:
 613            try:
 614                self._find_ms2_scans_for_mass_features(
 615                    mf_ids=None,  # Process all mass features
 616                    scan_filter=ms2_scan_filter,
 617                )
 618            except ValueError:
 619                # No MS2 scans found - this is okay, just skip
 620                pass
 621
 622        # Remove noisey mass features if designated in parameters
 623        if self.parameters.lc_ms.remove_redundant_mass_features and not targeted_search:
 624            self._remove_redundant_mass_features()
 625
 626    def integrate_mass_features(
 627        self,
 628        drop_if_fail=True,
 629        drop_duplicates=True,
 630        ms_level=1,
 631        induced_features=False,
 632    ):
 633        """Integrate mass features and extract EICs.
 634
 635        Populates the _eics attribute on the LCMSBase object for each unique mz in the mass_features dataframe and adds data (start_scan, final_scan, area) to the mass_features attribute.
 636
 637        Parameters
 638        ----------
 639        drop_if_fail : bool, optional
 640            Whether to drop mass features if the EIC limit calculations fail.
 641            Default is True.
 642        drop_duplicates : bool, optional
 643            Whether to mass features that appear to be duplicates
 644            (i.e., mz is similar to another mass feature and limits of the EIC are similar or encapsulating).
 645            Default is True.
 646        ms_level : int, optional
 647            The MS level to use. Default is 1.
 648        induced_features : bool, optional
 649            Whether the mass features to be intergrated were induced. Default is False.
 650
 651        Raises
 652        ------
 653        ValueError
 654            If no mass features are found.
 655            If no MS level data is found for the given MS level (either in data or in the scan data)
 656
 657        Returns
 658        -------
 659        None, but populates the eics attribute on the LCMSBase object and adds data (start_scan, final_scan, area) to the mass_features attribute.
 660
 661        Notes
 662        -----
 663        drop_if_fail is useful for discarding mass features that do not have good shapes, usually due to a detection on a shoulder of a peak or a noisy region (especially if minimal smoothing is used during mass feature detection).
 664        """
 665
 666        # Check if there is data
 667        if ms_level in self._ms_unprocessed.keys():
 668            raw_data = self._ms_unprocessed[ms_level].copy()
 669        else:
 670            raise ValueError("No MS level " + str(ms_level) + " data found")
 671
 672        # Check if mass_spectrum exists on each mass feature
 673        if induced_features:
 674            mf_dict = self.induced_mass_features
 675            if len(mf_dict) == 0:
 676                raise ValueError(
 677                    "No induced mass features found, did you run fill_missing_cluster_features() first?"
 678                )
 679
 680            ## remove not found induced mass features by mz <= 0 (-99 indicator)
 681            # also remove any where mz is nan
 682            mf_dict = {
 683                k: v for k, v in mf_dict.items() if v.mz > 0 and not np.isnan(v.mz)
 684            }
 685
 686        else:
 687            mf_dict = self.mass_features
 688            if len(mf_dict) == 0:
 689                raise ValueError(
 690                    "No mass features found, did you run find_mass_features() first?"
 691                )
 692
 693        # Subset scan data to only include correct ms_level
 694        scan_df_sub = self.scan_df[
 695            self.scan_df["ms_level"] == int(ms_level)
 696        ].reset_index(drop=True)
 697        if scan_df_sub.empty:
 698            raise ValueError("No MS level " + ms_level + " data found in scan data")
 699        scan_df_sub = scan_df_sub[["scan", "scan_time"]].copy()
 700
 701        mzs_to_extract = np.unique([mf.mz for mf in mf_dict.values()])
 702        mzs_to_extract.sort()
 703
 704        # Pre-sort raw_data by mz for faster filtering
 705        raw_data_sorted = raw_data.sort_values(["mz", "scan"]).reset_index(drop=True)
 706        raw_data_mz = raw_data_sorted["mz"].to_numpy()
 707        raw_scan = raw_data_sorted["scan"].to_numpy()
 708        raw_intensity = raw_data_sorted["intensity"].to_numpy()
 709        scan_ids = scan_df_sub["scan"].to_numpy()
 710        scan_time = scan_df_sub["scan_time"].to_numpy()
 711        n_scans = scan_ids.size
 712        # Sum into the scan-table row with the same float32 Kahan sum the
 713        # previous groupby used, so the chromatogram values stay identical.
 714        scan_order = np.argsort(scan_ids, kind="mergesort")
 715        sorted_scan_ids = scan_ids[scan_order]
 716        found = np.searchsorted(sorted_scan_ids, raw_scan)
 717        in_range = found < n_scans
 718        scan_col = np.empty(raw_scan.shape, dtype=np.intp)
 719        scan_col[in_range] = scan_order[found[in_range]]
 720        matched = np.zeros(raw_scan.shape, dtype=bool)
 721        matched[in_range] = sorted_scan_ids[found[in_range]] == raw_scan[in_range]
 722
 723        # Get EICs for each unique mz in mass features list
 724        for mz in mzs_to_extract:
 725            mz_max = mz + self.parameters.lc_ms.eic_tolerance_ppm * mz / 1e6
 726            mz_min = mz - self.parameters.lc_ms.eic_tolerance_ppm * mz / 1e6
 727
 728            # Use binary search for faster mz range filtering
 729            left_idx = int(np.searchsorted(raw_data_mz, mz_min, side="left"))
 730            right_idx = int(np.searchsorted(raw_data_mz, mz_max, side="right"))
 731            eic = np.zeros(n_scans, dtype=np.float32)
 732            if left_idx < right_idx:
 733                comp = np.zeros(n_scans, dtype=np.float32)
 734                cols = scan_col[left_idx:right_idx]
 735                vals = raw_intensity[left_idx:right_idx]
 736                ok = matched[left_idx:right_idx]
 737                for i in range(cols.size):
 738                    if not ok[i]:
 739                        continue
 740                    col = cols[i]
 741                    val = vals[i]
 742                    y = val - comp[col]
 743                    t = eic[col] + y
 744                    c = (t - eic[col]) - y
 745                    comp[col] = c if np.isfinite(c) else np.float32(0)
 746                    eic[col] = t
 747            myEIC = EIC_Data(
 748                scans=scan_ids.copy(),
 749                time=scan_time.copy(),
 750                eic=eic,
 751            )
 752            # Smooth EIC
 753            smoothed_eic = self.smooth_tic(myEIC.eic)
 754            smoothed_eic[smoothed_eic < 0] = 0
 755            myEIC.eic_smoothed = smoothed_eic
 756            self.eics[mz] = myEIC
 757
 758        # Get limits of mass features using EIC centroid detector and integrate
 759        for idx, mass_feature in list(mf_dict.items()):
 760            mz = mass_feature.mz
 761            apex_scan = mass_feature.apex_scan
 762
 763            # Pull EIC data and find apex scan index
 764            myEIC = self.eics[mz]
 765            mf_dict[idx]._eic_data = myEIC
 766            mf_dict[idx]._eic_mz = mz
 767            apex_index = np.searchsorted(myEIC.scans, apex_scan)
 768
 769            # Find left and right limits of peak using EIC centroid detector, add to EICData
 770            centroid_eics = self.eic_centroid_detector(
 771                myEIC.time,
 772                myEIC.eic_smoothed,
 773                mass_feature.intensity * 1.1,
 774                apex_indexes=[int(apex_index)],
 775            )
 776            l_a_r_scan_idx = [i for i in centroid_eics]
 777            if len(l_a_r_scan_idx) > 0:
 778                # Calculate number of consecutive scans with intensity > 0 and check if it is above the minimum consecutive scans
 779                # Find the number of consecutive non-zero values in the EIC segment
 780                mask = myEIC.eic[l_a_r_scan_idx[0][0] : l_a_r_scan_idx[0][2] + 1] > 0
 781                # Find the longest run of consecutive True values
 782                if np.any(mask):
 783                    # Find indices where mask changes value
 784                    diff = np.diff(np.concatenate(([0], mask.astype(int), [0])))
 785                    starts = np.where(diff == 1)[0]
 786                    ends = np.where(diff == -1)[0]
 787                    consecutive_scans = (ends - starts).max()
 788                else:
 789                    consecutive_scans = 0
 790                if consecutive_scans < self.parameters.lc_ms.consecutive_scan_min:
 791                    mf_dict.pop(idx)
 792                    continue
 793                # Add start and final scan to mass_features and EICData
 794                left_scan, right_scan = (
 795                    myEIC.scans[l_a_r_scan_idx[0][0]],
 796                    myEIC.scans[l_a_r_scan_idx[0][2]],
 797                )
 798                mf_scan_apex = [(left_scan, int(apex_scan), right_scan)]
 799                myEIC.apexes = myEIC.apexes + mf_scan_apex
 800                mf_dict[idx].start_scan = left_scan
 801                mf_dict[idx].final_scan = right_scan
 802
 803                # Find area under peak using limits from EIC centroid detector, add to mass_features and EICData
 804                area = np.trapezoid(
 805                    myEIC.eic_smoothed[l_a_r_scan_idx[0][0] : l_a_r_scan_idx[0][2] + 1],
 806                    myEIC.time[l_a_r_scan_idx[0][0] : l_a_r_scan_idx[0][2] + 1],
 807                )
 808                myEIC.areas = myEIC.areas + [area]
 809                self.eics[mz] = myEIC
 810                mf_dict[idx]._area = area
 811            else:
 812                if drop_if_fail is True:
 813                    mf_dict.pop(idx)
 814
 815        if drop_duplicates:
 816            # Prepare mass feature dataframe
 817            # Skip features without scan bounds (failed EIC limits when
 818            # drop_if_fail=False); they cannot participate in containment de-dup
 819            if induced_features:
 820                mf_df = self.mass_features_to_df(induced_features=True).copy()
 821                mf_df = mf_df[mf_df.start_scan.notna() & mf_df.final_scan.notna()]
 822            else:
 823                mf_df = self.mass_features_to_df(induced_features=False).copy()
 824                mf_df = mf_df[mf_df.start_scan.notna() & mf_df.final_scan.notna()]
 825
 826            # For each mass feature, find all mass features within the clustering tolerance ppm and drop if their start and end times are within another mass feature
 827            # Keep the first mass feature (highest persistence)
 828            for idx, mass_feature in mf_df.iterrows():
 829                mz = mass_feature.mz
 830                apex_scan = mass_feature.apex_scan
 831
 832                mf_df["mz_diff_ppm"] = np.abs(mf_df["mz"] - mz) / mz * 10**6
 833                mf_df_sub = mf_df[
 834                    mf_df["mz_diff_ppm"]
 835                    < self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel
 836                    * 10**6
 837                ].copy()
 838
 839                # For all mass features within the clustering tolerance, check if the start and end times are within the start and end times of the mass feature
 840                for idx2, mass_feature2 in mf_df_sub.iterrows():
 841                    if idx2 != idx:
 842                        s1, e1 = mass_feature.start_scan, mass_feature.final_scan
 843                        s2, e2 = mass_feature2.start_scan, mass_feature2.final_scan
 844                        # Guard incomplete bounds (None / NaN) so ordering
 845                        # comparisons never raise TypeError
 846                        if s1 is None or e1 is None or s2 is None or e2 is None:
 847                            continue
 848                        if pd.isna(s1) or pd.isna(e1) or pd.isna(s2) or pd.isna(e2):
 849                            continue
 850                        if s2 >= s1 and e2 <= e1:
 851                            if idx2 in self.mass_features.keys():
 852                                self.mass_features.pop(idx2)
 853
 854        # Filter MS2 scans to only include those within integration bounds
 855        # This ensures MS2 scans outside start_scan to final_scan are removed
 856        if induced_features:
 857            self._filter_ms2_scans_by_integration_bounds(
 858                mf_dict=self.induced_mass_features
 859            )
 860        else:
 861            self._filter_ms2_scans_by_integration_bounds(mf_dict=self.mass_features)
 862
 863    def find_c13_mass_features(self):
 864        """Mark likely C13 isotopes and connect to monoisoitopic mass features.
 865
 866        Returns
 867        -------
 868        None, but populates the monoisotopic_mf_id and isotopologue_type attributes to the indivual LCMSMassFeatures within the mass_features attribute of the LCMSBase object.
 869
 870        Raises
 871        ------
 872        ValueError
 873            If no mass features are found.
 874        """
 875        verbose = self.parameters.lc_ms.verbose_processing
 876        if verbose:
 877            print("evaluating mass features for C13 isotopes")
 878        if self.mass_features is None:
 879            raise ValueError("No mass features found, run find_mass_features() first")
 880
 881        # Data prep fo sparse distance matrix
 882        dims = ["mz", "scan_time"]
 883        mf_df = self.mass_features_to_df().copy()
 884        # Drop mass features that have no area (these are likely to be noise)
 885        mf_df = mf_df[mf_df["area"].notnull()]
 886        mf_df["mf_id"] = mf_df.index.values
 887        dims = ["mz", "scan_time"]
 888
 889        # Sort my ascending mz so we always get the monoisotopic mass first, regardless of the order/intensity of the mass features
 890        mf_df = mf_df.sort_values(by=["mz"]).reset_index(drop=True).copy()
 891
 892        mz_diff = 1.003355  # C13-C12 mass difference
 893        tol = [
 894            mf_df["mz"].median()
 895            * self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel,
 896            self.parameters.lc_ms.mass_feature_cluster_rt_tolerance * 0.5,
 897        ]  # mz, in relative; scan_time in minutes
 898
 899        # Compute inter-feature distances
 900        distances = None
 901        for i in range(len(dims)):
 902            # Construct k-d tree
 903            values = mf_df[dims[i]].values
 904            tree = KDTree(values.reshape(-1, 1))
 905
 906            max_tol = tol[i]
 907            if dims[i] == "mz":
 908                # Maximum absolute tolerance
 909                max_tol = mz_diff + tol[i]
 910
 911            # Compute sparse distance matrix
 912            # the larger the max_tol, the slower this operation is
 913            sdm = tree.sparse_distance_matrix(tree, max_tol, output_type="coo_matrix")
 914
 915            # Only consider forward case, exclude diagonal
 916            sdm = sparse.triu(sdm, k=1)
 917
 918            if dims[i] == "mz":
 919                min_tol = mz_diff - tol[i]
 920                # Get only the ones that are above the min tol
 921                idx = sdm.data > min_tol
 922
 923                # Reconstruct sparse distance matrix
 924                sdm = sparse.coo_matrix(
 925                    (sdm.data[idx], (sdm.row[idx], sdm.col[idx])),
 926                    shape=(len(values), len(values)),
 927                )
 928
 929            # Cast as binary matrix
 930            sdm.data = np.ones_like(sdm.data)
 931
 932            # Stack distances
 933            if distances is None:
 934                distances = sdm
 935            else:
 936                distances = distances.multiply(sdm)
 937
 938        # Extract indices of within-tolerance points
 939        distances = distances.tocoo()
 940        pairs = np.stack((distances.row, distances.col), axis=1)  # C12 to C13 pairs
 941
 942        # Turn pairs (which are index of mf_df) into mf_id and then into two dataframes to join to mf_df
 943        pairs_mf = pairs.copy()
 944        pairs_mf[:, 0] = mf_df.iloc[pairs[:, 0]].mf_id.values
 945        pairs_mf[:, 1] = mf_df.iloc[pairs[:, 1]].mf_id.values
 946
 947        # Connect monoisotopic masses with isotopologes within mass_features
 948        monos = np.setdiff1d(np.unique(pairs_mf[:, 0]), np.unique(pairs_mf[:, 1]))
 949        for mono in monos:
 950            self.mass_features[mono].monoisotopic_mf_id = mono
 951        pairs_iso_df = pd.DataFrame(pairs_mf, columns=["parent", "child"])
 952        while not pairs_iso_df.empty:
 953            pairs_iso_df = pairs_iso_df.set_index("parent", drop=False)
 954            m1_isos = pairs_iso_df.loc[monos, "child"].unique()
 955            for iso in m1_isos:
 956                # Set monoisotopic_mf_id and isotopologue_type for isotopologues
 957                parent = pairs_mf[pairs_mf[:, 1] == iso, 0]
 958                if len(parent) > 1:
 959                    # Choose the parent that is closest in time to the isotopologue
 960                    parent_time = [self.mass_features[p].retention_time for p in parent]
 961                    time_diff = [
 962                        np.abs(self.mass_features[iso].retention_time - x)
 963                        for x in parent_time
 964                    ]
 965                    parent = parent[np.argmin(time_diff)]
 966                else:
 967                    parent = parent[0]
 968                self.mass_features[iso].monoisotopic_mf_id = self.mass_features[
 969                    parent
 970                ].monoisotopic_mf_id
 971                if self.mass_features[iso].monoisotopic_mf_id is not None:
 972                    mass_diff = (
 973                        self.mass_features[iso].mz
 974                        - self.mass_features[
 975                            self.mass_features[iso].monoisotopic_mf_id
 976                        ].mz
 977                    )
 978                    self.mass_features[iso].isotopologue_type = "13C" + str(
 979                        int(round(mass_diff, 0))
 980                    )
 981
 982            # Drop the mono and iso from the pairs_iso_df
 983            pairs_iso_df = pairs_iso_df.drop(
 984                index=monos, errors="ignore"
 985            )  # Drop pairs where the parent is a child that is a child of a root
 986            pairs_iso_df = pairs_iso_df.set_index("child", drop=False)
 987            pairs_iso_df = pairs_iso_df.drop(index=m1_isos, errors="ignore")
 988
 989            if not pairs_iso_df.empty:
 990                # Get new monos, recognizing that these are just 13C isotopologues that are connected to other 13C isotopologues to repeat the process
 991                monos = np.setdiff1d(
 992                    np.unique(pairs_iso_df.parent), np.unique(pairs_iso_df.child)
 993                )
 994        if verbose:
 995            # Report fraction of compounds annotated with isotopes
 996            mf_df["c13_flag"] = np.where(
 997                np.logical_or(
 998                    np.isin(mf_df["mf_id"], pairs_mf[:, 0]),
 999                    np.isin(mf_df["mf_id"], pairs_mf[:, 1]),
1000                ),
1001                1,
1002                0,
1003            )
1004            print(
1005                str(round(len(mf_df[mf_df["c13_flag"] == 1]) / len(mf_df), ndigits=3))
1006                + " of mass features have or are C13 isotopes"
1007            )
1008
1009    def deconvolute_ms1_mass_features(self):
1010        """Deconvolute MS1 mass features
1011
1012        Deconvolute mass features ms1 spectrum based on the correlation of all masses within a spectrum over the EIC of the mass features
1013
1014        Parameters
1015        ----------
1016        None
1017
1018        Returns
1019        -------
1020        None, but assigns the _ms_deconvoluted_idx, mass_spectrum_deconvoluted_parent,
1021        and associated_mass_features_deconvoluted attributes to the mass features in the
1022        mass_features attribute of the LCMSBase object.
1023
1024        Raises
1025        ------
1026        ValueError
1027            If no mass features are found, must run find_mass_features() first.
1028            If no EICs are found, did you run integrate_mass_features() first?
1029
1030        """
1031        # Checks for set mass_features and eics
1032        if self.mass_features is None:
1033            raise ValueError(
1034                "No mass features found, did you run find_mass_features() first?"
1035            )
1036
1037        if self.eics == {}:
1038            raise ValueError(
1039                "No EICs found, did you run integrate_mass_features() first?"
1040            )
1041
1042        if 1 not in self._ms_unprocessed.keys():
1043            raise ValueError("No unprocessed MS1 spectra found.")
1044
1045        # Prep ms1 data
1046        ms1_data = self._ms_unprocessed[1].copy()
1047        ms1_data = ms1_data.set_index("scan")
1048
1049        # Prep mass feature summary
1050        mass_feature_df = self.mass_features_to_df()
1051
1052        # Loop through each mass feature
1053        for mf_id, mass_feature in self.mass_features.items():
1054            # Check that the mass_feature.mz attribute == the mz of the mass feature in the mass_feature_df
1055            if mass_feature.mz != mass_feature.ms1_peak.mz_exp:
1056                continue
1057
1058            # Get the left and right limits of the EIC of the mass feature
1059            l_scan, _, r_scan = mass_feature._eic_data.apexes[0]
1060
1061            # Pull from the _ms1_unprocessed data the scan range of interest and sort by mz
1062            ms1_data_sub = ms1_data.loc[l_scan:r_scan].copy()
1063            ms1_data_sub = ms1_data_sub.sort_values(by=["mz"]).reset_index(drop=False)
1064
1065            # Get the centroided masses of the mass feature
1066            mf_mspeak_mzs = mass_feature.mass_spectrum.mz_exp
1067
1068            # Find the closest mz in the ms1 data to the centroided masses of the mass feature
1069            ms1_data_sub["mass_feature_mz"] = mf_mspeak_mzs[
1070                find_closest(mf_mspeak_mzs, ms1_data_sub.mz.values)
1071            ]
1072
1073            # Drop rows with mz_diff > 0.01 between the mass feature mz and the ms1 data mz
1074            ms1_data_sub["mz_diff_rel"] = (
1075                np.abs(ms1_data_sub["mass_feature_mz"] - ms1_data_sub["mz"])
1076                / ms1_data_sub["mz"]
1077            )
1078            ms1_data_sub = ms1_data_sub[
1079                ms1_data_sub["mz_diff_rel"]
1080                < self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel
1081            ].reset_index(drop=True)
1082
1083            # Group by mass_feature_mz and scan and sum intensity
1084            ms1_data_sub_group = (
1085                ms1_data_sub.groupby(["mass_feature_mz", "scan"])["intensity"]
1086                .sum()
1087                .reset_index()
1088            )
1089
1090            # Calculate the correlation of the intensities of the mass feature and the ms1 data (set to 0 if no intensity)
1091            corr = (
1092                ms1_data_sub_group.pivot(
1093                    index="scan", columns="mass_feature_mz", values="intensity"
1094                )
1095                .fillna(0)
1096                .corr()
1097            )
1098
1099            # Subset the correlation matrix to only include the masses of the mass feature and those with a correlation > 0.8
1100            decon_corr_min = self.parameters.lc_ms.ms1_deconvolution_corr_min
1101
1102            # Try catch for KeyError in case the mass feature mz is not in the correlation matrix
1103            try:
1104                corr_subset = corr.loc[mass_feature.mz]
1105            except KeyError:
1106                # If the mass feature mz is not in the correlation matrix, skip to the next mass feature
1107                continue
1108
1109            corr_subset = corr_subset[corr_subset > decon_corr_min]
1110
1111            # Get the masses from the mass spectrum that are the result of the deconvolution
1112            mzs_decon = corr_subset.index.values
1113
1114            # Get the indices of the mzs_decon in mass_feature.mass_spectrum.mz_exp and assign to the mass feature
1115            mzs_decon_idx = [
1116                id
1117                for id, mz in enumerate(mass_feature.mass_spectrum.mz_exp)
1118                if mz in mzs_decon
1119            ]
1120            mass_feature._ms_deconvoluted_idx = mzs_decon_idx
1121
1122            # Check if the mass feature's ms1 peak is the largest in the deconvoluted mass spectrum
1123            if (
1124                mass_feature.ms1_peak.abundance
1125                == mass_feature.mass_spectrum.abundance[mzs_decon_idx].max()
1126            ):
1127                mass_feature.mass_spectrum_deconvoluted_parent = True
1128            else:
1129                mass_feature.mass_spectrum_deconvoluted_parent = False
1130
1131            # Check for other mass features that are in the deconvoluted mass spectrum and add the deconvoluted mass spectrum to the mass feature
1132            # Subset mass_feature_df to only include mass features that are within the clustering tolerance
1133            mass_feature_df_sub = mass_feature_df[
1134                abs(mass_feature.retention_time - mass_feature_df["scan_time"])
1135                < self.parameters.lc_ms.mass_feature_cluster_rt_tolerance
1136            ].copy()
1137            # Calculate the mz difference in ppm between the mass feature and the peaks in the deconvoluted mass spectrum
1138            mass_feature_df_sub["mz_diff_ppm"] = [
1139                np.abs(mzs_decon - mz).min() / mz * 10**6
1140                for mz in mass_feature_df_sub["mz"]
1141            ]
1142            # Subset mass_feature_df to only include mass features that are within 1 ppm of the deconvoluted masses
1143            mfs_associated_decon = mass_feature_df_sub[
1144                mass_feature_df_sub["mz_diff_ppm"]
1145                < self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel * 10**6
1146            ].index.values
1147
1148            mass_feature.associated_mass_features_deconvoluted = mfs_associated_decon
1149
1150    def _remove_redundant_mass_features(
1151        self,
1152    ) -> None:
1153        """
1154        Identify and remove redundant mass features that are likely contaminants based on their m/z values and scan frequency.
1155        Especially useful for HILIC data where signals do not return to baseline between peaks or for data with significant background noise.
1156
1157        Contaminants are characterized by:
1158        1. Similar m/z values (within ppm_tolerance)
1159        2. High frequency across scan numbers (ubiquitous presence)
1160
1161        Notes
1162        -----
1163        Depends on self.mass_features being populated, uses the parameters in self.parameters.lc_ms for tolerances (mass_feature_cluster_mz_tolerance_rel)
1164        """
1165        ppm_tolerance = (
1166            self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel * 1e6
1167        )
1168        min_scan_frequency = self.parameters.lc_ms.redundant_scan_frequency_min
1169        n_retain = self.parameters.lc_ms.redundant_feature_retain_n
1170
1171        df = self.mass_features_to_df()
1172
1173        if df.empty:
1174            return pd.DataFrame()
1175        # df index should be mf_id
1176        if "mf_id" not in df.columns:
1177            if "mf_id" in df.index.names:
1178                df = df.reset_index()
1179            else:
1180                raise ValueError("DataFrame must contain 'mf_id' column or index.")
1181
1182        # Sort by m/z for efficient grouping
1183        df_sorted = df.sort_values("mz").reset_index(drop=True)
1184
1185        # Calculate total number of unique scans for frequency calculation
1186        # Calculate total possible scans (check the cluster rt tolerance and the min rt and max rt of the data)
1187        total_time = self.scan_df["scan_time"].max() - self.scan_df["scan_time"].min()
1188        cluster_rt_tolerance = self.parameters.lc_ms.mass_feature_cluster_rt_tolerance
1189        # If the feature was detected in every possible scan (and then rolled up), it would be in this many scans
1190        total_scans = int(total_time / cluster_rt_tolerance) + 1
1191
1192        # Group similar m/z values using ppm tolerance
1193        mz_groups = []
1194        current_group = []
1195
1196        for i, row in df_sorted.iterrows():
1197            current_mz = row["mz"]
1198
1199            if not current_group:
1200                # Start first group
1201                current_group = [i]
1202            else:
1203                # Check if current m/z is within tolerance of group representative
1204                group_representative_mz = df_sorted.iloc[current_group[0]]["mz"]
1205                ppm_diff = (
1206                    abs(current_mz - group_representative_mz)
1207                    / group_representative_mz
1208                    * 1e6
1209                )
1210
1211                if ppm_diff <= ppm_tolerance:
1212                    # Add to current group
1213                    current_group.append(i)
1214                else:
1215                    # Start new group, but first process current group
1216                    if len(current_group) > 0:
1217                        mz_groups.append(current_group)
1218                    current_group = [i]
1219
1220        # Don't forget the last group
1221        if current_group:
1222            mz_groups.append(current_group)
1223
1224        # Analyze each m/z group for contaminant characteristics
1225
1226        for group_indices in mz_groups:
1227            group_data = df_sorted.iloc[group_indices]
1228
1229            # Calculate group statistics
1230            unique_scans = group_data["apex_scan"].nunique()
1231            scan_frequency = unique_scans / total_scans
1232
1233            # Check if this group meets contaminant criteria
1234            if scan_frequency >= min_scan_frequency:
1235                group_data = group_data.sort_values("intensity", ascending=False)
1236                non_representative_mf_id = group_data.iloc[n_retain:][
1237                    "mf_id"
1238                ].tolist()  # These will be removed
1239
1240                self.mass_features = {
1241                    k: v
1242                    for k, v in self.mass_features.items()
1243                    if k not in non_representative_mf_id
1244                }
1245
1246    def _remove_mass_features_by_peak_metrics(self, induced_features=False) -> None:
1247        """Remove mass features based on peak metrics defined in mass_feature_attribute_filter_dict.
1248
1249        This method filters mass features based on various peak shape metrics and quality indicators
1250        such as noise scores, Gaussian similarity, tailing factors, dispersity index, etc.
1251
1252        The filtering criteria are defined in the mass_feature_attribute_filter_dict parameter,
1253        which should contain attribute names as keys and filter specifications as values.
1254
1255        Filter specification format:
1256        {attribute_name: {'value': threshold, 'operator': comparison}}
1257
1258        Available operators:
1259        - '>' or 'greater': Keep features where attribute > threshold
1260        - '<' or 'less': Keep features where attribute < threshold
1261        - '>=' or 'greater_equal': Keep features where attribute >= threshold
1262        - '<=' or 'less_equal': Keep features where attribute <= threshold
1263
1264        Examples:
1265        - {'noise_score_max': {'value': 0.5, 'operator': '>='}} - Keep features with noise_score_max >= 0.5
1266        - {'dispersity_index': {'value': 0.1, 'operator': '<'}} - Keep features with dispersity_index < 0.1
1267        - {'gaussian_similarity': {'value': 0.7, 'operator': '>='}} - Keep features with gaussian_similarity >= 0.7
1268
1269        Parameters
1270        ----------
1271        induced_features : bool, optional
1272            If True, filter induced_mass_features instead of regular mass_features. Default is False.
1273
1274        Returns
1275        -------
1276        None
1277            Modifies self.mass_features or self.induced_mass_features in place by removing filtered features.
1278
1279        Raises
1280        ------
1281        ValueError
1282            If no mass features are found, if an invalid attribute is specified, or if filter specification is malformed.
1283        """
1284        # Select the appropriate mass features dictionary
1285        if induced_features:
1286            mf_dict = self.induced_mass_features
1287            mf_type = "induced mass features"
1288        else:
1289            mf_dict = self.mass_features
1290            mf_type = "mass features"
1291
1292        if mf_dict is None or len(mf_dict) == 0:
1293            raise ValueError(
1294                f"No {mf_type} found, run {'gap filling' if induced_features else 'find_mass_features()'} first"
1295            )
1296
1297        filter_dict = self.parameters.lc_ms.mass_feature_attribute_filter_dict
1298
1299        if not filter_dict:
1300            # No filtering criteria specified, return early
1301            return
1302
1303        verbose = self.parameters.lc_ms.verbose_processing
1304        initial_count = len(mf_dict)
1305
1306        if verbose:
1307            print(
1308                f"Filtering {mf_type} using peak metrics. Initial count: {initial_count}"
1309            )
1310
1311        # List to collect IDs of mass features to remove
1312        features_to_remove = []
1313
1314        for mf_id, mass_feature in mf_dict.items():
1315            should_remove = False
1316
1317            for attribute_name, filter_spec in filter_dict.items():
1318                # Validate filter specification structure
1319                if not isinstance(filter_spec, dict):
1320                    raise ValueError(
1321                        f"Filter specification for '{attribute_name}' must be a dictionary with 'value' and 'operator' keys"
1322                    )
1323
1324                if "value" not in filter_spec or "operator" not in filter_spec:
1325                    raise ValueError(
1326                        f"Filter specification for '{attribute_name}' must contain both 'value' and 'operator' keys"
1327                    )
1328
1329                threshold_value = filter_spec["value"]
1330                operator = filter_spec["operator"].lower().strip()
1331
1332                # Validate operator
1333                valid_operators = {
1334                    ">",
1335                    "<",
1336                    ">=",
1337                    "<=",
1338                    "greater",
1339                    "less",
1340                    "greater_equal",
1341                    "less_equal",
1342                }
1343                if operator not in valid_operators:
1344                    raise ValueError(
1345                        f"Invalid operator '{operator}' for attribute '{attribute_name}'. Valid operators: {valid_operators}"
1346                    )
1347
1348                # Normalize operator names
1349                operator_map = {
1350                    "greater": ">",
1351                    "less": "<",
1352                    "greater_equal": ">=",
1353                    "less_equal": "<=",
1354                }
1355                operator = operator_map.get(operator, operator)
1356
1357                # Get the attribute value from the mass feature
1358                try:
1359                    if hasattr(mass_feature, attribute_name):
1360                        attribute_value = getattr(mass_feature, attribute_name)
1361                    else:
1362                        raise ValueError(
1363                            f"Mass feature does not have attribute '{attribute_name}'"
1364                        )
1365
1366                    # Handle None values or attributes that haven't been calculated
1367                    if attribute_value is None:
1368                        if verbose:
1369                            print(
1370                                f"Warning: Mass feature {mf_id} has None value for '{attribute_name}'. Removing feature."
1371                            )
1372                        should_remove = True
1373                        break
1374
1375                    # Handle numpy arrays (like half_height_width which returns mean)
1376                    if hasattr(attribute_value, "__len__") and not isinstance(
1377                        attribute_value, str
1378                    ):
1379                        # For arrays, we use the mean or appropriate summary statistic
1380                        if attribute_name == "half_height_width":
1381                            # half_height_width property already returns the mean
1382                            pass
1383                        else:
1384                            attribute_value = float(np.mean(attribute_value))
1385
1386                    # Handle NaN values
1387                    if np.isnan(float(attribute_value)):
1388                        if verbose:
1389                            print(
1390                                f"Warning: Mass feature {mf_id} has NaN value for '{attribute_name}'. Removing feature."
1391                            )
1392                        should_remove = True
1393                        break
1394
1395                    # Apply the threshold comparison based on operator
1396                    attribute_value = float(attribute_value)
1397                    threshold_value = float(threshold_value)
1398
1399                    if operator == ">" and not (attribute_value > threshold_value):
1400                        should_remove = True
1401                        break
1402                    elif operator == "<" and not (attribute_value < threshold_value):
1403                        should_remove = True
1404                        break
1405                    elif operator == ">=" and not (attribute_value >= threshold_value):
1406                        should_remove = True
1407                        break
1408                    elif operator == "<=" and not (attribute_value <= threshold_value):
1409                        should_remove = True
1410                        break
1411
1412                except (AttributeError, ValueError, TypeError) as e:
1413                    if verbose:
1414                        print(
1415                            f"Error evaluating filter '{attribute_name}' for mass feature {mf_id}: {e}"
1416                        )
1417                    should_remove = True
1418                    break
1419
1420            if should_remove:
1421                features_to_remove.append(mf_id)
1422
1423        # Remove filtered mass features
1424        for mf_id in features_to_remove:
1425            del mf_dict[mf_id]
1426
1427        if verbose and len(features_to_remove) > 0:
1428            print(
1429                f"Removed {len(features_to_remove)} {mf_type} based on peak metrics. Remaining: {len(mf_dict)}"
1430            )
1431
1432        # Update the appropriate dictionary
1433        if induced_features:
1434            self.induced_mass_features = mf_dict
1435        else:
1436            self.mass_features = mf_dict
1437
1438        # Clean up unassociated EICs and ms1 data (only for regular features)
1439        self._remove_unassociated_eics()
1440        self._remove_unassociated_ms1_spectra()
1441
1442    def _remove_unassociated_eics(self) -> None:
1443        """Remove EICs that are not associated with any mass features.
1444
1445        This method cleans up the eics attribute by removing any EICs that do not correspond to
1446        any mass features currently stored in the mass_features attribute. This is useful for
1447        freeing up memory and ensuring that only relevant EICs are retained.
1448
1449        Returns
1450        -------
1451        None
1452            Modifies self.eics in place by removing unassociated EICs.
1453        """
1454        if self.mass_features is None or len(self.mass_features) == 0:
1455            self.eics = {}
1456            return
1457
1458        # Get the set of m/z values associated with current mass features
1459        associated_mzs = {mf.mz for mf in self.mass_features.values()}
1460
1461        # Remove EICs that are not associated with any mass features
1462        self.eics = {mz: eic for mz, eic in self.eics.items() if mz in associated_mzs}
1463
1464    def _remove_unassociated_ms1_spectra(self) -> None:
1465        """Remove MS1 spectra that are not associated with any mass features.
1466        This method cleans up the _ms_unprocessed attribute by removing any MS1 spectra that do not correspond to
1467        any mass features currently stored in the mass_features attribute. This is useful for freeing up memory
1468        and ensuring that only relevant MS1 spectra are retained.
1469
1470        Returns
1471        -------
1472        None
1473        """
1474        if self.mass_features is None or len(self.mass_features) == 0:
1475            self._ms_unprocessed = {}
1476            return
1477
1478        # Get the set of m/z values associated with current mass features
1479        associated_ms1_scans = {mf.apex_scan for mf in self.mass_features.values()}
1480        associated_ms1_scans = [int(scan) for scan in associated_ms1_scans]
1481
1482        # Get keys within the _ms attribute (these are individual MassSpectrum objects)
1483        current_stored_spectra = list(set(self._ms.keys()))
1484        if len(current_stored_spectra) == 0:
1485            return
1486        current_stored_spectra = [int(scan) for scan in current_stored_spectra]
1487
1488        # Filter the current_stored_spectra to only ms1 scans
1489        current_stored_spectra_ms1 = [
1490            scan for scan in current_stored_spectra if scan in self.ms1_scans
1491        ]
1492
1493        # Remove MS1 spectra that are not associated with any mass features
1494        scans_to_drop = [
1495            scan
1496            for scan in current_stored_spectra_ms1
1497            if scan not in associated_ms1_scans
1498        ]
1499        for scan in scans_to_drop:
1500            if scan in self._ms:
1501                del self._ms[scan]
1502
1503
1504class PHCalculations:
1505    """Methods for performing calculations related to 2D peak picking via persistent homology on LCMS data.
1506
1507    Notes
1508    -----
1509    This class is intended to be used as a mixin for the LCMSBase class.
1510
1511    Methods
1512    -------
1513    * sparse_mean_filter(idx, V, radius=[0, 1, 1]).
1514        Sparse implementation of a mean filter.
1515    * embed_unique_indices(a).
1516        Creates an array of indices, sorted by unique element.
1517    * sparse_upper_star(idx, V).
1518        Sparse implementation of an upper star filtration.
1519    * check_if_grid(data).
1520        Check if the data is gridded in mz space.
1521    * grid_data(data).
1522        Grid the data in the mz dimension.
1523    * find_mass_features_ph(ms_level=1, grid=True).
1524        Find mass features within an LCMSBase object using persistent homology.
1525    * cluster_mass_features(drop_children=True).
1526        Cluster regions of interest.
1527    """
1528
1529    @staticmethod
1530    def sparse_mean_filter(idx, V, radius=[0, 1, 1]):
1531        """Sparse implementation of a mean filter.
1532
1533        Parameters
1534        ----------
1535        idx : :obj:`~numpy.array`
1536            Edge indices for each dimension (MxN).
1537        V : :obj:`~numpy.array`
1538            Array of intensity data (Mx1).
1539        radius : float or list
1540            Radius of the sparse filter in each dimension. Values less than
1541            zero indicate no connectivity in that dimension.
1542
1543        Returns
1544        -------
1545        :obj:`~numpy.array`
1546            Filtered intensities (Mx1).
1547
1548        Notes
1549        -----
1550        This function has been adapted from the original implementation in the Deimos package: https://github.com/pnnl/deimos.
1551        This is a static method.
1552        """
1553
1554        # Copy indices
1555        idx = idx.copy().astype(V.dtype)
1556
1557        # Scale
1558        for i, r in enumerate(radius):
1559            # Increase inter-index distance
1560            if r < 1:
1561                idx[:, i] *= 2
1562
1563            # Do nothing
1564            elif r == 1:
1565                pass
1566
1567            # Decrease inter-index distance
1568            else:
1569                idx[:, i] /= r
1570
1571        # Connectivity matrix
1572        cmat = KDTree(idx)
1573        cmat = cmat.sparse_distance_matrix(cmat, 1, p=np.inf, output_type="coo_matrix")
1574        cmat.setdiag(1)
1575
1576        # Pair indices
1577        I, J = cmat.nonzero()
1578
1579        # Delete cmat
1580        cmat_shape = cmat.shape
1581        del cmat
1582
1583        # Sum over columns
1584        V_sum = sparse.bsr_matrix(
1585            (V[J], (I, I)), shape=cmat_shape, dtype=V.dtype
1586        ).diagonal(0)
1587
1588        # Count over columns
1589        V_count = sparse.bsr_matrix(
1590            (np.ones_like(J), (I, I)), shape=cmat_shape, dtype=V.dtype
1591        ).diagonal(0)
1592
1593        return V_sum / V_count
1594
1595    @staticmethod
1596    def embed_unique_indices(a):
1597        """Creates an array of indices, sorted by unique element.
1598
1599        Parameters
1600        ----------
1601        a : :obj:`~numpy.array`
1602            Array of unique elements (Mx1).
1603
1604        Returns
1605        -------
1606        :obj:`~numpy.array`
1607            Array of indices (Mx1).
1608
1609        Notes
1610        -----
1611        This function has been adapted from the original implementation in the Deimos package: https://github.com/pnnl/deimos
1612        This is a static method.
1613        """
1614
1615        def count_tens(n):
1616            # Count tens
1617            ntens = (n - 1) // 10
1618
1619            while True:
1620                ntens_test = (ntens + n - 1) // 10
1621
1622                if ntens_test == ntens:
1623                    return ntens
1624                else:
1625                    ntens = ntens_test
1626
1627        def arange_exclude_10s(n):
1628            # How many 10s will there be?
1629            ntens = count_tens(n)
1630
1631            # Base array
1632            arr = np.arange(0, n + ntens)
1633
1634            # Exclude 10s
1635            arr = arr[(arr == 0) | (arr % 10 != 0)][:n]
1636
1637            return arr
1638
1639        # Creates an array of indices, sorted by unique element
1640        idx_sort = np.argsort(a)
1641        idx_unsort = np.argsort(idx_sort)
1642
1643        # Sorts records array so all unique elements are together
1644        sorted_a = a[idx_sort]
1645
1646        # Returns the unique values, the index of the first occurrence,
1647        # and the count for each element
1648        vals, idx_start, count = np.unique(
1649            sorted_a, return_index=True, return_counts=True
1650        )
1651
1652        # Splits the indices into separate arrays
1653        splits = np.split(idx_sort, idx_start[1:])
1654
1655        # Creates unique indices for each split
1656        idx_unq = np.concatenate([arange_exclude_10s(len(x)) for x in splits])
1657
1658        # Reorders according to input array
1659        idx_unq = idx_unq[idx_unsort]
1660
1661        # Magnitude of each index
1662        exp = np.log10(
1663            idx_unq, where=idx_unq > 0, out=np.zeros_like(idx_unq, dtype=np.float64)
1664        )
1665        idx_unq_mag = np.power(10, np.floor(exp) + 1)
1666
1667        # Result
1668        return a + idx_unq / idx_unq_mag
1669
1670    @staticmethod
1671    def _greedy_keep_mask(all_values, tol, relative):
1672        """Rows to keep, highest sort_by first.
1673
1674        ``all_values`` is aligned to a frame already sorted by ``sort_by``
1675        descending. A later row is dropped when a kept row passes every
1676        dimension. Relative limits use the kept row. Values are compared
1677        as float32. No pair matrix is built.
1678        """
1679        values = [np.asarray(col, dtype=np.float32) for col in all_values]
1680        n = len(values[0]) if values else 0
1681        keep = np.ones(n, dtype=bool)
1682        if n == 0 or any(float(t) < 0 for t in tol):
1683            return keep
1684        n_dims = len(values)
1685        grid = {}
1686        loose = []
1687
1688        def bin_index(raw, tol_d, rel):
1689            fv = np.float32(raw)
1690            if float(tol_d) == 0.0:
1691                if fv == 0:
1692                    fv = np.float32(0.0)
1693                return int(np.asarray(fv).view(np.int32))
1694            if (not rel) and float(tol_d) > 0.0:
1695                return int(math.floor(float(fv) / float(tol_d)))
1696            if rel and 0.0 < float(tol_d) < 1.0 and float(fv) > 0.0:
1697                return int(
1698                    math.floor(math.log(float(fv)) / math.log(1.0 + float(tol_d)))
1699                )
1700            return None
1701
1702        def query_bins(raw, tol_d, rel):
1703            key = bin_index(raw, tol_d, rel)
1704            if key is None:
1705                return None
1706            if float(tol_d) == 0.0:
1707                return (key,)
1708            fv = float(np.float32(raw))
1709            if not rel:
1710                width = float(tol_d)
1711                lo = int(math.floor((fv - width) / width)) - 1
1712                hi = int(math.floor((fv + width) / width)) + 1
1713                return range(lo, hi + 1)
1714            log_width = math.log(1.0 + float(tol_d))
1715            lo_v = fv / (1.0 + float(tol_d))
1716            hi_v = fv / (1.0 - float(tol_d))
1717            if lo_v <= 0.0:
1718                return None
1719            lo = int(math.floor(math.log(lo_v) / log_width)) - 1
1720            hi = int(math.floor(math.log(hi_v) / log_width)) + 1
1721            return range(lo, hi + 1)
1722
1723        def claims(parent_i, child_i):
1724            for d in range(n_dims):
1725                pv = values[d][parent_i]
1726                cv = values[d][child_i]
1727                if not np.isfinite(pv) or not np.isfinite(cv):
1728                    return False
1729                diff = np.float32(np.abs(np.float32(cv) - np.float32(pv)))
1730                if relative[d]:
1731                    limit = np.float32(np.float32(tol[d]) * np.float32(pv))
1732                else:
1733                    limit = np.float32(tol[d])
1734                if diff > limit:
1735                    return False
1736            return True
1737
1738        for i in range(n):
1739            if not all(np.isfinite(values[d][i]) for d in range(n_dims)):
1740                continue
1741            bins = [
1742                query_bins(values[d][i], tol[d], relative[d]) for d in range(n_dims)
1743            ]
1744            if any(bucket is None for bucket in bins):
1745                search = list(loose)
1746                for bucket in grid.values():
1747                    search.extend(bucket)
1748            else:
1749                search = []
1750                for combo in product(*bins):
1751                    search.extend(grid.get(combo, ()))
1752                search.extend(loose)
1753            if any(claims(parent_i, i) for parent_i in search):
1754                keep[i] = False
1755                continue
1756            keys = [bin_index(values[d][i], tol[d], relative[d]) for d in range(n_dims)]
1757            if any(key is None for key in keys):
1758                loose.append(i)
1759            else:
1760                grid.setdefault(tuple(keys), []).append(i)
1761        return keep
1762
1763    @staticmethod
1764    def roll_up_dataframe(
1765        df: pd.DataFrame,
1766        sort_by: str,
1767        tol: list,
1768        relative: list,
1769        dims: list,
1770        memory_opt_threshold: int = 10000,
1771    ):
1772        """Subset data by rolling up into apex in appropriate dimensions.
1773
1774        Parameters
1775        ----------
1776        data : pd.DataFrame
1777            The input data containing "dims" columns and the "sort_by" column.
1778        sort_by : str
1779            The column to sort the data by, this will determine which mass features get rolled up into a parent mass feature
1780            (i.e., the mass feature with the highest value in the sort_by column).
1781        dims : list
1782            A list of dimension names (column names in the data DataFrame) to roll up the mass features by.
1783        tol : list
1784            A list of tolerances for each dimension. The length of the list must match the number of dimensions.
1785            The tolerances can be relative (as a fraction of the maximum value in that dimension) or absolute (in the units of that dimension).
1786            If relative is True, a kept row claims a later row when the absolute difference is at most tol times the kept row's value.
1787            If relative is False, the absolute difference must be at most tol.
1788        relative : list
1789            A list of booleans indicating whether the tolerance for each dimension is relative (True) or absolute (False).
1790        memory_opt_threshold : int, optional
1791            Accepted for compatibility. Ignored. Every frame uses the same greedy search. Default is 10000.
1792
1793        Returns
1794        -------
1795        pd.DataFrame
1796            A DataFrame with only the rolled up mass features, with the original index and columns.
1797
1798
1799        Raises
1800        ------
1801        ValueError
1802            If the input data is not a pandas DataFrame.
1803            If the input data does not have columns for each of the dimensions in "dims".
1804            If the length of "dims", "tol", and "relative" do not match.
1805        """
1806        og_columns = df.columns.copy()
1807
1808        # Unindex the data, but keep the original index
1809        if df.index.name is not None:
1810            og_index = df.index.name
1811        else:
1812            og_index = "index"
1813        df = df.reset_index(drop=False)
1814
1815        # Sort data by sort_by column, and reindex
1816        df = df.sort_values(by=sort_by, ascending=False).reset_index(drop=True)
1817
1818        # Check that data is a DataFrame and has columns for each of the dims
1819        if not isinstance(df, pd.DataFrame):
1820            raise ValueError("Data must be a pandas DataFrame")
1821        for dim in dims:
1822            if dim not in df.columns:
1823                raise ValueError(f"Data must have a column for {dim}")
1824        if len(dims) != len(tol) or len(dims) != len(relative):
1825            raise ValueError(
1826                "Dimensions, tolerances, and relative flags must be the same length"
1827            )
1828
1829        # Pre-compute all values arrays
1830        all_values = [df[dim].values for dim in dims]
1831
1832        keep = PHCalculations._greedy_keep_mask(all_values, tol, relative)
1833        df_sub = df.iloc[np.flatnonzero(keep)]
1834
1835        # Set index back to og_index and only keep original columns
1836        df_sub = df_sub.set_index(og_index).sort_index()[og_columns]
1837
1838        return df_sub
1839
1840    def sparse_upper_star(self, idx, V):
1841        """Sparse implementation of an upper star filtration.
1842
1843        Parameters
1844        ----------
1845        idx : :obj:`~numpy.array`
1846            Edge indices for each dimension (MxN).
1847        V : :obj:`~numpy.array`
1848            Array of intensity data (Mx1).
1849        Returns
1850        -------
1851        idx : :obj:`~numpy.array`
1852            Index of filtered points (Mx1).
1853        persistence : :obj:`~numpy.array`
1854            Persistence of each filtered point (Mx1).
1855
1856        Notes
1857        -----
1858        This function has been adapted from the original implementation in the Deimos package: https://github.com/pnnl/deimos
1859        """
1860
1861        # Invert
1862        V = -1 * V.copy().astype(int)
1863
1864        # Embed indices
1865        V = self.embed_unique_indices(V)
1866
1867        # Connectivity matrix
1868        cmat = KDTree(idx)
1869        cmat = cmat.sparse_distance_matrix(cmat, 1, p=np.inf, output_type="coo_matrix")
1870        cmat.setdiag(1)
1871        cmat = sparse.triu(cmat)
1872
1873        # Pairwise minimums
1874        I, J = cmat.nonzero()
1875        d = np.maximum(V[I], V[J])
1876
1877        # Delete connectiity matrix
1878        cmat_shape = cmat.shape
1879        del cmat
1880
1881        # Sparse distance matrix
1882        sdm = sparse.coo_matrix((d, (I, J)), shape=cmat_shape)
1883
1884        # Delete pairwise mins
1885        del d, I, J
1886
1887        # Persistence homology
1888        ph = ripser(sdm, distance_matrix=True, maxdim=0)["dgms"][0]
1889
1890        # Bound death values
1891        ph[ph[:, 1] == np.inf, 1] = np.max(V)
1892
1893        # Construct tree to query against
1894        tree = KDTree(V.reshape((-1, 1)))
1895
1896        # Get the indexes of the first nearest neighbor by birth
1897        _, nn = tree.query(ph[:, 0].reshape((-1, 1)), k=1, workers=-1)
1898
1899        return nn, -(ph[:, 0] // 1 - ph[:, 1] // 1)
1900
1901    def check_if_grid(self, data):
1902        """Check if the data are gridded in mz space.
1903
1904        Parameters
1905        ----------
1906        data : DataFrame
1907            DataFrame containing the mass spectrometry data.  Needs to have mz and scan columns.
1908
1909        Returns
1910        -------
1911        bool
1912            True if the data is gridded in the mz direction, False otherwise.
1913
1914        Notes
1915        -----
1916        This function is used within the grid_data function and the find_mass_features function and is not intended to be called directly.
1917        """
1918        # Calculate the difference between consecutive mz values in a single scan
1919        dat_check = data.copy().reset_index(drop=True)
1920        dat_check["mz_diff"] = np.abs(dat_check["mz"].diff())
1921        mz_diff_min = (
1922            dat_check.groupby("scan")["mz_diff"].min().min()
1923        )  # within each scan, what is the smallest mz difference between consecutive mz values
1924
1925        # Find the mininum mz difference between mz values in the data; regardless of scan
1926        dat_check_mz = dat_check[["mz"]].drop_duplicates().copy()
1927        dat_check_mz = dat_check_mz.sort_values(by=["mz"]).reset_index(drop=True)
1928        dat_check_mz["mz_diff"] = np.abs(dat_check_mz["mz"].diff())
1929
1930        # Get minimum mz_diff between mz values in the data
1931        mz_diff_min_raw = dat_check_mz["mz_diff"].min()
1932
1933        # If the minimum mz difference between mz values in the data is less than the minimum mz difference between mz values within a single scan, then the data is not gridded
1934        if mz_diff_min_raw < mz_diff_min:
1935            return False
1936        else:
1937            return True
1938
1939    def grid_data(self, data, attempts=5):
1940        """Grid the data in the mz dimension.
1941
1942        Data must be gridded prior to persistent homology calculations and computing average mass spectrum
1943
1944        Parameters
1945        ----------
1946        data : DataFrame
1947            The input data containing mz, scan, scan_time, and intensity columns.
1948        attempts : int, optional
1949            The number of attempts to grid the data. Default is 5.
1950
1951        Returns
1952        -------
1953        DataFrame
1954            The gridded data with mz, scan, scan_time, and intensity columns.
1955
1956        Raises
1957        ------
1958        ValueError
1959            If gridding fails after the specified number of attempts.
1960        """
1961        attempt_i = 0
1962        while attempt_i < attempts:
1963            attempt_i += 1
1964            data = self._grid_data(data)
1965
1966            if self.check_if_grid(data):
1967                return data
1968
1969        if not self.check_if_grid(data):
1970            raise ValueError(
1971                "Gridding failed after "
1972                + str(attempt_i)
1973                + " attempts. Please check the data."
1974            )
1975        else:
1976            return data
1977
1978    def _grid_data(self, data):
1979        """Internal method to grid the data in the mz dimension.
1980
1981        Notes
1982        -----
1983        This method is called by the grid_data method and should not be called directly.
1984        It will attempt to grid the data in the mz dimension by creating a grid of mz values based on the minimum mz difference within each scan,
1985        but it does not check if the data is already gridded or if the gridding is successful.
1986
1987        Parameters
1988        ----------
1989        data : pd.DataFrame or pl.DataFrame
1990            The input data to grid.
1991
1992        Returns
1993        -------
1994        pd.DataFrame or pl.DataFrame
1995            The data after attempting to grid it in the mz dimension.
1996        """
1997        # Calculate the difference between consecutive mz values in a single scan for grid spacing
1998        data_w = data.copy().reset_index(drop=True)
1999        data_w["mz_diff"] = np.abs(data_w["mz"].diff())
2000        mz_diff_min = data_w.groupby("scan")["mz_diff"].min().min() * 0.99999
2001
2002        # Need high intensity mz values first so they are parents in the output pairs stack
2003        dat_mz = data_w[["mz", "intensity"]].sort_values(
2004            by=["intensity"], ascending=False
2005        )
2006        dat_mz = dat_mz[["mz"]].drop_duplicates().reset_index(drop=True).copy()
2007        # Same kept centers as the distance-pair grid. A 1D index replaces the
2008        # pair matrix, which grew with the square of the scan count.
2009        mz_dat_np = np.sort(
2010            _mz_grid_centers(dat_mz["mz"].to_numpy(dtype=np.float64), mz_diff_min)
2011        )
2012
2013        # Sort data by mz and recast mz to nearest value in mz_dat_np
2014        data_w = data_w.sort_values(by=["mz"]).reset_index(drop=True).copy()
2015        data_w["mz_new"] = mz_dat_np[find_closest(mz_dat_np, data_w["mz"].values)]
2016        data_w["mz_diff"] = np.abs(data_w["mz"] - data_w["mz_new"])
2017
2018        # Rename mz_new as mz; drop mz_diff; groupby scan and mz and sum intensity
2019        new_data_w = data_w.rename(columns={"mz": "mz_orig", "mz_new": "mz"}).copy()
2020        new_data_w = (
2021            new_data_w.drop(columns=["mz_diff", "mz_orig"])
2022            .groupby(["scan", "mz"])["intensity"]
2023            .sum()
2024            .reset_index()
2025        )
2026        new_data_w = (
2027            new_data_w.sort_values(by=["scan", "mz"], ascending=[True, True])
2028            .reset_index(drop=True)
2029            .copy()
2030        )
2031
2032        return new_data_w
2033
2034    def _filter_data_by_targets(self, data, target_search_dict):
2035        """Filter MS data to only include m/z and RT windows around target values.
2036
2037        Parameters
2038        ----------
2039        data : pd.DataFrame
2040            MS data with 'mz' and 'scan_time' columns
2041        target_search_dict : dict
2042            Dictionary with target_mz_list, target_rt_list, mz_tolerance_ppm, rt_tolerance
2043
2044        Returns
2045        -------
2046        pd.DataFrame
2047            Filtered data containing only points within target windows
2048        """
2049        target_mz_list = target_search_dict["target_mz_list"]
2050        target_rt_list = target_search_dict["target_rt_list"]
2051        mz_tolerance_ppm = target_search_dict["mz_tolerance_ppm"]
2052        rt_tolerance = target_search_dict["rt_tolerance"]
2053
2054        # Create a mask for data points that fall within any target window
2055        mask = np.zeros(len(data), dtype=bool)
2056
2057        for target_mz, target_rt in zip(target_mz_list, target_rt_list):
2058            # Calculate m/z window
2059            mz_tol = target_mz * mz_tolerance_ppm / 1e6
2060            mz_min = target_mz - mz_tol
2061            mz_max = target_mz + mz_tol
2062
2063            # Calculate RT window
2064            rt_min = target_rt - rt_tolerance
2065            rt_max = target_rt + rt_tolerance
2066
2067            # Create mask for this target
2068            target_mask = (
2069                (data["mz"] >= mz_min)
2070                & (data["mz"] <= mz_max)
2071                & (data["scan_time"] >= rt_min)
2072                & (data["scan_time"] <= rt_max)
2073            )
2074
2075            # Combine with overall mask
2076            mask |= target_mask
2077
2078        return data[mask].reset_index(drop=True)
2079
2080    def find_mass_features_ph(
2081        self,
2082        ms_level=1,
2083        grid=True,
2084        targeted_search=False,
2085        target_search_dict=None,
2086        mf_type="untargeted",
2087        accumulate_features=False,
2088    ):
2089        """Find mass features within an LCMSBase object using persistent homology.
2090
2091        Assigns the mass_features attribute to the object (a dictionary of LCMSMassFeature objects, keyed by mass feature id)
2092
2093        Parameters
2094        ----------
2095        ms_level : int, optional
2096            The MS level to use. Default is 1.
2097        grid : bool, optional
2098            If True, will regrid the data before running the persistent homology calculations (after checking if the data is gridded). Default is True.
2099        targeted_search : bool, optional
2100            If True, perform targeted search mode. Default is False.
2101        target_search_dict : dict or None, optional
2102            Dictionary with target parameters for targeted search. Default is None.
2103        mf_type : str, optional
2104            Type label for the mass features. Default is "untargeted".
2105        accumulate_features : bool, optional
2106            If True, add to existing features rather than replacing them. Default is False.
2107
2108        Raises
2109        ------
2110        ValueError
2111            If no MS level data is found on the object.
2112            If data is not gridded and grid is False.
2113
2114        Returns
2115        -------
2116        None, but assigns the mass_features attribute to the object.
2117
2118        Notes
2119        -----
2120        This function has been adapted from the original implementation in the Deimos package: https://github.com/pnnl/deimos
2121        """
2122        # Check that ms_level is a key in self._ms_uprocessed
2123        if ms_level not in self._ms_unprocessed.keys():
2124            raise ValueError(
2125                "No MS level "
2126                + str(ms_level)
2127                + " data found, did you instantiate with parser specific to MS level?"
2128            )
2129
2130        # Get ms data
2131        data = self._ms_unprocessed[ms_level].copy()
2132
2133        # Drop rows with missing intensity values and reset index
2134        data = data.dropna(subset=["intensity"]).reset_index(drop=True)
2135
2136        # Add scan_time for filtering if in targeted mode
2137        if targeted_search:
2138            data = data.merge(
2139                self.scan_df[["scan", "scan_time"]], on="scan", how="left"
2140            )
2141
2142        # Threshold data (bypass thresholds in targeted mode)
2143        dims = ["mz", "scan_time"]
2144        if targeted_search:
2145            # In targeted mode, bypass intensity and persistence thresholds
2146            threshold = 0
2147            persistence_threshold = 0
2148            # Filter data to only target windows
2149            data_thres = self._filter_data_by_targets(data, target_search_dict)
2150            if len(data_thres) == 0:
2151                if self.parameters.lc_ms.verbose_processing:
2152                    print("No data found in target windows")
2153                self.mass_features = {}
2154                return
2155        else:
2156            threshold = max(
2157                self.parameters.lc_ms.ph_inten_min_rel * data.intensity.max(),
2158                self.parameters.lc_ms.ph_inten_min_abs,
2159            )
2160            persistence_threshold = (
2161                self.parameters.lc_ms.ph_persis_min_rel * data.intensity.max()
2162            )
2163            data_thres = (
2164                data[data["intensity"] > threshold].reset_index(drop=True).copy()
2165            )
2166
2167        # One m/z grid for the whole run. A separate grid in each scan window
2168        # moves the apex of the same ion.
2169        if not self.check_if_grid(data_thres):
2170            if grid is False:
2171                raise ValueError(
2172                    "Data are not gridded in mz dimension, try reprocessing with a different params or grid data before running this function"
2173                )
2174            data_thres = self.grid_data(data_thres)
2175        if "scan_time" not in data_thres.columns:
2176            data_thres = data_thres.merge(
2177                self.scan_df[["scan", "scan_time"]], on="scan"
2178            )
2179
2180        return self._find_mass_features_ph_partition(
2181            data_thres,
2182            dims,
2183            persistence_threshold,
2184            mf_type,
2185            accumulate_features,
2186        )
2187
2188    def _find_mass_features_ph_partition(
2189        self,
2190        data_thres,
2191        dims,
2192        persistence_threshold,
2193        mf_type="untargeted",
2194        accumulate_features=False,
2195    ):
2196        """Find mass features on overlapping scan partitions of one m/z grid.
2197
2198        Splits the already gridded input into overlapping scan partitions,
2199        detects mass features with persistent homology, rolls up duplicates
2200        within and across partitions, and populates mass_features.
2201
2202        Parameters
2203        ----------
2204        data_thres : pd.DataFrame
2205            The thresholded input data containing mass spectrometry information.
2206        dims : list
2207            List of dimension names (e.g., ["mz", "scan_time"]) used for feature detection.
2208        persistence_threshold : float
2209            Minimum persistence value required for a detected mass feature to be retained.
2210        mf_type : str, optional
2211            Type label for the mass features. Default is "untargeted".
2212        accumulate_features : bool, optional
2213            If True, add to existing features rather than replacing them. Default is False.
2214
2215        Returns
2216        -------
2217        None
2218            Populates the mass_features attribute of the object with detected mass features.
2219        """
2220        all_mass_features = []
2221
2222        # Split scans into partitions
2223        unique_scans = sorted(data_thres["scan"].unique())
2224        unique_scans_n = len(unique_scans)
2225        if unique_scans_n == 0:
2226            if not accumulate_features:
2227                self.mass_features = {}
2228            return
2229
2230        # Calculate partition size in scans based on goal
2231        partition_size_goal = 5000
2232        scans_per_partition = max(
2233            1, partition_size_goal // (len(data_thres) // unique_scans_n)
2234        )
2235
2236        # Make partitions based on scans, with overlapping in partitioned scans
2237        scan_overlap = 4
2238        partition_scans = []
2239        for i in range(0, unique_scans_n, scans_per_partition):
2240            start_idx = max(0, i - scan_overlap)
2241            end_idx = min(
2242                unique_scans_n - 1, i + scans_per_partition - 1 + scan_overlap
2243            )
2244            scans_group = [int(s) for s in unique_scans[start_idx : end_idx + 1]]
2245            partition_scans.append(scans_group)
2246
2247        # Set index to scan for faster filtering
2248        data_thres = data_thres.set_index("scan")
2249        for scans in partition_scans:
2250            # Determine start and end scan for partition, with 5 scans overlap
2251            partition_data = data_thres.loc[scans].reset_index(drop=False).copy()
2252
2253            if len(partition_data) == 0:
2254                continue
2255
2256            if "scan_time" not in partition_data.columns:
2257                partition_data = partition_data.merge(
2258                    self.scan_df[["scan", "scan_time"]], on="scan"
2259                )
2260
2261            # Build factors for this partition
2262            factors = {
2263                dim: pd.factorize(partition_data[dim], sort=True)[1].astype(np.float32)
2264                for dim in dims
2265            }
2266
2267            # Build indexes
2268            index = {
2269                dim: np.searchsorted(factors[dim], partition_data[dim]).astype(
2270                    np.float32
2271                )
2272                for dim in factors
2273            }
2274
2275            # Process partition
2276            partition_features = self._process_partition_ph(
2277                partition_data, index, dims, persistence_threshold
2278            )
2279
2280            if len(partition_features) == 0:
2281                continue
2282
2283            # Roll up within partition
2284            partition_features = self.roll_up_dataframe(
2285                df=partition_features,
2286                sort_by="persistence",
2287                dims=["mz", "scan_time"],
2288                tol=[
2289                    self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel,
2290                    self.parameters.lc_ms.mass_feature_cluster_rt_tolerance,
2291                ],
2292                relative=[True, False],
2293            )
2294            partition_features = partition_features.reset_index(drop=True)
2295
2296            if len(partition_features) > 0:
2297                all_mass_features.append(partition_features)
2298
2299        # Combine results from all partitions
2300        if all_mass_features:
2301            combined_features = pd.concat(all_mass_features, ignore_index=True)
2302
2303            # Sort by persistence
2304            combined_features = combined_features.sort_values(
2305                by="persistence", ascending=False
2306            ).reset_index(drop=True)
2307
2308            # Remove duplicates from overlapping regions
2309            combined_features = self.roll_up_dataframe(
2310                df=combined_features,
2311                sort_by="persistence",
2312                dims=["mz", "scan_time"],
2313                tol=[
2314                    self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel,
2315                    self.parameters.lc_ms.mass_feature_cluster_rt_tolerance,
2316                ],
2317                relative=[True, False],
2318            )
2319
2320            # resort by persistence and reset index
2321            combined_features = combined_features.reset_index(drop=True)
2322
2323            # Populate mass_features attribute
2324            self._populate_mass_features(
2325                combined_features, mf_type, accumulate_features
2326            )
2327        else:
2328            self.mass_features = {}
2329
2330    def _process_partition_ph(self, partition_data, index, dims, persistence_threshold):
2331        """Process a single partition with persistent homology."""
2332        # Smooth data
2333        iterations = self.parameters.lc_ms.ph_smooth_it
2334        smooth_radius = [
2335            self.parameters.lc_ms.ph_smooth_radius_mz,
2336            self.parameters.lc_ms.ph_smooth_radius_scan,
2337        ]
2338
2339        index_array = np.vstack([index[dim] for dim in dims]).T
2340        V = partition_data["intensity"].values
2341        resid = np.inf
2342
2343        for i in range(iterations):
2344            # Previous iteration
2345            V_prev = V.copy()
2346            resid_prev = resid
2347            V = self.sparse_mean_filter(index_array, V, radius=smooth_radius)
2348
2349            # Calculate residual with previous iteration
2350            resid = np.sqrt(np.mean(np.square(V - V_prev)))
2351
2352            # Evaluate convergence
2353            if i > 0:
2354                # Percent change in residual
2355                test = np.abs(resid - resid_prev) / resid_prev
2356
2357                # Exit criteria
2358                if test <= 0:
2359                    break
2360
2361        # Overwrite values
2362        partition_data = partition_data.copy()
2363        partition_data["intensity"] = V
2364
2365        # Use persistent homology to find regions of interest
2366        pidx, pers = self.sparse_upper_star(index_array, V)
2367        pidx = pidx[pers > 1]
2368        pers = pers[pers > 1]
2369
2370        if len(pidx) == 0:
2371            return pd.DataFrame()
2372
2373        # Get peaks
2374        peaks = partition_data.iloc[pidx, :].reset_index(drop=True)
2375
2376        # Add persistence column
2377        peaks["persistence"] = pers
2378        mass_features = peaks.sort_values(
2379            by="persistence", ascending=False
2380        ).reset_index(drop=True)
2381
2382        # Filter by persistence threshold
2383        mass_features = mass_features.loc[
2384            mass_features["persistence"] > persistence_threshold, :
2385        ].reset_index(drop=True)
2386
2387        return mass_features
2388
2389    def _populate_mass_features(
2390        self, mass_features_df, mf_type="untargeted", accumulate_features=False
2391    ):
2392        """Populate the mass_features attribute from a DataFrame.
2393
2394        Parameters
2395        ----------
2396        mass_features_df : pd.DataFrame
2397            DataFrame containing mass feature information.
2398            Note that the order of this DataFrame will determine the order of mass features in the mass_features attribute.
2399        mf_type : str, optional
2400            Type label for the mass features. Default is "untargeted".
2401        accumulate_features : bool, optional
2402            If True, new features will be added to existing features rather than replacing them.
2403            Mass feature IDs will be offset to avoid conflicts. Default is False.
2404
2405        Returns
2406        -------
2407        None, but assigns or updates the mass_features attribute to the object.
2408        """
2409        # Rename scan column to apex_scan
2410        mass_features_df = mass_features_df.rename(
2411            columns={"scan": "apex_scan", "scan_time": "retention_time"}
2412        )
2413
2414        # Initialize or preserve existing mass_features attribute
2415        if (
2416            accumulate_features
2417            and self.mass_features is not None
2418            and len(self.mass_features) > 0
2419        ):
2420            # Find the maximum existing ID to offset new IDs and avoid conflicts
2421            id_offset = max(self.mass_features.keys()) + 1
2422            initial_count = len(self.mass_features)
2423        else:
2424            # Replace mode (default/backwards compatible)
2425            self.mass_features = {}
2426            id_offset = 0
2427            initial_count = 0
2428
2429        # Add new mass features
2430        for idx, row in enumerate(mass_features_df.itertuples()):
2431            row_dict = mass_features_df.iloc[row.Index].to_dict()
2432            lcms_feature = LCMSMassFeature(self, **row_dict)
2433            lcms_feature.type = mf_type
2434            # Use sequential ID starting from id_offset to avoid conflicts with existing features
2435            new_id = idx + id_offset
2436            lcms_feature._id = new_id  # Update the internal ID
2437            self.mass_features[new_id] = lcms_feature
2438
2439        if self.parameters.lc_ms.verbose_processing:
2440            if accumulate_features and initial_count > 0:
2441                print(
2442                    f"Found {len(mass_features_df)} new mass features (total: {len(self.mass_features)})"
2443                )
2444            else:
2445                print("Found " + str(len(mass_features_df)) + " initial mass features")
2446
2447    def find_mass_features_ph_centroid(
2448        self,
2449        ms_level=1,
2450        targeted_search=False,
2451        target_search_dict=None,
2452        mf_type="untargeted",
2453        accumulate_features=False,
2454    ):
2455        """Find mass features within an LCMSBase object using persistent homology-type approach but with centroided data.
2456
2457        Parameters
2458        ----------
2459        ms_level : int, optional
2460            The MS level to use. Default is 1.
2461        targeted_search : bool, optional
2462            If True, perform targeted search mode. Default is False.
2463        target_search_dict : dict or None, optional
2464            Dictionary with target parameters for targeted search. Default is None.
2465        mf_type : str, optional
2466            Type label for the mass features. Default is "untargeted".
2467        accumulate_features : bool, optional
2468            If True, add to existing features rather than replacing them. Default is False.
2469
2470        Raises
2471        ------
2472        ValueError
2473            If no MS level data is found on the object.
2474
2475        Returns
2476        -------
2477        None, but assigns the mass_features attribute to the object.
2478        """
2479        # Check that ms_level is a key in self._ms_uprocessed
2480        if ms_level not in self._ms_unprocessed.keys():
2481            raise ValueError(
2482                "No MS level "
2483                + str(ms_level)
2484                + " data found, did you instantiate with parser specific to MS level?"
2485            )
2486
2487        # Work with reference instead of copy
2488        data = self._ms_unprocessed[ms_level]
2489
2490        # Merge with scan data first (needed for filtering in targeted mode)
2491        scan_subset = self.scan_df[["scan", "scan_time"]]
2492        data_with_time = data.merge(scan_subset, on="scan", how="inner")
2493
2494        # Calculate threshold and filter (bypass in targeted mode)
2495        if targeted_search:
2496            # In targeted mode, bypass intensity threshold
2497            threshold = 0
2498            valid_mask = data_with_time["intensity"].notna()
2499            required_cols = ["mz", "intensity", "scan", "scan_time"]
2500            data_thres = data_with_time.loc[valid_mask, required_cols].copy()
2501
2502            # Filter to target windows
2503            data_thres = self._filter_data_by_targets(data_thres, target_search_dict)
2504
2505            if len(data_thres) == 0:
2506                if self.parameters.lc_ms.verbose_processing:
2507                    print("No data found in target windows")
2508                self.mass_features = {}
2509                return
2510        else:
2511            # Normal mode with threshold
2512            max_intensity = data_with_time["intensity"].max()
2513            threshold = max(
2514                self.parameters.lc_ms.ph_inten_min_rel * max_intensity,
2515                self.parameters.lc_ms.ph_inten_min_abs,
2516            )
2517            valid_mask = data_with_time["intensity"].notna() & (
2518                data_with_time["intensity"] > threshold
2519            )
2520            required_cols = ["mz", "intensity", "scan", "scan_time"]
2521            data_thres = data_with_time.loc[valid_mask, required_cols].copy()
2522
2523        data_thres["persistence"] = data_thres["intensity"]
2524        mf_df = data_thres
2525        del data_thres, scan_subset, data_with_time
2526
2527        # Order by scan_time and then mz to ensure features near in rt are processed together
2528        # It's ok that different scans are in different partitions; we will roll up later
2529        mf_df = mf_df.sort_values(
2530            by=["scan_time", "mz"], ascending=[True, True]
2531        ).reset_index(drop=True)
2532        partition_size = 10000
2533        partitions = [
2534            mf_df.iloc[i : i + partition_size].reset_index(drop=True)
2535            for i in range(0, len(mf_df), partition_size)
2536        ]
2537        del mf_df
2538
2539        # Run roll_up_dataframe on each partition
2540        rolled_partitions = []
2541        for part in partitions:
2542            rolled = self.roll_up_dataframe(
2543                df=part,
2544                sort_by="persistence",
2545                dims=["mz", "scan_time"],
2546                tol=[
2547                    self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel,
2548                    self.parameters.lc_ms.mass_feature_cluster_rt_tolerance,
2549                ],
2550                relative=[True, False],
2551            )
2552            rolled_partitions.append(rolled)
2553        del partitions
2554
2555        # Run roll_up_dataframe on the rolled_up partitions to merge features near partition boundaries
2556
2557        # Combine results and run a final roll-up to merge features near partition boundaries
2558        mf_df_final = pd.concat(rolled_partitions, ignore_index=True)
2559        del rolled_partitions
2560
2561        # Reorder by persistence before final roll-up
2562        mf_df_final = mf_df_final.sort_values(
2563            by="persistence", ascending=False
2564        ).reset_index(drop=True)
2565
2566        mf_df_final = self.roll_up_dataframe(
2567            df=mf_df_final,
2568            sort_by="persistence",
2569            dims=["mz", "scan_time"],
2570            tol=[
2571                self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel,
2572                self.parameters.lc_ms.mass_feature_cluster_rt_tolerance,
2573            ],
2574            relative=[True, False],
2575        )
2576        # reset index
2577        mf_df_final = mf_df_final.reset_index(drop=True)
2578
2579        # Combine rename and sort operations
2580        mass_features = (
2581            mf_df_final.rename(
2582                columns={"scan": "apex_scan", "scan_time": "retention_time"}
2583            )
2584            .sort_values(by="persistence", ascending=False)
2585            .reset_index(drop=True)
2586        )
2587        del mf_df_final  # Free memory
2588
2589        # Order by persistence and reset index
2590        mass_features = mass_features.sort_values(
2591            by="persistence", ascending=False
2592        ).reset_index(drop=True)
2593
2594        self.mass_features = {}
2595        for idx, row in mass_features.iterrows():
2596            row_dict = row.to_dict()
2597            lcms_feature = LCMSMassFeature(self, **row_dict)
2598            lcms_feature.type = mf_type
2599            self.mass_features[lcms_feature.id] = lcms_feature
2600
2601        if self.parameters.lc_ms.verbose_processing:
2602            print("Found " + str(len(mass_features)) + " initial mass features")
2603
2604    def cluster_mass_features(self, drop_children=True, sort_by="persistence"):
2605        """Cluster mass features
2606
2607        Based on their proximity in the mz and scan_time dimensions, priorizies the mass features with the highest persistence.
2608
2609        Parameters
2610        ----------
2611        drop_children : bool, optional
2612            Whether to drop the mass features that are not cluster parents. Default is True.
2613        sort_by : str, optional
2614            The column to sort the mass features by, this will determine which mass features get rolled up into a parent mass feature. Default is "persistence".
2615
2616        Raises
2617        ------
2618        ValueError
2619            If no mass features are found.
2620            If too many mass features are found.
2621
2622        Returns
2623        -------
2624        None if drop_children is True, otherwise returns a list of mass feature ids that are not cluster parents.
2625        """
2626        if self.mass_features is None:
2627            raise ValueError("No mass features found, run find_mass_features() first")
2628        if len(self.mass_features) > 400000:
2629            raise ValueError(
2630                "Too many mass features of interest found, run find_mass_features() with a higher intensity threshold"
2631            )
2632        dims = ["mz", "scan_time"]
2633        mf_df_og = self.mass_features_to_df()
2634        mf_df = mf_df_og.copy()
2635
2636        tol = [
2637            self.parameters.lc_ms.mass_feature_cluster_mz_tolerance_rel,
2638            self.parameters.lc_ms.mass_feature_cluster_rt_tolerance,
2639        ]  # mz, in relative; scan_time in minutes
2640        relative = [True, False]
2641
2642        # Roll up mass features based on their proximity in the declared dimensions
2643        mf_df_new = self.roll_up_dataframe(
2644            df=mf_df, sort_by=sort_by, dims=dims, tol=tol, relative=relative
2645        )
2646
2647        mf_df["cluster_parent"] = np.where(
2648            np.isin(mf_df.index, mf_df_new.index), True, False
2649        )
2650
2651        # get mass feature ids of features that are not cluster parents
2652        cluster_daughters = mf_df[~mf_df["cluster_parent"]].index.values
2653        if drop_children is True:
2654            # Drop mass features that are not cluster parents from self
2655            self.mass_features = {
2656                k: v
2657                for k, v in self.mass_features.items()
2658                if k not in cluster_daughters
2659            }
2660        else:
2661            return cluster_daughters
2662
2663
2664class LCMSCollectionCalculations:
2665    """Methods for performing calculations related to LCMSCollection objects.
2666
2667    Notes
2668    -----
2669    This class is intended as a mixin for the LCMSCollection class.
2670    """
2671
2672    @staticmethod
2673    def _resolve_eic_query_mz(row):
2674        """Pick m/z used to look up ``sample.eics`` for a mass-feature row.
2675
2676        Prefer ``_eic_mz`` when present (set after integrate / gap-fill).
2677        Regular features in ``mass_features_dataframe`` often have NaN
2678        ``_eic_mz`` (never back-filled from objects on light load / reload),
2679        while induced features get ``_eic_mz`` written after gap-fill. Fall
2680        back to the feature ``mz`` so regular and representative traces are
2681        not silently dropped from multi-sample EIC plots.
2682        """
2683        if row is None:
2684            return None
2685        eic_mz = row.get("_eic_mz") if hasattr(row, "get") else None
2686        if eic_mz is not None and not pd.isna(eic_mz):
2687            return eic_mz
2688        mz = row.get("mz") if hasattr(row, "get") else None
2689        if mz is not None and not pd.isna(mz):
2690            return mz
2691        return None
2692
2693    @staticmethod
2694    def _get_eic_data_for_mz(sample, eic_mz, tolerance=0.0001):
2695        """Resolve EIC data for an m/z using exact key, then tolerance match.
2696
2697        Plotting previously used ``sample.eics.get(eic_mz)`` only. HDF5 EIC
2698        keys and feature ``_eic_mz`` / ``mz`` often differ by float noise, which
2699        silently dropped traces (including the representative). Elsewhere
2700        CoreMS uses the same default tolerance via
2701        ``get_eic_mz_for_mass_feature`` / ``associate_eics_with_mass_features``.
2702
2703        Parameters
2704        ----------
2705        sample : LCMSBase
2706            Sample that holds ``eics`` (dict keyed by m/z).
2707        eic_mz : float or None
2708            Target m/z (typically from the mass-feature ``_eic_mz`` or ``mz``).
2709        tolerance : float, optional
2710            Maximum |Δm/z| for fallback matching. Default is 0.0001 Da.
2711
2712        Returns
2713        -------
2714        EIC_Data or None
2715            Matching EIC data, or None if no key is within tolerance.
2716        """
2717        if eic_mz is None or pd.isna(eic_mz):
2718            return None
2719        if not hasattr(sample, "eics") or not sample.eics:
2720            return None
2721
2722        eic_data = sample.eics.get(eic_mz)
2723        if eic_data is not None:
2724            return eic_data
2725
2726        # Exact key miss: tolerance match (same default as associate_eics_with_mass_features)
2727        if hasattr(sample, "get_eic_mz_for_mass_feature"):
2728            matched_mz = sample.get_eic_mz_for_mass_feature(
2729                float(eic_mz), tolerance=tolerance
2730            )
2731            if matched_mz is not None:
2732                return sample.eics.get(matched_mz)
2733            return None
2734
2735        # Fallback if sample is a simple mock without the helper
2736        best_key = min(sample.eics, key=lambda k: abs(float(k) - float(eic_mz)))
2737        if abs(float(best_key) - float(eic_mz)) < tolerance:
2738            return sample.eics[best_key]
2739        return None
2740
2741    @staticmethod
2742    def _plot_multiple_eics(
2743        ax,
2744        cluster_mfs,
2745        induced_cluster_mfs,
2746        rep_sample_id,
2747        rep_mf_id,
2748        median_rt,
2749        eic_buffer_time,
2750        plot_smoothed=False,
2751        plot_datapoints=False,
2752        label_samples=False,
2753        lcms_collection=None,
2754    ):
2755        """Internal method to plot multiple EICs from different samples on a given axis.
2756
2757        Parameters
2758        ----------
2759        ax : matplotlib.axes.Axes
2760            The axis to plot on.
2761        cluster_mfs : pd.DataFrame
2762            DataFrame containing cluster mass features (non-induced).
2763        induced_cluster_mfs : pd.DataFrame or None
2764            DataFrame containing induced (gap-filled) mass features.
2765        rep_sample_id : int
2766            Sample ID of the representative mass feature.
2767        rep_mf_id : int
2768            Mass feature ID of the representative mass feature.
2769        median_rt : float
2770            Median retention time for the cluster.
2771        eic_buffer_time : float
2772            Time buffer around the peak (minutes).
2773        plot_smoothed : bool, optional
2774            If True, plot smoothed EICs. Default is False.
2775        plot_datapoints : bool, optional
2776            If True, plot EIC datapoints. Default is False.
2777        label_samples : bool, optional
2778            If True, label each sample individually. Default is False.
2779        lcms_collection : LCMSCollection, optional
2780            The parent collection object for accessing samples. Required.
2781        """
2782        ax.set_title("EICs from all samples", loc="left")
2783
2784        # Track if we've added labels for legend (to avoid duplicates)
2785        rep_labeled = False
2786        regular_labeled = False
2787        induced_labeled = False
2788
2789        # Plot regular (non-induced) mass features
2790        for _, row in cluster_mfs.iterrows():
2791            sample_id = int(row["sample_id"])
2792            mf_id = row["mf_id"]
2793            sample = lcms_collection[sample_id]
2794            sample_name = row["sample_name"]
2795
2796            # Prefer _eic_mz, else feature mz; then exact/tolerance key match (#257)
2797            eic_data = LCMSCollectionCalculations._get_eic_data_for_mz(
2798                sample, LCMSCollectionCalculations._resolve_eic_query_mz(row)
2799            )
2800
2801            if eic_data:
2802                # Determine line style and width
2803                if sample_id == rep_sample_id and mf_id == rep_mf_id:
2804                    # Representative feature - bold line
2805                    linewidth = 2.5
2806                    alpha = 1.0
2807                    color = "tab:blue"
2808                    if label_samples:
2809                        label = f"{sample_name} (representative)"
2810                    else:
2811                        label = "Representative" if not rep_labeled else None
2812                        rep_labeled = True
2813                else:
2814                    # Other features - thinner line
2815                    linewidth = 1.0
2816                    alpha = 0.5
2817                    color = "tab:blue"
2818                    if label_samples:
2819                        label = sample_name
2820                    else:
2821                        label = "Regular features" if not regular_labeled else None
2822                        regular_labeled = True
2823
2824                ax.plot(
2825                    eic_data.time,
2826                    eic_data.eic,
2827                    c=color,
2828                    linewidth=linewidth,
2829                    alpha=alpha,
2830                    linestyle="-",
2831                    label=label,
2832                )
2833
2834                if plot_datapoints:
2835                    ax.scatter(eic_data.time, eic_data.eic, c=color, alpha=alpha, s=10)
2836
2837                if plot_smoothed and hasattr(eic_data, "eic_smoothed"):
2838                    ax.plot(
2839                        eic_data.time,
2840                        eic_data.eic_smoothed,
2841                        c=color,
2842                        linestyle="--",
2843                        alpha=alpha * 0.8,
2844                        linewidth=linewidth * 0.8,
2845                    )
2846
2847        # Plot induced (gap-filled) mass features if available
2848        if induced_cluster_mfs is not None and not induced_cluster_mfs.empty:
2849            for _, row in induced_cluster_mfs.iterrows():
2850                sample_id = int(row["sample_id"])
2851                mf_id = row["mf_id"]
2852                sample = lcms_collection[sample_id]
2853                sample_name = row["sample_name"]
2854
2855                # Prefer _eic_mz, else feature mz; then exact/tolerance key match (#257)
2856                eic_data = LCMSCollectionCalculations._get_eic_data_for_mz(
2857                    sample, LCMSCollectionCalculations._resolve_eic_query_mz(row)
2858                )
2859
2860                if eic_data:
2861                    # Induced features - even thinner line
2862                    linewidth = 0.5
2863                    alpha = 0.4
2864                    color = "tab:orange"
2865
2866                    if label_samples:
2867                        label = f"{sample_name} (induced)"
2868                    else:
2869                        label = "Gap-filled features" if not induced_labeled else None
2870                        induced_labeled = True
2871
2872                    ax.plot(
2873                        eic_data.time,
2874                        eic_data.eic,
2875                        c=color,
2876                        linewidth=linewidth,
2877                        alpha=alpha,
2878                        linestyle="-",
2879                        label=label,
2880                    )
2881
2882                    if plot_datapoints:
2883                        ax.scatter(
2884                            eic_data.time, eic_data.eic, c=color, alpha=alpha, s=5
2885                        )
2886
2887                    if plot_smoothed and hasattr(eic_data, "eic_smoothed"):
2888                        ax.plot(
2889                            eic_data.time,
2890                            eic_data.eic_smoothed,
2891                            c=color,
2892                            linestyle="--",
2893                            alpha=alpha * 0.8,
2894                            linewidth=linewidth * 0.8,
2895                        )
2896
2897        # Add vertical line at median RT
2898        ax.axvline(x=median_rt, color="k", linestyle="--", alpha=0.7, label="Median RT")
2899
2900        ax.set_ylabel("Intensity")
2901        ax.set_xlabel("Time (minutes)")
2902        ax.set_xlim(
2903            median_rt - eic_buffer_time,
2904            median_rt + eic_buffer_time,
2905        )
2906        ax.legend(loc="upper left", fontsize=8)
2907        ax.yaxis.get_major_formatter().set_useOffset(False)
2908
2909    def clean_sparse_matrix(self, sparse_matrix):
2910        """Clean a sparse matrix by removing duplicates and sorting.
2911
2912        Parameters
2913        ----------
2914        sparse_matrix : :obj:`~numpy.array`
2915            A sparse matrix to clean.
2916
2917        Returns
2918        -------
2919        :obj:`~numpy.array`
2920            A cleaned sparse matrix.
2921        """
2922        for match in sparse_matrix:
2923            match.sort()
2924        sparse_matrix.sort()
2925        dereplicated_sparse_matrix = np.unique(sparse_matrix, axis=0)
2926        return dereplicated_sparse_matrix
2927
2928    def match_mfs(self, mf_c, mf_i):
2929        """Match mass features between two LCMS objects.
2930
2931        Parameters
2932        ----------
2933        mf_c : :obj:`~pandas.DataFrame`
2934            The mass features to match against.
2935        mf_i : :obj:`~pandas.DataFrame`
2936            The mass features to match.
2937
2938        Returns
2939        -------
2940        :obj:`~pandas.DataFrame`
2941            The matched mass features from mf_c.
2942        :obj:`~pandas.DataFrame`
2943            The matched mass features from mf_i.
2944
2945        Notes
2946        -----
2947        This function has been adapted from the original implementation in the Deimos package:
2948        https://github.com/pnnl/deimos
2949        """
2950        if mf_c is None or mf_i is None or len(mf_c.index) < 1 or len(mf_i.index) < 1:
2951            return None, None
2952
2953        # Prepare dataframes
2954        mf_c = mf_c.copy()
2955        mf_c["id_i"] = 0
2956        mf_i = mf_i.copy()
2957        mf_i["id_i"] = 1
2958
2959        # Set dimensions for matching
2960        dims = ["mz", "scan_time"]
2961        relative = [True, False]
2962        mz_tol = self.parameters.lcms_collection.alignment_mz_tol_ppm * 1e-6
2963        rt_tol = self.parameters.lcms_collection.alignment_rt_tol
2964        tol = [mz_tol, rt_tol]
2965
2966        # Compute inter-feature distances
2967        idx = []
2968        for i, f in enumerate(dims):
2969            # vectors
2970            v1 = mf_c[f].values.reshape(-1, 1)
2971            v2 = mf_i[f].values.reshape(-1, 1)
2972
2973            # Distances
2974            d = scipy.spatial.distance.cdist(v1, v2)
2975
2976            if relative[i] is True:
2977                # Divisor
2978                basis = np.repeat(v1, v2.shape[0], axis=1)
2979                fix = np.repeat(v2, v1.shape[0], axis=1).T
2980                basis = np.where(basis == 0, fix, basis)
2981
2982                # Divide
2983                d = np.divide(d, basis, out=np.zeros_like(basis), where=basis != 0)
2984
2985            # Check tol
2986            idx.append(d <= tol[i])
2987
2988        # Stack truth arrays
2989        idx = np.prod(np.dstack(idx), axis=-1, dtype=bool)
2990
2991        # Compute normalized 3d distance
2992        v1 = mf_c[dims].values / tol
2993        v2 = mf_i[dims].values / tol
2994        dist3d = scipy.spatial.distance.cdist(v1, v2, "cityblock")
2995
2996        # Separate features within tolerance from those outside
2997        # Features outside tolerance should be inf, features within tolerance keep their distance
2998        # Use idx mask: True for within tolerance, False for outside
2999        dist3d_within_tol = np.where(idx, dist3d, np.inf)
3000
3001        # Normalize to 0-1 (only affects within-tolerance distances)
3002        mx = np.max(dist3d_within_tol[idx]) if np.sum(idx) > 0 else 0
3003        if mx > 0:
3004            # Lower distance is better - normalize only the within-tolerance values
3005            dist3d_within_tol = np.where(idx, dist3d_within_tol / mx, np.inf)
3006        else:
3007            # All matches are perfect (distance=0), assign tiny value to within-tolerance pairs
3008            dist3d_within_tol = np.where(idx, 1e-10, np.inf)
3009
3010        # Use the masked distance matrix
3011        dist3d = dist3d_within_tol
3012
3013        # Min over dims
3014        mincols = np.min(dist3d, axis=0, keepdims=True)
3015
3016        # Zero out mincols over dims
3017        dist3d[dist3d != mincols] = np.inf
3018
3019        # Min over clusters
3020        minrows = np.min(dist3d, axis=1, keepdims=True)
3021
3022        # Where max and nonzero
3023        ii, jj = np.where((dist3d == minrows) & (dist3d < np.inf))
3024
3025        # Reorder
3026        mf_c = mf_c.iloc[ii]
3027        mf_i = mf_i.iloc[jj]
3028
3029        if len(mf_c.index) < 1 or len(mf_i.index) < 1:
3030            return None, None
3031
3032        return mf_c, mf_i
3033
3034    def fit_rts(self, a, b, align="scan_time", **kwargs):
3035        """
3036        Fit a support vector regressor to matched features.
3037
3038        Parameters
3039        ----------
3040        a : :obj:`~pandas.DataFrame`
3041            First set of input feature coordinates and intensities; the center object and the object to align to.
3042        b : :obj:`~pandas.DataFrame`
3043            Second set of input feature coordinates and intensities; the object to align to the center object.
3044        align : str
3045            Dimension to align.
3046        kwargs
3047            Keyword arguments for support vector regressor
3048            (:class:`sklearn.svm.SVR`).
3049
3050        Returns
3051        -------
3052        :obj:`~function`
3053            An interpolation function where one can input a retention time and get the predicted retention time.
3054
3055        Notes
3056        -----
3057        This function has been adapted from the original implementation in the Deimos package:
3058        https://github.com/pnnl/deimos
3059
3060        """
3061
3062        # Uniqueify
3063        x = a[align].values
3064        y = b[align].values
3065        arr = np.vstack((x, y)).T
3066        arr = np.unique(arr, axis=0)
3067
3068        # Safety check: ensure we have data to work with
3069        if len(arr) == 0:
3070            warnings.warn(
3071                "No data points available for retention time fitting. Returning identity function."
3072            )
3073            return lambda x: x
3074
3075        # Check kwargs
3076        if "kernel" in kwargs:
3077            kernel = kwargs.get("kernel")
3078        else:
3079            kernel = "linear"
3080
3081        # Construct interpolation axis
3082        newx = np.linspace(arr[:, 0].min(), arr[:, 0].max(), 1000)
3083
3084        # Linear kernel
3085        if kernel == "linear":
3086            reg = scipy.stats.linregress(x, y)
3087            newy = reg.slope * newx + reg.intercept
3088
3089        # Other kernels
3090        else:
3091            # Fit
3092            svr = SVR(**kwargs)
3093            svr.fit(arr[:, 1].reshape(-1, 1), arr[:, 0])
3094
3095            # Predict
3096            newy = svr.predict(newx.reshape(-1, 1))
3097
3098        # Pad x and y_pred with zeros to force interpolation to start at 0
3099        newx = np.concatenate(([0], newx))
3100        newy = np.concatenate(([0], newy))
3101
3102        # Pad x and y_pred with max time to force interpolation to end at max time to force interpolation to match at end max time
3103        max_time = self[0].scan_df["scan_time"].max()
3104        newx = np.concatenate((newx, [max_time]))
3105        newy = np.concatenate((newy, [max_time]))
3106
3107        # Return an interpolation function for the x and y_pred
3108        def interp(x):
3109            pred_y = np.interp(x, newx, newy)
3110            return pred_y
3111
3112        return interp
3113
3114    def get_anchor_mass_features(self, mf_df):
3115        """
3116        Get the anchor mass features from a DataFrame of mass features.
3117
3118        Parameters
3119        ----------
3120        mf_df : :obj:`~pandas.DataFrame`
3121            The mass features to filter to just the anchor mass features.
3122
3123        Returns
3124        -------
3125        :obj:`~pandas.DataFrame`
3126            The anchor mass features dataframe.
3127        """
3128        mf_df = mf_df.copy()
3129
3130        if (
3131            "deconvoluted_mass_spectra"
3132            in self.parameters.lcms_collection.mass_feature_anchor_technique
3133        ):
3134            # Drop features that are not mass_spectrum_deconvoluted_parent or are NA as mass_spectrum_deconvoluted_parent
3135            mf_df = mf_df.dropna(subset=["mass_spectrum_deconvoluted_parent"])
3136            mf_df = mf_df[mf_df["mass_spectrum_deconvoluted_parent"]]
3137
3138        if (
3139            "absolute_intensity"
3140            in self.parameters.lcms_collection.mass_feature_anchor_technique
3141        ):
3142            # Drop features that have an intensity lower than the threshold
3143            threshold = self.parameters.lcms_collection.mass_feature_anchor_absolute_intensity_threshold
3144            mf_df = mf_df[mf_df["intensity"] > threshold]
3145
3146        if (
3147            "relative_intensity"
3148            in self.parameters.lcms_collection.mass_feature_anchor_technique
3149        ):
3150            # Drop features in the lower fraction of intensities
3151            threshold_quantile = self.parameters.lcms_collection.mass_feature_anchor_relative_intensity_threshold
3152            intensity_threshold = mf_df["intensity"].quantile(threshold_quantile)
3153            mf_df = mf_df[mf_df["intensity"] >= intensity_threshold]
3154
3155        return mf_df
3156
3157    def attempt_alignment(self, matches_c, matches_i):
3158        """
3159        Check if alignment is needed for the LCMS objects in the collection.
3160        """
3161
3162        # Hold out a subset of matches_c and matches_i for spline fitting
3163        matches_c.reset_index(drop=False, inplace=True)
3164        matches_i.reset_index(drop=False, inplace=True)
3165
3166        # Check if there are enough matches to attempt alignment
3167        minimum_matches = self.parameters.lcms_collection.alignment_minimum_matches
3168        if len(matches_c) < minimum_matches:
3169            # Return False (no alignment) and identity function (returns original time)
3170            # which isn't used but is a placeholder to avoid errors in downstream code since
3171            # the function expects a callable to be returned
3172            return False, lambda x: x
3173
3174        # Rearrange matches_c and matches_i to be in the order of the scan_time of matches_c
3175        matches_c = matches_c.sort_values(by="scan_time")
3176        matches_i = matches_i.iloc[matches_c.index.values]
3177
3178        hold_out_fraction = self.parameters.lcms_collection.alignment_hold_out_fraction
3179        # starting with an array of length len(matches_c), select equally spaced indices to hold out
3180        idx_holdout = matches_c.index.values[
3181            np.arange(0, len(matches_c), int(1 / hold_out_fraction))
3182        ]
3183
3184        matches_c_holdout = matches_c.loc[idx_holdout].copy()
3185        matches_i_holdout = matches_i.loc[idx_holdout].copy()
3186
3187        # Remove the holdout matches from the matches_c and matches_i DataFrames and reset the index
3188        matches_c = matches_c.drop(index=idx_holdout).set_index("sample_name")
3189        matches_i = matches_i.drop(index=idx_holdout).set_index("sample_name")
3190
3191        # Reset the scan_time to the original scan_time
3192        matches_i = matches_i.copy()
3193        matches_i["scan_time"] = matches_i["scan_time_og"]
3194
3195        # Fit the retention times of the LCMS object to the center LCMS object using the matched mass features
3196        spl = self.fit_rts(matches_c, matches_i, kernel="rbf", C=1000)
3197
3198        # Check if the spline fitting improved the alignment for the holdout matches
3199        matches_i_holdout["scan_time_fit"] = spl(matches_i_holdout["scan_time"])
3200        og_diff = np.abs(
3201            matches_i_holdout["scan_time"] - matches_c_holdout["scan_time"]
3202        )
3203        fit_diff = np.abs(
3204            matches_i_holdout["scan_time_fit"] - matches_c_holdout["scan_time"]
3205        )
3206
3207        if (
3208            "fraction_improved"
3209            in self.parameters.lcms_collection.alignment_acceptance_technique
3210        ):
3211            fraction_improved = np.sum(fit_diff < og_diff) / len(og_diff)
3212            use_spline_alignment = (
3213                fraction_improved
3214                > self.parameters.lcms_collection.alignment_acceptance_fraction_improved_threshold
3215            )
3216        if (
3217            "mean_squared_error_improved"
3218            in self.parameters.lcms_collection.alignment_acceptance_technique
3219        ):
3220            mse_og = np.mean(og_diff**2)
3221            mse = np.mean(fit_diff**2)
3222            use_spline_alignment = mse < mse_og
3223            # Convert to boolean
3224            use_spline_alignment = bool(use_spline_alignment)
3225
3226        return use_spline_alignment, spl
3227
3228    def align_lcms_objects(self, overwrite=False):
3229        """
3230        Align LCMS objects in the collection.
3231
3232        Aligns the LCMS objects in the collection by aligning the retention times of the mass features in the LCMS objects.
3233        First, the mass features in the center LCMS object are matched to the mass features in the other LCMS objects,
3234        starting with the LCMS object immediately following the center LCMS object. The retention times of the LCMS objects
3235        are then fit to the center LCMS object using the matched mass features.
3236
3237        Returns
3238        -------
3239        None, but aligns the LCMS objects in the collection and sets the scan_time_aligned column in the scan_df attribute of each LCMS object.
3240
3241        Notes
3242        -----
3243        This function has been adapted from the original implementation in the Deimos package:
3244        https://github.com/pnnl/deimos
3245        """
3246
3247        # Prepare the center LCMS object
3248        center_obj_ids = self.manifest_dataframe[
3249            self.manifest_dataframe["center"]
3250        ].collection_id.values
3251
3252        full_mf_df = self.mass_features_dataframe
3253        # re-index to sample_name for faster lookups
3254        full_mf_df = full_mf_df.reset_index().set_index("sample_name")
3255        samples_with_features = set(full_mf_df.index.get_level_values("sample_name"))
3256
3257        if "scan_time_aligned" in full_mf_df.columns and not overwrite:
3258            raise ValueError("Mass features have already been aligned")
3259
3260        def _set_scan_time_alignment_for_sample(sample_idx, use_alignment, spline):
3261            """Set scan_time_aligned for one sample using spline or identity mapping."""
3262            if use_alignment and spline is not None:
3263                self[sample_idx]._scan_info["scan_time_aligned"] = {
3264                    k: spline(v)
3265                    for k, v in self[sample_idx]._scan_info["scan_time"].items()
3266                }
3267                return True
3268
3269            self[sample_idx]._scan_info["scan_time_aligned"] = (
3270                self[sample_idx]._scan_info["scan_time"].copy()
3271            )
3272            return False
3273
3274        def _get_feature_df_at_or_after(start_idx, index_step, use_alignment, spline):
3275            """Return next sample index/dataframe with features, aligning empty samples on the way."""
3276            i = start_idx
3277            while 0 <= i < len(self):
3278                sample_name = self.samples[i]
3279                if sample_name in samples_with_features:
3280                    mf_df_i = full_mf_df.loc[sample_name].copy()
3281                    mf_df_i["scan_time_og"] = mf_df_i["scan_time"]
3282                    mf_df_i = mf_df_i.reset_index(drop=False)
3283                    if use_alignment and spline is not None:
3284                        # Use previous step transform as a better matching starting point.
3285                        mf_df_i["scan_time"] = spline(mf_df_i["scan_time"])
3286                    return i, mf_df_i
3287
3288                _set_scan_time_alignment_for_sample(i, use_alignment, spline)
3289                self.rt_alignment_attempted = True
3290                i += index_step
3291
3292            return i, None
3293
3294        anchor_mf_dfs = []
3295        for center_obj_id in center_obj_ids:
3296            # Get the anchor mass features from the center LCMS object
3297            mf_df_c = full_mf_df.loc[self.samples[center_obj_id]]
3298            mf_df_c = self.get_anchor_mass_features(mf_df_c)
3299            anchor_mf_dfs.append(mf_df_c)
3300
3301            # Set scan_time_aligned to scan_time for the center LCMS object
3302            center_scan_df = self[center_obj_id].scan_df.copy()
3303            center_scan_df["scan_time_aligned"] = center_scan_df["scan_time"]
3304            self[center_obj_id].scan_df = center_scan_df
3305
3306            # Store alignment data for center object (identity mapping)
3307            center_sample_name = self.samples[center_obj_id]
3308
3309            index_steps = (1, -1)
3310            # Run this twice, once going forward (+1 indexing) and once going backward (-1 indexing)
3311            for index_step in index_steps:
3312                # Initialize spline for propagation to samples without features
3313                spl = None
3314                use_spline_alignment = False
3315
3316                # Loop through the other LCMS objects in this direction.
3317                i, mf_df_i = _get_feature_df_at_or_after(
3318                    center_obj_id + index_step,
3319                    index_step,
3320                    use_spline_alignment,
3321                    spl,
3322                )
3323
3324                while mf_df_i is not None:
3325                    mf_df_i = self.get_anchor_mass_features(mf_df_i)
3326
3327                    # Match the mass features in the LCMS object to the anchor mass features in the center LCMS object.
3328                    matches_c, matches_i = self.match_mfs(mf_df_c, mf_df_i)
3329
3330                    if matches_c is not None:
3331                        use_spline_alignment, spl = self.attempt_alignment(
3332                            matches_c, matches_i
3333                        )
3334
3335                        # Record if we used alignment for this sample
3336                        sample_name = self.samples[i]
3337                        self._manifest_dict[sample_name]["use_rt_alignment"] = (
3338                            use_spline_alignment
3339                        )
3340
3341                        if use_spline_alignment:
3342                            # Set new retention times on scan_df for lc_obj using the spline fitting
3343                            matches_i["scan_time_fit"] = spl(matches_i["scan_time"])
3344
3345                        self.rt_aligned = _set_scan_time_alignment_for_sample(
3346                            i, use_spline_alignment, spl
3347                        )
3348                        self.rt_alignment_attempted = True
3349
3350                        i, mf_df_i = _get_feature_df_at_or_after(
3351                            i + index_step,
3352                            index_step,
3353                            use_spline_alignment,
3354                            spl,
3355                        )
3356                    else:
3357                        # If no matches are found, propagate prior alignment from this index step.
3358                        sample_name = self.samples[i]
3359                        used_previous_alignment = (
3360                            use_spline_alignment and spl is not None
3361                        )
3362                        self._manifest_dict[sample_name]["use_rt_alignment"] = (
3363                            used_previous_alignment
3364                        )
3365
3366                        self.rt_aligned = _set_scan_time_alignment_for_sample(
3367                            i, used_previous_alignment, spl
3368                        )
3369                        self.rt_alignment_attempted = True
3370
3371                        i, mf_df_i = _get_feature_df_at_or_after(
3372                            i + index_step,
3373                            index_step,
3374                            used_previous_alignment,
3375                            spl,
3376                        )
3377
3378        # Now align each batch using the center objects as anchors with the other batches
3379        mf_df_c = anchor_mf_dfs[0]
3380        for i in center_obj_ids[1:]:
3381            mf_df_i = full_mf_df.loc[self.samples[i]].copy()
3382            mf_df_i["scan_time_og"] = mf_df_i["scan_time"]
3383            mf_df_i = self.get_anchor_mass_features(mf_df_i)
3384
3385            matches_c, matches_i = self.match_mfs(mf_df_c, mf_df_i)
3386            if matches_c is not None:
3387                use_spline_alignment, spl = self.attempt_alignment(matches_c, matches_i)
3388
3389                # Record if we used alignment for this sample
3390                sample_name = self.samples[i]
3391                self._manifest_dict[sample_name]["use_rt_alignment"] = (
3392                    use_spline_alignment
3393                )
3394
3395                if use_spline_alignment:
3396                    # Set new retention times on all this object's
3397                    new_times = spl(self[i].scan_df["scan_time"])
3398                    new_scan_info = self[i].scan_df.copy()
3399                    new_scan_info["scan_time_aligned"] = new_times
3400                    self[i].scan_df = new_scan_info
3401
3402                    # Get the batch that this object belongs to
3403                    batch = self.manifest[self.samples[i]]["batch"]
3404
3405                    for j in range(len(self)):
3406                        if self.manifest[self.samples[j]]["batch"] == batch:
3407                            if j != i:
3408                                sample_name_j = self.samples[j]
3409                                self._manifest_dict[sample_name_j][
3410                                    "use_rt_alignment"
3411                                ] = use_spline_alignment
3412                                new_scan_info = self[j].scan_df.copy()
3413                                aligned_times = spl(
3414                                    self[j].scan_df["scan_time_aligned"]
3415                                )
3416                                new_scan_info["scan_time_aligned"] = aligned_times
3417                                self[j].scan_df = new_scan_info
3418
3419        # Set final mass_features_dataframe with the aligned scan_time
3420        center_sample_name = self.samples[center_obj_ids[0]]
3421        self._manifest_dict[center_sample_name]["use_rt_alignment"] = False
3422        new_scan_info = self[center_obj_ids[0]].scan_df.copy()
3423        new_scan_info["scan_time_aligned"] = new_scan_info["scan_time"]
3424
3425    def add_consensus_mass_features(self):
3426        """
3427        Create consensus mass features by clustering aligned features across samples.
3428
3429        This method clusters mass features from all samples in the collection based on
3430        their m/z and aligned retention time proximity. Features that cluster together
3431        across samples are assigned a common cluster ID, creating consensus features
3432        that represent the same compound detected across multiple samples.
3433
3434        The clustering process:
3435        1. Partitions features by m/z to avoid large sparse matrices and enable parallelization
3436        2. Clusters features within each partition using hierarchical clustering
3437        3. Merges partition-boundary clusters that represent the same feature
3438        4. Filters out clusters not present in minimum fraction of samples
3439
3440        Must be run after align_lcms_objects(). Results are stored in the
3441        mass_features_dataframe with a 'cluster' column added.
3442
3443        Parameters
3444        ----------
3445        None
3446            Uses parameters from self.parameters.lcms_collection:
3447            - consensus_mz_tol_ppm: m/z tolerance for clustering (ppm)
3448            - consensus_rt_tol: retention time tolerance for clustering (minutes)
3449            - consensus_partition_size: target partition size for managing memory and parallelization
3450            - consensus_min_sample_fraction: minimum fraction of samples a cluster
3451              must appear in to be retained (0-1)
3452            - cores: number of CPU cores to use for parallel partition processing
3453
3454        Returns
3455        -------
3456        None
3457            Updates self.mass_features_dataframe in place by adding 'cluster' column
3458            and filtering to retain only clusters meeting minimum sample presence.
3459
3460        Raises
3461        ------
3462        ValueError
3463            If mass features have not been aligned (run align_lcms_objects() first).
3464
3465        Notes
3466        -----
3467        - Partitioning prevents memory issues with large sparse distance matrices
3468        - Each partition is processed in parallel (up to cores limit)
3469        - Clusters not meeting consensus_min_sample_fraction are automatically removed
3470        - Access cluster_summary_dataframe property for summary statistics
3471        - Use fill_missing_cluster_features() for gap-filling after clustering
3472
3473        See Also
3474        --------
3475        align_lcms_objects : Aligns retention times before consensus clustering
3476        cluster_summary_dataframe : Property that generates summary statistics for clusters
3477        fill_missing_cluster_features : Gap-fill missing features in clusters
3478        """
3479        # Get the combined mass features from all LCMS objects, keep the original index as a separate column
3480        combined_mfs = self.mass_features_dataframe.copy()
3481        combined_mfs["coll_mf_id"] = combined_mfs.index
3482
3483        # Check if the mass features have been aligned
3484        if "scan_time_aligned" not in combined_mfs.columns:
3485            raise ValueError(
3486                "Mass features have not been aligned, run align_lcms_objects() first"
3487            )
3488
3489        # Partition the mass features by mz so we can parallelize the matching before clustering
3490        from corems.chroma_peak.calc import subset as corems_subset
3491
3492        # get max mz from combined_mfs and calculate tolerance from ppm
3493        mz_tol = self.parameters.lcms_collection.consensus_mz_tol_ppm * 1e-6
3494        n_partition_size = self.parameters.lcms_collection.consensus_partition_size
3495        lazy_partitions = corems_subset.multi_sample_partition(
3496            combined_mfs,
3497            split_on="mz",
3498            size=n_partition_size,
3499            tol=mz_tol,
3500            relative=True,
3501        )
3502
3503        # If any of lazy_partitions._counts is 2xn_partition_size, issue a warning
3504        if np.array(lazy_partitions._counts).max() > 2 * n_partition_size:
3505            warnings.warn(
3506                "Some partitions are larger than 2x the goal partition size. Consider increasing the partition or decreasing the mz_tol."
3507            )
3508
3509        # Cluster the mass features within each partition
3510        if self.parameters.lcms_collection.cores > lazy_partitions.n_partitions:
3511            cores_to_use = lazy_partitions.n_partitions
3512        else:
3513            cores_to_use = self.parameters.lcms_collection.cores
3514        # mfs_with_clusters = lazy_partitions.map(self.cluster_mass_features, processes=cores_to_use)
3515        mfs_with_clusters = lazy_partitions.map(
3516            self.cluster_mass_features_agg_cluster, processes=cores_to_use
3517        )
3518
3519        # Clean up cluster id names after partitioning
3520        new_cluster_ids = (
3521            mfs_with_clusters[["cluster", "partition_idx"]]
3522            .drop_duplicates()
3523            .reset_index(drop=True)
3524        )
3525        new_cluster_ids["cluster_unqiue"] = new_cluster_ids.index
3526        mfs_with_clusters = mfs_with_clusters.merge(
3527            new_cluster_ids, on=["cluster", "partition_idx"]
3528        )
3529        mfs_with_clusters["cluster"] = mfs_with_clusters["cluster_unqiue"]
3530        mfs_with_clusters = mfs_with_clusters.drop(columns=["cluster_unqiue"])
3531
3532        # Embed a new cluster id into the mass features dataframe and set as index
3533        mfs_with_clusters["idx"] = mfs_with_clusters.index
3534
3535        try:
3536            # Check if any clusters can be merged into a single cluster
3537            eval_dict = self.evaluate_clusters_for_repeats(mfs_with_clusters)
3538
3539            # Merge clusters identified in eval_dict
3540            while len(eval_dict["merge_these_clusters"]) > 0:
3541                list_of_clusters_to_merge = [
3542                    [x[0], x[1]] for x in eval_dict["merge_these_clusters"]
3543                ]
3544                # Convert to a dataframe with columns "new_cluster" and "cluster"
3545                df = pd.DataFrame(
3546                    np.array(list_of_clusters_to_merge),
3547                    columns=["new_cluster", "cluster"],
3548                )
3549                # Drop duplicates of "child" clusters
3550                df = df.drop_duplicates("cluster", keep="first")
3551                df = df.drop_duplicates("new_cluster", keep="first")
3552                mfs_with_clusters = mfs_with_clusters.merge(
3553                    df, on="cluster", how="left"
3554                )
3555                mfs_with_clusters["cluster"] = mfs_with_clusters["new_cluster"].fillna(
3556                    mfs_with_clusters["cluster"]
3557                )
3558                mfs_with_clusters = mfs_with_clusters.drop(columns=["new_cluster"])
3559
3560                # Re-evaluate clusters for repeats
3561                eval_dict = self.evaluate_clusters_for_repeats(mfs_with_clusters)
3562                self.mass_features_dataframe = mfs_with_clusters
3563
3564        except:
3565            mfs_with_clusters.set_index("coll_mf_id", inplace=True)
3566            self.mass_features_dataframe = mfs_with_clusters
3567
3568        # Filter out clusters that don't meet minimum sample fraction
3569        self._filter_clusters_by_sample_presence()
3570
3571        # TODO KRH: Deal with isomers better? Pool them together and then split them out using samples with 2 as the template?
3572
3573    def _filter_clusters_by_sample_presence(self):
3574        """
3575        Filter out clusters that don't meet the minimum sample fraction threshold.
3576
3577        Removes clusters (and their associated mass features) from the mass_features_dataframe
3578        if they don't appear in at least consensus_min_sample_fraction of samples.
3579
3580        This is called automatically at the end of add_consensus_mass_features().
3581
3582        Returns
3583        -------
3584        None
3585            Updates self.mass_features_dataframe in place by removing clusters that don't
3586            meet the minimum sample presence threshold.
3587        """
3588        if (
3589            self.mass_features_dataframe is None
3590            or len(self.mass_features_dataframe) == 0
3591        ):
3592            return
3593
3594        min_sample_fraction = (
3595            self.parameters.lcms_collection.consensus_min_sample_fraction
3596        )
3597
3598        # Validate parameter
3599        if not 0 <= min_sample_fraction <= 1:
3600            raise ValueError("consensus_min_sample_fraction must be between 0 and 1")
3601
3602        # Calculate minimum number of samples required
3603        total_samples = len(self.samples)
3604        min_samples_required = min_sample_fraction * total_samples
3605
3606        # Count unique samples per cluster
3607        cluster_sample_counts = (
3608            self.mass_features_dataframe.groupby("cluster")["sample_id"]
3609            .nunique()
3610            .reset_index(name="sample_count")
3611        )
3612
3613        # Identify clusters to keep
3614        clusters_to_keep = cluster_sample_counts[
3615            cluster_sample_counts["sample_count"] > min_samples_required
3616        ]["cluster"].values
3617
3618        # Filter mass features dataframe
3619        self.mass_features_dataframe = self.mass_features_dataframe[
3620            self.mass_features_dataframe["cluster"].isin(clusters_to_keep)
3621        ]
3622
3623    def summarize_clusters(self):
3624        """
3625        Generate summary statistics for consensus mass feature clusters.
3626
3627        Computes aggregate statistics (median, mean, std, min, max) for each cluster
3628        across all samples. Combines both regular mass features and induced mass features
3629        (from gap-filling) when available to provide complete cluster statistics.
3630
3631        Must be run after add_consensus_mass_features() which creates the cluster assignments.
3632        Results are stored in cluster_summary_dataframe property and used by plotting methods.
3633
3634        Parameters
3635        ----------
3636        None
3637            Operates on self.mass_features_dataframe and self.induced_mass_features_dataframe.
3638            Both must contain 'cluster' column.
3639
3640        Returns
3641        -------
3642        :obj:`~pandas.DataFrame` or None
3643            DataFrame with one row per cluster containing summary statistics:
3644            - cluster: cluster ID
3645            - mz_{median,mean,std,max,min}: m/z statistics
3646            - scan_time_aligned_{median,mean,std,max,min}: aligned RT statistics
3647            - half_height_width_{median,mean,std,max,min}: peak width statistics
3648            - tailing_factor_{median,mean,std,max,min}: peak shape statistics
3649            - dispersity_index_{median,mean,std,max,min}: peak quality statistics
3650            - sample_id_nunique: number of unique samples containing the cluster
3651            - intensity_{max,median,mean,std,min}: intensity statistics
3652            - persistence_{max,median,mean,std,min}: persistence statistics
3653
3654            Returns None if mass_features_dataframe is empty.
3655
3656        Notes
3657        -----
3658        - Summary DataFrame is automatically stored in cluster_summary_dataframe property
3659        - Includes both regular and induced (gap-filled) mass features when available.
3660          ``group_consensus_features`` does not use this m/z median for ion-role
3661          geometry; that median is detected features only.
3662        - Used by plotting methods: plot_consensus_mz_features, plot_mz_features_per_cluster
3663        - Sample count (sample_id_nunique) indicates cluster prevalence across samples
3664        - Filters applied by consensus_min_sample_fraction affect which clusters appear
3665
3666        See Also
3667        --------
3668        add_consensus_mass_features : Creates clusters before summarization
3669        fill_missing_cluster_features : Creates induced mass features via gap-filling
3670        plot_consensus_mz_features : Visualizes cluster summaries
3671        plot_mz_features_per_cluster : Shows cluster size distribution
3672        """
3673        # First check if there are minimum columns in the features dataframe
3674        if (
3675            self.mass_features_dataframe is None
3676            or len(self.mass_features_dataframe.columns) < 1
3677        ):
3678            return None
3679
3680        # Combine regular and induced mass features
3681        mf_df = self.mass_features_dataframe.copy()
3682        mf_df = mf_df.reset_index(drop=False)
3683
3684        # Check if induced mass features are available and combine them
3685        if (
3686            self.induced_mass_features_dataframe is not None
3687            and len(self.induced_mass_features_dataframe) > 0
3688        ):
3689            imf_df = self.induced_mass_features_dataframe.copy()
3690            imf_df = imf_df.reset_index(drop=False)
3691            # Cluster column extracted from mf_id in _prepare_lcms_mass_features_for_combination
3692            # Combine regular and induced features
3693            mf_df = pd.concat([mf_df, imf_df], axis=0)
3694            mf_df = mf_df.reset_index(drop=True)
3695
3696        # Filter out any rows with NaN cluster values before converting to int
3697        if "cluster" in mf_df.columns:
3698            mf_df = mf_df.dropna(subset=["cluster"])
3699            mf_df["cluster"] = mf_df["cluster"].astype(int)
3700
3701        # Build aggregation dictionary based on available columns
3702        agg_dict = {
3703            "mz": ["median", "mean", "std", "max", "min"],
3704            "scan_time_aligned": ["median", "mean", "std", "max", "min"],
3705            "sample_id": ["nunique"],
3706            "intensity": ["max", "median", "mean", "std", "min"],
3707        }
3708
3709        # Add optional columns if they exist
3710        optional_columns = {
3711            "half_height_width": ["median", "mean", "std", "max", "min"],
3712            "tailing_factor": ["median", "mean", "std", "max", "min"],
3713            "dispersity_index": ["median", "mean", "std", "max", "min"],
3714            "persistence": ["max", "median", "mean", "std", "min"],
3715        }
3716
3717        for col, funcs in optional_columns.items():
3718            if col in mf_df.columns:
3719                agg_dict[col] = funcs
3720
3721        summary_df = mf_df.groupby("cluster").agg(agg_dict).reset_index()
3722
3723        # Fix the column names
3724        summary_df.columns = [
3725            "_".join(col).strip()
3726            for col in summary_df.columns.values
3727            if col != "cluster"
3728        ]
3729        summary_df = summary_df.rename(columns={"cluster_": "cluster"})
3730        # Set cluster as the index for easy lookup
3731        summary_df = summary_df.set_index("cluster")
3732
3733        # Merge consensus feature-group labels when present (isotopes + adducts).
3734        # Cluster summary exposes possible_ion_types (all
3735        # residual-tied candidates) and drops ion_type so a single type column
3736        # is not mistaken for unique assignment.
3737        fg = getattr(self, "feature_group_dataframe", None)
3738        if fg is not None and len(fg) > 0:
3739            from corems.mass_spectra.calc.feature_grouping import (
3740                GROUP_COLUMNS,
3741                merge_possible_ion_types,
3742            )
3743
3744            cols = [c for c in GROUP_COLUMNS if c in fg.columns]
3745            if cols:
3746                drop_cols = [c for c in cols if c in summary_df.columns]
3747                if drop_cols:
3748                    summary_df = summary_df.drop(columns=drop_cols)
3749                # Also drop stale ion_type from a prior join shape
3750                if "ion_type" in summary_df.columns:
3751                    summary_df = summary_df.drop(columns=["ion_type"])
3752                summary_df = summary_df.join(fg[cols], how="left")
3753
3754            if (
3755                "possible_ion_types" in summary_df.columns
3756                and "ion_type" in summary_df.columns
3757            ):
3758                for cid in summary_df.index:
3759                    it = summary_df.loc[cid, "ion_type"]
3760                    poss = summary_df.loc[cid, "possible_ion_types"]
3761                    filled = merge_possible_ion_types(it, poss)
3762                    if filled is not None:
3763                        summary_df.loc[cid, "possible_ion_types"] = filled
3764            if "ion_type" in summary_df.columns:
3765                summary_df = summary_df.drop(columns=["ion_type"])
3766
3767        return summary_df
3768
3769    def group_consensus_features(self):
3770        """
3771        Group consensus features into natural-abundance isotope and adduct families.
3772
3773        After consensus clustering and gap-filling, labels clusters that share
3774        similar retention time, an interpretable **singly-charged** mass offset,
3775        and correlated cross-sample **apex peak heights**.
3776
3777        Two stages
3778        ----------
3779        1. **Isotopes.** Unit ``Atoms`` edges (e.g. ¹²C/¹³C) plus roll-up to
3780           ¹³C₂, ¹³C₃, … when intermediates exist. Spacing is the raw Atoms
3781           Δm (``|z| = 1``). Multi-charge envelopes are out of scope.
3782           Chemical mono is the geometry root, not the tallest peak.
3783        2. **Adducts.** Pairwise ``ion_type_dict`` forms in
3784           ``feature_group_ion_types`` (must be ``|z| = 1``). Empty or a
3785           single configured type skips adduct linking (isotopes only).
3786
3787        **Not** for tracer, enriched, or isotopically labeled experiments.
3788
3789        Quant gate is fixed (no settings switch):
3790
3791        - Correlation: **Pearson** only, pairwise-complete on samples with both
3792          heights > 0 (threshold ``feature_group_corr_threshold``, default 0.80)
3793        - Abundance: mass-feature apex ``intensity`` only (not integrated area)
3794
3795        Coelution and m/z delta windows reuse collection
3796        ``alignment_rt_tol`` and ``alignment_mz_tol_ppm``. Other knobs
3797        (``feature_group_isotope_atoms``, shared-sample fraction,
3798        ``feature_group_ion_types``, …) live on ``parameters.lcms_collection``.
3799
3800        ``feature_group_ion_types`` is filtered to the collection polarity
3801        (all samples share one polarity; mixed collections cannot be
3802        instantiated). Only keys ending in ``+`` are kept on positive data
3803        and only keys ending in ``-`` on negative data.
3804
3805        Preferred for collection-level families. Optional per-file
3806        ``LCMSBase.find_c13_mass_features`` remains available for single-file
3807        workflows; do not treat both as authoritative in the same pipeline.
3808
3809        Must run after ``add_consensus_mass_features()``. The m/z median used for
3810        isotope and adduct spacing is calculated from detected features only.
3811        Gap-filled m/z values stay on the induced rows and do not enter that
3812        median. The height matrix still includes gap-filled apex intensities
3813        when present; missing entries are 0.
3814
3815        Returns
3816        -------
3817        pandas.DataFrame
3818            Cluster-indexed labels with columns ``feature_group_id``,
3819            ``ion_role``, ``ion_type``, ``possible_ion_types``,
3820            ``isotope_state``, ``mono_cluster_id`` (set only when the form
3821            has at least one isotope). Also stored on
3822            ``self.feature_group_dataframe`` and merged into
3823            ``mass_features_dataframe`` / ``induced_mass_features_dataframe``.
3824
3825        Raises
3826        ------
3827        ValueError
3828            If consensus features are missing or settings are invalid.
3829        """
3830        from corems.mass_spectra.calc.feature_grouping import (
3831            GROUP_COLUMNS,
3832            HEIGHT_COL,
3833            FeatureGroupParams,
3834            build_height_matrix_from_features,
3835            detected_cluster_mz_median,
3836            group_features_arrays,
3837            normalize_ms_polarity,
3838        )
3839
3840        if (
3841            self.mass_features_dataframe is None
3842            or len(self.mass_features_dataframe) == 0
3843        ):
3844            raise ValueError(
3845                "No mass features dataframe. Run add_consensus_mass_features() first."
3846            )
3847        if "cluster" not in self.mass_features_dataframe.columns:
3848            raise ValueError(
3849                "mass_features_dataframe has no 'cluster' column. "
3850                "Run add_consensus_mass_features() first."
3851            )
3852
3853        # Clear prior labels (idempotent re-run)
3854        self.feature_group_dataframe = None
3855        for col in GROUP_COLUMNS:
3856            if col in self.mass_features_dataframe.columns:
3857                self.mass_features_dataframe = self.mass_features_dataframe.drop(
3858                    columns=[col]
3859                )
3860            if (
3861                self.induced_mass_features_dataframe is not None
3862                and col in self.induced_mass_features_dataframe.columns
3863            ):
3864                self.induced_mass_features_dataframe = (
3865                    self.induced_mass_features_dataframe.drop(columns=[col])
3866                )
3867
3868        summary = self.summarize_clusters()
3869        if summary is None or len(summary) == 0:
3870            raise ValueError(
3871                "cluster summary is empty. Run add_consensus_mass_features() first."
3872            )
3873
3874        # Collection polarity: all samples share one polarity (enforced when
3875        # LCMSCollection is loaded). Filter opposite-sign adducts.
3876        collection_polarity = None
3877        if len(self):
3878            collection_polarity = normalize_ms_polarity(
3879                getattr(self[0], "polarity", None)
3880            )
3881
3882        params = FeatureGroupParams.from_lcms_collection_settings(
3883            self.parameters.lcms_collection,
3884            polarity=collection_polarity,
3885        )
3886
3887        cluster_ids = summary.index.to_numpy()
3888        # Gap-filled apex m/z is not part of ion-role geometry.
3889        mz = detected_cluster_mz_median(self.mass_features_dataframe, cluster_ids)
3890        if "scan_time_aligned_median" not in summary.columns:
3891            raise ValueError(
3892                "cluster summary missing scan_time_aligned_median; "
3893                "ensure features are RT-aligned before consensus."
3894            )
3895        rt = summary["scan_time_aligned_median"].to_numpy(dtype=float)
3896
3897        # Sample ids as used in mass_features_dataframe
3898        sample_ids = sorted(self.mass_features_dataframe["sample_id"].dropna().unique())
3899        if self.induced_mass_features_dataframe is not None and len(
3900            self.induced_mass_features_dataframe
3901        ):
3902            sample_ids = sorted(
3903                set(sample_ids)
3904                | set(
3905                    self.induced_mass_features_dataframe["sample_id"]
3906                    .dropna()
3907                    .unique()
3908                    .tolist()
3909                )
3910            )
3911
3912        t_height = time.perf_counter()
3913        heights = build_height_matrix_from_features(
3914            self.mass_features_dataframe,
3915            cluster_ids=cluster_ids,
3916            sample_ids=sample_ids,
3917            induced_df=self.induced_mass_features_dataframe,
3918            intensity_col=HEIGHT_COL,  # apex intensity only; not area
3919        )
3920        stage_timings = {
3921            "height_matrix": time.perf_counter() - t_height,
3922            "n_ion_types_configured": float(len(params.ion_types)),
3923            "collection_polarity": (
3924                str(collection_polarity) if collection_polarity else "unknown"
3925            ),
3926        }
3927
3928        algo_timings: dict = {}
3929        labels = group_features_arrays(
3930            cluster_ids, mz, rt, heights, params, timings_out=algo_timings
3931        )
3932        stage_timings.update(algo_timings)
3933        self._last_feature_group_timings = stage_timings
3934        self.feature_group_dataframe = labels
3935
3936        from corems.mass_spectra.calc.feature_grouping import (
3937            merge_feature_group_labels_into_frames,
3938        )
3939
3940        self.mass_features_dataframe, induced = merge_feature_group_labels_into_frames(
3941            labels,
3942            self.mass_features_dataframe,
3943            self.induced_mass_features_dataframe,
3944        )
3945        if self.induced_mass_features_dataframe is not None:
3946            self.induced_mass_features_dataframe = induced
3947
3948        return labels
3949
3950    def plot_mz_features_per_cluster(self, return_fig=False, path=None):
3951        """
3952        Plot the number of mass features in a cluster against how many clusters
3953        contain that number of mass features
3954
3955        Parameters
3956        -----------
3957        return_fig : bool, optional
3958            If True, return the open figure (caller owns lifecycle).
3959            Default is False.
3960        path : str or path-like, optional
3961            If set, save the figure to this path. When ``return_fig`` is False,
3962            the figure is closed after saving and ``plt.show()`` is not called.
3963
3964        Returns
3965        --------
3966        matplotlib.figure.Figure or None
3967            The figure if ``return_fig`` is True; otherwise None.
3968
3969        Raises
3970        ------
3971        ValueError
3972            If consensus features haven't been added to the object yet
3973        """
3974
3975        if not hasattr(self, "cluster_summary_dataframe"):
3976            raise ValueError(
3977                "cluster_summary_dataframe is not set, must run add_consensus_mass_features() first"
3978            )
3979        sum_data = self.cluster_summary_dataframe
3980        fig, ax = plt.subplots()
3981        sum_data.sample_id_nunique.value_counts().sort_index().plot(ax=ax, kind="bar")
3982        plt.xlabel("Number of mass features in a cluster")
3983        plt.ylabel("Number of clusters with this many mass features")
3984        return _finalize_plot(fig, return_fig=return_fig, path=path)
3985
3986    def plot_mz_features_across_samples(
3987        self, alpha=0.75, s=0.005, return_fig=False, path=None
3988    ):
3989        """
3990        Generate Scan Time vs m/z plot of all the mass features across all
3991        samples in collection where intensity of color on the plot indicates
3992        density of mass features, NOT INTENSITY
3993
3994        Parameters
3995        -----------
3996        alpha :  float
3997            Desired transparency for plotted m/z features.  Defaults to 0.75.
3998        s : float
3999            Desired size of plotted m/z features. Defaults to 0.005.
4000        return_fig : bool, optional
4001            If True, return the open figure (caller owns lifecycle).
4002            Default is False.
4003        path : str or path-like, optional
4004            If set, save the figure to this path. When ``return_fig`` is False,
4005            the figure is closed after saving and ``plt.show()`` is not called.
4006
4007        Returns
4008        --------
4009        matplotlib.figure.Figure or None
4010            A figure displaying a scan time vs m/z scatterplot of all the m/z
4011            features identified in the collection if ``return_fig`` is True;
4012            otherwise None. Parameters alpha (transparency) and s (marker size)
4013            allow the user to emphasize the density of features. Intensity of
4014            features is not represented.
4015        """
4016        df = self.mass_features_dataframe.copy()
4017        fig = plt.figure()
4018        plt.scatter(df.scan_time_aligned, df.mz, c="tab:gray", alpha=alpha, s=s)
4019
4020        plt.xlabel("Scan time")
4021        plt.ylabel("m/z")
4022        plt.ylim(0, np.ceil(np.max(df.mz)))
4023        plt.xlim(0, np.ceil(np.max(df.scan_time)))
4024        plt.title("All mass features, all samples")
4025
4026        return _finalize_plot(fig, return_fig=return_fig, path=path)
4027
4028    def plot_consensus_mz_features(
4029        self,
4030        xb="xb",
4031        xt="xt",
4032        yb="yb",
4033        yt="yt",
4034        show_all=True,
4035        return_fig=False,
4036        path=None,
4037    ):
4038        """
4039        Generate Scan Time vs m/z plot of the consensus features scaled by size
4040        with option ('show_all') of leaving the individual m/z features in the figure.
4041
4042        Parameters
4043        -----------
4044        xb :  float
4045            Desired starting scan time value for the x-axis. Defaults to 0.
4046        xt : float
4047            Desired ending scan time for the x-axis. Defaults to the maximum scan time value in the provided data.
4048        yb :  float
4049            Desired starting m/z value for the y-axis. Defaults to 0.
4050        yt : float
4051            Desired ending m/z for the y-axis. Defaults to the maximum m/z value in the provided data.
4052        show_all : boolean
4053            Indicates whether to display all identified m/z features (True) or just the consensus features (False). Defaults to True.
4054        return_fig : bool, optional
4055            If True, return the open figure (caller owns lifecycle).
4056            Default is False.
4057        path : str or path-like, optional
4058            If set, save the figure to this path. When ``return_fig`` is False,
4059            the figure is closed after saving and ``plt.show()`` is not called.
4060
4061        Returns
4062        --------
4063        matplotlib.figure.Figure or None
4064            A scalable figure that overlays the consensus features over all the
4065            m/z features identified in the collection if ``return_fig`` is True;
4066            otherwise None. Consensus features are scaled by how many m/z
4067            features are represented in the consensus. Figure can be scaled by
4068            inputting desired boundaries on the scan time (xb, xt) and m/z
4069            values (yb, yt).
4070        """
4071        df = self.cluster_summary_dataframe.copy()
4072        mfdf = self.mass_features_dataframe.copy()
4073
4074        fig = plt.figure()
4075        if show_all:
4076            plt.scatter(mfdf.scan_time_aligned, mfdf.mz, c="tab:gray", s=1)
4077
4078        m = plt.scatter(
4079            df.scan_time_aligned_median,
4080            df.mz_median,
4081            c="tab:orange",
4082            alpha=0.7,
4083            s=(df.sample_id_nunique**2) / 5,
4084        )
4085
4086        plt.xlabel("Scan time")
4087        plt.ylabel("m/z")
4088
4089        if xt == "xt":
4090            xt = np.ceil(np.max(mfdf.mz))
4091        if yt == "yt":
4092            yt = np.ceil(np.max(mfdf.scan_time))
4093        if xb == "xb":
4094            xb = 0
4095        if yb == "yb":
4096            yb = 0
4097        plt.ylim(xb, xt)
4098        plt.xlim(yb, yt)
4099
4100        kw = dict(
4101            prop="sizes",
4102            num=max(1, int(len(df.sample_id_nunique.unique()) / 3)),
4103            color="tab:orange",
4104            alpha=0.7,
4105            func=lambda s: np.sqrt(s * 5),
4106        )
4107
4108        plt.legend(
4109            *m.legend_elements(**kw),
4110            title="Features\nper cluster",
4111            bbox_to_anchor=(1.01, 0.4, 0.225, 0.5),
4112        )
4113        plt.tight_layout()
4114        plt.title("Consensus Features")
4115
4116        return _finalize_plot(fig, return_fig=return_fig, path=path)
4117
4118    def plot_cluster(
4119        self,
4120        cluster_id,
4121        to_plot=["EIC", "MS1", "MS2"],
4122        return_fig=False,
4123        path=None,
4124        plot_smoothed_eic=False,
4125        plot_eic_datapoints=False,
4126        eic_buffer_time=None,
4127        label_samples=False,
4128        molecular_metadata=None,
4129        spectral_library=None,
4130    ):
4131        """
4132        Plot a consensus mass feature cluster across all samples.
4133
4134        Similar to LCMSMassFeature.plot() but shows EICs from all samples in the cluster,
4135        highlighting the representative mass feature.
4136
4137        Parameters
4138        ----------
4139        cluster_id : int
4140            The cluster ID to plot
4141        to_plot : list, optional
4142            List of strings specifying what to plot: "EIC", "MS1", "MS2", "MS2_mirror".
4143            Default is ["EIC", "MS1", "MS2"].
4144        return_fig : bool, optional
4145            If True, return the open figure (caller owns lifecycle).
4146            Default is False.
4147        path : str or path-like, optional
4148            If set, save the figure to this path. When ``return_fig`` is False,
4149            the figure is closed after saving and ``plt.show()`` is not called.
4150        plot_smoothed_eic : bool, optional
4151            If True, plots smoothed EICs. Default is False.
4152        plot_eic_datapoints : bool, optional
4153            If True, plots EIC data points. Default is False.
4154        eic_buffer_time : float, optional
4155            Time buffer around the peak for EIC plotting (minutes).
4156            If None, uses parameter setting. Default is None.
4157        label_samples : bool, optional
4158            If True, labels each sample in the legend. Default is False.
4159        molecular_metadata : dict, optional
4160            Dictionary mapping molecular IDs to MetaboliteMetadata objects.
4161            Required for MS2_mirror plots. Default is None.
4162        spectral_library : FlashEntropySearch, optional
4163            FlashEntropy spectral library containing MS2 spectra.
4164            Required for MS2_mirror plots to retrieve library spectra. Default is None.
4165
4166        Returns
4167        -------
4168        matplotlib.figure.Figure or None
4169            The figure object if return_fig=True, otherwise None
4170
4171        Raises
4172        ------
4173        ValueError
4174            If cluster_id is not found or if required data is not loaded
4175        """
4176        import matplotlib.pyplot as plt
4177
4178        # Get cluster summary for median values
4179        if cluster_id not in self.cluster_summary_dataframe.index:
4180            raise ValueError(
4181                f"Cluster {cluster_id} not found in cluster_summary_dataframe. "
4182                f"Run add_consensus_mass_features() first."
4183            )
4184
4185        cluster_summary = self.cluster_summary_dataframe.loc[cluster_id]
4186
4187        # Get representative mass feature info
4188        rep_info = self.get_most_representative_sample_for_cluster(cluster_id)
4189        rep_sample_id = rep_info["sample_id"]
4190        rep_mf_id = rep_info["mf_id"]
4191        rep_sample = self[rep_sample_id]
4192
4193        # Check if representative mass feature is loaded
4194        if rep_mf_id not in rep_sample.mass_features:
4195            raise ValueError(
4196                f"Representative mass feature {rep_mf_id} not loaded in sample {rep_sample.sample_name}. "
4197                f"Run reload_representative_mass_features() or process_consensus_features() first."
4198            )
4199
4200        rep_mf = rep_sample.mass_features[rep_mf_id]
4201
4202        # Get eic buffer time
4203        if eic_buffer_time is None:
4204            eic_buffer_time = self[0].parameters.lc_ms.eic_buffer_time
4205
4206        # Adjust to_plot based on available data
4207        if rep_mf.mass_spectrum is None:
4208            to_plot = [x for x in to_plot if x != "MS1"]
4209        if len(rep_mf.ms2_mass_spectra) == 0:
4210            to_plot = [x for x in to_plot if x not in ["MS2", "MS2_mirror"]]
4211
4212        # Check if EICs are available
4213        cluster_mfs = self.mass_features_dataframe[
4214            self.mass_features_dataframe["cluster"] == cluster_id
4215        ]
4216
4217        has_eics = False
4218        # Check regular features
4219        for _, row in cluster_mfs.iterrows():
4220            sample_id = int(row["sample_id"])
4221            sample = self[sample_id]
4222            if hasattr(sample, "eics") and sample.eics:
4223                if len(sample.eics) > 0:
4224                    has_eics = True
4225                    break
4226
4227        # Also check induced features if available
4228        induced_cluster_mfs = None
4229        if not has_eics and self.induced_mass_features_dataframe is not None:
4230            induced_cluster_mfs = self.induced_mass_features_dataframe[
4231                self.induced_mass_features_dataframe["cluster"] == cluster_id
4232            ]
4233            for _, row in induced_cluster_mfs.iterrows():
4234                sample_id = int(row["sample_id"])
4235                sample = self[sample_id]
4236                if hasattr(sample, "eics") and sample.eics:
4237                    if len(sample.eics) > 0:
4238                        has_eics = True
4239                        break
4240
4241        if not has_eics:
4242            to_plot = [x for x in to_plot if x != "EIC"]
4243            if len(to_plot) == 0:
4244                raise ValueError(
4245                    f"No plottable data available for cluster {cluster_id}. "
4246                    f"Run process_consensus_features(gather_eics=True, add_ms1=True, add_ms2=True) first."
4247                )
4248
4249        # Get induced features if not already retrieved
4250        if (
4251            induced_cluster_mfs is None
4252            and self.induced_mass_features_dataframe is not None
4253        ):
4254            induced_cluster_mfs = self.induced_mass_features_dataframe[
4255                self.induced_mass_features_dataframe["cluster"] == cluster_id
4256            ]
4257
4258        # Check if MS1 is deconvoluted
4259        deconvoluted = rep_mf._ms_deconvoluted_idx is not None
4260
4261        # Create figure
4262        fig, axs = plt.subplots(
4263            len(to_plot), 1, figsize=(10, len(to_plot) * 4), squeeze=False
4264        )
4265
4266        fig.suptitle(
4267            f"Consensus Cluster {cluster_id}: "
4268            f"m/z = {cluster_summary['mz_median']:.4f} "
4269            f"(±{cluster_summary['mz_std']:.4f}); "
4270            f"RT = {cluster_summary['scan_time_aligned_median']:.2f} min "
4271            f"(±{cluster_summary['scan_time_aligned_std']:.2f}); "
4272            f"{int(cluster_summary['sample_id_nunique'])} samples"
4273        )
4274
4275        i = 0
4276
4277        # EIC plot - show all samples using helper method
4278        if "EIC" in to_plot:
4279            self._plot_multiple_eics(
4280                axs[i][0],
4281                cluster_mfs,
4282                induced_cluster_mfs,
4283                rep_sample_id,
4284                rep_mf_id,
4285                cluster_summary["scan_time_aligned_median"],
4286                eic_buffer_time,
4287                plot_smoothed=plot_smoothed_eic,
4288                plot_datapoints=plot_eic_datapoints,
4289                label_samples=label_samples,
4290                lcms_collection=self,
4291            )
4292            i += 1
4293
4294        # MS1 plot - from representative using helper method
4295        if "MS1" in to_plot:
4296            rep_mf._plot_ms1_spectrum(
4297                axs[i][0], deconvoluted=deconvoluted, sample_name=rep_sample.sample_name
4298            )
4299            i += 1
4300
4301        # MS2 plot - from representative using helper method
4302        if "MS2" in to_plot:
4303            rep_mf._plot_ms2_spectrum(axs[i][0], sample_name=rep_sample.sample_name)
4304            i += 1
4305
4306        # MS2 mirror plot - from representative using helper method
4307        if "MS2_mirror" in to_plot:
4308            rep_mf._plot_ms2_mirror(
4309                axs[i][0],
4310                molecular_metadata=molecular_metadata,
4311                spectral_library=spectral_library,
4312            )
4313            i += 1
4314
4315        plt.tight_layout()
4316
4317        return _finalize_plot(fig, return_fig=return_fig, path=path)
4318
4319    def get_representative_mass_features_for_all_clusters(
4320        self, representative_metric=None
4321    ):
4322        """
4323        Get the most representative mass feature for all clusters in bulk.
4324
4325        This is much more efficient than calling get_most_representative_sample_for_cluster
4326        in a loop, as it processes all clusters in a single pass over the dataframe.
4327
4328        Parameters
4329        ----------
4330        representative_metric : str, optional
4331            The metric to use to determine the most representative sample.
4332            If None, uses the value from self.parameters.lcms_collection.consensus_representative_metric.
4333            Options:
4334            - 'intensity': Selects the mass feature with the highest intensity
4335            - 'intensity_prefer_ms2': Selects the highest intensity feature that has MS2 scans,
4336              or the highest intensity overall if none have MS2
4337            Default is None (uses parameter setting).
4338
4339        Returns
4340        -------
4341        :obj:`~pandas.DataFrame`
4342            DataFrame with one row per cluster containing:
4343            - cluster: cluster ID
4344            - sample_id: The sample ID of the most representative sample
4345            - mf_id: The mass feature ID in the sample
4346            - coll_mf_id: The collection-level mass feature ID (index)
4347            - has_ms2: Whether this mass feature has MS2 scan numbers
4348            - intensity: The intensity value of the representative mass feature
4349        """
4350        # Use default from parameters if not specified
4351        if representative_metric is None:
4352            representative_metric = (
4353                self.parameters.lcms_collection.consensus_representative_metric
4354            )
4355
4356        mf_df = self.mass_features_dataframe.copy()
4357        # Reset index to make coll_mf_id a column we can work with
4358        mf_df = mf_df.reset_index(drop=False)
4359
4360        # Handle special metric 'intensity_prefer_ms2'
4361        if representative_metric == "intensity_prefer_ms2":
4362            if "intensity" not in mf_df.columns:
4363                raise ValueError(
4364                    f"'intensity' column not found in mass_features_dataframe. "
4365                    f"Available columns: {mf_df.columns.tolist()}"
4366                )
4367
4368            # Add has_ms2 flag if ms2_scan_numbers column exists
4369            if "ms2_scan_numbers" in mf_df.columns:
4370
4371                def has_ms2_scans(val):
4372                    if val is None:
4373                        return False
4374                    try:
4375                        return len(val) > 0
4376                    except (TypeError, ValueError):
4377                        return False
4378
4379                mf_df["has_ms2"] = mf_df["ms2_scan_numbers"].apply(has_ms2_scans)
4380
4381                # Sort by has_ms2 (descending) then intensity (descending)
4382                # This ensures features with MS2 are preferred when intensities are equal
4383                mf_df = mf_df.sort_values(
4384                    ["has_ms2", "intensity"], ascending=[False, False]
4385                )
4386            else:
4387                mf_df["has_ms2"] = False
4388                mf_df = mf_df.sort_values("intensity", ascending=False)
4389
4390            # Group by cluster and take the first (highest intensity, preferring MS2)
4391            representatives = mf_df.groupby("cluster").first().reset_index()
4392
4393        else:
4394            # Standard metric - check if it exists
4395            if representative_metric not in mf_df.columns:
4396                raise ValueError(
4397                    f"Metric '{representative_metric}' not found. Available columns: {mf_df.columns.tolist()}"
4398                )
4399
4400            # Add has_ms2 flag for consistency
4401            if "ms2_scan_numbers" in mf_df.columns:
4402
4403                def has_ms2_scans(val):
4404                    if val is None:
4405                        return False
4406                    try:
4407                        return len(val) > 0
4408                    except (TypeError, ValueError):
4409                        return False
4410
4411                mf_df["has_ms2"] = mf_df["ms2_scan_numbers"].apply(has_ms2_scans)
4412            else:
4413                mf_df["has_ms2"] = False
4414
4415            # Get the index of max value for each cluster
4416            idx = mf_df.groupby("cluster")[representative_metric].idxmax()
4417            representatives = mf_df.loc[idx].copy()
4418
4419        # Select only the columns we need
4420        result_cols = [
4421            "cluster",
4422            "sample_id",
4423            "mf_id",
4424            "coll_mf_id",
4425            "has_ms2",
4426            "intensity",
4427        ]
4428        representatives = representatives[result_cols]
4429
4430        return representatives
4431
4432    def get_sample_mf_map_for_representatives(
4433        self, representative_metric=None, include_cluster_id=True
4434    ):
4435        """
4436        Build a mapping of sample_id -> list of representative mass feature IDs to load.
4437
4438        This is a DRY helper method used by both process_consensus_features() and
4439        ReadSavedLCMSCollection to determine which mass features should be loaded
4440        for each sample when loading representatives.
4441
4442        Parameters
4443        ----------
4444        representative_metric : str, optional
4445            The metric to use to determine the most representative sample.
4446            If None, uses the value from self.parameters.lcms_collection.consensus_representative_metric.
4447            Default is None.
4448        include_cluster_id : bool, optional
4449            If True, returns tuples of (mf_id, cluster_id). If False, returns just mf_id.
4450            Default is True.
4451
4452        Returns
4453        -------
4454        dict
4455            Dictionary mapping sample_id (int) to list of mass feature identifiers.
4456            If include_cluster_id=True: list of tuples (mf_id, cluster_id)
4457            If include_cluster_id=False: list of mf_id integers
4458
4459        Examples
4460        --------
4461        >>> # Get map with cluster IDs for loading
4462        >>> sample_mf_map = collection.get_sample_mf_map_for_representatives()
4463        >>> # sample_mf_map = {0: [(123, 0), (456, 1)], 1: [(789, 2)], ...}
4464        >>>
4465        >>> # Get map without cluster IDs for pipeline
4466        >>> sample_mf_map = collection.get_sample_mf_map_for_representatives(include_cluster_id=False)
4467        >>> # sample_mf_map = {0: [123, 456], 1: [789], ...}
4468        """
4469        # Get all representative mass features in bulk (much faster than looping)
4470        representatives = self.get_representative_mass_features_for_all_clusters(
4471            representative_metric=representative_metric
4472        )
4473
4474        # Build sample_mf_map
4475        sample_mf_map = {}
4476        for _, row in representatives.iterrows():
4477            sample_id = row["sample_id"]
4478            mf_id = row["mf_id"]
4479            cluster_id = row["cluster"]
4480
4481            if sample_id not in sample_mf_map:
4482                sample_mf_map[sample_id] = []
4483
4484            if include_cluster_id:
4485                sample_mf_map[sample_id].append((mf_id, cluster_id))
4486            else:
4487                sample_mf_map[sample_id].append(mf_id)
4488
4489        return sample_mf_map
4490
4491    def get_most_representative_sample_for_cluster(
4492        self, cluster_id, representative_metric=None
4493    ):
4494        """
4495        Get the most representative sample for a given cluster based on a metric.
4496
4497        Parameters
4498        ----------
4499        cluster_id : int
4500            The cluster ID to find the representative sample for.
4501        representative_metric : str, optional
4502            The metric to use to determine the most representative sample.
4503            If None, uses the value from self.parameters.lcms_collection.consensus_representative_metric.
4504            Options:
4505            - 'intensity': Selects the mass feature with the highest intensity
4506            - 'intensity_prefer_ms2': Selects the highest intensity feature that has MS2 scans,
4507              or the highest intensity overall if none have MS2
4508            Default is None (uses parameter setting).
4509
4510        Returns
4511        -------
4512        dict
4513            Dictionary containing:
4514            - 'sample_id': The sample ID of the most representative sample
4515            - 'sample_name': The sample name of the most representative sample
4516            - 'mf_id': The mass feature ID in the sample
4517            - 'coll_mf_id': The collection-level mass feature ID (index)
4518            - 'has_ms2': Whether this mass feature has MS2 scan numbers
4519            - 'intensity': The intensity value of the representative mass feature
4520
4521        Raises
4522        ------
4523        ValueError
4524            If cluster_id is not found or if representative_metric is not a valid column.
4525        """
4526        # Use the bulk method to get all representatives, then filter to this cluster
4527        # This follows DRY principle and ensures consistency
4528        all_representatives = self.get_representative_mass_features_for_all_clusters(
4529            representative_metric=representative_metric
4530        )
4531
4532        # Filter to the requested cluster
4533        cluster_rep = all_representatives[all_representatives["cluster"] == cluster_id]
4534
4535        if len(cluster_rep) == 0:
4536            # Try to provide helpful error message
4537            available_clusters = self.mass_features_dataframe["cluster"].unique()
4538            raise ValueError(
4539                f"Cluster {cluster_id} not found in mass_features_dataframe. "
4540                f"Available clusters: {sorted(available_clusters[:10].tolist())}... "
4541                f"(showing first 10 of {len(available_clusters)} total clusters)"
4542            )
4543
4544        # Get the representative row (should only be one)
4545        rep_row = cluster_rep.iloc[0]
4546
4547        # Get sample name from sample_id (convert to int for list indexing)
4548        sample_id = int(rep_row["sample_id"])
4549        sample_name = self.samples[sample_id]
4550
4551        return {
4552            "sample_id": sample_id,
4553            "sample_name": sample_name,
4554            "mf_id": rep_row["mf_id"],
4555            "coll_mf_id": rep_row["coll_mf_id"],
4556            "has_ms2": rep_row["has_ms2"],
4557            "intensity": rep_row["intensity"],
4558        }
4559
4560    def reload_representative_mass_features(
4561        self,
4562        add_ms2=False,
4563        auto_process_ms2=True,
4564        ms2_spectrum_mode=None,
4565        ms2_scan_filter=None,
4566    ):
4567        """
4568        Reload mass features for all representative samples in the cluster summary.
4569
4570        This method is useful when the collection was loaded with load_light=True,
4571        which stores mass features only in the collection dataframe. This reloads
4572        the specific mass features that are representatives for each cluster,
4573        allowing them to be accessed as LCMSMassFeature objects.
4574
4575        Parameters
4576        ----------
4577        add_ms2 : bool, optional
4578            If True, also loads and associates MS2 spectra with mass features. Default is False.
4579        auto_process_ms2 : bool, optional
4580            If True and add_ms2=True, auto-processes MS2 spectra. Default is True.
4581        ms2_spectrum_mode : str or None, optional
4582            Spectrum mode for MS2 spectra. If None, determines from parser. Default is None.
4583        ms2_scan_filter : str or None, optional
4584            Filter string for MS2 scans (e.g., 'hcd'). Default is None.
4585
4586        Returns
4587        -------
4588        dict
4589            Dictionary mapping sample_id to list of reloaded mf_ids.
4590
4591        Raises
4592        ------
4593        ValueError
4594            If cluster_summary_dataframe is not set (run add_consensus_mass_features first).
4595
4596        Notes
4597        -----
4598        - Only reloads mass features that are cluster representatives
4599        - Uses get_most_representative_sample_for_cluster() to determine which to reload
4600        - More memory-efficient than reloading all mass features
4601        - Parallelized based on lcms_collection.cores parameter
4602        - MS2 association uses same logic as add_associated_ms2_dda()
4603
4604        See Also
4605        --------
4606        _reload_sample_mass_features : Low-level method to reload specific mass features
4607        get_most_representative_sample_for_cluster : Gets representative sample for cluster
4608        """
4609        # Validate prerequisites
4610        if (
4611            not hasattr(self, "cluster_summary_dataframe")
4612            or self.cluster_summary_dataframe is None
4613        ):
4614            raise ValueError(
4615                "cluster_summary_dataframe not found. Must run add_consensus_mass_features() first."
4616            )
4617
4618        # Get all representative mass features in bulk (much faster than looping)
4619        representatives = self.get_representative_mass_features_for_all_clusters()
4620
4621        # Build a dictionary of sample_id -> list of mf_ids that are representatives
4622        sample_mf_map = {}
4623        for _, row in representatives.iterrows():
4624            sample_id = row["sample_id"]
4625            mf_id = row["mf_id"]
4626
4627            if sample_id not in sample_mf_map:
4628                sample_mf_map[sample_id] = []
4629            sample_mf_map[sample_id].append(mf_id)
4630
4631        # Reload mass features for each sample (parallelized)
4632        if self.parameters.lcms_collection.cores == 1:
4633            # Serial processing
4634            from tqdm import tqdm
4635
4636            for sample_id in tqdm(
4637                sample_mf_map.keys(),
4638                desc="Reloading representative mass features",
4639                unit="sample",
4640            ):
4641                mf_ids = sample_mf_map[sample_id]
4642                self._reload_sample_mass_features(
4643                    sample_id,
4644                    mf_ids_to_load=mf_ids,
4645                    add_ms2=add_ms2,
4646                    auto_process_ms2=auto_process_ms2,
4647                    ms2_spectrum_mode=ms2_spectrum_mode,
4648                    ms2_scan_filter=ms2_scan_filter,
4649                )
4650        else:
4651            # Parallel processing
4652            from tqdm import tqdm
4653
4654            if self.parameters.lcms_collection.cores > len(sample_mf_map):
4655                ncores = len(sample_mf_map)
4656            else:
4657                ncores = self.parameters.lcms_collection.cores
4658
4659            # Build arguments list for starmap
4660            args_list = [
4661                (
4662                    sample_id,
4663                    sample_mf_map[sample_id],
4664                    add_ms2,
4665                    auto_process_ms2,
4666                    ms2_spectrum_mode,
4667                    ms2_scan_filter,
4668                    False,
4669                )
4670                for sample_id in sample_mf_map.keys()
4671            ]
4672
4673            # Execute in parallel
4674            with spawn_pool(ncores) as pool:
4675                mp_result = pool.starmap(self._reload_sample_mass_features, args_list)
4676
4677            # Collect results back into samples
4678            for i, sample_id in enumerate(
4679                tqdm(
4680                    sample_mf_map.keys(),
4681                    desc="Collecting reloaded mass features",
4682                    unit="sample",
4683                )
4684            ):
4685                self[sample_id].mass_features = mp_result[i]
4686
4687        return sample_mf_map
4688
4689    def _associate_ms2_with_mass_features(
4690        self,
4691        sample,
4692        local_mf_ids,
4693        auto_process=True,
4694        spectrum_mode=None,
4695        scan_filter=None,
4696    ):
4697        """
4698        Associate MS2 spectra with specific mass features in a sample.
4699
4700        Uses the LCMSBase helper method to find and load MS2 scans for the specified mass features.
4701
4702        Parameters
4703        ----------
4704        sample : LCMSBase
4705            The sample object containing mass features and scan data.
4706        local_mf_ids : list of int
4707            List of local (sample-level) mass feature IDs to find MS2 for.
4708        auto_process : bool, optional
4709            If True, auto-processes the MS2 spectra. Default is True.
4710        spectrum_mode : str or None, optional
4711            Spectrum mode for MS2 spectra. If None, determines from parser. Default is None.
4712        scan_filter : str or None, optional
4713            Filter string for MS2 scans (e.g., 'hcd'). Default is None.
4714
4715        Returns
4716        -------
4717        dict
4718            Dictionary of scan_number -> MassSpectrum objects for the loaded MS2 spectra.
4719        """
4720        # Check if we have scan data
4721        if not hasattr(sample, "scan_df") or sample.scan_df is None:
4722            return {}
4723
4724        # Separate mass features into those that need scan finding vs those that already have scans
4725        mfs_needing_scan_finding = []
4726        unique_dda_scans = set()
4727
4728        for mf_id in local_mf_ids:
4729            if mf_id not in sample.mass_features:
4730                continue
4731            mf = sample.mass_features[mf_id]
4732            # If this mass feature already has MS2 scans, add them to our set
4733            if mf.ms2_scan_numbers is not None and len(mf.ms2_scan_numbers) > 0:
4734                # Convert to integers in case they come from HDF5 as numpy types
4735                unique_dda_scans.update([int(scan) for scan in mf.ms2_scan_numbers])
4736            else:
4737                # Otherwise, we need to find scans for this mass feature
4738                mfs_needing_scan_finding.append(mf_id)
4739
4740        # Only run the scan finding for mass features that need it
4741        if mfs_needing_scan_finding:
4742            found_scans = sample._find_ms2_scans_for_mass_features(
4743                mf_ids=mfs_needing_scan_finding, scan_filter=scan_filter
4744            )
4745            unique_dda_scans.update(found_scans)
4746
4747        if len(unique_dda_scans) == 0:
4748            return {}
4749
4750        # Get ms2 parameters from sample
4751        # TODO KRH: deal with different ms2 scan types here (CID vs HCD), may need to add scan translator to the initializeion
4752        ms_params = sample.parameters.mass_spectrum["ms2"]
4753
4754        # Load MS2 spectra (convert set to list)
4755        sample.add_mass_spectra(
4756            scan_list=list(unique_dda_scans),
4757            auto_process=auto_process,
4758            spectrum_mode=spectrum_mode,
4759            ms_level=2,
4760            use_parser=True,
4761            ms_params=ms_params,
4762        )
4763
4764        # Associate MS2 spectra with mass features
4765        for mf_id in local_mf_ids:
4766            if mf_id not in sample.mass_features:
4767                continue
4768            if (
4769                sample.mass_features[mf_id].ms2_scan_numbers is not None
4770                and len(sample.mass_features[mf_id].ms2_scan_numbers) > 0
4771            ):
4772                for dda_scan in sample.mass_features[mf_id].ms2_scan_numbers:
4773                    if dda_scan in sample._ms:
4774                        sample.mass_features[mf_id].ms2_mass_spectra[dda_scan] = (
4775                            sample._ms[dda_scan]
4776                        )
4777
4778        # Return only the MS2 spectra we loaded (for parallel processing)
4779        return {
4780            scan: sample._ms[scan] for scan in unique_dda_scans if scan in sample._ms
4781        }
4782
4783    def _reload_sample_mass_features(
4784        self,
4785        sample_id,
4786        mf_ids_to_load=None,
4787        add_ms2=False,
4788        auto_process_ms2=True,
4789        ms2_spectrum_mode=None,
4790        ms2_scan_filter=None,
4791        inplace=True,
4792    ):
4793        """
4794        Reload specific mass features for a sample from HDF5.
4795
4796        This is useful when the collection was loaded with load_light=True,
4797        which stores mass features only in the collection dataframe and not
4798        as LCMSMassFeature objects in individual samples.
4799
4800        Parameters
4801        ----------
4802        sample_id : int
4803            The sample ID to reload mass features for.
4804        mf_ids_to_load : list of str, optional
4805            List of collection-level mf_ids (format: '{sample_id}_{local_mf_id}') to load.
4806            If None, loads all mass features for the sample.
4807        add_ms2 : bool, optional
4808            If True, also loads and associates MS2 spectra. Default is False.
4809        auto_process_ms2 : bool, optional
4810            If True, auto-processes MS2 spectra. Default is True.
4811        ms2_spectrum_mode : str or None, optional
4812            Spectrum mode for MS2 spectra. Default is None.
4813        ms2_scan_filter : str or None, optional
4814            Filter string for MS2 scans. Default is None.
4815        inplace : bool, optional
4816            If True, updates the sample's mass_features in place. If False, returns the
4817            mass_features dictionary (for multiprocessing). Default is True.
4818
4819        Returns
4820        -------
4821        dict or None
4822            If inplace=False, returns dictionary of mass features.
4823            Otherwise returns None and updates object in place.
4824        """
4825        sample = self[sample_id]
4826        sample_name = self.samples[sample_id]
4827
4828        # Check if we have a collection parser that can reload
4829        if not hasattr(self, "collection_parser") or self.collection_parser is None:
4830            print(
4831                "Warning: Cannot reload mass features - no collection_parser available"
4832            )
4833            if not inplace:
4834                return {}
4835            return
4836
4837        # Get the HDF5 file for this sample
4838        hdf5_file = (
4839            self.collection_parser.folder_location
4840            / f"{sample_name}.corems/{sample_name}.hdf5"
4841        )
4842
4843        if not hdf5_file.exists():
4844            print(f"Warning: HDF5 file not found for sample {sample_name}: {hdf5_file}")
4845            if not inplace:
4846                return {}
4847            return
4848
4849        # Import here to avoid circular imports
4850        from corems.mass_spectra.input.corems_hdf5 import ReadCoreMSHDFMassSpectra
4851
4852        # If specific mf_ids requested, use them directly
4853        local_mf_ids_to_load = None
4854        if mf_ids_to_load is not None:
4855            # mf_ids_to_load is already a list of sample-level mf_ids (integers)
4856            # No parsing needed - they come from the mf_id column in the dataframe
4857            local_mf_ids_to_load = set(mf_ids_to_load)
4858
4859        # Reload mass features from HDF5
4860        with ReadCoreMSHDFMassSpectra(hdf5_file) as parser:
4861            # Load mass features - if specific IDs requested, only load those
4862            parser.import_mass_features(sample, mf_ids=local_mf_ids_to_load)
4863
4864        # If add_ms2, associate MS2 spectra with the loaded mass features
4865        if add_ms2 and local_mf_ids_to_load is not None:
4866            self._associate_ms2_with_mass_features(
4867                sample,
4868                list(local_mf_ids_to_load),
4869                auto_process=auto_process_ms2,
4870                spectrum_mode=ms2_spectrum_mode,
4871                scan_filter=ms2_scan_filter,
4872            )
4873
4874        # Return mass features if not inplace (for multiprocessing)
4875        if not inplace:
4876            return sample.mass_features
4877
4878    def add_sparse_distance_matrix(self, features):
4879        if features is None:
4880            return None
4881        else:
4882            features = features.copy()
4883
4884        # Parameters for calculating distance between features
4885        dims = ["mz", "scan_time_aligned"]
4886        relative = [True, False]
4887        mz_tol_relative = self.parameters.lcms_collection.consensus_mz_tol_ppm * 1e-6
4888        tol = [mz_tol_relative, self.parameters.lcms_collection.consensus_rt_tol]
4889        dist_weight = [1, 1]
4890
4891        # Check that the dimensions and tolerances are the same length
4892        if (
4893            len(dims) != len(tol)
4894            or len(dims) != len(relative)
4895            or len(dims) != len(dist_weight)
4896        ):
4897            raise ValueError(
4898                "The dimensions, tolerances, relative, dist_weight, and na_allow lists must be the same length"
4899            )
4900
4901        # Make connectivity matrix for masking within sample mass features
4902        ## Masking matrix cmat will mark all features from the same sample as 0
4903        ## To mask, a matrix can be multiplied by cmat and features from same
4904        ## samples are multiplied by 0, while features from different samples
4905        ## are multiplied by 1
4906
4907        if "sample_id" not in features.columns:
4908            cmat = None
4909        else:
4910            vals = features["sample_id"].values.reshape(-1, 1)
4911            cmat = scipy.spatial.distance.cdist(vals, vals)
4912            # Convert to binary (0 if same sample, 1 if different)
4913            cmat = np.where(cmat == 0, 0, 1)
4914            # Convert to coordinate matrix for sparse operations later
4915            cmat = sparse.coo_matrix(cmat)
4916
4917        # Compute inter-feature distances using sparse matrix approach
4918        distances = None  # clear the distances object before starting
4919        for i in range(len(dims)):  # iterate through all dimensions to be considered
4920            # Construct k-d tree
4921            values = features[dims[i]].values
4922
4923            tree = KDTree(values.reshape(-1, 1))
4924
4925            max_tol = tol[i]
4926            if relative[i] is True:
4927                # Maximum absolute tolerance
4928                max_tol = tol[i] * values.max()
4929
4930            # Compute sparse distance matrix
4931            # the larger the max_tol, the slower this operation is
4932            sdm = tree.sparse_distance_matrix(tree, max_tol, output_type="coo_matrix")
4933
4934            # Only consider forward case, exclude diagonal
4935            sdm = sparse.triu(sdm, k=1)
4936
4937            # Filter relative distances
4938            if relative[i] is True:
4939                # Compute relative distances
4940                rel_dists = sdm.data / values[sdm.row]
4941
4942                # Indices of relative distances less than tolerance
4943                idx = rel_dists <= tol[i]
4944
4945                # Reconstruct sparse distance matrix
4946                sdm = sparse.coo_matrix(
4947                    (rel_dists[idx], (sdm.row[idx], sdm.col[idx])),
4948                    shape=(len(values), len(values)),
4949                )
4950
4951            # Scaled distances wrt the maximum tolerance for the dimension
4952            sdm.data = sdm.data / tol[i]
4953
4954            # Stack distances for dimensions where na_allow is False
4955            if distances is None:
4956                sdm.data = sdm.data * dist_weight[i]
4957                # Replace zeros with epsilon to handle perfect matches
4958                sdm.data[sdm.data == 0] = 1e-10
4959                distances = sdm
4960            else:
4961                # Prepare sdm to match shape of existing distances
4962                distances_truth = distances.copy()
4963                # make new sparse matrix with same positions as previous
4964                # distance matrix but all ones for values
4965                distances_truth.data = np.ones_like(distances_truth.data)
4966
4967                # Replace zeros with epsilon BEFORE multiply to prevent sparse matrix from dropping them
4968                sdm.data[sdm.data == 0] = 1e-10
4969
4970                # multiply the new sparse matrix (sdm) by this mask to remove
4971                # data that doesn't exist in original sparse matrix
4972                sdm = distances_truth.multiply(sdm)
4973
4974                sdm.data = sdm.data * dist_weight[i]
4975                # Replace zeros with epsilon to handle perfect matches
4976                sdm.data[sdm.data == 0] = 1e-10
4977
4978                # use same process as before to remove data from previous
4979                # distances matrix that isn't in new distances matrix
4980                sdm_truth = sdm.copy()
4981                sdm_truth.data = np.ones_like(sdm_truth.data)
4982
4983                # remove the distances that are not sdm
4984                distances = distances.multiply(sdm_truth)
4985
4986                # Sum the new distances
4987                distances = distances + sdm
4988
4989        # Multiply by connectivity matrix for more masking
4990        distances = distances.multiply(cmat)
4991
4992        # Set attribute holding distance matrix
4993        self._sparse_distance_matrix = distances
4994
4995    def evaluate_clusters_for_repeats(self, features):
4996        raise NotImplementedError("evaluate_clusters_for_repeats not implemented yet")
4997        summary_df = self.cluster_summary_dataframe.copy()
4998
4999        # Arrange by decreasing median intensity
5000        summary_df = summary_df.sort_values(
5001            by="intensity_median", ascending=False
5002        ).reset_index(drop=True)
5003
5004        # Find clusters that are within the mz_tol and rt_tol of each other (on the medians)
5005        # Create a distance matrix
5006        # Define how to calculate the distance between features
5007        dims = ["mz_median", "scan_time_aligned_median"]
5008        relative = [True, False]
5009        mz_tol_relative = self.parameters.lcms_collection.consensus_mz_tol_ppm * 1e-6
5010        tol = [mz_tol_relative, self.parameters.lcms_collection.consensus_rt_tol]
5011
5012        # Compute inter-feature distances
5013        distances = None
5014        for i in range(len(dims)):
5015            # Construct k-d tree
5016            values = summary_df[dims[i]].values
5017            tree = KDTree(values.reshape(-1, 1))
5018
5019            max_tol = tol[i]
5020            if relative[i] is True:
5021                # Maximum absolute tolerance
5022                max_tol = tol[i] * values.max()
5023
5024            # Compute sparse distance matrix
5025            # the larger the max_tol, the slower this operation is
5026            sdm = tree.sparse_distance_matrix(tree, max_tol, output_type="coo_matrix")
5027
5028            # Only consider forward case, exclude diagonal
5029            sdm = sparse.triu(sdm, k=1)
5030
5031            # Filter relative distances
5032            if relative[i] is True:
5033                # Compute relative distances
5034                rel_dists = sdm.data / values[sdm.row]  # or col?
5035
5036                # Indices of relative distances less than tolerance
5037                idx = rel_dists <= tol[i]
5038
5039                # Reconstruct sparse distance matrix
5040                sdm = sparse.coo_matrix(
5041                    (rel_dists[idx], (sdm.row[idx], sdm.col[idx])),
5042                    shape=(len(values), len(values)),
5043                )
5044
5045            # Cast as binary matrix
5046            sdm.data = np.ones_like(sdm.data)
5047
5048            # Stack distances
5049            if distances is None:
5050                distances = sdm
5051            else:
5052                distances = distances.multiply(sdm)
5053
5054        # Roll up features
5055        # Extract indices of within-tolerance points
5056        distances = distances.tocoo()
5057        pairs = np.stack(
5058            (distances.row, distances.col), axis=1
5059        )  # These are the index values of the clusters, not the cluster ids
5060        # Conver to cluster ids
5061        pairs_df = pd.DataFrame(pairs, columns=["parent", "child"])
5062        pairs_df["parent"] = summary_df.loc[pairs[:, 0]]["cluster"].values
5063        pairs_df["child"] = summary_df.loc[pairs[:, 1]]["cluster"].values
5064        pairs_df = pairs_df.set_index("parent")
5065
5066        merge_these_clusters = []
5067        possible_overlaps = []
5068        root_parents = np.setdiff1d(
5069            np.unique(pairs_df.index.values), np.unique(pairs_df.child.values)
5070        )
5071        for parent in root_parents:
5072            parent_features = features[features["cluster"] == parent]
5073            children = pairs_df.loc[[parent], "child"].tolist()
5074            for child in children:
5075                overlap = self.check_merge(parent_features, child, features)
5076                if len(overlap) == 0:
5077                    merge_these_clusters.append((parent, child, len(overlap)))
5078                else:
5079                    possible_overlaps.append((parent, child, len(overlap)))
5080
5081        result_dict = {}
5082        result_dict["merge_these_clusters"] = merge_these_clusters
5083        result_dict["possible_overlaps"] = possible_overlaps
5084
5085        return result_dict
5086
5087    def check_merge(self, parent_features, child, features):
5088        # Grab the features of the parent and children
5089        child_features = features[features["cluster"] == child]
5090
5091        # Check if there is an overlap between mf_coll_id in the parent and child clusters
5092        overlap = np.intersect1d(
5093            parent_features["sample_id"].values, child_features["sample_id"].values
5094        )
5095
5096        return overlap
5097
5098    def cluster_mass_features_agg_cluster(self, features):
5099        if features is None:
5100            return None
5101
5102        features = features.copy()
5103
5104        self.add_sparse_distance_matrix(features)
5105
5106        distances = self._sparse_distance_matrix
5107
5108        # Convert to full matrix
5109        distances = distances.todense()
5110
5111        # Cast all 0s to 1s for a distance matrix
5112        distances[distances == 0] = 1
5113        distances = np.asarray(distances)
5114
5115        # Perform clustering
5116        try:
5117            clustering = AgglomerativeClustering(
5118                n_clusters=None,
5119                linkage="complete",
5120                # using complete linkage will prevent one sample from being assigned to multiple clusters
5121                metric="precomputed",
5122                distance_threshold=1,
5123            ).fit(distances)
5124            features["cluster"] = clustering.labels_
5125
5126        # All data points are singleton clusters
5127        except:
5128            features["cluster"] = np.arange(len(features.index))
5129
5130        return features
5131
5132    def cluster_inspection_plot(self, clu, return_fig=False, path=None):
5133        """
5134        Generate Scan Time vs m/z plot for a narrow range around a given
5135        cluster. This tool is meant to support the user in fine tuning the
5136        tolerances used for the clustering algorithm. The user-provided cluster
5137        ID is highlighted in larger, magenta marker and the ten largest of the
5138        remaining clusters are idenfitied with different colors while the
5139        smallest clusters are light gray.
5140
5141        Parameters
5142        -----------
5143        clu :  integer
5144            A cluster ID that exists in self.mass_features_dataframe
5145        return_fig : bool, optional
5146            If True, return the open figure (caller owns lifecycle).
5147            Default is False.
5148        path : str or path-like, optional
5149            If set, save the figure to this path. When ``return_fig`` is False,
5150            the figure is closed after saving and ``plt.show()`` is not called.
5151
5152        Returns
5153        --------
5154        matplotlib.figure.Figure or None
5155            A figure displaying a scan time vs m/z scatterplot of small region
5156            around a given cluster with the ten largest clusters in the region
5157            distinctly identified if ``return_fig`` is True; otherwise None.
5158
5159        Raises
5160        ------
5161        ValueError
5162            If cluster data haven't been added to the object yet
5163        """
5164
5165        if "cluster" not in self.mass_features_dataframe.columns:
5166            raise ValueError(
5167                "Cluster information is not yet added to mass_features_dataframe, must run add_consensus_mass_features() first"
5168            )
5169
5170        else:
5171            mztol = self.parameters.lcms_collection.consensus_mz_tol_ppm * 1e-6
5172            rttol = self.parameters.lcms_collection.consensus_rt_tol
5173            clu_features = self.mass_features_dataframe.copy()
5174
5175            inclu = clu_features[clu_features.cluster == clu]
5176            exclu = clu_features[clu_features.cluster != clu]
5177
5178            dt_ymin = np.floor(min(inclu.mz)) - 1
5179            dt_ymax = np.ceil(max(inclu.mz)) + 1
5180            dt_xmin = np.floor(min(inclu.scan_time_aligned)) - 1
5181            dt_xmax = np.ceil(max(inclu.scan_time_aligned)) + 1
5182
5183            exclu = exclu[
5184                (exclu.mz.between(dt_ymin, dt_ymax, inclusive="both"))
5185                & (exclu.scan_time_aligned.between(dt_xmin, dt_xmax, inclusive="both"))
5186            ]
5187
5188            bigclulist = list(exclu.cluster.value_counts()[:10].index)
5189            bigclu = exclu[exclu.cluster.isin(bigclulist)]
5190            smclu = exclu[~exclu.cluster.isin(bigclulist)]
5191
5192            colors = np.arange(0, 10)
5193            colordict = dict(zip(bigclulist, colors))
5194            bigclu["color"] = bigclu.cluster.apply(lambda x: colordict[x])
5195
5196            fig = plt.figure(figsize=(7.5, 5))
5197
5198            plt.scatter(
5199                inclu.scan_time_aligned,
5200                inclu.mz,
5201                c="m",
5202                s=3,
5203                label="Cluster " + str(clu),
5204            )
5205
5206            plt.scatter(
5207                bigclu.scan_time_aligned, bigclu.mz, c=bigclu.color, cmap="tab10", s=1.5
5208            )
5209
5210            plt.scatter(
5211                smclu.scan_time_aligned,
5212                smclu.mz,
5213                c="silver",
5214                s=2,
5215                label="Small clusters",
5216            )
5217
5218            plt.ylim(dt_ymin, dt_ymax)
5219            plt.xlim(dt_xmin, dt_xmax)
5220            plt.legend(ncol=2, bbox_to_anchor=(0.8, -0.1))
5221            plt.xlabel("Scan time")
5222            plt.ylabel("m/z")
5223            title_str = "Cluster " + str(clu)
5224            title_str += ": representing " + str(len(inclu.sample_id.unique()))
5225            title_str += " of " + str(len(clu_features.sample_id.unique()))
5226            title_str += " samples\n"
5227            title_str += "M/Z tolerance: " + str(mztol) + "\n"
5228            title_str += "Scan Time tolerance: " + str(rttol)
5229            plt.title(title_str, fontsize=10)
5230
5231            return _finalize_plot(fig, return_fig=return_fig, path=path)
5232
5233    def plot_cluster_outlier_frequency(
5234        self,
5235        dim_list=["mz", "scan_time_aligned"],
5236        clu_size_thresh=0.5,
5237        return_fig=False,
5238        path=None,
5239    ):
5240        """
5241        Generate histogram showing the frequency of outlier occurrences by
5242        clustering dimension across all clusters
5243
5244        Parameters
5245        -----------
5246        dim_list :  list
5247            List of strings describing dimensions that can be used in
5248            clustering. Available list items:
5249                - 'mz'
5250                - 'scan_time_aligned'
5251                - 'half_height_width'
5252                - 'tailing_factor'
5253                - 'dispersity_index'
5254                - 'intensity'
5255                - 'persistence'
5256        clu_size_thresh : float
5257            Value between 0 and 1 that indicates what percentage of samples
5258            need to be present in a cluster before it's evaluated for outliers.
5259            Defaults to 0.5.
5260        return_fig : bool, optional
5261            If True, return the open figure (caller owns lifecycle).
5262            Default is False.
5263        path : str or path-like, optional
5264            If set, save the figure to this path. When ``return_fig`` is False,
5265            the figure is closed after saving and ``plt.show()`` is not called.
5266
5267        Returns
5268        --------
5269        matplotlib.figure.Figure or None
5270            A figure displaying the frequency of outlier occurrences across all
5271            clusters in the provided measurement dimensions if ``return_fig`` is
5272            True; otherwise None.
5273
5274        Raises
5275        ------
5276        ValueError
5277            If cluster data haven't been added to the object yet
5278        """
5279
5280        if not hasattr(self, "cluster_summary_dataframe"):
5281            raise ValueError(
5282                "cluster_summary_dataframe is not yet added, must run add_consensus_mass_features() first"
5283            )
5284
5285        mfdf = self.mass_features_dataframe.copy()
5286        summarydf = self.cluster_summary_dataframe
5287
5288        numsamples = len(self)
5289        sumdf = (
5290            summarydf[summarydf.sample_id_nunique > numsamples * clu_size_thresh]
5291            .reset_index(drop=True)
5292            .copy()
5293        )
5294
5295        ## find the ranges for non-outlier values and add them to sumdf
5296        mergelist = ["cluster"]
5297        for dim in dim_list:
5298            maxtag = dim + "_outmax"
5299            mintag = dim + "_outmin"
5300            mergelist.append(maxtag)
5301            mergelist.append(mintag)
5302            # Calculate outlier thresholds using vectorized operations
5303            sumdf[mintag] = sumdf[dim + "_mean"] - 3 * sumdf[dim + "_std"]
5304            sumdf[maxtag] = sumdf[dim + "_mean"] + 3 * sumdf[dim + "_std"]
5305            ## If NaN shows up anywhere in dim_min, dim_max calculations, value is set to NaN and it's
5306            ## not flagged. This happens when there's not enough values to compute median/std for that
5307            ## dimension therefore can't have outliers
5308
5309        ## add ranges to mfdf and identify mass features that fall outside the ranges
5310        # Merge without dropping NaN - we'll handle it per-dimension
5311        outdf = pd.merge(mfdf, sumdf[mergelist], on="cluster")
5312
5313        outtags = ["cluster"]
5314        for dim in dim_list:
5315            dimtag = dim + "_outlier"
5316            outtags.append(dimtag)
5317            maxtag = dim + "_outmax"
5318            mintag = dim + "_outmin"
5319            # Only flag as outlier if thresholds are valid (not NaN)
5320            outdf[dimtag] = np.where(
5321                (outdf[maxtag].notna() & outdf[mintag].notna())
5322                & ((outdf[dim] > outdf[maxtag]) | (outdf[dim] < outdf[mintag])),
5323                True,
5324                False,
5325            )
5326
5327        ## identify number of outliers in each cluster
5328        outliers = outdf[outtags]
5329        outliers = outliers.groupby(["cluster"]).sum()
5330
5331        ## plot number of clusters that contain any outliers
5332        fig = plt.figure()
5333        plt.bar(dim_list, outliers.sum().values, width=0.5)
5334        plt.xticks(rotation=90)
5335        plt.title("Frequency of outliers across all clusters by category")
5336
5337        return _finalize_plot(fig, return_fig=return_fig, path=path)
5338
5339    def _search_for_targeted_mass_features_in_sample(
5340        self, obj_idx, missingdf, cluster_dict, expand_on_miss=False, inplace=True
5341    ):
5342        """
5343        Helper method to search for missing mass features in a single sample.
5344
5345        Internal method called by fill_missing_cluster_features() to perform
5346        gap-filling for one sample in the collection.
5347
5348        Parameters
5349        ----------
5350        obj_idx : int
5351            Index of the sample being processed
5352        missingdf : pd.DataFrame
5353            DataFrame containing cluster information with columns:
5354            'cluster', 'sample_id_nunique', 'mz_min', 'mz_max',
5355            'scan_time_aligned_min', 'scan_time_aligned_max', 'mz_min_allowed',
5356            'mz_max_allowed', 'scan_time_aligned_min_allowed',
5357            'scan_time_aligned_max_allowed', 'missing_samples'
5358        cluster_dict : dict
5359            Pre-computed cluster feature dictionary to avoid recomputation
5360        expand_on_miss : bool
5361            If True, expands search window when no peak found initially
5362        inplace : bool
5363            If True, assigns induced_mass_features in place. If False, returns the
5364            induced features dictionary (for multiprocessing)
5365
5366        Returns
5367        -------
5368        dict or None
5369            If inplace=False, returns dictionary of induced mass features.
5370            Otherwise returns None and updates object in place.
5371        """
5372        ## Use the pre-computed cluster dictionary passed as parameter
5373
5374        ## to get clusters missing data based on sample index:
5375        sampledf = (
5376            missingdf[missingdf.missing_samples.apply(lambda x: obj_idx in x)]
5377            .reset_index(drop=True)
5378            .copy()
5379        )
5380
5381        # Skip if no missing features for this sample
5382        if len(sampledf) == 0:
5383            if not inplace:
5384                return {}
5385            return
5386
5387        self.load_raw_data(obj_idx, 1)
5388
5389        ## this is the line that bugs due to _ms_unprocessed not having key 1
5390        ms1df = self[obj_idx]._ms_unprocessed[1].copy()
5391        scan_df = self[obj_idx].scan_df[["scan", "scan_time_aligned"]]
5392        ms1df = pd.merge(ms1df, scan_df, on="scan")
5393
5394        # Pre-extract all values from sampledf to avoid repeated .iloc calls
5395        clusters = sampledf.cluster.values
5396        mz_mins = sampledf.mz_min.values
5397        mz_maxs = sampledf.mz_max.values
5398        st_mins = sampledf.scan_time_aligned_min.values
5399        st_maxs = sampledf.scan_time_aligned_max.values
5400
5401        if expand_on_miss:
5402            mz_mins_allowed = sampledf.mz_min_allowed.values
5403            mz_maxs_allowed = sampledf.mz_max_allowed.values
5404            st_mins_allowed = sampledf.sta_min_allowed.values
5405            st_maxs_allowed = sampledf.sta_max_allowed.values
5406
5407        # Pre-filter ms1df to reduce search space
5408        mz_global_min = mz_mins.min()
5409        mz_global_max = mz_maxs.max()
5410        st_global_min = st_mins.min()
5411        st_global_max = st_maxs.max()
5412
5413        if expand_on_miss:
5414            mz_global_min = min(mz_global_min, mz_mins_allowed.min())
5415            mz_global_max = max(mz_global_max, mz_maxs_allowed.max())
5416            st_global_min = min(st_global_min, st_mins_allowed.min())
5417            st_global_max = max(st_global_max, st_maxs_allowed.max())
5418
5419        ms1df_filtered = ms1df[
5420            (ms1df.mz >= mz_global_min)
5421            & (ms1df.mz <= mz_global_max)
5422            & (ms1df.scan_time_aligned >= st_global_min)
5423            & (ms1df.scan_time_aligned <= st_global_max)
5424        ].copy()
5425
5426        # Generate set_ids for all features
5427        set_ids = [f"c{clusters[i]}_{i}_i" for i in range(len(sampledf))]
5428
5429        # Use batch method to process all features at once
5430        if expand_on_miss:
5431            # First try with normal bounds
5432            peaks_dict = self[obj_idx].search_for_targeted_mass_features_batch(
5433                ms1df_filtered,
5434                mz_mins,
5435                mz_maxs,
5436                st_mins,
5437                st_maxs,
5438                set_ids,
5439                obj_idx=obj_idx,
5440                st_aligned=True,
5441            )
5442
5443            # Retry failed features with expanded bounds
5444            failed_indices = [
5445                i for i, sid in enumerate(set_ids) if peaks_dict[sid].apex_scan == -99
5446            ]
5447            if failed_indices:
5448                failed_ids = [set_ids[i] for i in failed_indices]
5449                retry_peaks = self[obj_idx].search_for_targeted_mass_features_batch(
5450                    ms1df_filtered,
5451                    mz_mins_allowed[failed_indices],
5452                    mz_maxs_allowed[failed_indices],
5453                    st_mins_allowed[failed_indices],
5454                    st_maxs_allowed[failed_indices],
5455                    failed_ids,
5456                    obj_idx=obj_idx,
5457                    st_aligned=True,
5458                )
5459                peaks_dict.update(retry_peaks)
5460        else:
5461            peaks_dict = self[obj_idx].search_for_targeted_mass_features_batch(
5462                ms1df_filtered,
5463                mz_mins,
5464                mz_maxs,
5465                st_mins,
5466                st_maxs,
5467                set_ids,
5468                obj_idx=obj_idx,
5469                st_aligned=True,
5470            )
5471
5472        # Assign peaks to induced_mass_features and cluster_dict
5473        for i in range(len(sampledf)):
5474            peak = peaks_dict[set_ids[i]]
5475            self[obj_idx].induced_mass_features[peak.id] = peak
5476            cluster_dict[clusters[i]] += [set_ids[i]]
5477
5478        # TODO KRH: Let's try to avoid these steps unless asked for by parameters to pick up speed
5479        if False:
5480            self[obj_idx].add_associated_ms1(induced_features=True)
5481            # need to set drop_if_fail to false for induced features as they will fail
5482            self[obj_idx].add_peak_metrics(induced_features=True)
5483
5484        # The duplicate scan only deletes keys from mass_features. Induced ids
5485        # are strings such as c12_0_i, so that scan does not remove them.
5486        self[obj_idx].integrate_mass_features(
5487            drop_if_fail=False, induced_features=True, drop_duplicates=False
5488        )
5489
5490        if not inplace:
5491            return self[obj_idx].induced_mass_features
5492
5493    def _missing_sample_ids_by_cluster(self, mfdf, cluster_ids):
5494        """Sample ids with no feature row for each cluster, in sample-id order.
5495
5496        Includes samples that have no mass features. Lookup matches numpy
5497        integers and integral floats to the cluster key.
5498        """
5499        all_sample_ids = list(range(len(self.samples)))
5500        if mfdf is None or len(mfdf) == 0:
5501            present_by_cluster = pd.Series(dtype=object)
5502        else:
5503            present_by_cluster = mfdf.groupby("cluster", sort=False)["sample_id"].agg(
5504                lambda ids: set(ids.to_numpy())
5505            )
5506        missing = []
5507        for cluster_id in cluster_ids:
5508            present = present_by_cluster.get(cluster_id, ())
5509            missing.append(
5510                [sample_id for sample_id in all_sample_ids if sample_id not in present]
5511            )
5512        return missing
5513
5514    def fill_missing_cluster_features(self):
5515        """
5516        Gap-filling for consensus mass features across collection samples.
5517
5518        For clusters present in multiple samples but missing from others, searches
5519        raw MS1 data to find peaks in expected m/z and retention time windows. This
5520        creates "induced" mass features for peaks that exist in the data but weren't
5521        detected in the initial peak detection.
5522
5523        Must be run after add_consensus_mass_features(). Results are accessible via
5524        induced_mass_features_dataframe property and included in collection_pivot_table
5525        and collection_consensus_report outputs.
5526
5527        Parameters
5528        ----------
5529        None
5530            Uses parameters from self.parameters.lcms_collection:
5531            - consensus_min_sample_fraction: Minimum fraction of samples (0-1) that must contain
5532              a cluster before gap-filling is attempted
5533            - gap_fill_expand_on_miss: If True, expands search window when no peak is found
5534
5535        Returns
5536        -------
5537        None
5538            Updates induced_mass_features attribute for each LCMSBase object and
5539            combines them into induced_mass_features_dataframe.
5540
5541        Raises
5542        ------
5543        ValueError
5544            If cluster_summary_dataframe is not set (must run add_consensus_mass_features first).
5545
5546        Notes
5547        -----
5548        - Loads raw MS1 data for each sample, which may be memory intensive
5549        - Induced features are integrated and metrics calculated automatically
5550        - Processing can be parallelized using parameters.lcms_collection.cores
5551
5552        See Also
5553        --------
5554        add_consensus_mass_features : Creates consensus features before gap-filling
5555        collection_pivot_table : Includes both regular and induced features
5556        collection_consensus_report : Reports on complete feature matrix
5557        """
5558
5559        # Validate prerequisites
5560        if (
5561            not hasattr(self, "cluster_summary_dataframe")
5562            or self.cluster_summary_dataframe is None
5563        ):
5564            raise ValueError(
5565                "cluster_summary_dataframe not found. Must run add_consensus_mass_features() first."
5566            )
5567
5568        # Get parameters from settings
5569        min_cluster_presence = (
5570            self.parameters.lcms_collection.consensus_min_sample_fraction
5571        )
5572        expand_on_miss = self.parameters.lcms_collection.gap_fill_expand_on_miss
5573
5574        # Validate parameters
5575        if not 0 <= min_cluster_presence <= 1:
5576            raise ValueError("consensus_min_sample_fraction must be between 0 and 1")
5577
5578        summarydf = self.cluster_summary_dataframe
5579        mfdf = self.mass_features_dataframe
5580
5581        sample_ct = len(self.samples)
5582
5583        # Identify clusters present in sufficient samples but not all samples
5584        missingdf = summarydf[
5585            [
5586                "cluster",
5587                "sample_id_nunique",
5588                "mz_min",
5589                "mz_max",
5590                "scan_time_aligned_min",
5591                "scan_time_aligned_max",
5592            ]
5593        ]
5594        missingdf = missingdf[
5595            missingdf.sample_id_nunique > min_cluster_presence * sample_ct
5596        ]
5597        missingdf = missingdf[missingdf.sample_id_nunique != sample_ct]
5598
5599        # Check if there are any clusters to gap-fill
5600        if len(missingdf) == 0:
5601            return
5602
5603        # One groupby. Sample ids stay in collection order, including samples
5604        # that have no mass features.
5605        missingdf["missing_samples"] = self._missing_sample_ids_by_cluster(
5606            mfdf, missingdf.cluster.to_numpy()
5607        )
5608
5609        # Calculate expanded search windows for expand_on_miss option
5610        mz_clu_tol = self.parameters.lcms_collection.consensus_mz_tol_ppm * 1e-6
5611        rt_clu_tol = self.parameters.lcms_collection.consensus_rt_tol
5612        missingdf["mz_max_allowed"] = missingdf.mz_max + mz_clu_tol * missingdf.mz_max
5613        missingdf["mz_min_allowed"] = missingdf.mz_min - mz_clu_tol * missingdf.mz_min
5614        # consensus_rt_tol is an absolute tolerance in minutes, same as clustering.
5615        missingdf["sta_max_allowed"] = missingdf.scan_time_aligned_max + rt_clu_tol
5616        missingdf["sta_min_allowed"] = missingdf.scan_time_aligned_min - rt_clu_tol
5617
5618        # Compute cluster dictionary once to avoid recomputing for each sample
5619        cluster_dict = self.cluster_feature_dictionary
5620
5621        # Process each sample to search for missing features
5622        if self.parameters.lcms_collection.cores == 1:
5623            for i in tqdm(range(sample_ct), desc="Gap-filling samples", unit="sample"):
5624                self._search_for_targeted_mass_features_in_sample(
5625                    i, missingdf, cluster_dict, expand_on_miss
5626                )
5627
5628        if self.parameters.lcms_collection.cores > 1:
5629            if self.parameters.lcms_collection.cores > len(self):
5630                ncores = len(self)
5631            else:
5632                ncores = self.parameters.lcms_collection.cores
5633            with spawn_pool(ncores) as pool:
5634                mp_result = pool.starmap(
5635                    self._search_for_targeted_mass_features_in_sample,
5636                    [
5637                        (x, missingdf, cluster_dict, expand_on_miss, False)
5638                        for x in range(sample_ct)
5639                    ],
5640                )
5641
5642            for i in tqdm(
5643                range(sample_ct), desc="Collecting gap-filled features", unit="sample"
5644            ):
5645                self[i].induced_mass_features = mp_result[i]
5646
5647        self._combine_mass_features(induced_features=True)
5648
5649        # Mark that gap-filling has been performed
5650        self.missing_mass_features_searched = True
5651
5652        for sample_name in self.samples:
5653            self._lcms[sample_name].mass_features = {}
5654
5655    def process_samples_pipeline(
5656        self, operations, description=None, keep_raw_data=False, show_progress=True
5657    ):
5658        """
5659        Execute a pipeline of operations on all samples in parallel.
5660
5661        This method provides a flexible framework for performing multiple
5662        sample-level operations in a single parallelized pass, which is more
5663        efficient than calling separate methods sequentially.
5664
5665        Parameters
5666        ----------
5667        operations : list of SampleOperation
5668            List of operations to perform on each sample, in order.
5669            Each operation should be an instance of a class derived from
5670            SampleOperation (see lc_calc_operations module).
5671        description : str or None, optional
5672            Progress bar description. If None, automatically generates description
5673            from operation descriptions (e.g., "gap-filling, reloading features").
5674            Default is None.
5675        keep_raw_data : bool, optional
5676            If True, keeps raw MS data loaded in memory after pipeline completes.
5677            If False, cleans up raw data to free memory. Default is False.
5678        show_progress : bool, optional
5679            If True, displays progress bars during processing. If False, runs silently.
5680            Default is True.
5681
5682        Returns
5683        -------
5684        dict
5685            Dictionary with results from pipeline execution, keyed by operation name.
5686            Structure: {operation_name: {sample_id: result, ...}, ...}
5687
5688        Raises
5689        ------
5690        ValueError
5691            If operations list is empty or contains invalid operations.
5692
5693        Notes
5694        -----
5695        - Operations are executed sequentially within each sample
5696        - Samples are processed in parallel based on parameters.lcms_collection.cores
5697        - Each operation can have conditional execution via can_execute()
5698        - Results are collected back via collect_results() method of each operation
5699        - Failed operations for a sample are logged but don't halt processing
5700        - Raw MS data loaded by operations is automatically cleaned up unless keep_raw_data=True
5701
5702        Examples
5703        --------
5704        >>> from corems.mass_spectra.calc.lc_calc_operations import (
5705        ...     GapFillOperation, ReloadFeaturesOperation
5706        ... )
5707        >>> ops = [
5708        ...     GapFillOperation('gap_fill', expand_on_miss=True),
5709        ...     ReloadFeaturesOperation('reload', add_ms2=True)
5710        ... ]
5711        >>> results = lcms_collection.process_samples_pipeline(ops)
5712
5713        See Also
5714        --------
5715        lc_calc_operations : Module containing built-in operation classes
5716        fill_and_process_features : Convenience method combining common operations
5717        """
5718        from corems.mass_spectra.calc.lc_calc_operations import SampleOperation
5719
5720        # Validate operations
5721        if not operations or len(operations) == 0:
5722            raise ValueError("operations list cannot be empty")
5723
5724        for op in operations:
5725            if not isinstance(op, SampleOperation):
5726                raise ValueError(
5727                    f"All operations must be SampleOperation instances, got {type(op)}"
5728                )
5729
5730        # Generate description from operations if not provided
5731        if description is None:
5732            operation_descriptions = [op.description for op in operations]
5733            description = ", ".join(operation_descriptions).capitalize()
5734
5735        # Prepare runtime parameters for each operation
5736        # This is where we gather collection-level data that operations need
5737        runtime_params = self._prepare_pipeline_runtime_params(operations)
5738        runtime_params["keep_raw_data"] = keep_raw_data
5739
5740        # Execute pipeline
5741        sample_ct = len(self.samples)
5742
5743        if self.parameters.lcms_collection.cores == 1:
5744            # Serial processing
5745            results_by_operation = {op.name: {} for op in operations}
5746
5747            if show_progress:
5748                from tqdm import tqdm
5749
5750                # Print description on its own line before progress bar
5751                print(f"\n{description.capitalize()}:")
5752                iterator = tqdm(range(sample_ct), unit="sample", ncols=80)
5753            else:
5754                iterator = range(sample_ct)
5755
5756            for sample_id in iterator:
5757                sample_results = self._execute_sample_pipeline(
5758                    sample_id, operations, runtime_params, inplace=True
5759                )
5760                # Collect results (collect_results already called in _execute_sample_pipeline when inplace=True)
5761                # Skip 'sample_id' key which is added for tracking
5762                for op_name, result in sample_results.items():
5763                    if op_name != "sample_id":
5764                        results_by_operation[op_name][sample_id] = result
5765        else:
5766            # Parallel processing
5767            if self.parameters.lcms_collection.cores > sample_ct:
5768                ncores = sample_ct
5769            else:
5770                ncores = self.parameters.lcms_collection.cores
5771
5772            # Build arguments for each sample
5773            args_list = [
5774                (sample_id, operations, runtime_params, False)
5775                for sample_id in range(sample_ct)
5776            ]
5777
5778            # Execute in parallel with progress tracking
5779            results_by_operation = {op.name: {} for op in operations}
5780
5781            with spawn_pool(ncores) as pool:
5782                if show_progress:
5783                    from tqdm import tqdm
5784                    import time
5785
5786                    # Use starmap_async for parallel execution with progress tracking
5787                    async_result = pool.starmap_async(
5788                        self._execute_sample_pipeline, args_list
5789                    )
5790
5791                    # Poll for completion and update progress bar
5792                    print(description)
5793                    pbar = tqdm(
5794                        total=sample_ct,
5795                        desc="",
5796                        unit="sample",
5797                        position=0,
5798                        leave=True,
5799                        dynamic_ncols=True,
5800                    )
5801                    prev_completed = 0
5802                    while not async_result.ready():
5803                        # Get number of completed tasks by checking remaining
5804                        completed = sample_ct - async_result._number_left
5805                        if completed > prev_completed:
5806                            pbar.update(completed - prev_completed)
5807                            prev_completed = completed
5808                        time.sleep(0.5)  # Poll every 500ms to avoid spam
5809
5810                    # Final update to 100%
5811                    if prev_completed < sample_ct:
5812                        pbar.update(sample_ct - prev_completed)
5813                    pbar.close()
5814
5815                    # Get all results
5816                    mp_results = async_result.get()
5817                else:
5818                    # Execute without progress
5819                    mp_results = pool.starmap(self._execute_sample_pipeline, args_list)
5820
5821            # Collect results back into collection
5822            for result in mp_results:
5823                sample_id = result.get("sample_id")
5824                for op in operations:
5825                    op_result = result.get(op.name)
5826                    if op_result is not None:
5827                        op.collect_results(sample_id, op_result, self)
5828                        results_by_operation[op.name][sample_id] = op_result
5829
5830        return results_by_operation
5831
5832    def _prepare_pipeline_runtime_params(self, operations):
5833        """
5834        Prepare runtime parameters needed by operations in the pipeline.
5835
5836        This method gathers collection-level data that operations need,
5837        such as cluster information for gap-filling or mf_ids for reloading.
5838
5839        Parameters
5840        ----------
5841        operations : list of SampleOperation
5842            List of operations that will be executed
5843
5844        Returns
5845        -------
5846        dict
5847            Dictionary of runtime parameters for operations
5848        """
5849        from corems.mass_spectra.calc.lc_calc_operations import (
5850            GapFillOperation,
5851            ReloadFeaturesOperation,
5852            MS2SpectralSearchOperation,
5853            LoadEICsOperation,
5854        )
5855
5856        runtime_params = {}
5857
5858        # Check if any operation needs gap-fill parameters
5859        needs_gap_fill = any(isinstance(op, GapFillOperation) for op in operations)
5860        if needs_gap_fill:
5861            # Prepare gap-fill parameters (same as fill_missing_cluster_features)
5862            min_cluster_presence = (
5863                self.parameters.lcms_collection.consensus_min_sample_fraction
5864            )
5865            expand_on_miss = self.parameters.lcms_collection.gap_fill_expand_on_miss
5866
5867            summarydf = self.cluster_summary_dataframe
5868            mfdf = self.mass_features_dataframe
5869            sample_ct = len(self.samples)
5870
5871            # Identify clusters needing gap-filling
5872            # Note: cluster_summary_dataframe has 'cluster' as index, need to reset it
5873            missingdf = summarydf.reset_index()[
5874                [
5875                    "cluster",
5876                    "sample_id_nunique",
5877                    "mz_min",
5878                    "mz_max",
5879                    "scan_time_aligned_min",
5880                    "scan_time_aligned_max",
5881                ]
5882            ].copy()
5883            missingdf = missingdf[
5884                missingdf.sample_id_nunique > min_cluster_presence * sample_ct
5885            ]
5886            missingdf = missingdf[missingdf.sample_id_nunique != sample_ct]
5887
5888            if len(missingdf) > 0:
5889                # One groupby. Sample ids stay in collection order, including samples
5890                # that have no mass features.
5891                missingdf["missing_samples"] = self._missing_sample_ids_by_cluster(
5892                    mfdf, missingdf.cluster.to_numpy()
5893                )
5894
5895                # Calculate expanded search windows
5896                mz_clu_tol = self.parameters.lcms_collection.consensus_mz_tol_ppm * 1e-6
5897                rt_clu_tol = self.parameters.lcms_collection.consensus_rt_tol
5898                missingdf["mz_max_allowed"] = (
5899                    missingdf.mz_max + mz_clu_tol * missingdf.mz_max
5900                )
5901                missingdf["mz_min_allowed"] = (
5902                    missingdf.mz_min - mz_clu_tol * missingdf.mz_min
5903                )
5904                # consensus_rt_tol is an absolute tolerance in minutes, same as clustering.
5905                missingdf["sta_max_allowed"] = (
5906                    missingdf.scan_time_aligned_max + rt_clu_tol
5907                )
5908                missingdf["sta_min_allowed"] = (
5909                    missingdf.scan_time_aligned_min - rt_clu_tol
5910                )
5911
5912                runtime_params["missingdf"] = missingdf
5913                runtime_params["cluster_dict"] = self.cluster_feature_dictionary
5914                runtime_params["expand_on_miss"] = expand_on_miss
5915
5916        # Check if any operation needs reload parameters
5917        needs_reload = any(isinstance(op, ReloadFeaturesOperation) for op in operations)
5918        if needs_reload:
5919            # Use DRY helper method to build sample_mf_map
5920            sample_mf_map = self.get_sample_mf_map_for_representatives(
5921                include_cluster_id=False
5922            )
5923            runtime_params["sample_mf_map"] = sample_mf_map
5924
5925        # Check if any operation needs MS2 spectral search parameters
5926        needs_ms2_search = any(
5927            isinstance(op, MS2SpectralSearchOperation) for op in operations
5928        )
5929        if needs_ms2_search:
5930            # Pass through pre-prepared spectral library
5931            if hasattr(self, "_spectral_lib") and self._spectral_lib is not None:
5932                runtime_params["fe_lib"] = self._spectral_lib
5933            if hasattr(self, "_spectral_search_molecular_metadata"):
5934                runtime_params["molecular_metadata"] = (
5935                    self._spectral_search_molecular_metadata
5936                )
5937
5938        # Check if any operation needs EIC loading parameters
5939        needs_eic_loading = any(isinstance(op, LoadEICsOperation) for op in operations)
5940        if needs_eic_loading:
5941            # Map sample_id -> m/z list for LoadEICsOperation.
5942            # Must use collection *dataframes*, not only in-memory mass_features:
5943            # after load_representatives=True each sample holds only a sparse set
5944            # of representative objects, so object-only walks miss most cluster
5945            # members and plot_cluster shows incomplete multi-sample EICs (#258).
5946            runtime_params["cluster_mz_dict"] = (
5947                self._build_cluster_mz_dict_for_eic_loading()
5948            )
5949
5950        return runtime_params
5951
5952    def _build_cluster_mz_dict_for_eic_loading(self):
5953        """Build sample_id -> list of EIC m/z targets for all clustered features.
5954
5955        Uses ``mass_features_dataframe`` (and induced dataframe if present).
5956        Prefers ``_eic_mz``, falls back to ``mz``. Does not require feature
5957        objects to be loaded in ``sample.mass_features``.
5958
5959        Returns
5960        -------
5961        dict
5962            Mapping of sample_id (int) to unique float m/z values to load.
5963        """
5964        cluster_mz_dict = {}
5965
5966        def _add_from_df(df):
5967            if df is None or len(df) == 0:
5968                return
5969            if "cluster" in df.columns:
5970                df = df[df["cluster"].notna()]
5971            if len(df) == 0 or "sample_id" not in df.columns:
5972                return
5973            for sample_id, sample_df in df.groupby("sample_id"):
5974                mzs = []
5975                if "_eic_mz" in sample_df.columns:
5976                    mzs.extend(sample_df["_eic_mz"].dropna().tolist())
5977                if "mz" in sample_df.columns:
5978                    # Include mz for rows with missing _eic_mz
5979                    if "_eic_mz" in sample_df.columns:
5980                        missing = sample_df["_eic_mz"].isna()
5981                        mzs.extend(sample_df.loc[missing, "mz"].dropna().tolist())
5982                    else:
5983                        mzs.extend(sample_df["mz"].dropna().tolist())
5984                sid = int(sample_id)
5985                existing = cluster_mz_dict.get(sid, [])
5986                cluster_mz_dict[sid] = list(set(existing + [float(m) for m in mzs]))
5987
5988        # Regular clustered features only: these are what were exported to per-sample
5989        # HDF5 EICs. Induced (gap-filled) features get EICs at gap-fill time on the
5990        # sample object; they are typically not present in the original HDF5 eics group.
5991        _add_from_df(self.mass_features_dataframe)
5992
5993        return cluster_mz_dict
5994
5995    def _execute_sample_pipeline(
5996        self, sample_id, operations, runtime_params, inplace=True
5997    ):
5998        """
5999        Execute a pipeline of operations on a single sample.
6000
6001        This is the worker function called (potentially in parallel) for each sample.
6002
6003        Parameters
6004        ----------
6005        sample_id : int
6006            Sample ID to process
6007        operations : list of SampleOperation
6008            Operations to execute in order
6009        runtime_params : dict
6010            Runtime parameters prepared by _prepare_pipeline_runtime_params
6011        inplace : bool, optional
6012            If True, updates sample in place. If False, returns results for
6013            multiprocessing. Default is True.
6014
6015        Returns
6016        -------
6017        dict
6018            Dictionary with results from each operation, keyed by operation name.
6019            If inplace=True, returns results that need to be collected.
6020            If inplace=False, returns all results for multiprocessing collection.
6021        """
6022        results = {}
6023
6024        # Check if any operations need raw MS data
6025        needs_raw_data = {}  # {ms_level: True/False}
6026        for op in operations:
6027            needs_raw, ms_level = op.needs_raw_ms_data()
6028            if needs_raw and ms_level:
6029                needs_raw_data[ms_level] = True
6030
6031        # Load raw data once if any operations need it
6032        # Note: For gap-filling, it loads data internally, so we just track it here
6033        for ms_level in needs_raw_data.keys():
6034            # Gap-filling loads its own data, but we want to keep track that it's loaded
6035            # Other operations can then use the loaded data
6036            pass
6037
6038        for op in operations:
6039            # Check if operation can execute on this sample
6040            sample = self[sample_id]
6041            if not op.can_execute(sample, self):
6042                # Skip this operation for this sample if prerequisites aren't met
6043                # This allows processing to continue for samples that don't have
6044                # all required data (e.g., MS2 spectra)
6045                results[op.name] = None
6046                continue
6047
6048            # Prepare operation-specific runtime params
6049            op_runtime_params = {}
6050
6051            # Add gap-fill params if this is a gap-fill operation
6052            from corems.mass_spectra.calc.lc_calc_operations import (
6053                GapFillOperation,
6054                ReloadFeaturesOperation,
6055                MS2SpectralSearchOperation,
6056                LoadEICsOperation,
6057            )
6058
6059            if isinstance(op, GapFillOperation):
6060                if "missingdf" in runtime_params:
6061                    op_runtime_params["missingdf"] = runtime_params["missingdf"]
6062                    op_runtime_params["cluster_dict"] = runtime_params["cluster_dict"]
6063                    op_runtime_params["expand_on_miss"] = runtime_params[
6064                        "expand_on_miss"
6065                    ]
6066
6067            elif isinstance(op, ReloadFeaturesOperation):
6068                if "sample_mf_map" in runtime_params:
6069                    sample_mf_map = runtime_params["sample_mf_map"]
6070                    # Always pass mf_ids_to_load to ensure we only load what's needed
6071                    # If sample not in map, it has no representatives - pass empty list
6072                    op_runtime_params["mf_ids_to_load"] = sample_mf_map.get(
6073                        sample_id, []
6074                    )
6075
6076            elif isinstance(op, MS2SpectralSearchOperation):
6077                # Add MS2 spectral search parameters
6078                if "fe_lib" in runtime_params:
6079                    op_runtime_params["fe_lib"] = runtime_params["fe_lib"]
6080                if "molecular_metadata" in runtime_params:
6081                    op_runtime_params["molecular_metadata"] = runtime_params[
6082                        "molecular_metadata"
6083                    ]
6084
6085            elif isinstance(op, LoadEICsOperation):
6086                # Add EIC loading parameters
6087                if "cluster_mz_dict" in runtime_params:
6088                    op_runtime_params["cluster_mz_dict"] = runtime_params[
6089                        "cluster_mz_dict"
6090                    ]
6091
6092            # Execute the operation
6093            result = op.execute(sample_id, self, **op_runtime_params)
6094            results[op.name] = result
6095
6096            # If inplace, collect immediately
6097            if inplace and result is not None:
6098                op.collect_results(sample_id, result, self)
6099
6100        # Clean up raw data if requested
6101        keep_raw_data = runtime_params.get("keep_raw_data", False)
6102        if not keep_raw_data:
6103            for ms_level in needs_raw_data.keys():
6104                if ms_level in self[sample_id]._ms_unprocessed:
6105                    del self[sample_id]._ms_unprocessed[ms_level]
6106
6107        # Include sample_id in results for tracking (especially important for imap_unordered)
6108        results["sample_id"] = sample_id
6109        return results
6110
6111    def process_consensus_features(
6112        self,
6113        load_representatives=True,
6114        perform_gap_filling=True,
6115        add_ms1=False,
6116        add_ms2=False,
6117        ms2_scan_filter=None,
6118        molecular_formula_search=False,
6119        ms2_spectral_search=False,
6120        spectral_lib=None,
6121        molecular_metadata=None,
6122        gather_eics=False,
6123        group_features=False,
6124        keep_raw_data=False,
6125        show_progress=True,
6126    ):
6127        """
6128        Process consensus mass features across the collection in a single parallelized pass.
6129
6130        This method provides a convenient interface to the sample processing pipeline,
6131        allowing multiple operations (gap-filling, feature reloading, MS1/MS2 association,
6132        molecular formula search, and MS2 spectral search) to be performed efficiently in
6133        a single pass through all samples.
6134
6135        Parameters
6136        ----------
6137        load_representatives : bool, optional
6138            If True, loads representative mass features from HDF5. Default is True.
6139        perform_gap_filling : bool, optional
6140            If True, performs gap-filling for missing cluster features. Default is True.
6141            This operation loads raw MS1 data which can be reused by subsequent operations.
6142        group_features : bool, optional
6143            If True, run ``group_consensus_features()`` after gap-fill (and before
6144            molecular formula / MS2 search when those are also enabled).
6145            Collection-level natural-abundance **isotope and adduct** family
6146            labeling (singly charged only; not tracer/enriched). Default is False.
6147        add_ms1 : bool, optional
6148            If True and load_representatives=True, associates MS1 spectra with
6149            loaded features. Automatically uses raw data from gap-filling if available,
6150            otherwise uses parser. Spectrum mode is auto-detected. Default is False.
6151        add_ms2 : bool, optional
6152            If True and load_representatives=True, associates MS2 spectra with
6153            loaded features and automatically processes them. Spectrum mode is auto-detected. Default is False.
6154        ms2_scan_filter : str or None, optional
6155            Filter string for MS2 scans (e.g., 'hcd'). Default is None.
6156        molecular_formula_search : bool, optional
6157            If True, performs molecular formula search on mass features using
6158            associated MS1 spectra. Requires add_ms1=True or that MS1 spectra
6159            are already associated. Uses parameters from
6160            parameters.mass_spectrum["ms1"].molecular_search. Default is False.
6161        ms2_spectral_search : bool, optional
6162            If True, performs MS2 spectral library search using FlashEntropy.
6163            Requires add_ms2=True and spectral_lib to be provided. Default is False.
6164        spectral_lib : FlashEntropy library, optional
6165            Pre-prepared FlashEntropy spectral library for MS2 search.
6166            Create using MSPInterface.get_metabolomics_spectra_library().
6167            Required if ms2_spectral_search=True. Default is None.
6168        molecular_metadata : pd.DataFrame, optional
6169            Molecular metadata corresponding to spectral_lib.
6170            Returned from MSPInterface.get_metabolomics_spectra_library().
6171            Stored as self.spectral_search_molecular_metadata for later export.
6172            Default is None.
6173        gather_eics : bool, optional
6174            If True, loads extracted ion chromatograms (EICs) from HDF5 for all
6175            mass features with assigned cluster_index (including gap-filled features).
6176            Enables access to EICs via get_eics_for_cluster(cluster_id) method.
6177            Requires that EICs were previously exported with export_eics=True.
6178            Default is False.
6179        keep_raw_data : bool, optional
6180            If True, keeps raw MS data loaded in memory after pipeline completes.
6181            If False, cleans up raw data to free memory. Default is False.
6182        show_progress : bool, optional
6183            If True, displays progress bars during processing. If False, runs silently.
6184            Default is True.
6185
6186        Returns
6187        -------
6188        dict
6189            Dictionary with pipeline results. Keys include:
6190            - 'gap_fill': dict mapping sample_id to induced mass features (if gap-filling)
6191            - 'reload': dict mapping sample_id to reloaded mass features (if reloading)
6192            - 'mf_search': dict mapping sample_id to number of features searched (if molecular formula search)
6193            - 'ms2_search': dict mapping sample_id to number of spectra searched (if MS2 spectral search)
6194
6195        Raises
6196        ------
6197        ValueError
6198            If neither operation is enabled, or if required parameters are missing.
6199
6200        Notes
6201        -----
6202        - Must run add_consensus_mass_features() before calling this method
6203        - Processes samples in parallel based on parameters.lcms_collection.cores
6204        - Raw MS1 data loaded by gap-filling is automatically reused by MS1 association
6205        - MS2 spectral search requires add_ms2=True and msp_file_path
6206        - FlashEntropy library is created once and reused across all samples
6207        - More efficient than calling individual methods separately
6208        - After gap-filling, sets missing_mass_features_searched = True
6209        - Mass features remain loaded in memory for downstream processing
6210        - For more advanced workflows, use process_samples_pipeline() directly
6211
6212        Examples
6213        --------
6214        >>> # Prepare spectral library for MS2 search
6215        >>> from corems.molecular_id.search.database_interfaces import MSPInterface
6216        >>> my_msp = MSPInterface(file_path='path/to/library.msp')
6217        >>> spectral_lib, molecular_metadata = my_msp.get_metabolomics_spectra_library(
6218        ...     polarity='negative',
6219        ...     format='flashentropy',
6220        ...     normalize=True,
6221        ...     fe_kwargs={
6222        ...         'normalize_intensity': True,
6223        ...         'min_ms2_difference_in_da': 0.02,
6224        ...         'max_ms2_tolerance_in_da': 0.01,
6225        ...         'max_indexed_mz': 3000,
6226        ...         'precursor_ions_removal_da': None,
6227        ...         'noise_threshold': 0,
6228        ...     }
6229        ... )
6230        >>>
6231        >>> # Gap-fill, reload with MS1/MS2, perform molecular formula and spectral search
6232        >>> results = lcms_collection.process_consensus_features(
6233        ...     load_representatives=True,
6234        ...     perform_gap_filling=True,
6235        ...     add_ms1=True,
6236        ...     add_ms2=True,
6237        ...     molecular_formula_search=True,
6238        ...     ms2_spectral_search=True,
6239        ...     spectral_lib=spectral_lib,
6240        ...     molecular_metadata=molecular_metadata
6241        ... )
6242
6243        See Also
6244        --------
6245        process_samples_pipeline : Generic pipeline executor for custom workflows
6246        fill_missing_cluster_features : Original gap-filling method
6247        reload_representative_mass_features : Original reload method
6248        """
6249        from corems.mass_spectra.calc.lc_calc_operations import (
6250            GapFillOperation,
6251            ReloadFeaturesOperation,
6252            MolecularFormulaSearchOperation,
6253            MS2SpectralSearchOperation,
6254            LoadEICsOperation,
6255        )
6256
6257        # Validate that at least one meaningful operation is enabled
6258        has_operations = (
6259            perform_gap_filling
6260            or load_representatives
6261            or molecular_formula_search
6262            or ms2_spectral_search
6263            or gather_eics
6264            or add_ms1
6265            or add_ms2
6266            or group_features
6267        )
6268
6269        if not has_operations:
6270            raise ValueError(
6271                "At least one operation must be enabled: perform_gap_filling, load_representatives, "
6272                "molecular_formula_search, ms2_spectral_search, gather_eics, add_ms1, add_ms2, "
6273                "or group_features"
6274            )
6275
6276        # Validate prerequisites for gap-filling
6277        if perform_gap_filling:
6278            if (
6279                not hasattr(self, "cluster_summary_dataframe")
6280                or self.cluster_summary_dataframe is None
6281            ):
6282                raise ValueError(
6283                    "Cannot perform gap-filling: cluster_summary_dataframe not set. "
6284                    "You must run add_consensus_mass_features() before calling process_consensus_features()."
6285                )
6286
6287        # Validate prerequisites for MS2 spectral search
6288        if ms2_spectral_search:
6289            if spectral_lib is None:
6290                raise ValueError(
6291                    "MS2 spectral search requires spectral_lib to be provided. "
6292                    "Create it using MSPInterface.get_metabolomics_spectra_library() before calling this method."
6293                )
6294            # Check if mass features will be loaded OR are already loaded
6295            # (The operation's can_execute will check if MS2 spectra are actually present)
6296            if not load_representatives and not perform_gap_filling:
6297                # Check if at least one sample has mass features loaded
6298                # This allows MS2 search on already-loaded features
6299                has_loaded_features = any(
6300                    len(self[i].mass_features) > 0
6301                    if hasattr(self[i], "mass_features")
6302                    and self[i].mass_features is not None
6303                    else False
6304                    for i in range(len(self.samples))
6305                )
6306                if not has_loaded_features:
6307                    raise ValueError(
6308                        "MS2 spectral search requires mass features to be loaded. "
6309                        "Either set load_representatives=True or perform_gap_filling=True to load them, "
6310                        "or load them in a previous call to process_consensus_features() before calling "
6311                        "with ms2_spectral_search=True."
6312                    )
6313
6314        # Build pipeline in two phases so group_features runs after gap-fill
6315        # and before formula / MS2 search when those are requested in the same call.
6316        pre_operations = []
6317        post_operations = []
6318
6319        if perform_gap_filling:
6320            expand_on_miss = self.parameters.lcms_collection.gap_fill_expand_on_miss
6321            pre_operations.append(
6322                GapFillOperation("gap_fill", expand_on_miss=expand_on_miss)
6323            )
6324
6325        if load_representatives:
6326            # MS1/MS2 association for representatives can happen in pre-pass;
6327            # formula/MS2 library search stay in post-pass after grouping.
6328            pre_operations.append(
6329                ReloadFeaturesOperation(
6330                    "reload",
6331                    add_ms1=add_ms1,
6332                    add_ms2=add_ms2,
6333                    auto_process_ms2=add_ms2,
6334                    ms2_scan_filter=ms2_scan_filter,
6335                )
6336            )
6337
6338        if molecular_formula_search:
6339            post_operations.append(MolecularFormulaSearchOperation("mf_search"))
6340
6341        if ms2_spectral_search:
6342            post_operations.append(
6343                MS2SpectralSearchOperation(
6344                    "ms2_search", ms2_scan_filter=ms2_scan_filter
6345                )
6346            )
6347            self._spectral_lib = spectral_lib
6348            self._spectral_search_molecular_metadata = molecular_metadata
6349
6350        if gather_eics:
6351            # After grouping when grouping is requested; otherwise with pre-pass
6352            # (or post-pass if only annotation ops are enabled).
6353            if group_features or not pre_operations:
6354                post_operations.append(LoadEICsOperation("load_eics"))
6355            else:
6356                pre_operations.append(LoadEICsOperation("load_eics"))
6357
6358        results = {}
6359
6360        def _run_ops(operations):
6361            if not operations:
6362                return {}
6363            return self.process_samples_pipeline(
6364                operations, keep_raw_data=keep_raw_data, show_progress=show_progress
6365            )
6366
6367        # Phase 1: gap-fill / reload (and eics if not grouping)
6368        if pre_operations:
6369            results.update(_run_ops(pre_operations))
6370
6371        # Post-processing for gap-fill must happen before grouping so heights include induced
6372        if perform_gap_filling:
6373            self._combine_mass_features(induced_features=True)
6374            self.missing_mass_features_searched = True
6375
6376            if (
6377                self.induced_mass_features_dataframe is not None
6378                and len(self.induced_mass_features_dataframe) > 0
6379            ):
6380                eics_mz = []
6381                for i, row in self.induced_mass_features_dataframe.iterrows():
6382                    sample_id = row["sample_id"]
6383                    sample = self[sample_id]
6384                    if row["mf_id"] in sample.induced_mass_features.keys():
6385                        eic_mz = sample.induced_mass_features[row["mf_id"]]._eic_mz
6386                        eics_mz.append(eic_mz)
6387                    else:
6388                        eics_mz.append(None)
6389                self.induced_mass_features_dataframe["_eic_mz"] = eics_mz
6390
6391        # Collection-level isotope + adduct grouping
6392        if group_features:
6393            t_group = time.perf_counter()
6394            labels = self.group_consensus_features()
6395            group_wall = time.perf_counter() - t_group
6396            results["group_features"] = {
6397                "n_clusters": len(labels),
6398                "n_grouped": int(labels["feature_group_id"].notna().sum()),
6399                "n_groups": int(labels["feature_group_id"].nunique(dropna=True)),
6400                "wall_seconds": group_wall,
6401                "stage_timings": getattr(self, "_last_feature_group_timings", None),
6402            }
6403
6404        # Phase 2: annotation / optional EICs after grouping
6405        if post_operations:
6406            # If only post_operations and no pre, still need a pipeline run
6407            post_results = _run_ops(post_operations)
6408            for k, v in post_results.items():
6409                results[k] = v
6410
6411        # Store molecular metadata if spectral search was performed
6412        if ms2_spectral_search and hasattr(self, "_spectral_search_molecular_metadata"):
6413            self.spectral_search_molecular_metadata = (
6414                self._spectral_search_molecular_metadata
6415            )
6416
6417        # Associate EICs while induced feature objects still exist (before any clear).
6418        if gather_eics:
6419            print("\nAssociating EICs with mass features:")
6420            from tqdm import tqdm
6421
6422            for sample_id in tqdm(range(len(self.samples)), unit="sample", ncols=80):
6423                sample = self[sample_id]
6424                if sample.eics:
6425                    sample.associate_eics_with_mass_features(induced=False)
6426                    sample.associate_eics_with_mass_features(induced=True)
6427
6428        # Drop induced feature objects to free memory after grouping + post-ops
6429        if perform_gap_filling:
6430            for sample_name in self.samples:
6431                self._lcms[sample_name].induced_mass_features = {}
6432
6433        return results