Results

AnimalOrganizer

Bases: AnimalFeatureParser

Source code in pythoneeg/visualization/results.py

class AnimalOrganizer(AnimalFeatureParser):
    def __init__(
        self,
        base_folder_path,
        anim_id: str,
        day_sep: str | None = None,
        mode: Literal["nest", "concat", "base", "noday"] = "concat",
        assume_from_number=False,
        skip_days: list[str] = [],
        truncate: bool | int = False,
        lro_kwargs: dict = {},
    ) -> None:
        """
        AnimalOrganizer is used to organize data from a single animal into a format that can be used for analysis.
        It is used to organize data from a single animal into a format that can be used for analysis.

        Args:
            base_folder_path (str): The path to the base folder of the animal data.
            anim_id (str): The ID of the animal. This should correspond to only one animal.
            day_sep (str, optional): Separator for day in folder name. Set to None or empty string to get all folders. Defaults to None.
            mode (Literal["nest", "concat", "base", "noday"], optional): The mode of the AnimalOrganizer. Defaults to "concat".
                "nest": base_folder_path / animal_id / date_format
                "concat": base_folder_path / animal_id*date_format
                "base": base_folder_path
                "noday": base_folder_path / animal_id
            assume_from_number (bool, optional): Whether to assume the animal ID is a number. Defaults to False.
            skip_days (list[str], optional): The days to skip. Defaults to [].
            truncate (bool|int, optional): Whether to truncate the data. Defaults to False.
            lro_kwargs (dict, optional): Keyword arguments for LongRecordingOrganizer. Defaults to {}.
        """

        self.base_folder_path = Path(base_folder_path)
        self.anim_id = anim_id
        self.animal_param = [anim_id]
        self.day_sep = day_sep
        self.read_mode = mode
        self.assume_from_number = assume_from_number

        match mode:
            case "nest":
                self.bin_folder_pattern = self.base_folder_path / f"*{self.anim_id}*" / "*"
            case "concat" | "noday":
                self.bin_folder_pattern = self.base_folder_path / f"*{self.anim_id}*"
                # self.bin_folder_pat = self.base_folder_path / f"*{self.anim_id}*{self.date_format}*"
            case "base":
                self.bin_folder_pattern = self.base_folder_path
            # case 'noday':
            #     self.bin_folder_pat = self.base_folder_path / f"*{self.anim_id}*"
            case _:
                raise ValueError(f"Invalid mode: {mode}")

        self._bin_folders = glob.glob(str(self.bin_folder_pattern))
        # if mode != 'noday':
        #     self.__bin_folders = [x for x in self.__bin_folders if datetime.strptime(Path(x).name, self.date_format)]
        truncate = core.utils.parse_truncate(truncate)
        if truncate:
            warnings.warn(f"AnimalOrganizer will be truncated to the first {truncate} LongRecordings")
            self._bin_folders = self._bin_folders[:truncate]
        self._bin_folders = [x for x in self._bin_folders if not any(y in x for y in skip_days)]
        self.bin_folder_names = [Path(x).name for x in self._bin_folders]
        logging.info(f"bin_folder_pattern: {self.bin_folder_pattern}")
        logging.info(f"self._bin_folders: {self._bin_folders}")
        logging.info(f"self.bin_folder_names: {self.bin_folder_names}")

        if mode == "noday" and len(self._bin_folders) > 1:
            raise ValueError(f"Animal ID '{self.anim_id}' is not unique, found: {', '.join(self._bin_folders)}")
        elif len(self._bin_folders) == 0:
            raise ValueError(f"No files found for animal ID {self.anim_id}")

        animalday_dicts = [
            core.parse_path_to_animalday(e, animal_param=self.animal_param, day_sep=self.day_sep, mode=self.read_mode)
            for e in self._bin_folders
        ]
        self.animaldays = [x["animalday"] for x in animalday_dicts]
        logging.info(f"self.animaldays: {self.animaldays}")

        genotypes = [x["genotype"] for x in animalday_dicts]
        if len(set(genotypes)) > 1:
            warnings.warn(f"Inconsistent genotypes in {genotypes}")
        self.genotype = genotypes[0]
        logging.info(f"self.genotype: {self.genotype}")

        self.long_analyzers: list[core.LongRecordingAnalyzer] = []
        logging.debug(f"Creating {len(self._bin_folders)} LongRecordings")
        self.long_recordings = [core.LongRecordingOrganizer(e, **lro_kwargs) for e in self._bin_folders]

        channel_names = [x.channel_names for x in self.long_recordings]
        if len(set([" ".join(x) for x in channel_names])) > 1:
            warnings.warn(f"Inconsistent channel names in long_recordings: {channel_names}")
        self.channel_names = channel_names[0]
        self.bad_channels_dict = {}

        animal_ids = [x["animal"] for x in animalday_dicts]
        if len(set(animal_ids)) > 1:
            warnings.warn(f"Inconsistent animal IDs in {animal_ids}")
        self.animal_id = animal_ids[0]

        self.features_df: pd.DataFrame = pd.DataFrame()
        self.features_avg_df: pd.DataFrame = pd.DataFrame()

    def convert_colbins_to_rowbins(self, overwrite=False, multiprocess_mode: Literal["dask", "serial"] = "serial"):
        for lrec in tqdm(self.long_recordings, desc="Converting column bins to row bins"):
            lrec.convert_colbins_to_rowbins(overwrite=overwrite, multiprocess_mode=multiprocess_mode)

    def convert_rowbins_to_rec(self, multiprocess_mode: Literal["dask", "serial"] = "serial"):
        for lrec in tqdm(self.long_recordings, desc="Converting row bins to recs"):
            lrec.convert_rowbins_to_rec(multiprocess_mode=multiprocess_mode)

    def cleanup_rec(self):
        for lrec in self.long_recordings:
            lrec.cleanup_rec()

    def compute_bad_channels(self, lof_threshold: float = 1.5):
        for lrec in self.long_recordings:
            lrec.compute_bad_channels(lof_threshold=lof_threshold)
        self.bad_channels_dict = {
            animalday: lrec.bad_channel_names for animalday, lrec in zip(self.animaldays, self.long_recordings)
        }

    def compute_windowed_analysis(
        self,
        features: list[str],
        exclude: list[str] = ["nspike", "lognspike"],
        window_s=4,
        multiprocess_mode: Literal["dask", "serial"] = "serial",
        suppress_short_interval_error=False,
        **kwargs,
    ):
        """Computes windowed analysis of animal recordings. The data is divided into windows (time bins), then features are extracted from each window. The result is
        formatted to a Dataframe and wrapped into a WindowAnalysisResult object.

        Args:
            features (list[str]): List of features to compute. See individual compute_...() functions for output format
            exclude (list[str], optional): List of features to ignore. Will override the features parameter. Defaults to [].
            window_s (int, optional): Length of each window in seconds. Note that some features break with very short window times. Defaults to 4.
            suppress_short_interval_error (bool, optional): If True, suppress ValueError for short intervals between timestamps in resulting WindowAnalysisResult. Useful for aggregated WARs. Defaults to False.

        Raises:
            AttributeError: If a feature's compute_...() function was not implemented, this error will be raised.

        Returns:
            window_analysis_result: a WindowAnalysisResult object
        """
        features = _sanitize_feature_request(features, exclude)

        dataframes = []
        for lrec in self.long_recordings:  # Iterate over all long recordings
            logging.info(f"Computing windowed analysis for {lrec.base_folder_path}")
            lan = core.LongRecordingAnalyzer(lrec, fragment_len_s=window_s)
            if lan.n_fragments == 0:
                logging.warning(f"No fragments found for {lrec.base_folder_path}. Skipping.")
                continue

            logging.debug(f"Processing {lan.n_fragments} fragments")
            miniters = int(lan.n_fragments / 100)
            match multiprocess_mode:
                case "dask":
                    # The last fragment is not included because it makes the dask array ragged
                    logging.debug("Converting LongRecording to numpy array")

                    n_fragments_war = max(lan.n_fragments - 1, 1)
                    first_fragment = lan.get_fragment_np(0)
                    np_fragments = np.empty((n_fragments_war,) + first_fragment.shape, dtype=first_fragment.dtype)
                    logging.debug(f"np_fragments.shape: {np_fragments.shape}")
                    for idx in range(n_fragments_war):
                        np_fragments[idx] = lan.get_fragment_np(idx)

                    # Cache fragments to zarr
                    tmppath, _ = core.cache_fragments_to_zarr(np_fragments, n_fragments_war)
                    del np_fragments

                    logging.debug("Processing metadata serially")
                    metadatas = [self._process_fragment_metadata(idx, lan, window_s) for idx in range(n_fragments_war)]
                    meta_df = pd.DataFrame(metadatas)

                    logging.debug("Processing features in parallel")
                    np_fragments_reconstruct = da.from_zarr(tmppath, chunks=("auto", -1, -1))
                    logging.debug(f"Dask array shape: {np_fragments_reconstruct.shape}")
                    logging.debug(f"Dask array chunks: {np_fragments_reconstruct.chunks}")

                    # Create delayed tasks for each fragment using efficient dependency resolution
                    feature_values = [
                        delayed(FragmentAnalyzer.process_fragment_with_dependencies)(
                            np_fragments_reconstruct[idx], lan.f_s, features, kwargs
                        )
                        for idx in range(n_fragments_war)
                    ]

                    # Compute features in parallel
                    feature_values = dask.compute(*feature_values)

                    # Clean up temp directory after processing
                    logging.debug("Cleaning up temp directory")
                    try:
                        import shutil

                        shutil.rmtree(tmppath)
                    except (OSError, FileNotFoundError) as e:
                        logging.warning(f"Failed to remove temporary directory {tmppath}: {e}")

                    logging.debug("Combining metadata and feature values")
                    feat_df = pd.DataFrame(feature_values)
                    lan_df = pd.concat([meta_df, feat_df], axis=1)

                case _:
                    logging.debug("Processing serially")
                    lan_df = []
                    for idx in tqdm(range(lan.n_fragments), desc="Processing rows", miniters=miniters):
                        lan_df.append(self._process_fragment_serial(idx, features, lan, window_s, kwargs))

            lan_df = pd.DataFrame(lan_df)

            logging.debug("Validating timestamps")
            core.validate_timestamps(lan_df["timestamp"].tolist())
            lan_df = lan_df.sort_values("timestamp").reset_index(drop=True)

            self.long_analyzers.append(lan)
            dataframes.append(lan_df)

        self.features_df = pd.concat(dataframes)
        self.features_df = self.features_df

        self.window_analysis_result = WindowAnalysisResult(
            self.features_df,
            self.animal_id,
            self.genotype,
            self.channel_names,
            self.assume_from_number,
            self.bad_channels_dict,
            suppress_short_interval_error,
        )

        return self.window_analysis_result

    def compute_spike_analysis(self, multiprocess_mode: Literal["dask", "serial"] = "serial"):
        """Compute spike sorting on all long recordings and return a list of SpikeAnalysisResult objects

        Args:
            multiprocess_mode (Literal['dask', 'serial']): Whether to use Dask for parallel processing. Defaults to 'serial'.

        Returns:
            spike_analysis_results: list[SpikeAnalysisResult]. Each SpikeAnalysisResult object corresponds to a LongRecording object,
            typically a different day or recording session.

        Raises:
            ImportError: If mountainsort5 is not available.
        """
        # Check if mountainsort5 is available
        from ..core.analyze_sort import MOUNTAINSORT_AVAILABLE
        if not MOUNTAINSORT_AVAILABLE:
            raise ImportError(
                "Spike analysis requires mountainsort5. Install it with: pip install mountainsort5"
            )
        sars = []
        lrec_sorts = []
        lrec_recs = []
        recs = [lrec.LongRecording for lrec in self.long_recordings]
        logging.info(f"Sorting {len(recs)} recordings")
        for rec in recs:
            if rec.get_total_samples() == 0:
                logging.warning(f"Skipping {rec.__str__()} because it has no samples")
                sortings, recordings = [], []
            else:
                sortings, recordings = core.MountainSortAnalyzer.sort_recording(
                    rec, multiprocess_mode=multiprocess_mode
                )
            lrec_sorts.append(sortings)
            lrec_recs.append(recordings)

        if multiprocess_mode == "dask":
            lrec_sorts = dask.compute(*lrec_sorts)

        lrec_sas = [
            [
                si.create_sorting_analyzer(sorting, recording, sparse=False)
                for sorting, recording in zip(sortings, recordings)
            ]
            for sortings, recordings in zip(lrec_sorts, lrec_recs)
        ]
        sars = [
            SpikeAnalysisResult(
                result_sas=sas,
                result_mne=None,
                animal_id=self.animal_id,
                genotype=self.genotype,
                animal_day=self.animaldays[i],
                bin_folder_name=self.bin_folder_names[i],
                metadata=self.long_recordings[i].meta,
                channel_names=self.channel_names,
                assume_from_number=self.assume_from_number,
            )
            for i, sas in enumerate(lrec_sas)
        ]

        self.spike_analysis_results = sars
        return self.spike_analysis_results

    def _process_fragment_serial(self, idx, features, lan: core.LongRecordingAnalyzer, window_s, kwargs: dict):
        row = self._process_fragment_metadata(idx, lan, window_s)
        row.update(self._process_fragment_features(idx, features, lan, kwargs))
        return row

    def _process_fragment_metadata(self, idx, lan: core.LongRecordingAnalyzer, window_s):
        row = {}

        lan_folder = lan.LongRecording.base_folder_path
        animalday_dict = core.parse_path_to_animalday(
            lan_folder, animal_param=self.animal_param, day_sep=self.day_sep, mode=self.read_mode
        )
        row["animalday"] = animalday_dict["animalday"]
        row["animal"] = animalday_dict["animal"]
        row["day"] = animalday_dict["day"]
        row["genotype"] = animalday_dict["genotype"]
        row["duration"] = lan.LongRecording.get_dur_fragment(window_s, idx)
        row["endfile"] = lan.get_file_end(idx)

        frag_dt = lan.LongRecording.get_datetime_fragment(window_s, idx)
        row["timestamp"] = frag_dt
        row["isday"] = core.is_day(frag_dt)

        return row

    def _process_fragment_features(self, idx, features, lan: core.LongRecordingAnalyzer, kwargs: dict):
        row = {}
        for feat in features:
            func = getattr(lan, f"compute_{feat}")
            if callable(func):
                row[feat] = func(idx, **kwargs)
            else:
                raise AttributeError(f"Invalid function {func}")
        return row

__init__(base_folder_path, anim_id, day_sep=None, mode='concat', assume_from_number=False, skip_days=[], truncate=False, lro_kwargs={})

AnimalOrganizer is used to organize data from a single animal into a format that can be used for analysis. It is used to organize data from a single animal into a format that can be used for analysis.

Parameters:

Name	Type	Description	Default
base_folder_path	str	The path to the base folder of the animal data.	required
anim_id	str	The ID of the animal. This should correspond to only one animal.	required
day_sep	str	Separator for day in folder name. Set to None or empty string to get all folders. Defaults to None.	None
mode	Literal['nest', 'concat', 'base', 'noday']	The mode of the AnimalOrganizer. Defaults to "concat". "nest": base_folder_path / animal_id / date_format "concat": base_folder_path / animal_id*date_format "base": base_folder_path "noday": base_folder_path / animal_id	'concat'
assume_from_number	bool	Whether to assume the animal ID is a number. Defaults to False.	False
skip_days	list[str]	The days to skip. Defaults to [].	[]
truncate	bool | int	Whether to truncate the data. Defaults to False.	False
lro_kwargs	dict	Keyword arguments for LongRecordingOrganizer. Defaults to {}.	{}

Source code in pythoneeg/visualization/results.py

def __init__(
    self,
    base_folder_path,
    anim_id: str,
    day_sep: str | None = None,
    mode: Literal["nest", "concat", "base", "noday"] = "concat",
    assume_from_number=False,
    skip_days: list[str] = [],
    truncate: bool | int = False,
    lro_kwargs: dict = {},
) -> None:
    """
    AnimalOrganizer is used to organize data from a single animal into a format that can be used for analysis.
    It is used to organize data from a single animal into a format that can be used for analysis.

    Args:
        base_folder_path (str): The path to the base folder of the animal data.
        anim_id (str): The ID of the animal. This should correspond to only one animal.
        day_sep (str, optional): Separator for day in folder name. Set to None or empty string to get all folders. Defaults to None.
        mode (Literal["nest", "concat", "base", "noday"], optional): The mode of the AnimalOrganizer. Defaults to "concat".
            "nest": base_folder_path / animal_id / date_format
            "concat": base_folder_path / animal_id*date_format
            "base": base_folder_path
            "noday": base_folder_path / animal_id
        assume_from_number (bool, optional): Whether to assume the animal ID is a number. Defaults to False.
        skip_days (list[str], optional): The days to skip. Defaults to [].
        truncate (bool|int, optional): Whether to truncate the data. Defaults to False.
        lro_kwargs (dict, optional): Keyword arguments for LongRecordingOrganizer. Defaults to {}.
    """

    self.base_folder_path = Path(base_folder_path)
    self.anim_id = anim_id
    self.animal_param = [anim_id]
    self.day_sep = day_sep
    self.read_mode = mode
    self.assume_from_number = assume_from_number

    match mode:
        case "nest":
            self.bin_folder_pattern = self.base_folder_path / f"*{self.anim_id}*" / "*"
        case "concat" | "noday":
            self.bin_folder_pattern = self.base_folder_path / f"*{self.anim_id}*"
            # self.bin_folder_pat = self.base_folder_path / f"*{self.anim_id}*{self.date_format}*"
        case "base":
            self.bin_folder_pattern = self.base_folder_path
        # case 'noday':
        #     self.bin_folder_pat = self.base_folder_path / f"*{self.anim_id}*"
        case _:
            raise ValueError(f"Invalid mode: {mode}")

    self._bin_folders = glob.glob(str(self.bin_folder_pattern))
    # if mode != 'noday':
    #     self.__bin_folders = [x for x in self.__bin_folders if datetime.strptime(Path(x).name, self.date_format)]
    truncate = core.utils.parse_truncate(truncate)
    if truncate:
        warnings.warn(f"AnimalOrganizer will be truncated to the first {truncate} LongRecordings")
        self._bin_folders = self._bin_folders[:truncate]
    self._bin_folders = [x for x in self._bin_folders if not any(y in x for y in skip_days)]
    self.bin_folder_names = [Path(x).name for x in self._bin_folders]
    logging.info(f"bin_folder_pattern: {self.bin_folder_pattern}")
    logging.info(f"self._bin_folders: {self._bin_folders}")
    logging.info(f"self.bin_folder_names: {self.bin_folder_names}")

    if mode == "noday" and len(self._bin_folders) > 1:
        raise ValueError(f"Animal ID '{self.anim_id}' is not unique, found: {', '.join(self._bin_folders)}")
    elif len(self._bin_folders) == 0:
        raise ValueError(f"No files found for animal ID {self.anim_id}")

    animalday_dicts = [
        core.parse_path_to_animalday(e, animal_param=self.animal_param, day_sep=self.day_sep, mode=self.read_mode)
        for e in self._bin_folders
    ]
    self.animaldays = [x["animalday"] for x in animalday_dicts]
    logging.info(f"self.animaldays: {self.animaldays}")

    genotypes = [x["genotype"] for x in animalday_dicts]
    if len(set(genotypes)) > 1:
        warnings.warn(f"Inconsistent genotypes in {genotypes}")
    self.genotype = genotypes[0]
    logging.info(f"self.genotype: {self.genotype}")

    self.long_analyzers: list[core.LongRecordingAnalyzer] = []
    logging.debug(f"Creating {len(self._bin_folders)} LongRecordings")
    self.long_recordings = [core.LongRecordingOrganizer(e, **lro_kwargs) for e in self._bin_folders]

    channel_names = [x.channel_names for x in self.long_recordings]
    if len(set([" ".join(x) for x in channel_names])) > 1:
        warnings.warn(f"Inconsistent channel names in long_recordings: {channel_names}")
    self.channel_names = channel_names[0]
    self.bad_channels_dict = {}

    animal_ids = [x["animal"] for x in animalday_dicts]
    if len(set(animal_ids)) > 1:
        warnings.warn(f"Inconsistent animal IDs in {animal_ids}")
    self.animal_id = animal_ids[0]

    self.features_df: pd.DataFrame = pd.DataFrame()
    self.features_avg_df: pd.DataFrame = pd.DataFrame()

compute_spike_analysis(multiprocess_mode='serial')

Compute spike sorting on all long recordings and return a list of SpikeAnalysisResult objects

Parameters:

Name	Type	Description	Default
multiprocess_mode	Literal['dask', 'serial']	Whether to use Dask for parallel processing. Defaults to 'serial'.	'serial'

Returns:

Name	Type	Description
spike_analysis_results		list[SpikeAnalysisResult]. Each SpikeAnalysisResult object corresponds to a LongRecording object,
		typically a different day or recording session.

Raises:

Type	Description
ImportError	If mountainsort5 is not available.

Source code in pythoneeg/visualization/results.py

def compute_spike_analysis(self, multiprocess_mode: Literal["dask", "serial"] = "serial"):
    """Compute spike sorting on all long recordings and return a list of SpikeAnalysisResult objects

    Args:
        multiprocess_mode (Literal['dask', 'serial']): Whether to use Dask for parallel processing. Defaults to 'serial'.

    Returns:
        spike_analysis_results: list[SpikeAnalysisResult]. Each SpikeAnalysisResult object corresponds to a LongRecording object,
        typically a different day or recording session.

    Raises:
        ImportError: If mountainsort5 is not available.
    """
    # Check if mountainsort5 is available
    from ..core.analyze_sort import MOUNTAINSORT_AVAILABLE
    if not MOUNTAINSORT_AVAILABLE:
        raise ImportError(
            "Spike analysis requires mountainsort5. Install it with: pip install mountainsort5"
        )
    sars = []
    lrec_sorts = []
    lrec_recs = []
    recs = [lrec.LongRecording for lrec in self.long_recordings]
    logging.info(f"Sorting {len(recs)} recordings")
    for rec in recs:
        if rec.get_total_samples() == 0:
            logging.warning(f"Skipping {rec.__str__()} because it has no samples")
            sortings, recordings = [], []
        else:
            sortings, recordings = core.MountainSortAnalyzer.sort_recording(
                rec, multiprocess_mode=multiprocess_mode
            )
        lrec_sorts.append(sortings)
        lrec_recs.append(recordings)

    if multiprocess_mode == "dask":
        lrec_sorts = dask.compute(*lrec_sorts)

    lrec_sas = [
        [
            si.create_sorting_analyzer(sorting, recording, sparse=False)
            for sorting, recording in zip(sortings, recordings)
        ]
        for sortings, recordings in zip(lrec_sorts, lrec_recs)
    ]
    sars = [
        SpikeAnalysisResult(
            result_sas=sas,
            result_mne=None,
            animal_id=self.animal_id,
            genotype=self.genotype,
            animal_day=self.animaldays[i],
            bin_folder_name=self.bin_folder_names[i],
            metadata=self.long_recordings[i].meta,
            channel_names=self.channel_names,
            assume_from_number=self.assume_from_number,
        )
        for i, sas in enumerate(lrec_sas)
    ]

    self.spike_analysis_results = sars
    return self.spike_analysis_results

compute_windowed_analysis(features, exclude=['nspike', 'lognspike'], window_s=4, multiprocess_mode='serial', suppress_short_interval_error=False, **kwargs)

Computes windowed analysis of animal recordings. The data is divided into windows (time bins), then features are extracted from each window. The result is formatted to a Dataframe and wrapped into a WindowAnalysisResult object.

Parameters:

Name	Type	Description	Default
features	list[str]	List of features to compute. See individual compute_...() functions for output format	required
exclude	list[str]	List of features to ignore. Will override the features parameter. Defaults to [].	['nspike', 'lognspike']
window_s	int	Length of each window in seconds. Note that some features break with very short window times. Defaults to 4.	4
suppress_short_interval_error	bool	If True, suppress ValueError for short intervals between timestamps in resulting WindowAnalysisResult. Useful for aggregated WARs. Defaults to False.	False

Raises:

Type	Description
AttributeError	If a feature's compute_...() function was not implemented, this error will be raised.

Returns:

Name	Type	Description
window_analysis_result		a WindowAnalysisResult object

Source code in pythoneeg/visualization/results.py

def compute_windowed_analysis(
    self,
    features: list[str],
    exclude: list[str] = ["nspike", "lognspike"],
    window_s=4,
    multiprocess_mode: Literal["dask", "serial"] = "serial",
    suppress_short_interval_error=False,
    **kwargs,
):
    """Computes windowed analysis of animal recordings. The data is divided into windows (time bins), then features are extracted from each window. The result is
    formatted to a Dataframe and wrapped into a WindowAnalysisResult object.

    Args:
        features (list[str]): List of features to compute. See individual compute_...() functions for output format
        exclude (list[str], optional): List of features to ignore. Will override the features parameter. Defaults to [].
        window_s (int, optional): Length of each window in seconds. Note that some features break with very short window times. Defaults to 4.
        suppress_short_interval_error (bool, optional): If True, suppress ValueError for short intervals between timestamps in resulting WindowAnalysisResult. Useful for aggregated WARs. Defaults to False.

    Raises:
        AttributeError: If a feature's compute_...() function was not implemented, this error will be raised.

    Returns:
        window_analysis_result: a WindowAnalysisResult object
    """
    features = _sanitize_feature_request(features, exclude)

    dataframes = []
    for lrec in self.long_recordings:  # Iterate over all long recordings
        logging.info(f"Computing windowed analysis for {lrec.base_folder_path}")
        lan = core.LongRecordingAnalyzer(lrec, fragment_len_s=window_s)
        if lan.n_fragments == 0:
            logging.warning(f"No fragments found for {lrec.base_folder_path}. Skipping.")
            continue

        logging.debug(f"Processing {lan.n_fragments} fragments")
        miniters = int(lan.n_fragments / 100)
        match multiprocess_mode:
            case "dask":
                # The last fragment is not included because it makes the dask array ragged
                logging.debug("Converting LongRecording to numpy array")

                n_fragments_war = max(lan.n_fragments - 1, 1)
                first_fragment = lan.get_fragment_np(0)
                np_fragments = np.empty((n_fragments_war,) + first_fragment.shape, dtype=first_fragment.dtype)
                logging.debug(f"np_fragments.shape: {np_fragments.shape}")
                for idx in range(n_fragments_war):
                    np_fragments[idx] = lan.get_fragment_np(idx)

                # Cache fragments to zarr
                tmppath, _ = core.cache_fragments_to_zarr(np_fragments, n_fragments_war)
                del np_fragments

                logging.debug("Processing metadata serially")
                metadatas = [self._process_fragment_metadata(idx, lan, window_s) for idx in range(n_fragments_war)]
                meta_df = pd.DataFrame(metadatas)

                logging.debug("Processing features in parallel")
                np_fragments_reconstruct = da.from_zarr(tmppath, chunks=("auto", -1, -1))
                logging.debug(f"Dask array shape: {np_fragments_reconstruct.shape}")
                logging.debug(f"Dask array chunks: {np_fragments_reconstruct.chunks}")

                # Create delayed tasks for each fragment using efficient dependency resolution
                feature_values = [
                    delayed(FragmentAnalyzer.process_fragment_with_dependencies)(
                        np_fragments_reconstruct[idx], lan.f_s, features, kwargs
                    )
                    for idx in range(n_fragments_war)
                ]

                # Compute features in parallel
                feature_values = dask.compute(*feature_values)

                # Clean up temp directory after processing
                logging.debug("Cleaning up temp directory")
                try:
                    import shutil

                    shutil.rmtree(tmppath)
                except (OSError, FileNotFoundError) as e:
                    logging.warning(f"Failed to remove temporary directory {tmppath}: {e}")

                logging.debug("Combining metadata and feature values")
                feat_df = pd.DataFrame(feature_values)
                lan_df = pd.concat([meta_df, feat_df], axis=1)

            case _:
                logging.debug("Processing serially")
                lan_df = []
                for idx in tqdm(range(lan.n_fragments), desc="Processing rows", miniters=miniters):
                    lan_df.append(self._process_fragment_serial(idx, features, lan, window_s, kwargs))

        lan_df = pd.DataFrame(lan_df)

        logging.debug("Validating timestamps")
        core.validate_timestamps(lan_df["timestamp"].tolist())
        lan_df = lan_df.sort_values("timestamp").reset_index(drop=True)

        self.long_analyzers.append(lan)
        dataframes.append(lan_df)

    self.features_df = pd.concat(dataframes)
    self.features_df = self.features_df

    self.window_analysis_result = WindowAnalysisResult(
        self.features_df,
        self.animal_id,
        self.genotype,
        self.channel_names,
        self.assume_from_number,
        self.bad_channels_dict,
        suppress_short_interval_error,
    )

    return self.window_analysis_result

SpikeAnalysisResult

Bases: AnimalFeatureParser

Source code in pythoneeg/visualization/results.py

class SpikeAnalysisResult(AnimalFeatureParser):
    def __init__(
        self,
        result_sas: list[si.SortingAnalyzer],
        result_mne: mne.io.RawArray = None,
        animal_id: str = None,
        genotype: str = None,
        animal_day: str = None,
        bin_folder_name: str = None,
        metadata: core.DDFBinaryMetadata = None,
        channel_names: list[str] = None,
        assume_from_number=False,
    ) -> None:
        """
        Args:
            result (list[si.SortingAnalyzer]): Result comes from AnimalOrganizer.compute_spike_analysis(). Each SortingAnalyzer is a single channel.
            animal_id (str, optional): Identifier for the animal where result was computed from. Defaults to None.
            genotype (str, optional): Genotype of animal. Defaults to None.
            channel_names (list[str], optional): List of channel names. Defaults to None.
            assume_channels (bool, optional): If true, assumes channel names according to AnimalFeatureParser.DEFAULT_CHNUM_TO_NAME. Defaults to False.
        """
        self.result_sas = result_sas
        self.result_mne = result_mne
        if (result_mne is None) == (result_sas is None):
            raise ValueError("Exactly one of result_sas or result_mne must be provided")
        self.animal_id = animal_id
        self.genotype = genotype
        self.animal_day = animal_day
        self.bin_folder_name = bin_folder_name
        self.metadata = metadata
        self.channel_names = channel_names
        self.assume_from_number = assume_from_number
        self.channel_abbrevs = [
            core.parse_chname_to_abbrev(x, assume_from_number=assume_from_number) for x in self.channel_names
        ]

        logging.info(f"Channel names: \t{self.channel_names}")
        logging.info(f"Channel abbreviations: \t{self.channel_abbrevs}")

    def convert_to_mne(self, chunk_len: float = 60, save_raw=True) -> mne.io.RawArray:
        if self.result_mne is None:
            result_mne = SpikeAnalysisResult.convert_sas_to_mne(self.result_sas, chunk_len)
            if save_raw:
                self.result_mne = result_mne
            else:
                return result_mne
        return self.result_mne

    def save_fif_and_json(
        self,
        folder: str | Path,
        convert_to_mne=True,
        make_folder=True,
        slugify_filebase=True,
        save_abbrevs_as_chnames=False,
        overwrite=False,
    ):
        """Archive spike analysis result into the folder specified, as a fif and json file.

        Args:
            folder (str | Path): Destination folder to save results to
            convert_to_mne (bool, optional): If True, convert the SortingAnalyzers to a MNE RawArray if self.result_mne is None. Defaults to True.
            make_folder (bool, optional): If True, create the folder if it doesn't exist. Defaults to True.
            slugify_filebase (bool, optional): If True, slugify the filebase (replace special characters). Defaults to True.
            save_abbrevs_as_chnames (bool, optional): If True, save the channel abbreviations as the channel names in the json file. Defaults to False.
            overwrite (bool, optional): If True, overwrite the existing files. Defaults to False.
        """
        if self.result_mne is None:
            if convert_to_mne:
                result_mne = self.convert_to_mne(save_raw=True)
                if result_mne is None:
                    warnings.warn("No SortingAnalyzers found, skipping saving")
                    return
            else:
                raise ValueError("No MNE RawArray found, and convert_to_mne is False. Run convert_to_mne() first.")
        else:
            result_mne = self.result_mne

        folder = Path(folder)
        if make_folder:
            folder.mkdir(parents=True, exist_ok=True)

        if slugify_filebase:
            filebase = folder / slugify(f"{self.animal_id}-{self.genotype}-{self.animal_day}")
        else:
            filebase = folder / f"{self.animal_id}-{self.genotype}-{self.animal_day}"
        filebase = str(filebase)

        if not overwrite:
            if filebase + ".json" in folder.glob("*.json"):
                raise FileExistsError(f"File {filebase}.json already exists")
            if filebase + ".fif" in folder.glob("*.fif"):
                raise FileExistsError(f"File {filebase}.fif already exists")
        else:
            for f in folder.glob("*"):
                f.unlink()
        result_mne.save(filebase + "-raw.fif", overwrite=overwrite)
        del result_mne

        json_dict = {
            "animal_id": self.animal_id,
            "genotype": self.genotype,
            "animal_day": self.animal_day,
            "bin_folder_name": self.bin_folder_name,
            "metadata": self.metadata.metadata_path,
            "channel_names": self.channel_abbrevs if save_abbrevs_as_chnames else self.channel_names,
            "assume_from_number": False if save_abbrevs_as_chnames else self.assume_from_number,
        }
        with open(filebase + ".json", "w") as f:
            json.dump(json_dict, f, indent=2)

    @classmethod
    def load_fif_and_json(cls, folder: str | Path):
        folder = Path(folder)
        if not folder.exists():
            raise ValueError(f"Folder {folder} does not exist")

        fif_files = list(folder.glob("*.fif"))  # there may be more than 1 fif file
        json_files = list(folder.glob("*.json"))

        if len(json_files) != 1:
            raise ValueError(f"Expected exactly one json file in {folder}")

        fif_path = fif_files[0]
        json_path = json_files[0]

        with open(json_path, "r") as f:
            data = json.load(f)
        # data['metadata'] = core.DDFBinaryMetadata(data['metadata'])
        data["result_mne"] = mne.io.read_raw_fif(fif_path)
        data["result_sas"] = None
        return cls(**data)

    @staticmethod
    def convert_sas_to_mne(sas: list[si.SortingAnalyzer], chunk_len: float = 60) -> mne.io.RawArray:
        """Convert a list of SortingAnalyzers to a MNE RawArray.

        Args:
            sas (list[si.SortingAnalyzer]): The list of SortingAnalyzers to convert
            chunk_len (float, optional): The length of the chunks to use for the conversion. Defaults to 60.

        Returns:
            mne.io.RawArray: The converted RawArray, with spikes labeled as annotations
        """
        if len(sas) == 0:
            return None

        # Check that all SortingAnalyzers have the same sampling frequency
        sfreqs = [sa.recording.get_sampling_frequency() for sa in sas]
        if not all(sf == sfreqs[0] for sf in sfreqs):
            raise ValueError(f"All SortingAnalyzers must have the same sampling frequency. Got frequencies: {sfreqs}")

        # Preallocate data array
        total_frames = int(sas[0].recording.get_duration() * sfreqs[0])
        n_channels = len(sas)
        data = np.empty((n_channels, total_frames))
        logging.debug(f"Data shape: {data.shape}")

        # Fill data array one channel at a time
        for i, sa in enumerate(sas):
            logging.debug(f"Converting channel {i + 1} of {n_channels}")
            data[i, :] = SpikeAnalysisResult.convert_sa_to_np(sa, chunk_len)

        channel_names = [str(sa.recording.get_channel_ids().item()) for sa in sas]
        logging.debug(f"Channel names: {channel_names}")
        sfreq = sfreqs[0]

        # Extract spike times for each unit and create annotations
        onset = []
        description = []
        for sa in sas:
            for unit_id in sa.sorting.get_unit_ids():
                spike_train = sa.sorting.get_unit_spike_train(unit_id)
                # Convert to seconds and filter to recording duration
                spike_times = spike_train / sa.sorting.get_sampling_frequency()
                mask = spike_times < sa.recording.get_duration()
                spike_times = spike_times[mask]

                # Create annotation for each spike
                onset.extend(spike_times)
                description.extend(
                    [sa.recording.get_channel_ids().item()] * len(spike_times)
                )  # collapse all units into 1 spike train
        annotations = mne.Annotations(onset, duration=0, description=description)

        info = mne.create_info(ch_names=channel_names, sfreq=sfreq, ch_types="eeg")
        raw = mne.io.RawArray(data=data, info=info)
        raw = raw.set_annotations(annotations)
        return raw

    @staticmethod
    def convert_sa_to_np(sa: si.SortingAnalyzer, chunk_len: float = 60) -> np.ndarray:
        """Convert a SortingAnalyzer to an MNE RawArray.

        Args:
            sa (si.SortingAnalyzer): The SortingAnalyzer to convert. Must have only 1 channel.
            chunk_len (float, optional): The length of the chunks to use for the conversion. Defaults to 60.
        Returns:
            np.ndarray: The converted traces
        """
        # Check that SortingAnalyzer only has 1 channel
        if len(sa.recording.get_channel_ids()) != 1:
            raise ValueError(
                f"Expected SortingAnalyzer to have 1 channel, but got {len(sa.recording.get_channel_ids())} channels"
            )

        rec = sa.recording
        logging.debug(f"Recording info: {rec}")

        # Calculate total number of frames and chunks
        total_frames = int(rec.get_duration() * rec.get_sampling_frequency())
        frames_per_chunk = round(chunk_len * rec.get_sampling_frequency())
        n_chunks = total_frames // frames_per_chunk

        traces = np.empty(total_frames)

        for j in range(n_chunks):
            start_frame = j * frames_per_chunk
            if j == n_chunks - 1:
                end_frame = total_frames
            else:
                end_frame = (j + 1) * frames_per_chunk
            traces[start_frame:end_frame] = rec.get_traces(
                start_frame=start_frame, end_frame=end_frame, return_scaled=True
            ).flatten()
        traces *= 1e-6  # convert from uV to V
        return traces

__init__(result_sas, result_mne=None, animal_id=None, genotype=None, animal_day=None, bin_folder_name=None, metadata=None, channel_names=None, assume_from_number=False)

Parameters:

Name	Type	Description	Default
result	list[SortingAnalyzer]	Result comes from AnimalOrganizer.compute_spike_analysis(). Each SortingAnalyzer is a single channel.	required
animal_id	str	Identifier for the animal where result was computed from. Defaults to None.	None
genotype	str	Genotype of animal. Defaults to None.	None
channel_names	list[str]	List of channel names. Defaults to None.	None
assume_channels	bool	If true, assumes channel names according to AnimalFeatureParser.DEFAULT_CHNUM_TO_NAME. Defaults to False.	required

Source code in pythoneeg/visualization/results.py

def __init__(
    self,
    result_sas: list[si.SortingAnalyzer],
    result_mne: mne.io.RawArray = None,
    animal_id: str = None,
    genotype: str = None,
    animal_day: str = None,
    bin_folder_name: str = None,
    metadata: core.DDFBinaryMetadata = None,
    channel_names: list[str] = None,
    assume_from_number=False,
) -> None:
    """
    Args:
        result (list[si.SortingAnalyzer]): Result comes from AnimalOrganizer.compute_spike_analysis(). Each SortingAnalyzer is a single channel.
        animal_id (str, optional): Identifier for the animal where result was computed from. Defaults to None.
        genotype (str, optional): Genotype of animal. Defaults to None.
        channel_names (list[str], optional): List of channel names. Defaults to None.
        assume_channels (bool, optional): If true, assumes channel names according to AnimalFeatureParser.DEFAULT_CHNUM_TO_NAME. Defaults to False.
    """
    self.result_sas = result_sas
    self.result_mne = result_mne
    if (result_mne is None) == (result_sas is None):
        raise ValueError("Exactly one of result_sas or result_mne must be provided")
    self.animal_id = animal_id
    self.genotype = genotype
    self.animal_day = animal_day
    self.bin_folder_name = bin_folder_name
    self.metadata = metadata
    self.channel_names = channel_names
    self.assume_from_number = assume_from_number
    self.channel_abbrevs = [
        core.parse_chname_to_abbrev(x, assume_from_number=assume_from_number) for x in self.channel_names
    ]

    logging.info(f"Channel names: \t{self.channel_names}")
    logging.info(f"Channel abbreviations: \t{self.channel_abbrevs}")

convert_sa_to_np(sa, chunk_len=60) staticmethod

Convert a SortingAnalyzer to an MNE RawArray.

Parameters:

Name	Type	Description	Default
sa	SortingAnalyzer	The SortingAnalyzer to convert. Must have only 1 channel.	required
chunk_len	float	The length of the chunks to use for the conversion. Defaults to 60.	60

Returns: np.ndarray: The converted traces

Source code in pythoneeg/visualization/results.py

@staticmethod
def convert_sa_to_np(sa: si.SortingAnalyzer, chunk_len: float = 60) -> np.ndarray:
    """Convert a SortingAnalyzer to an MNE RawArray.

    Args:
        sa (si.SortingAnalyzer): The SortingAnalyzer to convert. Must have only 1 channel.
        chunk_len (float, optional): The length of the chunks to use for the conversion. Defaults to 60.
    Returns:
        np.ndarray: The converted traces
    """
    # Check that SortingAnalyzer only has 1 channel
    if len(sa.recording.get_channel_ids()) != 1:
        raise ValueError(
            f"Expected SortingAnalyzer to have 1 channel, but got {len(sa.recording.get_channel_ids())} channels"
        )

    rec = sa.recording
    logging.debug(f"Recording info: {rec}")

    # Calculate total number of frames and chunks
    total_frames = int(rec.get_duration() * rec.get_sampling_frequency())
    frames_per_chunk = round(chunk_len * rec.get_sampling_frequency())
    n_chunks = total_frames // frames_per_chunk

    traces = np.empty(total_frames)

    for j in range(n_chunks):
        start_frame = j * frames_per_chunk
        if j == n_chunks - 1:
            end_frame = total_frames
        else:
            end_frame = (j + 1) * frames_per_chunk
        traces[start_frame:end_frame] = rec.get_traces(
            start_frame=start_frame, end_frame=end_frame, return_scaled=True
        ).flatten()
    traces *= 1e-6  # convert from uV to V
    return traces

convert_sas_to_mne(sas, chunk_len=60) staticmethod

Convert a list of SortingAnalyzers to a MNE RawArray.

Parameters:

Name	Type	Description	Default
sas	list[SortingAnalyzer]	The list of SortingAnalyzers to convert	required
chunk_len	float	The length of the chunks to use for the conversion. Defaults to 60.	60

Returns:

Type	Description
RawArray	mne.io.RawArray: The converted RawArray, with spikes labeled as annotations

Source code in pythoneeg/visualization/results.py

@staticmethod
def convert_sas_to_mne(sas: list[si.SortingAnalyzer], chunk_len: float = 60) -> mne.io.RawArray:
    """Convert a list of SortingAnalyzers to a MNE RawArray.

    Args:
        sas (list[si.SortingAnalyzer]): The list of SortingAnalyzers to convert
        chunk_len (float, optional): The length of the chunks to use for the conversion. Defaults to 60.

    Returns:
        mne.io.RawArray: The converted RawArray, with spikes labeled as annotations
    """
    if len(sas) == 0:
        return None

    # Check that all SortingAnalyzers have the same sampling frequency
    sfreqs = [sa.recording.get_sampling_frequency() for sa in sas]
    if not all(sf == sfreqs[0] for sf in sfreqs):
        raise ValueError(f"All SortingAnalyzers must have the same sampling frequency. Got frequencies: {sfreqs}")

    # Preallocate data array
    total_frames = int(sas[0].recording.get_duration() * sfreqs[0])
    n_channels = len(sas)
    data = np.empty((n_channels, total_frames))
    logging.debug(f"Data shape: {data.shape}")

    # Fill data array one channel at a time
    for i, sa in enumerate(sas):
        logging.debug(f"Converting channel {i + 1} of {n_channels}")
        data[i, :] = SpikeAnalysisResult.convert_sa_to_np(sa, chunk_len)

    channel_names = [str(sa.recording.get_channel_ids().item()) for sa in sas]
    logging.debug(f"Channel names: {channel_names}")
    sfreq = sfreqs[0]

    # Extract spike times for each unit and create annotations
    onset = []
    description = []
    for sa in sas:
        for unit_id in sa.sorting.get_unit_ids():
            spike_train = sa.sorting.get_unit_spike_train(unit_id)
            # Convert to seconds and filter to recording duration
            spike_times = spike_train / sa.sorting.get_sampling_frequency()
            mask = spike_times < sa.recording.get_duration()
            spike_times = spike_times[mask]

            # Create annotation for each spike
            onset.extend(spike_times)
            description.extend(
                [sa.recording.get_channel_ids().item()] * len(spike_times)
            )  # collapse all units into 1 spike train
    annotations = mne.Annotations(onset, duration=0, description=description)

    info = mne.create_info(ch_names=channel_names, sfreq=sfreq, ch_types="eeg")
    raw = mne.io.RawArray(data=data, info=info)
    raw = raw.set_annotations(annotations)
    return raw

save_fif_and_json(folder, convert_to_mne=True, make_folder=True, slugify_filebase=True, save_abbrevs_as_chnames=False, overwrite=False)

Archive spike analysis result into the folder specified, as a fif and json file.

Parameters:

Name	Type	Description	Default
folder	str | Path	Destination folder to save results to	required
convert_to_mne	bool	If True, convert the SortingAnalyzers to a MNE RawArray if self.result_mne is None. Defaults to True.	True
make_folder	bool	If True, create the folder if it doesn't exist. Defaults to True.	True
slugify_filebase	bool	If True, slugify the filebase (replace special characters). Defaults to True.	True
save_abbrevs_as_chnames	bool	If True, save the channel abbreviations as the channel names in the json file. Defaults to False.	False
overwrite	bool	If True, overwrite the existing files. Defaults to False.	False

Source code in pythoneeg/visualization/results.py

def save_fif_and_json(
    self,
    folder: str | Path,
    convert_to_mne=True,
    make_folder=True,
    slugify_filebase=True,
    save_abbrevs_as_chnames=False,
    overwrite=False,
):
    """Archive spike analysis result into the folder specified, as a fif and json file.

    Args:
        folder (str | Path): Destination folder to save results to
        convert_to_mne (bool, optional): If True, convert the SortingAnalyzers to a MNE RawArray if self.result_mne is None. Defaults to True.
        make_folder (bool, optional): If True, create the folder if it doesn't exist. Defaults to True.
        slugify_filebase (bool, optional): If True, slugify the filebase (replace special characters). Defaults to True.
        save_abbrevs_as_chnames (bool, optional): If True, save the channel abbreviations as the channel names in the json file. Defaults to False.
        overwrite (bool, optional): If True, overwrite the existing files. Defaults to False.
    """
    if self.result_mne is None:
        if convert_to_mne:
            result_mne = self.convert_to_mne(save_raw=True)
            if result_mne is None:
                warnings.warn("No SortingAnalyzers found, skipping saving")
                return
        else:
            raise ValueError("No MNE RawArray found, and convert_to_mne is False. Run convert_to_mne() first.")
    else:
        result_mne = self.result_mne

    folder = Path(folder)
    if make_folder:
        folder.mkdir(parents=True, exist_ok=True)

    if slugify_filebase:
        filebase = folder / slugify(f"{self.animal_id}-{self.genotype}-{self.animal_day}")
    else:
        filebase = folder / f"{self.animal_id}-{self.genotype}-{self.animal_day}"
    filebase = str(filebase)

    if not overwrite:
        if filebase + ".json" in folder.glob("*.json"):
            raise FileExistsError(f"File {filebase}.json already exists")
        if filebase + ".fif" in folder.glob("*.fif"):
            raise FileExistsError(f"File {filebase}.fif already exists")
    else:
        for f in folder.glob("*"):
            f.unlink()
    result_mne.save(filebase + "-raw.fif", overwrite=overwrite)
    del result_mne

    json_dict = {
        "animal_id": self.animal_id,
        "genotype": self.genotype,
        "animal_day": self.animal_day,
        "bin_folder_name": self.bin_folder_name,
        "metadata": self.metadata.metadata_path,
        "channel_names": self.channel_abbrevs if save_abbrevs_as_chnames else self.channel_names,
        "assume_from_number": False if save_abbrevs_as_chnames else self.assume_from_number,
    }
    with open(filebase + ".json", "w") as f:
        json.dump(json_dict, f, indent=2)

WindowAnalysisResult

Bases: AnimalFeatureParser

Wrapper for output of windowed analysis. Has useful functions like group-wise and global averaging, filtering, and saving

Source code in pythoneeg/visualization/results.py

class WindowAnalysisResult(AnimalFeatureParser):
    """
    Wrapper for output of windowed analysis. Has useful functions like group-wise and global averaging, filtering, and saving
    """

    def __init__(
        self,
        result: pd.DataFrame,
        animal_id: str = None,
        genotype: str = None,
        channel_names: list[str] = None,
        assume_from_number=False,
        bad_channels_dict: dict[str, list[str]] = {},
        suppress_short_interval_error=False,
    ) -> None:
        """
        Args:
            result (pd.DataFrame): Result comes from AnimalOrganizer.compute_windowed_analysis()
            animal_id (str, optional): Identifier for the animal where result was computed from. Defaults to None.
            genotype (str, optional): Genotype of animal. Defaults to None.
            channel_names (list[str], optional): List of channel names. Defaults to None.
            assume_channels (bool, optional): If true, assumes channel names according to AnimalFeatureParser.DEFAULT_CHNUM_TO_NAME. Defaults to False.
            bad_channels_dict (dict[str, list[str]], optional): Dictionary of channels to reject for each recording session. Defaults to {}.
            suppress_short_interval_error (bool, optional): If True, suppress ValueError for short intervals between timestamps. Useful for aggregated WARs with large window sizes. Defaults to False.
        """
        self.result = result
        self.animal_id = animal_id
        self.genotype = genotype
        self.channel_names = channel_names
        self.assume_from_number = assume_from_number
        self.bad_channels_dict = bad_channels_dict
        self.suppress_short_interval_error = suppress_short_interval_error

        self.__update_instance_vars()

        print(f"Channel names: \t{self.channel_names}")
        print(f"Channel abbreviations: \t{self.channel_abbrevs}")

    def __str__(self) -> str:
        return f"{self.animaldays}"

    def __update_instance_vars(self):
        """Run after updating self.result, or other init values"""
        if "index" in self.result.columns:
            warnings.warn("Dropping column 'index'")
            self.result = self.result.drop(columns=["index"])

        # Check if timestamps are sorted and sort if needed
        if "timestamp" in self.result.columns:
            if not self.result["timestamp"].is_monotonic_increasing:
                warnings.warn("Timestamps are not sorted. Sorting result DataFrame by timestamp.")
                self.result = self.result.sort_values("timestamp")

        # Check for unusually short intervals between timestamps
        if "timestamp" in self.result.columns and "duration" in self.result.columns:
            median_duration = self.result["duration"].median()
            timestamp_diffs = self.result["timestamp"].diff()
            short_intervals = timestamp_diffs < pd.Timedelta(seconds=median_duration)

            # Skip first row since diff() produces NaT
            short_intervals = short_intervals[1:]

            if short_intervals.any():
                n_short = short_intervals.sum()
                pct_short = (n_short / len(short_intervals)) * 100

                warning_msg = (
                    f"Found {n_short} intervals ({pct_short:.1f}%) between timestamps "
                    f"that are shorter than the median duration of {median_duration:.1f}s"
                )

                if pct_short > 1.0 and not self.suppress_short_interval_error:  # More than 1% of intervals are short
                    raise ValueError(warning_msg)
                elif not self.suppress_short_interval_error:
                    warnings.warn(warning_msg)

        if "animal" in self.result.columns:
            unique_animals = self.result["animal"].unique()
            if len(unique_animals) > 1:
                raise ValueError(f"Multiple animals found in result: {unique_animals}")
            if unique_animals[0] != self.animal_id:
                raise ValueError(
                    f"Animal ID mismatch: result has {unique_animals[0]}, but self.animal_id is {self.animal_id}"
                )

        self._feature_columns = [x for x in self.result.columns if x in constants.FEATURES]
        self._nonfeature_columns = [x for x in self.result.columns if x not in constants.FEATURES]
        self.animaldays = self.result.loc[:, "animalday"].unique()

        self.channel_abbrevs = [
            core.parse_chname_to_abbrev(x, assume_from_number=self.assume_from_number) for x in self.channel_names
        ]

    def reorder_and_pad_channels(
        self, target_channels: list[str], use_abbrevs: bool = True, inplace: bool = True
    ) -> pd.DataFrame:
        """Reorder and pad channels to match a target channel list.

        This method ensures that the data has a consistent channel order and structure
        by reordering existing channels and padding missing channels with NaNs.

        Args:
            target_channels (list[str]): List of target channel names to match
            use_abbrevs (bool, optional): If True, target channel names are read as channel abbreviations instead of channel names. Defaults to True.
            inplace (bool, optional): If True, modify the result in place. Defaults to True.
        Returns:
            pd.DataFrame: DataFrame with reordered and padded channels
        """
        duplicates = [ch for ch in target_channels if target_channels.count(ch) > 1]
        if duplicates:
            raise ValueError(f"Target channels must be unique. Found duplicates: {duplicates}")

        result = self.result.copy()

        channel_map = {ch: i for i, ch in enumerate(target_channels)}
        channel_names = self.channel_names if not use_abbrevs else self.channel_abbrevs

        valid_channels = [ch for ch in channel_names if ch in channel_map]
        if not valid_channels:
            warnings.warn(
                f"None of the channel names {channel_names} were found in target channels {target_channels}. Is use_abbrevs correctly set?"
            )

        for feature in self._feature_columns:
            match feature:
                case _ if feature in constants.LINEAR_FEATURES + constants.BAND_FEATURES:
                    if feature in constants.BAND_FEATURES:
                        df_bands = pd.DataFrame(result[feature].tolist())
                        vals = np.array(df_bands.values.tolist())
                        vals = vals.transpose((0, 2, 1))
                        keys = df_bands.keys()
                    else:
                        vals = np.array(result[feature].tolist())

                    new_vals = np.full((vals.shape[0], len(target_channels), *vals.shape[2:]), np.nan)  # dubious

                    for i, ch in enumerate(channel_names):
                        if ch in channel_map:
                            new_vals[:, channel_map[ch]] = vals[:, i]

                    if feature in constants.BAND_FEATURES:
                        new_vals = new_vals.transpose((0, 2, 1))
                        result[feature] = [dict(zip(keys, vals)) for vals in new_vals]
                    else:
                        result[feature] = [list(x) for x in new_vals]

                case _ if feature in constants.MATRIX_FEATURES:
                    if feature in ["cohere", "zcohere", "imcoh", "zimcoh"]:
                        df_bands = pd.DataFrame(result[feature].tolist())
                        vals = np.array(df_bands.values.tolist())
                        keys = df_bands.keys()
                    else:
                        vals = np.array(result[feature].tolist())

                    logging.debug(f"vals.shape: {vals.shape}")
                    new_shape = list(vals.shape[:-2]) + [len(target_channels), len(target_channels)]
                    new_vals = np.full(new_shape, np.nan)

                    ch1_valid = np.array([ch in channel_map for ch in channel_names])
                    ch2_valid = ch1_valid.copy()
                    valid_pairs = np.logical_and(ch1_valid[:, None], ch2_valid[None, :])  # 2D boolean mask

                    for i, j in zip(*np.where(valid_pairs)):
                        ch1, ch2 = channel_names[i], channel_names[j]
                        new_vals[..., channel_map[ch1], channel_map[ch2]] = vals[..., i, j]

                    triu_mask = np.triu_indices(len(target_channels), k=0)
                    new_vals += new_vals.transpose((*range(new_vals.ndim - 2), -1, -2))
                    new_vals[..., triu_mask[0], triu_mask[1]] = 0

                    if feature in ["cohere", "zcohere", "imcoh", "zimcoh"]:
                        result[feature] = [dict(zip(keys, vals)) for vals in new_vals]
                    else:
                        result[feature] = [list(x) for x in new_vals]

                case _ if feature in constants.HIST_FEATURES:
                    coords = np.array([x[0] for x in result[feature].tolist()])
                    vals = np.array([x[1] for x in result[feature].tolist()])
                    new_vals = np.full((*vals.shape[0:-1], len(target_channels)), np.nan)

                    for i, ch in enumerate(channel_names):
                        if ch in channel_map:
                            new_vals[:, ..., channel_map[ch]] = vals[:, ..., i]

                    result[feature] = [(coords[i], new_vals[i]) for i in range(len(coords))]

                case _:
                    raise ValueError(f"Invalid feature: {feature}")

        if inplace:
            self.result = result

            logging.debug(f"Old channel names: {self.channel_names}")
            self.channel_names = target_channels
            logging.debug(f"New channel names: {self.channel_names}")

            logging.debug(f"Old channel abbreviations: {self.channel_abbrevs}")
            self.__update_instance_vars()
            logging.debug(f"New channel abbreviations: {self.channel_abbrevs}")

        return result

    def read_sars_spikes(self, sars: list["SpikeAnalysisResult"], read_mode: Literal["sa", "mne"] = "sa", inplace=True):
        match read_mode:
            case "sa":
                spikes_all = []
                for sar in sars:  # for each continuous recording session
                    spikes_channel = []
                    for i, sa in enumerate(sar.result_sas):  # for each channel
                        spike_times = []
                        for unit in sa.sorting.get_unit_ids():  # Flatten units
                            spike_times.extend(sa.sorting.get_unit_spike_train(unit_id=unit).tolist())
                        spike_times = np.array(spike_times) / sa.sorting.get_sampling_frequency()
                        spikes_channel.append(spike_times)
                    spikes_all.append(spikes_channel)
                return self._read_from_spikes_all(spikes_all, inplace=inplace)
            case "mne":
                raws = [sar.result_mne for sar in sars]
                return self.read_mnes_spikes(raws, inplace=inplace)
            case _:
                raise ValueError(f"Invalid read_mode: {read_mode}")

    def read_mnes_spikes(self, raws: list[mne.io.RawArray], inplace=True):
        spikes_all = []
        for raw in raws:
            # each mne is a contiguous recording session
            events, event_id = mne.events_from_annotations(raw)
            event_id = {k.item(): v for k, v in event_id.items()}

            spikes_channel = []
            for channel in raw.ch_names:
                if channel not in event_id.keys():
                    logging.warning(f"Channel {channel} not found in event_id")
                    spikes_channel.append([])
                    continue
                event_id_channel = event_id[channel]
                spike_times = events[events[:, 2] == event_id_channel, 0]
                spike_times = spike_times / raw.info["sfreq"]
                spikes_channel.append(spike_times)
            spikes_all.append(spikes_channel)
        return self._read_from_spikes_all(spikes_all, inplace=inplace)

    def _read_from_spikes_all(self, spikes_all: list[list[list[float]]], inplace=True):
        # Each groupby animalday is a recording session
        grouped = self.result.groupby("animalday")
        animaldays = grouped.groups.keys()
        logging.debug(f"Animal days: {animaldays}")
        spike_counts = dict(zip(animaldays, spikes_all))
        spike_counts = grouped.apply(lambda x: _bin_spike_df(x, spikes_channel=spike_counts[x.name]))
        spike_counts: pd.Series = spike_counts.explode()

        if spike_counts.size != self.result.shape[0]:
            logging.warning(f"Spike counts size {spike_counts.size} does not match result size {self.result.shape[0]}")

        result = self.result.copy()
        result["nspike"] = spike_counts.tolist()
        result["lognspike"] = list(core.log_transform(np.stack(result["nspike"].tolist(), axis=0)))
        if inplace:
            self.result = result
        return result

    def get_info(self):
        """Returns a formatted string with basic information about the WindowAnalysisResult object"""
        info = []
        info.append(f"feature names: {', '.join(self._feature_columns)}")
        info.append(f"animaldays: {', '.join(self.result['animalday'].unique())}")
        info.append(
            f"animal_id: {self.result['animal'].unique()[0] if 'animal' in self.result.columns else self.animal_id}"
        )
        info.append(
            f"genotype: {self.result['genotype'].unique()[0] if 'genotype' in self.result.columns else self.genotype}"
        )
        info.append(f"channel_names: {', '.join(self.channel_names) if self.channel_names else 'None'}")

        return "\n".join(info)

    def get_result(self, features: list[str], exclude: list[str] = [], allow_missing=False):
        """Get windowed analysis result dataframe, with helpful filters

        Args:
            features (list[str]): List of features to get from result
            exclude (list[str], optional): List of features to exclude from result; will override the features parameter. Defaults to [].
            allow_missing (bool, optional): If True, will return all requested features as columns regardless if they exist in result. Defaults to False.

        Returns:
            result: pd.DataFrame object with features in columns and windows in rows
        """
        features = _sanitize_feature_request(features, exclude)
        if not allow_missing:
            return self.result.loc[:, self._nonfeature_columns + features]
        else:
            return self.result.reindex(columns=self._nonfeature_columns + features)

    def get_groupavg_result(
        self, features: list[str], exclude: list[str] = [], df: pd.DataFrame = None, groupby="animalday"
    ):
        """Group result and average within groups. Preserves data structure and shape for each feature.

        Args:
            features (list[str]): List of features to get from result
            exclude (list[str], optional): List of features to exclude from result. Will override the features parameter. Defaults to [].
            df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
            groupby (str, optional): Feature or list of features to group by before averaging. Passed to the `by` parameter in pd.DataFrame.groupby(). Defaults to "animalday".

        Returns:
            grouped_result: result grouped by `groupby` and averaged for each group.
        """
        result_grouped, result_validcols = self.__get_groups(features=features, exclude=exclude, df=df, groupby=groupby)
        features = _sanitize_feature_request(features, exclude)

        avg_results = []
        for f in features:
            if f in result_validcols:
                avg_result_col = result_grouped.apply(self._average_feature, f, "duration", include_groups=False)
                avg_result_col.name = f
                avg_results.append(avg_result_col)
            else:
                logging.warning(f"{f} not calculated, skipping")

        return pd.concat(avg_results, axis=1)

    def __get_groups(self, features: list[str], exclude: list[str] = [], df: pd.DataFrame = None, groupby="animalday"):
        features = _sanitize_feature_request(features, exclude)
        result_win = self.result if df is None else df
        return result_win.groupby(groupby), result_win.columns

    def get_grouprows_result(
        self,
        features: list[str],
        exclude: list[str] = [],
        df: pd.DataFrame = None,
        multiindex=["animalday", "animal", "genotype"],
        include=["duration", "endfile"],
    ):
        features = _sanitize_feature_request(features, exclude)
        result_win = self.result if df is None else df
        result_win = result_win.filter(features + multiindex + include)
        return result_win.set_index(multiindex)

    def get_filter_logrms_range(self, df: pd.DataFrame = None, z_range=3, **kwargs):
        """Filter windows based on log(rms).

        Args:
            df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
            z_range (float, optional): The z-score range to filter by. Values outside this range will be set to NaN.

        Returns:
            out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
        """
        result = df.copy() if df is not None else self.result.copy()
        z_range = abs(z_range)
        np_rms = np.array(result["rms"].tolist())
        np_logrms = np.log(np_rms)
        del np_rms
        np_logrmsz = zscore(np_logrms, axis=0, nan_policy="omit")
        np_logrms[(np_logrmsz > z_range) | (np_logrmsz < -z_range)] = np.nan

        out = np.full(np_logrms.shape, True)
        out[(np_logrmsz > z_range) | (np_logrmsz < -z_range)] = False
        return out

    def get_filter_high_rms(self, df: pd.DataFrame = None, max_rms=500, **kwargs):
        """Filter windows based on rms.

        Args:
            df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
            max_rms (float, optional): The maximum rms value to filter by. Values above this will be set to NaN.

        Returns:
            out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
        """
        result = df.copy() if df is not None else self.result.copy()
        np_rms = np.array(result["rms"].tolist())
        np_rmsnan = np_rms.copy()
        # Convert to float to allow NaN assignment for integer arrays
        if np_rmsnan.dtype.kind in ("i", "u"):  # integer types
            np_rmsnan = np_rmsnan.astype(float)
        np_rmsnan[np_rms > max_rms] = np.nan
        result["rms"] = np_rmsnan.tolist()

        out = np.full(np_rms.shape, True)
        out[np_rms > max_rms] = False
        return out

    def get_filter_low_rms(self, df: pd.DataFrame = None, min_rms=30, **kwargs):
        """Filter windows based on rms.

        Args:
            df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
            min_rms (float, optional): The minimum rms value to filter by. Values below this will be set to NaN.

        Returns:
            out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
        """
        result = df.copy() if df is not None else self.result.copy()
        np_rms = np.array(result["rms"].tolist())
        np_rmsnan = np_rms.copy()
        np_rmsnan[np_rms < min_rms] = np.nan
        result["rms"] = np_rmsnan.tolist()

        out = np.full(np_rms.shape, True)
        out[np_rms < min_rms] = False
        return out

    def get_filter_high_beta(self, df: pd.DataFrame = None, max_beta_prop=0.4, **kwargs):
        """Filter windows based on beta power.

        Args:
            df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
            max_beta_prop (float, optional): The maximum beta power to filter by. Values above this will be set to NaN. Defaults to 0.4.

        Returns:
            out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
        """
        result = df.copy() if df is not None else self.result.copy()
        if "psdfrac" in result.columns:
            df_psdfrac = pd.DataFrame(result["psdfrac"].tolist())
            np_prop = np.array(df_psdfrac["beta"].tolist())
        elif "psdband" in result.columns and "psdtotal" in result.columns:
            df_psdband = pd.DataFrame(result["psdband"].tolist())
            np_beta = np.array(df_psdband["beta"].tolist())
            np_total = np.array(result["psdtotal"].tolist())
            np_prop = np_beta / np_total
        else:
            raise ValueError("psdfrac or psdband+psdtotal required for beta power filtering")

        out = np.full(np_prop.shape, True)
        out[np_prop > max_beta_prop] = False
        out = np.broadcast_to(np.all(out, axis=-1)[:, np.newaxis], out.shape)
        return out

    def get_filter_reject_channels(self, bad_channels: list[str], use_abbrevs: bool = None, **kwargs):
        """Filter channels to reject.

        Args:
            bad_channels (list[str]): List of channels to reject. Can be either full channel names or abbreviations.
                The method will automatically detect which format is being used.
            use_abbrevs (bool, optional): Override automatic detection. If True, channels are assumed to be channel abbreviations. If False, channels are assumed to be channel names.
                If None, channels are parsed to abbreviations and matched against self.channel_abbrevs.

        Returns:
            out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
        """
        channel_targets = (
            self.channel_abbrevs if use_abbrevs or use_abbrevs is None else self.channel_names
        )  # Match to appropriate target
        if use_abbrevs is None:  # Match channels as abbreviations
            bad_channels = [
                core.utils.parse_chname_to_abbrev(ch, assume_from_number=self.assume_from_number) for ch in bad_channels
            ]

        n_samples = len(self.result)
        n_channels = len(channel_targets)
        mask = np.ones((n_samples, n_channels), dtype=bool)

        # Match channels to channel_targets
        for ch in bad_channels:
            if ch in channel_targets:
                mask[:, channel_targets.index(ch)] = False
            else:
                warnings.warn(f"Channel {ch} not found in {channel_targets}")
        return mask

    def get_filter_reject_channels_by_recording_session(
        self, bad_channels_dict: dict[str, list[str]] = None, use_abbrevs: bool = None
    ):
        """Filter channels to reject for each recording session

        Args:
            bad_channels_dict (dict[str, list[str]]): Dictionary of list of channels to reject for each recording session.
                Can be either full channel names or abbreviations. The method will automatically detect which format is being used.
                If None, the method will use the bad_channels_dict passed to the constructor.
            use_abbrevs (bool, optional): Override automatic detection. If True, channels are assumed to be channel abbreviations. If False, channels are assumed to be channel names.
                If None, channels are parsed to abbreviations and matched against self.channel_abbrevs.

        Returns:
            out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
        """
        if bad_channels_dict is None:
            bad_channels_dict = self.bad_channels_dict

        n_samples = len(self.result)
        n_channels = len(self.channel_names)
        mask = np.ones((n_samples, n_channels), dtype=bool)

        # Group by animalday to apply filters per recording session
        for animalday, group in self.result.groupby("animalday"):
            if animalday not in bad_channels_dict:
                logging.warning(f"No bad channels specified for recording session {animalday}")
                continue

            bad_channels = bad_channels_dict[animalday]
            channel_targets = self.channel_abbrevs if use_abbrevs or use_abbrevs is None else self.channel_names
            if use_abbrevs is None:
                bad_channels = [
                    core.parse_chname_to_abbrev(ch, assume_from_number=self.assume_from_number) for ch in bad_channels
                ]

            # Get indices for this recording session
            session_indices = group.index

            # Apply channel filtering for this session
            for ch in bad_channels:
                if ch in channel_targets:
                    ch_idx = channel_targets.index(ch)
                    mask[session_indices, ch_idx] = False
                else:
                    logging.warning(f"Channel {ch} not found in {channel_targets} for session {animalday}")

        return mask

    def get_filter_morphological_smoothing(
        self, filter_mask: np.ndarray, smoothing_seconds: float, **kwargs
    ) -> np.ndarray:
        """Apply morphological smoothing to a filter mask.

        Args:
            filter_mask (np.ndarray): Input boolean mask of shape (n_windows, n_channels)
            smoothing_seconds (float): Time window in seconds for morphological operations

        Returns:
            np.ndarray: Smoothed boolean mask
        """
        if "duration" not in self.result.columns:
            raise ValueError("Cannot calculate window duration - 'duration' column missing")

        window_duration = self.result["duration"].median()
        structure_size = max(1, int(smoothing_seconds / window_duration))

        if structure_size <= 1:
            return filter_mask

        smoothed_mask = filter_mask.copy()
        for ch_idx in range(filter_mask.shape[1]):
            channel_mask = filter_mask[:, ch_idx]
            # Opening removes small isolated artifacts
            channel_mask = binary_opening(channel_mask, structure=np.ones(structure_size))
            # Closing fills small gaps in valid data
            channel_mask = binary_closing(channel_mask, structure=np.ones(structure_size))
            smoothed_mask[:, ch_idx] = channel_mask

        return smoothed_mask

    def filter_morphological_smoothing(self, smoothing_seconds: float) -> "WindowAnalysisResult":
        """Apply morphological smoothing to all data.

        Args:
            smoothing_seconds (float): Time window in seconds for morphological operations

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        # Start with all-True mask and smooth it
        base_mask = np.ones((len(self.result), len(self.channel_names)), dtype=bool)
        smoothed_mask = self.get_filter_morphological_smoothing(base_mask, smoothing_seconds)
        return self._create_filtered_copy(smoothed_mask)

    def filter_all(
        self,
        df: pd.DataFrame = None,
        inplace=True,
        bad_channels: list[str] = None,
        min_valid_channels=3,
        filters: list[callable] = None,
        morphological_smoothing_seconds: float = None,
        **kwargs,
    ):
        """Apply a list of filters to the data. Filtering should be performed before aggregation.

        Args:
            df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
            inplace (bool, optional): If True, modify the result in place. Defaults to True.
            bad_channels (list[str], optional): List of channels to reject. Defaults to None.
            min_valid_channels (int, optional): Minimum number of valid channels required per window. Defaults to 3.
            filters (list[callable], optional): List of filter functions to apply. Each function should return a boolean mask.
                If None, uses default filters: [get_filter_logrms_range, get_filter_high_rms, get_filter_low_rms, get_filter_high_beta].
                Defaults to None.
            morphological_smoothing_seconds (float, optional): If provided, apply morphological opening/closing to smooth the filter mask.
                This removes isolated false positives/negatives along the time axis for each channel independently.
                The value specifies the time window in seconds for the morphological operations. Defaults to None.
            **kwargs: Additional keyword arguments to pass to filter functions.

        Returns:
            WindowAnalysisResult: Filtered result
        """
        if filters is None:
            # TODO refactor these into standalone functions, which take in a war as the first parameter, then pass
            # filt_bool = filt(self, df, **kwargs) as needed
            filters = [
                self.get_filter_logrms_range,
                self.get_filter_high_rms,
                self.get_filter_low_rms,
                self.get_filter_high_beta,
                self.get_filter_reject_channels_by_recording_session,
            ]

        filt_bools = []
        # Apply each filter function
        for filter_function in filters:
            filt_bool = filter_function(df, **kwargs)
            filt_bools.append(filt_bool)
            logging.info(
                f"{filter_function.__name__}:\tfiltered {filt_bool.size - np.count_nonzero(filt_bool)}/{filt_bool.size}"
            )

        # Filter channels manually
        # REVIEW somehow add this to the main list of filters, but I'm not sure how to do this.
        if bad_channels is not None:
            filt_bools.append(self.get_filter_reject_channels(bad_channels))
            logging.debug(f"Reject channels: {filt_bools[-1]}")

        # Apply all filters
        filt_bool_all = np.prod(np.stack(filt_bools, axis=-1), axis=-1).astype(bool)
        logging.debug(f"filt_bool_all.shape: {filt_bool_all.shape}")  # (windows, channels)

        # Apply morphological smoothing if requested
        if morphological_smoothing_seconds is not None:
            if "duration" not in self.result.columns:
                raise ValueError("Cannot calculate window duration - 'duration' column missing from result dataframe")
            window_duration = self.result["duration"].median()

            # Calculate number of windows for the smoothing
            structure_size = max(1, int(morphological_smoothing_seconds / window_duration))

            if structure_size > 1:
                logging.info(
                    f"Applying morphological smoothing with {structure_size} windows ({morphological_smoothing_seconds}s / {window_duration}s per window)"
                )
                # Apply channel-wise temporal smoothing (each channel processed independently)
                # This avoids spatial assumptions while smoothing temporal artifacts
                for ch_idx in range(filt_bool_all.shape[1]):
                    channel_mask = filt_bool_all[:, ch_idx]
                    # Opening removes small isolated artifacts
                    channel_mask = binary_opening(channel_mask, structure=np.ones(structure_size))
                    # Closing fills small gaps in valid data
                    channel_mask = binary_closing(channel_mask, structure=np.ones(structure_size))
                    filt_bool_all[:, ch_idx] = channel_mask
            else:
                logging.info("Skipping morphological smoothing - structure size would be 1 (no effect)")

        # Filter windows based on number of valid channels
        valid_channels_per_window = np.sum(filt_bool_all, axis=1)  # axis 1 = channel
        window_mask = valid_channels_per_window >= min_valid_channels  # True if window has enough valid channels
        filt_bool_all = filt_bool_all & window_mask[:, np.newaxis]  # Apply window mask to all channels

        filtered_result = self._apply_filter(filt_bool_all)
        if inplace:
            del self.result
            self.result = filtered_result
        return WindowAnalysisResult(
            filtered_result,
            self.animal_id,
            self.genotype,
            self.channel_names,
            self.assume_from_number,
            self.bad_channels_dict,
            self.suppress_short_interval_error,
        )

    def _create_filtered_copy(self, filter_mask: np.ndarray) -> "WindowAnalysisResult":
        """Create a new WindowAnalysisResult with the filter applied.

        Args:
            filter_mask (np.ndarray): Boolean mask of shape (n_windows, n_channels)

        Returns:
            WindowAnalysisResult: New instance with filter applied
        """
        filtered_result = self._apply_filter(filter_mask)
        return WindowAnalysisResult(
            filtered_result,
            self.animal_id,
            self.genotype,
            self.channel_names,
            self.assume_from_number,
            self.bad_channels_dict,
            self.suppress_short_interval_error,
        )

    def filter_logrms_range(self, z_range: float = 3) -> "WindowAnalysisResult":
        """Filter based on log(rms) z-score range.

        Args:
            z_range (float): Z-score range threshold. Defaults to 3.

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        mask = self.get_filter_logrms_range(z_range=z_range)
        return self._create_filtered_copy(mask)

    def filter_high_rms(self, max_rms: float = 500) -> "WindowAnalysisResult":
        """Filter out windows with RMS above threshold.

        Args:
            max_rms (float): Maximum RMS threshold. Defaults to 500.

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        mask = self.get_filter_high_rms(max_rms=max_rms)
        return self._create_filtered_copy(mask)

    def filter_low_rms(self, min_rms: float = 50) -> "WindowAnalysisResult":
        """Filter out windows with RMS below threshold.

        Args:
            min_rms (float): Minimum RMS threshold. Defaults to 50.

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        mask = self.get_filter_low_rms(min_rms=min_rms)
        return self._create_filtered_copy(mask)

    def filter_high_beta(self, max_beta_prop: float = 0.4) -> "WindowAnalysisResult":
        """Filter out windows with high beta power.

        Args:
            max_beta_prop (float): Maximum beta power proportion. Defaults to 0.4.

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        mask = self.get_filter_high_beta(max_beta_prop=max_beta_prop)
        return self._create_filtered_copy(mask)

    def filter_reject_channels(self, bad_channels: list[str], use_abbrevs: bool = None) -> "WindowAnalysisResult":
        """Filter out specified bad channels.

        Args:
            bad_channels (list[str]): List of channel names to reject
            use_abbrevs (bool, optional): Whether to use abbreviations. Defaults to None.

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        mask = self.get_filter_reject_channels(bad_channels, use_abbrevs=use_abbrevs)
        return self._create_filtered_copy(mask)

    def filter_reject_channels_by_session(
        self, bad_channels_dict: dict[str, list[str]] = None, use_abbrevs: bool = None
    ) -> "WindowAnalysisResult":
        """Filter out bad channels by recording session.

        Args:
            bad_channels_dict (dict[str, list[str]], optional): Dict of bad channels per session
            use_abbrevs (bool, optional): Whether to use abbreviations. Defaults to None.

        Returns:
            WindowAnalysisResult: New filtered instance
        """
        mask = self.get_filter_reject_channels_by_recording_session(bad_channels_dict, use_abbrevs=use_abbrevs)
        return self._create_filtered_copy(mask)

    def apply_filters(
        self, filter_config: dict = None, min_valid_channels: int = 3, morphological_smoothing_seconds: float = None
    ) -> "WindowAnalysisResult":
        """Apply multiple filters using configuration.

        Args:
            filter_config (dict, optional): Dictionary of filter names and parameters.
                Available filters: 'logrms_range', 'high_rms', 'low_rms', 'high_beta',
                'reject_channels', 'reject_channels_by_session', 'morphological_smoothing'
            min_valid_channels (int): Minimum valid channels per window. Defaults to 3.
            morphological_smoothing_seconds (float, optional): Temporal smoothing window (deprecated, use config instead)

        Returns:
            WindowAnalysisResult: New filtered instance

        Examples:
            >>> config = {
            ...     'logrms_range': {'z_range': 3},
            ...     'high_rms': {'max_rms': 500},
            ...     'reject_channels': {'bad_channels': ['LMot', 'RMot']},
            ...     'morphological_smoothing': {'smoothing_seconds': 8.0}
            ... }
            >>> filtered_war = war.apply_filters(config)
        """
        if filter_config is None:
            filter_config = {
                "logrms_range": {"z_range": 3},
                "high_rms": {"max_rms": 500},
                "low_rms": {"min_rms": 50},
                "high_beta": {"max_beta_prop": 0.4},
                "reject_channels_by_session": {},
            }

        filter_methods = {
            "logrms_range": self.get_filter_logrms_range,
            "high_rms": self.get_filter_high_rms,
            "low_rms": self.get_filter_low_rms,
            "high_beta": self.get_filter_high_beta,
            "reject_channels": self.get_filter_reject_channels,
            "reject_channels_by_session": self.get_filter_reject_channels_by_recording_session,
        }

        filt_bools = []
        morphological_params = None

        for filter_name, filter_params in filter_config.items():
            if filter_name == "morphological_smoothing":
                morphological_params = filter_params
                continue

            if filter_name not in filter_methods:
                raise ValueError(
                    f"Unknown filter: {filter_name}. Available: {list(filter_methods.keys()) + ['morphological_smoothing']}"
                )

            filter_func = filter_methods[filter_name]
            filt_bool = filter_func(**filter_params)
            filt_bools.append(filt_bool)
            logging.info(f"{filter_name}: filtered {filt_bool.size - np.count_nonzero(filt_bool)}/{filt_bool.size}")

        # Combine all filter masks
        if filt_bools:
            filt_bool_all = np.prod(np.stack(filt_bools, axis=-1), axis=-1).astype(bool)
        else:
            filt_bool_all = np.ones((len(self.result), len(self.channel_names)), dtype=bool)

        # Apply morphological smoothing if requested (either from config or parameter)
        if morphological_params or morphological_smoothing_seconds is not None:
            if morphological_params:
                smoothing_seconds = morphological_params["smoothing_seconds"]
            else:
                smoothing_seconds = morphological_smoothing_seconds

            filt_bool_all = self.get_filter_morphological_smoothing(filt_bool_all, smoothing_seconds)
            logging.info(f"Applied morphological smoothing: {smoothing_seconds}s")

        # Filter windows based on minimum valid channels
        valid_channels_per_window = np.sum(filt_bool_all, axis=1)
        window_mask = valid_channels_per_window >= min_valid_channels
        filt_bool_all = filt_bool_all & window_mask[:, np.newaxis]

        return self._create_filtered_copy(filt_bool_all)

    def _apply_filter(self, filter_tfs: np.ndarray):
        result = self.result.copy()
        filter_tfs = np.array(filter_tfs, dtype=bool)  # (M fragments, N channels)
        for feat in constants.FEATURES:
            if feat not in result.columns:
                logging.info(f"Skipping {feat} because it is not in result")
                continue
            logging.info(f"Filtering {feat}")
            match feat:  # NOTE refactor this to use constants
                case "rms" | "ampvar" | "psdtotal" | "nspike" | "logrms" | "logampvar" | "logpsdtotal" | "lognspike":
                    vals = np.array(result[feat].tolist())
                    # Convert to float to allow NaN assignment for integer features
                    if vals.dtype.kind in ("i", "u"):  # integer types
                        vals = vals.astype(float)
                    vals[~filter_tfs] = np.nan
                    result[feat] = vals.tolist()
                case "psd":
                    # FIXME The sampling rates have changed between computation passes so WARs have different shapes.
                    # Add a check for same sampling frequency, other war-relevant properties etc.
                    # The logging lines below should be removed at some point, but I'll keep it this way for now
                    logging.info(
                        f"set([x[0].shape for x in result[feat].tolist()]) = {list(set([x[0].shape for x in result[feat].tolist()]))}"
                    )
                    logging.info(
                        f"set([x[1].shape for x in result[feat].tolist()]) = {list(set([x[1].shape for x in result[feat].tolist()]))}"
                    )
                    coords = np.array([x[0] for x in result[feat].tolist()])
                    vals = np.array([x[1] for x in result[feat].tolist()])
                    mask = np.broadcast_to(filter_tfs[:, np.newaxis, :], vals.shape)
                    vals[~mask] = np.nan
                    outs = [(c, vals[i, :, :]) for i, c in enumerate(coords)]
                    result[feat] = outs
                case "psdband" | "psdfrac" | "logpsdband" | "logpsdfrac":
                    vals = pd.DataFrame(result[feat].tolist())
                    for colname in vals.columns:
                        v = np.array(vals[colname].tolist())
                        v[~filter_tfs] = np.nan
                        vals[colname] = v.tolist()
                    result[feat] = vals.to_dict("records")
                case "psdslope":
                    vals = np.array(result[feat].tolist())
                    mask = np.broadcast_to(filter_tfs[:, :, np.newaxis], vals.shape)
                    vals[~mask] = np.nan
                    # vals = [list(map(tuple, x)) for x in vals.tolist()]
                    result[feat] = vals.tolist()
                case "cohere" | "zcohere" | "imcoh" | "zimcoh":
                    vals = pd.DataFrame(result[feat].tolist())
                    shape = np.array(vals.iloc[:, 0].tolist()).shape
                    mask = np.broadcast_to(filter_tfs[:, :, np.newaxis], shape)
                    for colname in vals.columns:
                        v = np.array(vals[colname].tolist())
                        v[~mask] = np.nan
                        v[~mask.transpose(0, 2, 1)] = np.nan
                        vals[colname] = v.tolist()
                    result[feat] = vals.to_dict("records")
                case "pcorr" | "zpcorr":
                    vals = np.array(result[feat].tolist())
                    mask = np.broadcast_to(filter_tfs[:, :, np.newaxis], vals.shape)
                    vals[~mask] = np.nan
                    vals[~mask.transpose(0, 2, 1)] = np.nan
                    result[feat] = vals.tolist()
                case _:
                    raise ValueError(f"Unknown feature to filter {feat}")
        return result

    def save_pickle_and_json(
        self, folder: str | Path, make_folder=True, slugify_filebase=True, save_abbrevs_as_chnames=False
    ):
        """Archive window analysis result into the folder specified, as a pickle and json file.

        Args:
            folder (str | Path): Destination folder to save results to
            make_folder (bool, optional): If True, create the folder if it doesn't exist. Defaults to True.
            slugify_filebase (bool, optional): If True, slugify the filebase (replace special characters). Defaults to True.
            save_abbrevs_as_chnames (bool, optional): If True, save the channel abbreviations as the channel names in the json file. Defaults to False.
        """
        folder = Path(folder)
        if make_folder:
            folder.mkdir(parents=True, exist_ok=True)

        if slugify_filebase:
            filebase = folder / slugify(f"{self.animal_id}-{self.genotype}")
        else:
            filebase = folder / f"{self.animal_id}-{self.genotype}"
        filebase = str(filebase)

        self.result.to_pickle(filebase + ".pkl")
        json_dict = {
            "animal_id": self.animal_id,
            "genotype": self.genotype,
            "channel_names": self.channel_abbrevs if save_abbrevs_as_chnames else self.channel_names,
            "assume_from_number": False if save_abbrevs_as_chnames else self.assume_from_number,
            "bad_channels_dict": self.bad_channels_dict,
            "suppress_short_interval_error": self.suppress_short_interval_error,
        }

        with open(filebase + ".json", "w") as f:
            json.dump(json_dict, f, indent=2)

    @classmethod
    def load_pickle_and_json(cls, folder_path=None):
        """Load WindowAnalysisResult from folder

        Args:
            folder_path (str, optional): Path of folder containing one .pkl and .json file each. Defaults to None.
            df_pickle_path (str, optional): Path of .pkl file. If this and folder_path are not None, raises an error. Defaults to None.
            json_path (str, optional): Path of .json file. If this and folder_path are not None, raises an error. Defaults to None.

        Raises:
            ValueError: Both df_pickle_path and json_path must be None if folder_path is provided
            ValueError: Expected exactly one pickle and one json file in folder_path

        Returns:
            result: WindowAnalysisResult object
        """
        if folder_path is not None:
            folder_path = Path(folder_path)
            if not folder_path.exists():
                raise ValueError(f"Folder path {folder_path} does not exist")

            pkl_files = list(folder_path.glob("*.pkl"))
            json_files = list(folder_path.glob("*.json"))

            if len(pkl_files) != 1 or len(json_files) != 1:
                raise ValueError(f"Expected exactly one pickle and one json file in {folder_path}")

            df_pickle_path = pkl_files[0]
            json_path = json_files[0]

        with open(df_pickle_path, "rb") as f:
            data = pd.read_pickle(f)
        with open(json_path, "r") as f:
            metadata = json.load(f)
        return cls(data, **metadata)

    def aggregate_time_windows(self, groupby: list[str] | str = ["animalday", "isday"]) -> None:
        """Aggregate time windows into a single data point per groupby by averaging features. This reduces the number of rows in the result.

        Args:
            groupby (list[str] | str, optional): Columns to group by. Defaults to ['animalday', 'isday'], which groups by animalday (recording session) and isday (day/night).

        Raises:
            ValueError: groupby must be from ['animalday', 'isday']
            ValueError: Columns in groupby not found in result
            ValueError: Columns in groupby are not constant in groups
        """
        if isinstance(groupby, str):
            groupby = [groupby]
        if not all(col in ["animalday", "isday"] for col in groupby):
            raise ValueError(f"groupby must be from ['animalday', 'isday']. Got {groupby}")
        if not all(col in self.result.columns for col in groupby):
            raise ValueError(f"Columns {groupby} not found in result. Columns: {self.result.columns.tolist()}")

        features = [f for f in constants.FEATURES if f in self.result.columns]
        logging.debug(f"Aggregating {features}")
        result_grouped = self.result.groupby(groupby)

        agg_dict = {}

        if "animalday" not in groupby:
            agg_dict["animalday"] = lambda df: None
        if "isday" not in groupby:
            agg_dict["isday"] = lambda df: None

        constant_cols = ["animal", "day", "genotype"]
        for col in constant_cols:
            if col in self.result.columns:
                is_constant = result_grouped[col].nunique() == 1
                if not is_constant.all():
                    non_constant_groups = is_constant[~is_constant].index.tolist()
                    raise ValueError(f"Column {col} is not constant in groups: {non_constant_groups}")
                agg_dict[col] = lambda df, col=col: df[col].iloc[0]

        if "duration" in self.result.columns:
            agg_dict["duration"] = lambda df: np.sum(df["duration"])

        if "endfile" in self.result.columns:
            agg_dict["endfile"] = lambda df: df["endfile"].iloc[-1]

        if "timestamp" in self.result.columns:
            agg_dict["timestamp"] = lambda df: df["timestamp"].iloc[0]

        for feat in features:
            agg_dict[feat] = lambda df, feat=feat: self._average_feature(df, feat, "duration")

        aggregated_df = result_grouped.apply(
            lambda df: pd.Series({col: agg_dict[col](df) for col in self.result.columns if col not in groupby})
        )

        self.result = aggregated_df.reset_index(drop=False)  # Keep animalday/isday as a column

        self.suppress_short_interval_error = True
        logging.info("Setting suppress_short_interval_error to True")
        self.__update_instance_vars()

    def add_unique_hash(self, nbytes: int | None = None):
        """Adds a hex hash to the animal ID to ensure uniqueness. This prevents collisions when, for example, multiple animals in ExperimentPlotter have the same animal ID.

        Args:
            nbytes (int, optional): Number of bytes to generate. This is passed directly to secrets.token_hex(). Defaults to None, which generates 16 hex characters (8 bytes).
        """
        import secrets

        hash_suffix = secrets.token_hex(nbytes)
        new_animal_id = f"{self.animal_id}_{hash_suffix}"

        if "animal" in self.result.columns:
            self.result["animal"] = new_animal_id
        if "animalday" in self.result.columns:
            self.result["animalday"] = self.result["animalday"].str.replace(self.animal_id, new_animal_id)
        self.animal_id = new_animal_id

        self.__update_instance_vars()

__init__(result, animal_id=None, genotype=None, channel_names=None, assume_from_number=False, bad_channels_dict={}, suppress_short_interval_error=False)

Parameters:

Name	Type	Description	Default
result	DataFrame	Result comes from AnimalOrganizer.compute_windowed_analysis()	required
animal_id	str	Identifier for the animal where result was computed from. Defaults to None.	None
genotype	str	Genotype of animal. Defaults to None.	None
channel_names	list[str]	List of channel names. Defaults to None.	None
assume_channels	bool	If true, assumes channel names according to AnimalFeatureParser.DEFAULT_CHNUM_TO_NAME. Defaults to False.	required
bad_channels_dict	dict[str, list[str]]	Dictionary of channels to reject for each recording session. Defaults to {}.	{}
suppress_short_interval_error	bool	If True, suppress ValueError for short intervals between timestamps. Useful for aggregated WARs with large window sizes. Defaults to False.	False

Source code in pythoneeg/visualization/results.py

def __init__(
    self,
    result: pd.DataFrame,
    animal_id: str = None,
    genotype: str = None,
    channel_names: list[str] = None,
    assume_from_number=False,
    bad_channels_dict: dict[str, list[str]] = {},
    suppress_short_interval_error=False,
) -> None:
    """
    Args:
        result (pd.DataFrame): Result comes from AnimalOrganizer.compute_windowed_analysis()
        animal_id (str, optional): Identifier for the animal where result was computed from. Defaults to None.
        genotype (str, optional): Genotype of animal. Defaults to None.
        channel_names (list[str], optional): List of channel names. Defaults to None.
        assume_channels (bool, optional): If true, assumes channel names according to AnimalFeatureParser.DEFAULT_CHNUM_TO_NAME. Defaults to False.
        bad_channels_dict (dict[str, list[str]], optional): Dictionary of channels to reject for each recording session. Defaults to {}.
        suppress_short_interval_error (bool, optional): If True, suppress ValueError for short intervals between timestamps. Useful for aggregated WARs with large window sizes. Defaults to False.
    """
    self.result = result
    self.animal_id = animal_id
    self.genotype = genotype
    self.channel_names = channel_names
    self.assume_from_number = assume_from_number
    self.bad_channels_dict = bad_channels_dict
    self.suppress_short_interval_error = suppress_short_interval_error

    self.__update_instance_vars()

    print(f"Channel names: \t{self.channel_names}")
    print(f"Channel abbreviations: \t{self.channel_abbrevs}")

__update_instance_vars()

Run after updating self.result, or other init values

Source code in pythoneeg/visualization/results.py

def __update_instance_vars(self):
    """Run after updating self.result, or other init values"""
    if "index" in self.result.columns:
        warnings.warn("Dropping column 'index'")
        self.result = self.result.drop(columns=["index"])

    # Check if timestamps are sorted and sort if needed
    if "timestamp" in self.result.columns:
        if not self.result["timestamp"].is_monotonic_increasing:
            warnings.warn("Timestamps are not sorted. Sorting result DataFrame by timestamp.")
            self.result = self.result.sort_values("timestamp")

    # Check for unusually short intervals between timestamps
    if "timestamp" in self.result.columns and "duration" in self.result.columns:
        median_duration = self.result["duration"].median()
        timestamp_diffs = self.result["timestamp"].diff()
        short_intervals = timestamp_diffs < pd.Timedelta(seconds=median_duration)

        # Skip first row since diff() produces NaT
        short_intervals = short_intervals[1:]

        if short_intervals.any():
            n_short = short_intervals.sum()
            pct_short = (n_short / len(short_intervals)) * 100

            warning_msg = (
                f"Found {n_short} intervals ({pct_short:.1f}%) between timestamps "
                f"that are shorter than the median duration of {median_duration:.1f}s"
            )

            if pct_short > 1.0 and not self.suppress_short_interval_error:  # More than 1% of intervals are short
                raise ValueError(warning_msg)
            elif not self.suppress_short_interval_error:
                warnings.warn(warning_msg)

    if "animal" in self.result.columns:
        unique_animals = self.result["animal"].unique()
        if len(unique_animals) > 1:
            raise ValueError(f"Multiple animals found in result: {unique_animals}")
        if unique_animals[0] != self.animal_id:
            raise ValueError(
                f"Animal ID mismatch: result has {unique_animals[0]}, but self.animal_id is {self.animal_id}"
            )

    self._feature_columns = [x for x in self.result.columns if x in constants.FEATURES]
    self._nonfeature_columns = [x for x in self.result.columns if x not in constants.FEATURES]
    self.animaldays = self.result.loc[:, "animalday"].unique()

    self.channel_abbrevs = [
        core.parse_chname_to_abbrev(x, assume_from_number=self.assume_from_number) for x in self.channel_names
    ]

add_unique_hash(nbytes=None)

Adds a hex hash to the animal ID to ensure uniqueness. This prevents collisions when, for example, multiple animals in ExperimentPlotter have the same animal ID.

Parameters:

Name	Type	Description	Default
nbytes	int	Number of bytes to generate. This is passed directly to secrets.token_hex(). Defaults to None, which generates 16 hex characters (8 bytes).	None

Source code in pythoneeg/visualization/results.py

def add_unique_hash(self, nbytes: int | None = None):
    """Adds a hex hash to the animal ID to ensure uniqueness. This prevents collisions when, for example, multiple animals in ExperimentPlotter have the same animal ID.

    Args:
        nbytes (int, optional): Number of bytes to generate. This is passed directly to secrets.token_hex(). Defaults to None, which generates 16 hex characters (8 bytes).
    """
    import secrets

    hash_suffix = secrets.token_hex(nbytes)
    new_animal_id = f"{self.animal_id}_{hash_suffix}"

    if "animal" in self.result.columns:
        self.result["animal"] = new_animal_id
    if "animalday" in self.result.columns:
        self.result["animalday"] = self.result["animalday"].str.replace(self.animal_id, new_animal_id)
    self.animal_id = new_animal_id

    self.__update_instance_vars()

aggregate_time_windows(groupby=['animalday', 'isday'])

Aggregate time windows into a single data point per groupby by averaging features. This reduces the number of rows in the result.

Parameters:

Name	Type	Description	Default
groupby	list[str] | str	Columns to group by. Defaults to ['animalday', 'isday'], which groups by animalday (recording session) and isday (day/night).	['animalday', 'isday']

Raises:

Type	Description
ValueError	groupby must be from ['animalday', 'isday']
ValueError	Columns in groupby not found in result
ValueError	Columns in groupby are not constant in groups

Source code in pythoneeg/visualization/results.py

def aggregate_time_windows(self, groupby: list[str] | str = ["animalday", "isday"]) -> None:
    """Aggregate time windows into a single data point per groupby by averaging features. This reduces the number of rows in the result.

    Args:
        groupby (list[str] | str, optional): Columns to group by. Defaults to ['animalday', 'isday'], which groups by animalday (recording session) and isday (day/night).

    Raises:
        ValueError: groupby must be from ['animalday', 'isday']
        ValueError: Columns in groupby not found in result
        ValueError: Columns in groupby are not constant in groups
    """
    if isinstance(groupby, str):
        groupby = [groupby]
    if not all(col in ["animalday", "isday"] for col in groupby):
        raise ValueError(f"groupby must be from ['animalday', 'isday']. Got {groupby}")
    if not all(col in self.result.columns for col in groupby):
        raise ValueError(f"Columns {groupby} not found in result. Columns: {self.result.columns.tolist()}")

    features = [f for f in constants.FEATURES if f in self.result.columns]
    logging.debug(f"Aggregating {features}")
    result_grouped = self.result.groupby(groupby)

    agg_dict = {}

    if "animalday" not in groupby:
        agg_dict["animalday"] = lambda df: None
    if "isday" not in groupby:
        agg_dict["isday"] = lambda df: None

    constant_cols = ["animal", "day", "genotype"]
    for col in constant_cols:
        if col in self.result.columns:
            is_constant = result_grouped[col].nunique() == 1
            if not is_constant.all():
                non_constant_groups = is_constant[~is_constant].index.tolist()
                raise ValueError(f"Column {col} is not constant in groups: {non_constant_groups}")
            agg_dict[col] = lambda df, col=col: df[col].iloc[0]

    if "duration" in self.result.columns:
        agg_dict["duration"] = lambda df: np.sum(df["duration"])

    if "endfile" in self.result.columns:
        agg_dict["endfile"] = lambda df: df["endfile"].iloc[-1]

    if "timestamp" in self.result.columns:
        agg_dict["timestamp"] = lambda df: df["timestamp"].iloc[0]

    for feat in features:
        agg_dict[feat] = lambda df, feat=feat: self._average_feature(df, feat, "duration")

    aggregated_df = result_grouped.apply(
        lambda df: pd.Series({col: agg_dict[col](df) for col in self.result.columns if col not in groupby})
    )

    self.result = aggregated_df.reset_index(drop=False)  # Keep animalday/isday as a column

    self.suppress_short_interval_error = True
    logging.info("Setting suppress_short_interval_error to True")
    self.__update_instance_vars()

apply_filters(filter_config=None, min_valid_channels=3, morphological_smoothing_seconds=None)

Apply multiple filters using configuration.

Parameters:

Name	Type	Description	Default
filter_config	dict	Dictionary of filter names and parameters. Available filters: 'logrms_range', 'high_rms', 'low_rms', 'high_beta', 'reject_channels', 'reject_channels_by_session', 'morphological_smoothing'	None
min_valid_channels	int	Minimum valid channels per window. Defaults to 3.	3
morphological_smoothing_seconds	float	Temporal smoothing window (deprecated, use config instead)	None

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Examples:

>>> config = {
...     'logrms_range': {'z_range': 3},
...     'high_rms': {'max_rms': 500},
...     'reject_channels': {'bad_channels': ['LMot', 'RMot']},
...     'morphological_smoothing': {'smoothing_seconds': 8.0}
... }
>>> filtered_war = war.apply_filters(config)

Source code in pythoneeg/visualization/results.py

def apply_filters(
    self, filter_config: dict = None, min_valid_channels: int = 3, morphological_smoothing_seconds: float = None
) -> "WindowAnalysisResult":
    """Apply multiple filters using configuration.

    Args:
        filter_config (dict, optional): Dictionary of filter names and parameters.
            Available filters: 'logrms_range', 'high_rms', 'low_rms', 'high_beta',
            'reject_channels', 'reject_channels_by_session', 'morphological_smoothing'
        min_valid_channels (int): Minimum valid channels per window. Defaults to 3.
        morphological_smoothing_seconds (float, optional): Temporal smoothing window (deprecated, use config instead)

    Returns:
        WindowAnalysisResult: New filtered instance

    Examples:
        >>> config = {
        ...     'logrms_range': {'z_range': 3},
        ...     'high_rms': {'max_rms': 500},
        ...     'reject_channels': {'bad_channels': ['LMot', 'RMot']},
        ...     'morphological_smoothing': {'smoothing_seconds': 8.0}
        ... }
        >>> filtered_war = war.apply_filters(config)
    """
    if filter_config is None:
        filter_config = {
            "logrms_range": {"z_range": 3},
            "high_rms": {"max_rms": 500},
            "low_rms": {"min_rms": 50},
            "high_beta": {"max_beta_prop": 0.4},
            "reject_channels_by_session": {},
        }

    filter_methods = {
        "logrms_range": self.get_filter_logrms_range,
        "high_rms": self.get_filter_high_rms,
        "low_rms": self.get_filter_low_rms,
        "high_beta": self.get_filter_high_beta,
        "reject_channels": self.get_filter_reject_channels,
        "reject_channels_by_session": self.get_filter_reject_channels_by_recording_session,
    }

    filt_bools = []
    morphological_params = None

    for filter_name, filter_params in filter_config.items():
        if filter_name == "morphological_smoothing":
            morphological_params = filter_params
            continue

        if filter_name not in filter_methods:
            raise ValueError(
                f"Unknown filter: {filter_name}. Available: {list(filter_methods.keys()) + ['morphological_smoothing']}"
            )

        filter_func = filter_methods[filter_name]
        filt_bool = filter_func(**filter_params)
        filt_bools.append(filt_bool)
        logging.info(f"{filter_name}: filtered {filt_bool.size - np.count_nonzero(filt_bool)}/{filt_bool.size}")

    # Combine all filter masks
    if filt_bools:
        filt_bool_all = np.prod(np.stack(filt_bools, axis=-1), axis=-1).astype(bool)
    else:
        filt_bool_all = np.ones((len(self.result), len(self.channel_names)), dtype=bool)

    # Apply morphological smoothing if requested (either from config or parameter)
    if morphological_params or morphological_smoothing_seconds is not None:
        if morphological_params:
            smoothing_seconds = morphological_params["smoothing_seconds"]
        else:
            smoothing_seconds = morphological_smoothing_seconds

        filt_bool_all = self.get_filter_morphological_smoothing(filt_bool_all, smoothing_seconds)
        logging.info(f"Applied morphological smoothing: {smoothing_seconds}s")

    # Filter windows based on minimum valid channels
    valid_channels_per_window = np.sum(filt_bool_all, axis=1)
    window_mask = valid_channels_per_window >= min_valid_channels
    filt_bool_all = filt_bool_all & window_mask[:, np.newaxis]

    return self._create_filtered_copy(filt_bool_all)

filter_all(df=None, inplace=True, bad_channels=None, min_valid_channels=3, filters=None, morphological_smoothing_seconds=None, **kwargs)

Apply a list of filters to the data. Filtering should be performed before aggregation.

Parameters:

Name	Type	Description	Default
df	DataFrame	If not None, this function will use this dataframe instead of self.result. Defaults to None.	None
inplace	bool	If True, modify the result in place. Defaults to True.	True
bad_channels	list[str]	List of channels to reject. Defaults to None.	None
min_valid_channels	int	Minimum number of valid channels required per window. Defaults to 3.	3
filters	list[callable]	List of filter functions to apply. Each function should return a boolean mask. If None, uses default filters: [get_filter_logrms_range, get_filter_high_rms, get_filter_low_rms, get_filter_high_beta]. Defaults to None.	None
morphological_smoothing_seconds	float	If provided, apply morphological opening/closing to smooth the filter mask. This removes isolated false positives/negatives along the time axis for each channel independently. The value specifies the time window in seconds for the morphological operations. Defaults to None.	None
**kwargs		Additional keyword arguments to pass to filter functions.	{}

Returns:

Name	Type	Description
WindowAnalysisResult		Filtered result

Source code in pythoneeg/visualization/results.py

def filter_all(
    self,
    df: pd.DataFrame = None,
    inplace=True,
    bad_channels: list[str] = None,
    min_valid_channels=3,
    filters: list[callable] = None,
    morphological_smoothing_seconds: float = None,
    **kwargs,
):
    """Apply a list of filters to the data. Filtering should be performed before aggregation.

    Args:
        df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
        inplace (bool, optional): If True, modify the result in place. Defaults to True.
        bad_channels (list[str], optional): List of channels to reject. Defaults to None.
        min_valid_channels (int, optional): Minimum number of valid channels required per window. Defaults to 3.
        filters (list[callable], optional): List of filter functions to apply. Each function should return a boolean mask.
            If None, uses default filters: [get_filter_logrms_range, get_filter_high_rms, get_filter_low_rms, get_filter_high_beta].
            Defaults to None.
        morphological_smoothing_seconds (float, optional): If provided, apply morphological opening/closing to smooth the filter mask.
            This removes isolated false positives/negatives along the time axis for each channel independently.
            The value specifies the time window in seconds for the morphological operations. Defaults to None.
        **kwargs: Additional keyword arguments to pass to filter functions.

    Returns:
        WindowAnalysisResult: Filtered result
    """
    if filters is None:
        # TODO refactor these into standalone functions, which take in a war as the first parameter, then pass
        # filt_bool = filt(self, df, **kwargs) as needed
        filters = [
            self.get_filter_logrms_range,
            self.get_filter_high_rms,
            self.get_filter_low_rms,
            self.get_filter_high_beta,
            self.get_filter_reject_channels_by_recording_session,
        ]

    filt_bools = []
    # Apply each filter function
    for filter_function in filters:
        filt_bool = filter_function(df, **kwargs)
        filt_bools.append(filt_bool)
        logging.info(
            f"{filter_function.__name__}:\tfiltered {filt_bool.size - np.count_nonzero(filt_bool)}/{filt_bool.size}"
        )

    # Filter channels manually
    # REVIEW somehow add this to the main list of filters, but I'm not sure how to do this.
    if bad_channels is not None:
        filt_bools.append(self.get_filter_reject_channels(bad_channels))
        logging.debug(f"Reject channels: {filt_bools[-1]}")

    # Apply all filters
    filt_bool_all = np.prod(np.stack(filt_bools, axis=-1), axis=-1).astype(bool)
    logging.debug(f"filt_bool_all.shape: {filt_bool_all.shape}")  # (windows, channels)

    # Apply morphological smoothing if requested
    if morphological_smoothing_seconds is not None:
        if "duration" not in self.result.columns:
            raise ValueError("Cannot calculate window duration - 'duration' column missing from result dataframe")
        window_duration = self.result["duration"].median()

        # Calculate number of windows for the smoothing
        structure_size = max(1, int(morphological_smoothing_seconds / window_duration))

        if structure_size > 1:
            logging.info(
                f"Applying morphological smoothing with {structure_size} windows ({morphological_smoothing_seconds}s / {window_duration}s per window)"
            )
            # Apply channel-wise temporal smoothing (each channel processed independently)
            # This avoids spatial assumptions while smoothing temporal artifacts
            for ch_idx in range(filt_bool_all.shape[1]):
                channel_mask = filt_bool_all[:, ch_idx]
                # Opening removes small isolated artifacts
                channel_mask = binary_opening(channel_mask, structure=np.ones(structure_size))
                # Closing fills small gaps in valid data
                channel_mask = binary_closing(channel_mask, structure=np.ones(structure_size))
                filt_bool_all[:, ch_idx] = channel_mask
        else:
            logging.info("Skipping morphological smoothing - structure size would be 1 (no effect)")

    # Filter windows based on number of valid channels
    valid_channels_per_window = np.sum(filt_bool_all, axis=1)  # axis 1 = channel
    window_mask = valid_channels_per_window >= min_valid_channels  # True if window has enough valid channels
    filt_bool_all = filt_bool_all & window_mask[:, np.newaxis]  # Apply window mask to all channels

    filtered_result = self._apply_filter(filt_bool_all)
    if inplace:
        del self.result
        self.result = filtered_result
    return WindowAnalysisResult(
        filtered_result,
        self.animal_id,
        self.genotype,
        self.channel_names,
        self.assume_from_number,
        self.bad_channels_dict,
        self.suppress_short_interval_error,
    )

filter_high_beta(max_beta_prop=0.4)

Filter out windows with high beta power.

Parameters:

Name	Type	Description	Default
max_beta_prop	float	Maximum beta power proportion. Defaults to 0.4.	0.4

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_high_beta(self, max_beta_prop: float = 0.4) -> "WindowAnalysisResult":
    """Filter out windows with high beta power.

    Args:
        max_beta_prop (float): Maximum beta power proportion. Defaults to 0.4.

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    mask = self.get_filter_high_beta(max_beta_prop=max_beta_prop)
    return self._create_filtered_copy(mask)

filter_high_rms(max_rms=500)

Filter out windows with RMS above threshold.

Parameters:

Name	Type	Description	Default
max_rms	float	Maximum RMS threshold. Defaults to 500.	500

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_high_rms(self, max_rms: float = 500) -> "WindowAnalysisResult":
    """Filter out windows with RMS above threshold.

    Args:
        max_rms (float): Maximum RMS threshold. Defaults to 500.

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    mask = self.get_filter_high_rms(max_rms=max_rms)
    return self._create_filtered_copy(mask)

filter_logrms_range(z_range=3)

Filter based on log(rms) z-score range.

Parameters:

Name	Type	Description	Default
z_range	float	Z-score range threshold. Defaults to 3.	3

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_logrms_range(self, z_range: float = 3) -> "WindowAnalysisResult":
    """Filter based on log(rms) z-score range.

    Args:
        z_range (float): Z-score range threshold. Defaults to 3.

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    mask = self.get_filter_logrms_range(z_range=z_range)
    return self._create_filtered_copy(mask)

filter_low_rms(min_rms=50)

Filter out windows with RMS below threshold.

Parameters:

Name	Type	Description	Default
min_rms	float	Minimum RMS threshold. Defaults to 50.	50

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_low_rms(self, min_rms: float = 50) -> "WindowAnalysisResult":
    """Filter out windows with RMS below threshold.

    Args:
        min_rms (float): Minimum RMS threshold. Defaults to 50.

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    mask = self.get_filter_low_rms(min_rms=min_rms)
    return self._create_filtered_copy(mask)

filter_morphological_smoothing(smoothing_seconds)

Apply morphological smoothing to all data.

Parameters:

Name	Type	Description	Default
smoothing_seconds	float	Time window in seconds for morphological operations	required

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_morphological_smoothing(self, smoothing_seconds: float) -> "WindowAnalysisResult":
    """Apply morphological smoothing to all data.

    Args:
        smoothing_seconds (float): Time window in seconds for morphological operations

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    # Start with all-True mask and smooth it
    base_mask = np.ones((len(self.result), len(self.channel_names)), dtype=bool)
    smoothed_mask = self.get_filter_morphological_smoothing(base_mask, smoothing_seconds)
    return self._create_filtered_copy(smoothed_mask)

filter_reject_channels(bad_channels, use_abbrevs=None)

Filter out specified bad channels.

Parameters:

Name	Type	Description	Default
bad_channels	list[str]	List of channel names to reject	required
use_abbrevs	bool	Whether to use abbreviations. Defaults to None.	None

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_reject_channels(self, bad_channels: list[str], use_abbrevs: bool = None) -> "WindowAnalysisResult":
    """Filter out specified bad channels.

    Args:
        bad_channels (list[str]): List of channel names to reject
        use_abbrevs (bool, optional): Whether to use abbreviations. Defaults to None.

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    mask = self.get_filter_reject_channels(bad_channels, use_abbrevs=use_abbrevs)
    return self._create_filtered_copy(mask)

filter_reject_channels_by_session(bad_channels_dict=None, use_abbrevs=None)

Filter out bad channels by recording session.

Parameters:

Name	Type	Description	Default
bad_channels_dict	dict[str, list[str]]	Dict of bad channels per session	None
use_abbrevs	bool	Whether to use abbreviations. Defaults to None.	None

Returns:

Name	Type	Description
WindowAnalysisResult	WindowAnalysisResult	New filtered instance

Source code in pythoneeg/visualization/results.py

def filter_reject_channels_by_session(
    self, bad_channels_dict: dict[str, list[str]] = None, use_abbrevs: bool = None
) -> "WindowAnalysisResult":
    """Filter out bad channels by recording session.

    Args:
        bad_channels_dict (dict[str, list[str]], optional): Dict of bad channels per session
        use_abbrevs (bool, optional): Whether to use abbreviations. Defaults to None.

    Returns:
        WindowAnalysisResult: New filtered instance
    """
    mask = self.get_filter_reject_channels_by_recording_session(bad_channels_dict, use_abbrevs=use_abbrevs)
    return self._create_filtered_copy(mask)

get_filter_high_beta(df=None, max_beta_prop=0.4, **kwargs)

Filter windows based on beta power.

Parameters:

Name	Type	Description	Default
df	DataFrame	If not None, this function will use this dataframe instead of self.result. Defaults to None.	None
max_beta_prop	float	The maximum beta power to filter by. Values above this will be set to NaN. Defaults to 0.4.	0.4

Returns:

Name	Type	Description
out		np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window

Source code in pythoneeg/visualization/results.py

def get_filter_high_beta(self, df: pd.DataFrame = None, max_beta_prop=0.4, **kwargs):
    """Filter windows based on beta power.

    Args:
        df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
        max_beta_prop (float, optional): The maximum beta power to filter by. Values above this will be set to NaN. Defaults to 0.4.

    Returns:
        out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
    """
    result = df.copy() if df is not None else self.result.copy()
    if "psdfrac" in result.columns:
        df_psdfrac = pd.DataFrame(result["psdfrac"].tolist())
        np_prop = np.array(df_psdfrac["beta"].tolist())
    elif "psdband" in result.columns and "psdtotal" in result.columns:
        df_psdband = pd.DataFrame(result["psdband"].tolist())
        np_beta = np.array(df_psdband["beta"].tolist())
        np_total = np.array(result["psdtotal"].tolist())
        np_prop = np_beta / np_total
    else:
        raise ValueError("psdfrac or psdband+psdtotal required for beta power filtering")

    out = np.full(np_prop.shape, True)
    out[np_prop > max_beta_prop] = False
    out = np.broadcast_to(np.all(out, axis=-1)[:, np.newaxis], out.shape)
    return out

get_filter_high_rms(df=None, max_rms=500, **kwargs)

Filter windows based on rms.

Parameters:

Name	Type	Description	Default
df	DataFrame	If not None, this function will use this dataframe instead of self.result. Defaults to None.	None
max_rms	float	The maximum rms value to filter by. Values above this will be set to NaN.	500

Returns:

Name	Type	Description
out		np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window

Source code in pythoneeg/visualization/results.py

def get_filter_high_rms(self, df: pd.DataFrame = None, max_rms=500, **kwargs):
    """Filter windows based on rms.

    Args:
        df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
        max_rms (float, optional): The maximum rms value to filter by. Values above this will be set to NaN.

    Returns:
        out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
    """
    result = df.copy() if df is not None else self.result.copy()
    np_rms = np.array(result["rms"].tolist())
    np_rmsnan = np_rms.copy()
    # Convert to float to allow NaN assignment for integer arrays
    if np_rmsnan.dtype.kind in ("i", "u"):  # integer types
        np_rmsnan = np_rmsnan.astype(float)
    np_rmsnan[np_rms > max_rms] = np.nan
    result["rms"] = np_rmsnan.tolist()

    out = np.full(np_rms.shape, True)
    out[np_rms > max_rms] = False
    return out

get_filter_logrms_range(df=None, z_range=3, **kwargs)

Filter windows based on log(rms).

Parameters:

Name	Type	Description	Default
df	DataFrame	If not None, this function will use this dataframe instead of self.result. Defaults to None.	None
z_range	float	The z-score range to filter by. Values outside this range will be set to NaN.	3

Returns:

Name	Type	Description
out		np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window

Source code in pythoneeg/visualization/results.py

def get_filter_logrms_range(self, df: pd.DataFrame = None, z_range=3, **kwargs):
    """Filter windows based on log(rms).

    Args:
        df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
        z_range (float, optional): The z-score range to filter by. Values outside this range will be set to NaN.

    Returns:
        out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
    """
    result = df.copy() if df is not None else self.result.copy()
    z_range = abs(z_range)
    np_rms = np.array(result["rms"].tolist())
    np_logrms = np.log(np_rms)
    del np_rms
    np_logrmsz = zscore(np_logrms, axis=0, nan_policy="omit")
    np_logrms[(np_logrmsz > z_range) | (np_logrmsz < -z_range)] = np.nan

    out = np.full(np_logrms.shape, True)
    out[(np_logrmsz > z_range) | (np_logrmsz < -z_range)] = False
    return out

get_filter_low_rms(df=None, min_rms=30, **kwargs)

Filter windows based on rms.

Parameters:

Name	Type	Description	Default
df	DataFrame	If not None, this function will use this dataframe instead of self.result. Defaults to None.	None
min_rms	float	The minimum rms value to filter by. Values below this will be set to NaN.	30

Returns:

Name	Type	Description
out		np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window

Source code in pythoneeg/visualization/results.py

def get_filter_low_rms(self, df: pd.DataFrame = None, min_rms=30, **kwargs):
    """Filter windows based on rms.

    Args:
        df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
        min_rms (float, optional): The minimum rms value to filter by. Values below this will be set to NaN.

    Returns:
        out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
    """
    result = df.copy() if df is not None else self.result.copy()
    np_rms = np.array(result["rms"].tolist())
    np_rmsnan = np_rms.copy()
    np_rmsnan[np_rms < min_rms] = np.nan
    result["rms"] = np_rmsnan.tolist()

    out = np.full(np_rms.shape, True)
    out[np_rms < min_rms] = False
    return out

get_filter_morphological_smoothing(filter_mask, smoothing_seconds, **kwargs)

Apply morphological smoothing to a filter mask.

Parameters:

Name	Type	Description	Default
filter_mask	ndarray	Input boolean mask of shape (n_windows, n_channels)	required
smoothing_seconds	float	Time window in seconds for morphological operations	required

Returns:

Type	Description
ndarray	np.ndarray: Smoothed boolean mask

Source code in pythoneeg/visualization/results.py

def get_filter_morphological_smoothing(
    self, filter_mask: np.ndarray, smoothing_seconds: float, **kwargs
) -> np.ndarray:
    """Apply morphological smoothing to a filter mask.

    Args:
        filter_mask (np.ndarray): Input boolean mask of shape (n_windows, n_channels)
        smoothing_seconds (float): Time window in seconds for morphological operations

    Returns:
        np.ndarray: Smoothed boolean mask
    """
    if "duration" not in self.result.columns:
        raise ValueError("Cannot calculate window duration - 'duration' column missing")

    window_duration = self.result["duration"].median()
    structure_size = max(1, int(smoothing_seconds / window_duration))

    if structure_size <= 1:
        return filter_mask

    smoothed_mask = filter_mask.copy()
    for ch_idx in range(filter_mask.shape[1]):
        channel_mask = filter_mask[:, ch_idx]
        # Opening removes small isolated artifacts
        channel_mask = binary_opening(channel_mask, structure=np.ones(structure_size))
        # Closing fills small gaps in valid data
        channel_mask = binary_closing(channel_mask, structure=np.ones(structure_size))
        smoothed_mask[:, ch_idx] = channel_mask

    return smoothed_mask

get_filter_reject_channels(bad_channels, use_abbrevs=None, **kwargs)

Filter channels to reject.

Parameters:

Name	Type	Description	Default
bad_channels	list[str]	List of channels to reject. Can be either full channel names or abbreviations. The method will automatically detect which format is being used.	required
use_abbrevs	bool	Override automatic detection. If True, channels are assumed to be channel abbreviations. If False, channels are assumed to be channel names. If None, channels are parsed to abbreviations and matched against self.channel_abbrevs.	None

Returns:

Name	Type	Description
out		np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window

Source code in pythoneeg/visualization/results.py

def get_filter_reject_channels(self, bad_channels: list[str], use_abbrevs: bool = None, **kwargs):
    """Filter channels to reject.

    Args:
        bad_channels (list[str]): List of channels to reject. Can be either full channel names or abbreviations.
            The method will automatically detect which format is being used.
        use_abbrevs (bool, optional): Override automatic detection. If True, channels are assumed to be channel abbreviations. If False, channels are assumed to be channel names.
            If None, channels are parsed to abbreviations and matched against self.channel_abbrevs.

    Returns:
        out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
    """
    channel_targets = (
        self.channel_abbrevs if use_abbrevs or use_abbrevs is None else self.channel_names
    )  # Match to appropriate target
    if use_abbrevs is None:  # Match channels as abbreviations
        bad_channels = [
            core.utils.parse_chname_to_abbrev(ch, assume_from_number=self.assume_from_number) for ch in bad_channels
        ]

    n_samples = len(self.result)
    n_channels = len(channel_targets)
    mask = np.ones((n_samples, n_channels), dtype=bool)

    # Match channels to channel_targets
    for ch in bad_channels:
        if ch in channel_targets:
            mask[:, channel_targets.index(ch)] = False
        else:
            warnings.warn(f"Channel {ch} not found in {channel_targets}")
    return mask

get_filter_reject_channels_by_recording_session(bad_channels_dict=None, use_abbrevs=None)

Filter channels to reject for each recording session

Parameters:

Name	Type	Description	Default
bad_channels_dict	dict[str, list[str]]	Dictionary of list of channels to reject for each recording session. Can be either full channel names or abbreviations. The method will automatically detect which format is being used. If None, the method will use the bad_channels_dict passed to the constructor.	None
use_abbrevs	bool	Override automatic detection. If True, channels are assumed to be channel abbreviations. If False, channels are assumed to be channel names. If None, channels are parsed to abbreviations and matched against self.channel_abbrevs.	None

Returns:

Name	Type	Description
out		np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window

Source code in pythoneeg/visualization/results.py

def get_filter_reject_channels_by_recording_session(
    self, bad_channels_dict: dict[str, list[str]] = None, use_abbrevs: bool = None
):
    """Filter channels to reject for each recording session

    Args:
        bad_channels_dict (dict[str, list[str]]): Dictionary of list of channels to reject for each recording session.
            Can be either full channel names or abbreviations. The method will automatically detect which format is being used.
            If None, the method will use the bad_channels_dict passed to the constructor.
        use_abbrevs (bool, optional): Override automatic detection. If True, channels are assumed to be channel abbreviations. If False, channels are assumed to be channel names.
            If None, channels are parsed to abbreviations and matched against self.channel_abbrevs.

    Returns:
        out: np.ndarray of bool, (M fragments, N channels). True = keep window, False = remove window
    """
    if bad_channels_dict is None:
        bad_channels_dict = self.bad_channels_dict

    n_samples = len(self.result)
    n_channels = len(self.channel_names)
    mask = np.ones((n_samples, n_channels), dtype=bool)

    # Group by animalday to apply filters per recording session
    for animalday, group in self.result.groupby("animalday"):
        if animalday not in bad_channels_dict:
            logging.warning(f"No bad channels specified for recording session {animalday}")
            continue

        bad_channels = bad_channels_dict[animalday]
        channel_targets = self.channel_abbrevs if use_abbrevs or use_abbrevs is None else self.channel_names
        if use_abbrevs is None:
            bad_channels = [
                core.parse_chname_to_abbrev(ch, assume_from_number=self.assume_from_number) for ch in bad_channels
            ]

        # Get indices for this recording session
        session_indices = group.index

        # Apply channel filtering for this session
        for ch in bad_channels:
            if ch in channel_targets:
                ch_idx = channel_targets.index(ch)
                mask[session_indices, ch_idx] = False
            else:
                logging.warning(f"Channel {ch} not found in {channel_targets} for session {animalday}")

    return mask

get_groupavg_result(features, exclude=[], df=None, groupby='animalday')

Group result and average within groups. Preserves data structure and shape for each feature.

Parameters:

Name	Type	Description	Default
features	list[str]	List of features to get from result	required
exclude	list[str]	List of features to exclude from result. Will override the features parameter. Defaults to [].	[]
df	DataFrame	If not None, this function will use this dataframe instead of self.result. Defaults to None.	None
groupby	str	Feature or list of features to group by before averaging. Passed to the by parameter in pd.DataFrame.groupby(). Defaults to "animalday".	'animalday'

Returns:

Name	Type	Description
grouped_result		result grouped by groupby and averaged for each group.

Source code in pythoneeg/visualization/results.py

def get_groupavg_result(
    self, features: list[str], exclude: list[str] = [], df: pd.DataFrame = None, groupby="animalday"
):
    """Group result and average within groups. Preserves data structure and shape for each feature.

    Args:
        features (list[str]): List of features to get from result
        exclude (list[str], optional): List of features to exclude from result. Will override the features parameter. Defaults to [].
        df (pd.DataFrame, optional): If not None, this function will use this dataframe instead of self.result. Defaults to None.
        groupby (str, optional): Feature or list of features to group by before averaging. Passed to the `by` parameter in pd.DataFrame.groupby(). Defaults to "animalday".

    Returns:
        grouped_result: result grouped by `groupby` and averaged for each group.
    """
    result_grouped, result_validcols = self.__get_groups(features=features, exclude=exclude, df=df, groupby=groupby)
    features = _sanitize_feature_request(features, exclude)

    avg_results = []
    for f in features:
        if f in result_validcols:
            avg_result_col = result_grouped.apply(self._average_feature, f, "duration", include_groups=False)
            avg_result_col.name = f
            avg_results.append(avg_result_col)
        else:
            logging.warning(f"{f} not calculated, skipping")

    return pd.concat(avg_results, axis=1)

get_info()

Returns a formatted string with basic information about the WindowAnalysisResult object

Source code in pythoneeg/visualization/results.py

def get_info(self):
    """Returns a formatted string with basic information about the WindowAnalysisResult object"""
    info = []
    info.append(f"feature names: {', '.join(self._feature_columns)}")
    info.append(f"animaldays: {', '.join(self.result['animalday'].unique())}")
    info.append(
        f"animal_id: {self.result['animal'].unique()[0] if 'animal' in self.result.columns else self.animal_id}"
    )
    info.append(
        f"genotype: {self.result['genotype'].unique()[0] if 'genotype' in self.result.columns else self.genotype}"
    )
    info.append(f"channel_names: {', '.join(self.channel_names) if self.channel_names else 'None'}")

    return "\n".join(info)

get_result(features, exclude=[], allow_missing=False)

Get windowed analysis result dataframe, with helpful filters

Parameters:

Name	Type	Description	Default
features	list[str]	List of features to get from result	required
exclude	list[str]	List of features to exclude from result; will override the features parameter. Defaults to [].	[]
allow_missing	bool	If True, will return all requested features as columns regardless if they exist in result. Defaults to False.	False

Returns:

Name	Type	Description
result		pd.DataFrame object with features in columns and windows in rows

Source code in pythoneeg/visualization/results.py

def get_result(self, features: list[str], exclude: list[str] = [], allow_missing=False):
    """Get windowed analysis result dataframe, with helpful filters

    Args:
        features (list[str]): List of features to get from result
        exclude (list[str], optional): List of features to exclude from result; will override the features parameter. Defaults to [].
        allow_missing (bool, optional): If True, will return all requested features as columns regardless if they exist in result. Defaults to False.

    Returns:
        result: pd.DataFrame object with features in columns and windows in rows
    """
    features = _sanitize_feature_request(features, exclude)
    if not allow_missing:
        return self.result.loc[:, self._nonfeature_columns + features]
    else:
        return self.result.reindex(columns=self._nonfeature_columns + features)

load_pickle_and_json(folder_path=None) classmethod

Load WindowAnalysisResult from folder

Parameters:

Name	Type	Description	Default
folder_path	str	Path of folder containing one .pkl and .json file each. Defaults to None.	None
df_pickle_path	str	Path of .pkl file. If this and folder_path are not None, raises an error. Defaults to None.	required
json_path	str	Path of .json file. If this and folder_path are not None, raises an error. Defaults to None.	required

Raises:

Type	Description
ValueError	Both df_pickle_path and json_path must be None if folder_path is provided
ValueError	Expected exactly one pickle and one json file in folder_path

Returns:

Name	Type	Description
result		WindowAnalysisResult object

Source code in pythoneeg/visualization/results.py

@classmethod
def load_pickle_and_json(cls, folder_path=None):
    """Load WindowAnalysisResult from folder

    Args:
        folder_path (str, optional): Path of folder containing one .pkl and .json file each. Defaults to None.
        df_pickle_path (str, optional): Path of .pkl file. If this and folder_path are not None, raises an error. Defaults to None.
        json_path (str, optional): Path of .json file. If this and folder_path are not None, raises an error. Defaults to None.

    Raises:
        ValueError: Both df_pickle_path and json_path must be None if folder_path is provided
        ValueError: Expected exactly one pickle and one json file in folder_path

    Returns:
        result: WindowAnalysisResult object
    """
    if folder_path is not None:
        folder_path = Path(folder_path)
        if not folder_path.exists():
            raise ValueError(f"Folder path {folder_path} does not exist")

        pkl_files = list(folder_path.glob("*.pkl"))
        json_files = list(folder_path.glob("*.json"))

        if len(pkl_files) != 1 or len(json_files) != 1:
            raise ValueError(f"Expected exactly one pickle and one json file in {folder_path}")

        df_pickle_path = pkl_files[0]
        json_path = json_files[0]

    with open(df_pickle_path, "rb") as f:
        data = pd.read_pickle(f)
    with open(json_path, "r") as f:
        metadata = json.load(f)
    return cls(data, **metadata)

reorder_and_pad_channels(target_channels, use_abbrevs=True, inplace=True)

Reorder and pad channels to match a target channel list.

This method ensures that the data has a consistent channel order and structure by reordering existing channels and padding missing channels with NaNs.

Parameters:

Name	Type	Description	Default
target_channels	list[str]	List of target channel names to match	required
use_abbrevs	bool	If True, target channel names are read as channel abbreviations instead of channel names. Defaults to True.	True
inplace	bool	If True, modify the result in place. Defaults to True.	True

Returns: pd.DataFrame: DataFrame with reordered and padded channels

Source code in pythoneeg/visualization/results.py

def reorder_and_pad_channels(
    self, target_channels: list[str], use_abbrevs: bool = True, inplace: bool = True
) -> pd.DataFrame:
    """Reorder and pad channels to match a target channel list.

    This method ensures that the data has a consistent channel order and structure
    by reordering existing channels and padding missing channels with NaNs.

    Args:
        target_channels (list[str]): List of target channel names to match
        use_abbrevs (bool, optional): If True, target channel names are read as channel abbreviations instead of channel names. Defaults to True.
        inplace (bool, optional): If True, modify the result in place. Defaults to True.
    Returns:
        pd.DataFrame: DataFrame with reordered and padded channels
    """
    duplicates = [ch for ch in target_channels if target_channels.count(ch) > 1]
    if duplicates:
        raise ValueError(f"Target channels must be unique. Found duplicates: {duplicates}")

    result = self.result.copy()

    channel_map = {ch: i for i, ch in enumerate(target_channels)}
    channel_names = self.channel_names if not use_abbrevs else self.channel_abbrevs

    valid_channels = [ch for ch in channel_names if ch in channel_map]
    if not valid_channels:
        warnings.warn(
            f"None of the channel names {channel_names} were found in target channels {target_channels}. Is use_abbrevs correctly set?"
        )

    for feature in self._feature_columns:
        match feature:
            case _ if feature in constants.LINEAR_FEATURES + constants.BAND_FEATURES:
                if feature in constants.BAND_FEATURES:
                    df_bands = pd.DataFrame(result[feature].tolist())
                    vals = np.array(df_bands.values.tolist())
                    vals = vals.transpose((0, 2, 1))
                    keys = df_bands.keys()
                else:
                    vals = np.array(result[feature].tolist())

                new_vals = np.full((vals.shape[0], len(target_channels), *vals.shape[2:]), np.nan)  # dubious

                for i, ch in enumerate(channel_names):
                    if ch in channel_map:
                        new_vals[:, channel_map[ch]] = vals[:, i]

                if feature in constants.BAND_FEATURES:
                    new_vals = new_vals.transpose((0, 2, 1))
                    result[feature] = [dict(zip(keys, vals)) for vals in new_vals]
                else:
                    result[feature] = [list(x) for x in new_vals]

            case _ if feature in constants.MATRIX_FEATURES:
                if feature in ["cohere", "zcohere", "imcoh", "zimcoh"]:
                    df_bands = pd.DataFrame(result[feature].tolist())
                    vals = np.array(df_bands.values.tolist())
                    keys = df_bands.keys()
                else:
                    vals = np.array(result[feature].tolist())

                logging.debug(f"vals.shape: {vals.shape}")
                new_shape = list(vals.shape[:-2]) + [len(target_channels), len(target_channels)]
                new_vals = np.full(new_shape, np.nan)

                ch1_valid = np.array([ch in channel_map for ch in channel_names])
                ch2_valid = ch1_valid.copy()
                valid_pairs = np.logical_and(ch1_valid[:, None], ch2_valid[None, :])  # 2D boolean mask

                for i, j in zip(*np.where(valid_pairs)):
                    ch1, ch2 = channel_names[i], channel_names[j]
                    new_vals[..., channel_map[ch1], channel_map[ch2]] = vals[..., i, j]

                triu_mask = np.triu_indices(len(target_channels), k=0)
                new_vals += new_vals.transpose((*range(new_vals.ndim - 2), -1, -2))
                new_vals[..., triu_mask[0], triu_mask[1]] = 0

                if feature in ["cohere", "zcohere", "imcoh", "zimcoh"]:
                    result[feature] = [dict(zip(keys, vals)) for vals in new_vals]
                else:
                    result[feature] = [list(x) for x in new_vals]

            case _ if feature in constants.HIST_FEATURES:
                coords = np.array([x[0] for x in result[feature].tolist()])
                vals = np.array([x[1] for x in result[feature].tolist()])
                new_vals = np.full((*vals.shape[0:-1], len(target_channels)), np.nan)

                for i, ch in enumerate(channel_names):
                    if ch in channel_map:
                        new_vals[:, ..., channel_map[ch]] = vals[:, ..., i]

                result[feature] = [(coords[i], new_vals[i]) for i in range(len(coords))]

            case _:
                raise ValueError(f"Invalid feature: {feature}")

    if inplace:
        self.result = result

        logging.debug(f"Old channel names: {self.channel_names}")
        self.channel_names = target_channels
        logging.debug(f"New channel names: {self.channel_names}")

        logging.debug(f"Old channel abbreviations: {self.channel_abbrevs}")
        self.__update_instance_vars()
        logging.debug(f"New channel abbreviations: {self.channel_abbrevs}")

    return result

save_pickle_and_json(folder, make_folder=True, slugify_filebase=True, save_abbrevs_as_chnames=False)

Archive window analysis result into the folder specified, as a pickle and json file.

Parameters:

Name	Type	Description	Default
folder	str | Path	Destination folder to save results to	required
make_folder	bool	If True, create the folder if it doesn't exist. Defaults to True.	True
slugify_filebase	bool	If True, slugify the filebase (replace special characters). Defaults to True.	True
save_abbrevs_as_chnames	bool	If True, save the channel abbreviations as the channel names in the json file. Defaults to False.	False

Source code in pythoneeg/visualization/results.py

def save_pickle_and_json(
    self, folder: str | Path, make_folder=True, slugify_filebase=True, save_abbrevs_as_chnames=False
):
    """Archive window analysis result into the folder specified, as a pickle and json file.

    Args:
        folder (str | Path): Destination folder to save results to
        make_folder (bool, optional): If True, create the folder if it doesn't exist. Defaults to True.
        slugify_filebase (bool, optional): If True, slugify the filebase (replace special characters). Defaults to True.
        save_abbrevs_as_chnames (bool, optional): If True, save the channel abbreviations as the channel names in the json file. Defaults to False.
    """
    folder = Path(folder)
    if make_folder:
        folder.mkdir(parents=True, exist_ok=True)

    if slugify_filebase:
        filebase = folder / slugify(f"{self.animal_id}-{self.genotype}")
    else:
        filebase = folder / f"{self.animal_id}-{self.genotype}"
    filebase = str(filebase)

    self.result.to_pickle(filebase + ".pkl")
    json_dict = {
        "animal_id": self.animal_id,
        "genotype": self.genotype,
        "channel_names": self.channel_abbrevs if save_abbrevs_as_chnames else self.channel_names,
        "assume_from_number": False if save_abbrevs_as_chnames else self.assume_from_number,
        "bad_channels_dict": self.bad_channels_dict,
        "suppress_short_interval_error": self.suppress_short_interval_error,
    }

    with open(filebase + ".json", "w") as f:
        json.dump(json_dict, f, indent=2)

bin_spike_times(spike_times, fragment_durations)

Bin spike times into counts based on fragment durations.

Parameters:

Name	Type	Description	Default
spike_times	list[float]	List of spike timestamps in seconds	required
fragment_durations	list[float]	List of fragment durations in seconds	required

Returns:

Type	Description
list[int]	list[int]: List of spike counts per fragment

Source code in pythoneeg/visualization/results.py

def bin_spike_times(spike_times: list[float], fragment_durations: list[float]) -> list[int]:
    """Bin spike times into counts based on fragment durations.

    Args:
        spike_times (list[float]): List of spike timestamps in seconds
        fragment_durations (list[float]): List of fragment durations in seconds

    Returns:
        list[int]: List of spike counts per fragment
    """
    # Convert fragment durations to bin edges
    bin_edges = np.cumsum([0] + fragment_durations)

    # Use numpy's histogram function to count spikes in each bin
    counts, _ = np.histogram(spike_times, bins=bin_edges)

    return counts.tolist()