API Reference

This page documents all public Pyxations modules and functions. Docstrings are rendered automatically from the source.

Top-level entry points

The most common entry points are re-exported from the package root:

from pyxations import (
    dataset_to_bids,
    compute_derivatives_for_dataset,
    PreProcessing,
    RemodnavDetection,
    EngbertDetection,
    Visualization,
    SampleVisualization,
    Experiment,
    VisualSearchExperiment,
    get_ordered_trials_from_psycopy_logs,
)

bids_formatting

BIDS conversion and dataset-level derivatives computation.

dataset_to_bids(target_folder_path, files_folder_path, dataset_name, session_substrings=1, format_name='eyelink')

Convert a dataset to BIDS format.

Args: target_folder_path (str): Path to the folder where the BIDS dataset will be created. files_folder_path (str): Path to the folder containing the EDF files. The EDF files are assumed to have the ID of the subject at the beginning of the file name, separated by an underscore. dataset_name (str): Name of the BIDS dataset. session_substrings (int): Number of substrings to use for the session ID. Default is 1.

Returns: None

Source code in pyxations/bids_formatting.py

def dataset_to_bids(target_folder_path, files_folder_path, dataset_name, session_substrings=1, format_name='eyelink'):
    """
    Convert a dataset to BIDS format.

    Args:
        target_folder_path (str): Path to the folder where the BIDS dataset will be created.
        files_folder_path (str): Path to the folder containing the EDF files.
        The EDF files are assumed to have the ID of the subject at the beginning of the file name, separated by an underscore.
        dataset_name (str): Name of the BIDS dataset.
        session_substrings (int): Number of substrings to use for the session ID. Default is 1.

    Returns:
        None
    """
    converter = get_converter(format_name)

    # Create a metadata tsv file
    metadata = pd.DataFrame(columns=['subject_id', 'old_subject_id'])
    files_folder_path = Path(files_folder_path)
    # List all file paths in the folder
    file_paths = []
    for file_path in files_folder_path.rglob('*'):  # Recursively go through all files
        if file_path.is_file():
            file_paths.append(file_path)

    file_paths = [file for file in file_paths if file.suffix.lower() in converter.relevant_extensions()]

    bids_folder_path = Path(target_folder_path) / dataset_name
    bids_folder_path.mkdir(parents=True, exist_ok=True)

    subj_ids = converter.get_subject_ids(file_paths)

    # If all of the subjects have numerical IDs, sort them numerically, else sort them alphabetically
    if all(subject_id.isdigit() for subject_id in subj_ids):
        subj_ids.sort(key=int)
    else:
        subj_ids.sort()
    new_subj_ids = [str(subject_index).zfill(4) for subject_index in range(1, len(subj_ids) + 1)]

    # Create subfolders for each session for each subject
    for subject_id in new_subj_ids:
        old_subject_id = subj_ids[int(subject_id) - 1]
        for file in file_paths:
            file_name = Path(file).name
            session_id = "_".join("".join(file_name.split(".")[:-1]).split("_")[1:session_substrings + 1])
            converter.move_file_to_bids_folder(file, bids_folder_path, subject_id, old_subject_id, session_id)

        metadata.loc[len(metadata.index)] = [subject_id, old_subject_id]
    # Save metadata to tsv file
    metadata.to_csv(bids_folder_path / "participants.tsv", sep="\t", index=False)
    return bids_folder_path

pre_processing

Per-recording parsing and trial segmentation.

PreProcessing

Pyxations preprocessing: trial segmentation, quality flags, and saccade direction. All mutating functions are safe (copy-aware) and validate required columns.

Tables (pd.DataFrame) expected: samples: typically contains 'tSample' (ms), gaze columns (e.g., 'LX','LY','RX','RY' or 'X','Y') fixations: typically contains 'tStart','tEnd' and optional 'xAvg','yAvg' saccades: typically contains 'tStart','tEnd','xStart','yStart','xEnd','yEnd' blinks: typically contains 'tStart','tEnd' (optional) user_messages: must contain 'timestamp','message'

New columns created: - All tables after trialing: 'phase', 'trial_number', 'trial_label' (optional) - samples/fixations/saccades: 'bad' (bool) after bad_samples() - saccades: 'deg' (float degrees), 'dir' (str) after saccades_direction()

Source code in pyxations/pre_processing.py

class PreProcessing:
    """
    Pyxations preprocessing: trial segmentation, quality flags, and saccade direction.
    All mutating functions are safe (copy-aware) and validate required columns.

    Tables (pd.DataFrame) expected:
        samples:   typically contains 'tSample' (ms), gaze columns (e.g., 'LX','LY','RX','RY' or 'X','Y')
        fixations: typically contains 'tStart','tEnd' and optional 'xAvg','yAvg'
        saccades:  typically contains 'tStart','tEnd','xStart','yStart','xEnd','yEnd'
        blinks:    typically contains 'tStart','tEnd' (optional)
        user_messages: must contain 'timestamp','message'

    New columns created:
        - All tables after trialing: 'phase', 'trial_number', 'trial_label' (optional)
        - samples/fixations/saccades: 'bad' (bool) after bad_samples()
        - saccades: 'deg' (float degrees), 'dir' (str) after saccades_direction()
    """

    VERSION = "0.2.0"

    def __init__(
        self,
        samples: pd.DataFrame,
        fixations: pd.DataFrame,
        saccades: pd.DataFrame,
        blinks: pd.DataFrame,
        user_messages: pd.DataFrame,
        session_path: PathLike,
        metadata: Optional[SessionMetadata] = None,
    ):
        self.samples = samples.copy()
        self.fixations = fixations.copy()
        self.saccades = saccades.copy()
        self.blinks = blinks.copy()
        self.user_messages = user_messages.copy()
        self.session_path = Path(session_path)
        self.metadata = metadata or SessionMetadata()

        # Normalize dtypes where possible (strings for messages)
        if "message" in self.user_messages.columns:
            self.user_messages["message"] = self.user_messages["message"].astype(str)

    # ------------------------------- Utilities ------------------------------- #

    @staticmethod
    def _require_columns(
        df: pd.DataFrame, cols: Sequence[str], context: str
    ) -> None:
        missing = [c for c in cols if c not in df.columns]
        if missing:
            raise ValueError(
                f"[{context}] Missing required columns: {missing}. "
                f"Available: {list(df.columns)}"
            )

    @staticmethod
    def _assert_nonoverlap(starts: Sequence[int], ends: Sequence[int], key: str, session: Path) -> None:
        if len(starts) != len(ends):
            raise ValueError(
                f"[{key}] start_times and end_times must have the same length, "
                f"got {len(starts)} vs {len(ends)} in session: {session}"
            )
        for i, (s, e) in enumerate(zip(starts, ends)):
            if not (s < e):
                raise ValueError(
                    f"[{key}] Non-positive interval at trial {i}: start={s}, end={e} "
                    f"in session: {session}"
                )
            if i < len(starts) - 1:
                if e > starts[i + 1]:
                    raise ValueError(
                        f"[{key}] Overlapping trials {i}–{i+1}: end[i]={e} > start[i+1]={starts[i+1]} "
                        f"in session: {session}"
                    )

    @staticmethod
    def _ensure_columns_exist(df: pd.DataFrame, cols: Sequence[str]) -> List[str]:
        """Return the subset of 'cols' that actually exist in df."""
        return [c for c in cols if c in df.columns]

    def _save_json_sidecar(self, obj: dict, filename: str) -> None:
        outdir = self.session_path
        outdir.mkdir(parents=True, exist_ok=True)
        with open(outdir / filename, "w", encoding="utf-8") as f:
            json.dump(obj, f, indent=2, ensure_ascii=False)

    # ---------------------------- Public API: Meta ---------------------------- #

    def set_metadata(
        self,
        coords_unit: Optional[str] = None,
        time_unit: Optional[str] = None,
        pupil_unit: Optional[str] = None,
        screen_width: Optional[int] = None,
        screen_height: Optional[int] = None,
        **extra,
    ) -> None:
        """Update session-level metadata used in bounds checks and documentation."""
        if coords_unit is not None:
            self.metadata.coords_unit = coords_unit
        if time_unit is not None:
            self.metadata.time_unit = time_unit
        if pupil_unit is not None:
            self.metadata.pupil_unit = pupil_unit
        if screen_width is not None:
            self.metadata.screen_width = screen_width
        if screen_height is not None:
            self.metadata.screen_height = screen_height
        self.metadata.extra.update(extra)

    def save_metadata(self, filename: str = "metadata.json") -> None:
        """Persist metadata next to derivatives for reproducibility."""
        self._save_json_sidecar(self.metadata.to_dict(), filename)

    # ----------------------- Public API: Message Parsing ---------------------- #

    def get_timestamps_from_messages(
        self,
        messages_dict: Dict[str, List[str]],
        *,
        case_insensitive: bool = True,
        use_regex: bool = True,
        return_match_token: bool = False,
    ) -> Dict[str, List[int]]:
        """
        Extract ordered timestamps per phase by matching message substrings/patterns.

        Parameters
        ----------
        messages_dict : dict
            e.g., {'trial': ['TRIAL_START', 'BEGIN_TRIAL'], 'stim': ['STIM_ONSET']}
        case_insensitive : bool
            If True, ignore case during matching.
        use_regex : bool
            If True, treat entries as regex patterns joined by '|'; otherwise escape literals.
        return_match_token : bool
            If True, also creates/updates a 'matched_token' column with the first matched pattern.

        Returns
        -------
        Dict[str, List[int]]
            Ordered timestamps in ms for each key.
        """
        df = self.user_messages
        self._require_columns(df, ["timestamp", "message"], "get_timestamps_from_messages")

        timestamps_dict: Dict[str, List[int]] = {}
        flags = re.I if case_insensitive else 0

        # Prepare an optional matched_token column for traceability
        if return_match_token and "matched_token" not in df.columns:
            df = df.copy()
            df["matched_token"] = pd.Series([None] * len(df), index=df.index)

        for key, tokens in messages_dict.items():
            if not tokens:
                raise ValueError(f"[{key}] Empty token list passed to get_timestamps_from_messages.")
            parts = tokens if use_regex else [re.escape(t) for t in tokens]
            pat = re.compile("|".join(parts), flags=flags)

            hits = df[df["message"].str.contains(pat, regex=True, na=False)].copy()
            hits.sort_values(by="timestamp", inplace=True)

            if return_match_token and not hits.empty:
                # record which token matched first for each hit
                def _which(m: str) -> Optional[str]:
                    for t in tokens:
                        if (re.search(t, m, flags=flags) if use_regex else re.search(re.escape(t), m, flags=flags)):
                            return t
                    return None

                hits["matched_token"] = hits["message"].apply(_which)
                # write back those rows (optional traceability)
                df.loc[hits.index, "matched_token"] = hits["matched_token"]

            stamps = hits["timestamp"].astype(int).tolist()
            if len(stamps) == 0:
                raise ValueError(
                    f"[{key}] No timestamps found for messages {tokens} "
                    f"in session: {self.session_path}"
                )
            timestamps_dict[key] = stamps

        # Persist updated matched_token if requested
        if return_match_token:
            self.user_messages = df

        return timestamps_dict

    # ---------------------- Public API: Trial Segmentation -------------------- #

    def split_all_into_trials(
        self,
        start_times: Dict[str, List[int]],
        end_times: Dict[str, List[int]],
        trial_labels: Optional[Dict[str, List[str]]] = None,
        *,
        allow_open_last: bool = True,
        require_nonoverlap: bool = True,
    ) -> None:
        """Segment samples/fixations/saccades/blinks using explicit times."""
        for df in (self.samples, self.fixations, self.saccades, self.blinks):
            self._split_into_trials_df(
                df, start_times, end_times, trial_labels,
                allow_open_last=allow_open_last,
                require_nonoverlap=require_nonoverlap,
            )

    def split_all_into_trials_by_msgs(
        self,
        start_msgs: Dict[str, List[str]],
        end_msgs: Dict[str, List[str]],
        trial_labels: Optional[Dict[str, List[str]]] = None,
        **msg_kwargs,
    ) -> None:
        """Segment tables using start and end message patterns."""
        starts = self.get_timestamps_from_messages(start_msgs, **msg_kwargs)
        ends = self.get_timestamps_from_messages(end_msgs, **msg_kwargs)
        self.split_all_into_trials(starts, ends, trial_labels)

    def split_all_into_trials_by_durations(
        self,
        start_msgs: Dict[str, List[str]],
        durations: Dict[str, List[int]],
        trial_labels: Optional[Dict[str, List[str]]] = None,
        **msg_kwargs,
    ) -> None:
        """Segment using start message patterns and per-trial durations (ms)."""
        starts = self.get_timestamps_from_messages(start_msgs, **msg_kwargs)
        end_times: Dict[str, List[int]] = {}
        for key, durs in durations.items():
            s = starts.get(key, [])
            if len(durs) < len(s):
                raise ValueError(
                    f"[{key}] Provided {len(durs)} durations but found {len(s)} start times "
                    f"in session: {self.session_path}"
                )
            end_times[key] = [st + du for st, du in zip(s, durs)]
        self.split_all_into_trials(starts, end_times, trial_labels)

    def _split_into_trials_df(
        self,
        data: pd.DataFrame,
        start_times: Dict[str, List[int]],
        end_times: Dict[str, List[int]],
        trial_labels: Optional[Dict[str, List[str]]] = None,
        *,
        allow_open_last: bool = True,
        require_nonoverlap: bool = True,
    ) -> None:
        """
        Core segmentation for a single table. Works with 'tSample' OR ('tStart','tEnd').
        Adds 'phase', 'trial_number', 'trial_label'.
        """
        if data is self.samples:
            time_mode = "sample"
            self._require_columns(data, ["tSample"], "split_into_trials(samples)")
        else:
            # events (fixations/saccades/blinks)
            time_mode = "event"
            self._require_columns(data, ["tStart", "tEnd"], "split_into_trials(events)")

        df = data.copy()
        # Initialize columns deterministically
        df["phase"] = ""
        df["trial_number"] = -1
        df["trial_label"] = ""

        for key in start_times.keys():
            start_list = list(start_times[key])
            end_list = list(end_times[key])

            # Discard starts after last end (common partial last-trial artifact)
            if allow_open_last and end_list:
                last_end = end_list[-1]
                start_list = [st for st in start_list if st < last_end]

            # Sanity checks
            if require_nonoverlap:
                self._assert_nonoverlap(start_list, end_list, key, self.session_path)
            elif len(start_list) != len(end_list):
                raise ValueError(
                    f"[{key}] start_times and end_times length mismatch: {len(start_list)} vs {len(end_list)} "
                    f"in session: {self.session_path}"
                )

            labels = trial_labels.get(key) if (trial_labels and key in trial_labels) else None
            if labels is not None and len(labels) != len(start_list):
                raise ValueError(
                    f"[{key}] Computed {len(start_list)} trials but got {len(labels)} trial labels "
                    f"in session: {self.session_path}"
                )

            # Apply segmentation
            if time_mode == "sample":
                t = df["tSample"].values
                for i, (st, en) in enumerate(zip(start_list, end_list)):
                    mask = (t >= st) & (t <= en)
                    if not np.any(mask):
                        continue
                    df.loc[mask, "trial_number"] = i
                    df.loc[mask, "phase"] = str(key)
                    if labels is not None:
                        df.loc[mask, "trial_label"] = labels[i]
            else:
                t0 = df["tStart"].values
                t1 = df["tEnd"].values
                for i, (st, en) in enumerate(zip(start_list, end_list)):
                    mask = (t0 >= st) & (t1 <= en)
                    if not np.any(mask):
                        continue
                    df.loc[mask, "trial_number"] = i
                    df.loc[mask, "phase"] = str(key)
                    if labels is not None:
                        df.loc[mask, "trial_label"] = labels[i]

        # Commit
        if data is self.samples:
            self.samples = df
        elif data is self.fixations:
            self.fixations = df
        elif data is self.saccades:
            self.saccades = df
        elif data is self.blinks:
            self.blinks = df

    # ------------------------- Public API: QC / Flags ------------------------- #

    def bad_samples(
        self,
        screen_height: Optional[int] = None,
        screen_width: Optional[int] = None,
        *,
        mark_nan_as_bad: bool = True,
        inclusive_bounds: bool = True,
    ) -> None:
        """
        Mark rows as 'bad' if any available coordinate falls outside screen bounds.
        Applies to samples, fixations, saccades. (Blinks unaffected.)

        If width/height not provided, will use metadata.screen_* if available.
        """
        H = screen_height if screen_height is not None else self.metadata.screen_height
        W = screen_width if screen_width is not None else self.metadata.screen_width
        if H is None or W is None:
            raise ValueError(
                "bad_samples requires screen_height and screen_width (either passed "
                "or set via set_metadata())."
            )

        def _mark(df: pd.DataFrame) -> pd.DataFrame:
            d = df.copy()

            # Gather candidate coordinate columns if present
            coord_cols = self._ensure_columns_exist(
                d,
                [
                    "LX", "LY", "RX", "RY", "X", "Y",
                    "xStart", "xEnd", "yStart", "yEnd", "xAvg", "yAvg",
                ],
            )
            if not coord_cols:
                # If no coords present, default to 'not bad'
                if "bad" not in d.columns:
                    d["bad"] = False
                return d

            xcols = [c for c in coord_cols if c.lower().startswith("x")]
            ycols = [c for c in coord_cols if c.lower().startswith("y")]

            # Validity masks for each axis
            if inclusive_bounds:
                valid_w = np.logical_and.reduce([d[c].ge(0) & d[c].le(W) for c in xcols]) if xcols else True
                valid_h = np.logical_and.reduce([d[c].ge(0) & d[c].le(H) for c in ycols]) if ycols else True
            else:
                valid_w = np.logical_and.reduce([d[c].gt(0) & d[c].lt(W) for c in xcols]) if xcols else True
                valid_h = np.logical_and.reduce([d[c].gt(0) & d[c].lt(H) for c in ycols]) if ycols else True

            bad = ~(valid_w & valid_h)
            if mark_nan_as_bad:
                bad |= d[coord_cols].isna().any(axis=1)

            d["bad"] = bad.values
            return d

        self.samples = _mark(self.samples)
        self.fixations = _mark(self.fixations)
        self.saccades = _mark(self.saccades)

    # ---------------------- Public API: Saccade Direction --------------------- #

    def saccades_direction(self, tol_deg: float = 15.0) -> None:
        """
        Compute saccade angle (deg) and cardinal direction with tolerance bands.

        Parameters
        ----------
        tol_deg : float
            Half-width of the acceptance band around 0°, ±90°, and ±180°
            for classifying right/left/up/down.
        """
        df = self.saccades.copy()
        self._require_columns(
            df, ["xStart", "xEnd", "yStart", "yEnd"], "saccades_direction"
        )

        x_dif = df["xEnd"].astype(float) - df["xStart"].astype(float)
        y_dif = df["yEnd"].astype(float) - df["yStart"].astype(float)
        deg = np.degrees(np.arctan2(y_dif.to_numpy(), x_dif.to_numpy()))
        df["deg"] = deg.astype(float)

        # Tolerant direction bins
        right = (-tol_deg < df["deg"]) & (df["deg"] < tol_deg)
        left = (df["deg"] > 180 - tol_deg) | (df["deg"] < -180 + tol_deg)
        down = ((90 - tol_deg) < df["deg"]) & (df["deg"] < (90 + tol_deg))
        up = ((-90 - tol_deg) < df["deg"]) & (df["deg"] < (-90 + tol_deg))

        df["dir"] = ""
        df.loc[right, "dir"] = "right"
        df.loc[left, "dir"] = "left"
        df.loc[down, "dir"] = "down"
        df.loc[up, "dir"] = "up"

        self.saccades = df

    # -------------------------- Public API: Orchestrator ---------------------- #

    def process(
        self,
        functions_and_params: Dict[str, Dict],
        *,
        log_recipe: bool = True,
        recipe_filename: str = "preprocessing_recipe.json",
        provenance_filename: str = "preprocessing_provenance.json",
    ) -> None:
        """
        Run a declarative preprocessing recipe, e.g.:
            pp.process({
                "split_all_into_trials_by_msgs": {
                    "start_msgs": {"trial": ["TRIAL_START"]},
                    "end_msgs": {"trial": ["TRIAL_END"]},
                },
                "bad_samples": {"screen_height": 1080, "screen_width": 1920},
                "saccades_direction": {"tol_deg": 15},
            })

        Unknown function names raise a helpful error.
        """
        # Optional: save the declared recipe for exact reproducibility
        if log_recipe:
            recipe_obj = {
                "declared_recipe": functions_and_params,
                "tool_version": self.VERSION,
                "timestamp_utc": datetime.now(timezone.utc).isoformat(),
                "session_path": str(self.session_path),
            }
            self._save_json_sidecar(recipe_obj, recipe_filename)

        for func_name, params in functions_and_params.items():
            if not hasattr(self, func_name):
                raise AttributeError(
                    f"Unknown preprocessing function '{func_name}'. "
                    f"Available: {[m for m in dir(self) if not m.startswith('_')]}"
                )
            fn = getattr(self, func_name)
            if not isinstance(params, dict):
                raise TypeError(
                    f"Parameters for '{func_name}' must be a dict, got {type(params)}"
                )
            fn(**params)

        # Save lightweight provenance after successful run
        if log_recipe:
            prov = {
                "completed_recipe": list(functions_and_params.keys()),
                "tool_version": self.VERSION,
                "timestamp_utc": datetime.now(timezone.utc).isoformat(),
                "metadata": self.metadata.to_dict(),
            }
            self._save_json_sidecar(prov, provenance_filename)

    # -------------------- Public API: Behavioral Metadata -------------------- #

    def add_trial_metadata(
        self,
        metadata_df: pd.DataFrame,
        columns: List[str],
    ) -> None:
        """
        Propagate experiment-specific columns from behavioral data to gaze samples.

        Joins by ``trial_index`` (one row per trial in *metadata_df*, many rows
        per trial in *samples*).  After this call ``self.samples`` will contain
        the requested columns, so downstream analysis never needs to re-read
        the original behavioral file.

        Parameters
        ----------
        metadata_df : pd.DataFrame
            Behavioral data table with one row per trial.
            Must contain a ``trial_index`` column.
        columns : list of str
            Column names from *metadata_df* to propagate to samples.

        Raises
        ------
        ValueError
            If ``trial_index`` is missing from *metadata_df* or *self.samples*.
        """
        if "trial_index" not in metadata_df.columns:
            raise ValueError(
                "[add_trial_metadata] metadata_df must contain a 'trial_index' column."
            )
        if "trial_index" not in self.samples.columns:
            raise ValueError(
                "[add_trial_metadata] samples must contain a 'trial_index' column. "
                "Make sure the parser preserves it."
            )

        available = [c for c in columns if c in metadata_df.columns]
        skipped = [c for c in columns if c not in metadata_df.columns]
        if skipped:
            print(f"[add_trial_metadata] Columns not found in metadata_df, skipping: {skipped}")
        if not available:
            return

        trial_meta = (
            metadata_df[["trial_index"] + available]
            .drop_duplicates("trial_index")
        )
        self.samples = self.samples.merge(trial_meta, on="trial_index", how="left")

add_trial_metadata(metadata_df, columns)

Propagate experiment-specific columns from behavioral data to gaze samples.

Joins by trial_index (one row per trial in metadata_df, many rows per trial in samples). After this call self.samples will contain the requested columns, so downstream analysis never needs to re-read the original behavioral file.

Parameters:

Name	Type	Description	Default
metadata_df	DataFrame	Behavioral data table with one row per trial. Must contain a trial_index column.	required
columns	list of str	Column names from metadata_df to propagate to samples.	required

Raises:

Type	Description
ValueError	If trial_index is missing from metadata_df or self.samples.

Source code in pyxations/pre_processing.py

def add_trial_metadata(
    self,
    metadata_df: pd.DataFrame,
    columns: List[str],
) -> None:
    """
    Propagate experiment-specific columns from behavioral data to gaze samples.

    Joins by ``trial_index`` (one row per trial in *metadata_df*, many rows
    per trial in *samples*).  After this call ``self.samples`` will contain
    the requested columns, so downstream analysis never needs to re-read
    the original behavioral file.

    Parameters
    ----------
    metadata_df : pd.DataFrame
        Behavioral data table with one row per trial.
        Must contain a ``trial_index`` column.
    columns : list of str
        Column names from *metadata_df* to propagate to samples.

    Raises
    ------
    ValueError
        If ``trial_index`` is missing from *metadata_df* or *self.samples*.
    """
    if "trial_index" not in metadata_df.columns:
        raise ValueError(
            "[add_trial_metadata] metadata_df must contain a 'trial_index' column."
        )
    if "trial_index" not in self.samples.columns:
        raise ValueError(
            "[add_trial_metadata] samples must contain a 'trial_index' column. "
            "Make sure the parser preserves it."
        )

    available = [c for c in columns if c in metadata_df.columns]
    skipped = [c for c in columns if c not in metadata_df.columns]
    if skipped:
        print(f"[add_trial_metadata] Columns not found in metadata_df, skipping: {skipped}")
    if not available:
        return

    trial_meta = (
        metadata_df[["trial_index"] + available]
        .drop_duplicates("trial_index")
    )
    self.samples = self.samples.merge(trial_meta, on="trial_index", how="left")

bad_samples(screen_height=None, screen_width=None, *, mark_nan_as_bad=True, inclusive_bounds=True)

Mark rows as 'bad' if any available coordinate falls outside screen bounds. Applies to samples, fixations, saccades. (Blinks unaffected.)

If width/height not provided, will use metadata.screen_* if available.

Source code in pyxations/pre_processing.py

def bad_samples(
    self,
    screen_height: Optional[int] = None,
    screen_width: Optional[int] = None,
    *,
    mark_nan_as_bad: bool = True,
    inclusive_bounds: bool = True,
) -> None:
    """
    Mark rows as 'bad' if any available coordinate falls outside screen bounds.
    Applies to samples, fixations, saccades. (Blinks unaffected.)

    If width/height not provided, will use metadata.screen_* if available.
    """
    H = screen_height if screen_height is not None else self.metadata.screen_height
    W = screen_width if screen_width is not None else self.metadata.screen_width
    if H is None or W is None:
        raise ValueError(
            "bad_samples requires screen_height and screen_width (either passed "
            "or set via set_metadata())."
        )

    def _mark(df: pd.DataFrame) -> pd.DataFrame:
        d = df.copy()

        # Gather candidate coordinate columns if present
        coord_cols = self._ensure_columns_exist(
            d,
            [
                "LX", "LY", "RX", "RY", "X", "Y",
                "xStart", "xEnd", "yStart", "yEnd", "xAvg", "yAvg",
            ],
        )
        if not coord_cols:
            # If no coords present, default to 'not bad'
            if "bad" not in d.columns:
                d["bad"] = False
            return d

        xcols = [c for c in coord_cols if c.lower().startswith("x")]
        ycols = [c for c in coord_cols if c.lower().startswith("y")]

        # Validity masks for each axis
        if inclusive_bounds:
            valid_w = np.logical_and.reduce([d[c].ge(0) & d[c].le(W) for c in xcols]) if xcols else True
            valid_h = np.logical_and.reduce([d[c].ge(0) & d[c].le(H) for c in ycols]) if ycols else True
        else:
            valid_w = np.logical_and.reduce([d[c].gt(0) & d[c].lt(W) for c in xcols]) if xcols else True
            valid_h = np.logical_and.reduce([d[c].gt(0) & d[c].lt(H) for c in ycols]) if ycols else True

        bad = ~(valid_w & valid_h)
        if mark_nan_as_bad:
            bad |= d[coord_cols].isna().any(axis=1)

        d["bad"] = bad.values
        return d

    self.samples = _mark(self.samples)
    self.fixations = _mark(self.fixations)
    self.saccades = _mark(self.saccades)

get_timestamps_from_messages(messages_dict, *, case_insensitive=True, use_regex=True, return_match_token=False)

Extract ordered timestamps per phase by matching message substrings/patterns.

Parameters:

Name	Type	Description	Default
messages_dict	dict	e.g., {'trial': ['TRIAL_START', 'BEGIN_TRIAL'], 'stim': ['STIM_ONSET']}	required
case_insensitive	bool	If True, ignore case during matching.	True
use_regex	bool	If True, treat entries as regex patterns joined by '|'; otherwise escape literals.	True
return_match_token	bool	If True, also creates/updates a 'matched_token' column with the first matched pattern.	False

Returns:

Type	Description
Dict[str, List[int]]	Ordered timestamps in ms for each key.

Source code in pyxations/pre_processing.py

def get_timestamps_from_messages(
    self,
    messages_dict: Dict[str, List[str]],
    *,
    case_insensitive: bool = True,
    use_regex: bool = True,
    return_match_token: bool = False,
) -> Dict[str, List[int]]:
    """
    Extract ordered timestamps per phase by matching message substrings/patterns.

    Parameters
    ----------
    messages_dict : dict
        e.g., {'trial': ['TRIAL_START', 'BEGIN_TRIAL'], 'stim': ['STIM_ONSET']}
    case_insensitive : bool
        If True, ignore case during matching.
    use_regex : bool
        If True, treat entries as regex patterns joined by '|'; otherwise escape literals.
    return_match_token : bool
        If True, also creates/updates a 'matched_token' column with the first matched pattern.

    Returns
    -------
    Dict[str, List[int]]
        Ordered timestamps in ms for each key.
    """
    df = self.user_messages
    self._require_columns(df, ["timestamp", "message"], "get_timestamps_from_messages")

    timestamps_dict: Dict[str, List[int]] = {}
    flags = re.I if case_insensitive else 0

    # Prepare an optional matched_token column for traceability
    if return_match_token and "matched_token" not in df.columns:
        df = df.copy()
        df["matched_token"] = pd.Series([None] * len(df), index=df.index)

    for key, tokens in messages_dict.items():
        if not tokens:
            raise ValueError(f"[{key}] Empty token list passed to get_timestamps_from_messages.")
        parts = tokens if use_regex else [re.escape(t) for t in tokens]
        pat = re.compile("|".join(parts), flags=flags)

        hits = df[df["message"].str.contains(pat, regex=True, na=False)].copy()
        hits.sort_values(by="timestamp", inplace=True)

        if return_match_token and not hits.empty:
            # record which token matched first for each hit
            def _which(m: str) -> Optional[str]:
                for t in tokens:
                    if (re.search(t, m, flags=flags) if use_regex else re.search(re.escape(t), m, flags=flags)):
                        return t
                return None

            hits["matched_token"] = hits["message"].apply(_which)
            # write back those rows (optional traceability)
            df.loc[hits.index, "matched_token"] = hits["matched_token"]

        stamps = hits["timestamp"].astype(int).tolist()
        if len(stamps) == 0:
            raise ValueError(
                f"[{key}] No timestamps found for messages {tokens} "
                f"in session: {self.session_path}"
            )
        timestamps_dict[key] = stamps

    # Persist updated matched_token if requested
    if return_match_token:
        self.user_messages = df

    return timestamps_dict

process(functions_and_params, *, log_recipe=True, recipe_filename='preprocessing_recipe.json', provenance_filename='preprocessing_provenance.json')

Run a declarative preprocessing recipe, e.g.: pp.process({ "split_all_into_trials_by_msgs": { "start_msgs": {"trial": ["TRIAL_START"]}, "end_msgs": {"trial": ["TRIAL_END"]}, }, "bad_samples": {"screen_height": 1080, "screen_width": 1920}, "saccades_direction": {"tol_deg": 15}, })

Unknown function names raise a helpful error.

Source code in pyxations/pre_processing.py

def process(
    self,
    functions_and_params: Dict[str, Dict],
    *,
    log_recipe: bool = True,
    recipe_filename: str = "preprocessing_recipe.json",
    provenance_filename: str = "preprocessing_provenance.json",
) -> None:
    """
    Run a declarative preprocessing recipe, e.g.:
        pp.process({
            "split_all_into_trials_by_msgs": {
                "start_msgs": {"trial": ["TRIAL_START"]},
                "end_msgs": {"trial": ["TRIAL_END"]},
            },
            "bad_samples": {"screen_height": 1080, "screen_width": 1920},
            "saccades_direction": {"tol_deg": 15},
        })

    Unknown function names raise a helpful error.
    """
    # Optional: save the declared recipe for exact reproducibility
    if log_recipe:
        recipe_obj = {
            "declared_recipe": functions_and_params,
            "tool_version": self.VERSION,
            "timestamp_utc": datetime.now(timezone.utc).isoformat(),
            "session_path": str(self.session_path),
        }
        self._save_json_sidecar(recipe_obj, recipe_filename)

    for func_name, params in functions_and_params.items():
        if not hasattr(self, func_name):
            raise AttributeError(
                f"Unknown preprocessing function '{func_name}'. "
                f"Available: {[m for m in dir(self) if not m.startswith('_')]}"
            )
        fn = getattr(self, func_name)
        if not isinstance(params, dict):
            raise TypeError(
                f"Parameters for '{func_name}' must be a dict, got {type(params)}"
            )
        fn(**params)

    # Save lightweight provenance after successful run
    if log_recipe:
        prov = {
            "completed_recipe": list(functions_and_params.keys()),
            "tool_version": self.VERSION,
            "timestamp_utc": datetime.now(timezone.utc).isoformat(),
            "metadata": self.metadata.to_dict(),
        }
        self._save_json_sidecar(prov, provenance_filename)

saccades_direction(tol_deg=15.0)

Compute saccade angle (deg) and cardinal direction with tolerance bands.

Parameters:

Name	Type	Description	Default
tol_deg	float	Half-width of the acceptance band around 0°, ±90°, and ±180° for classifying right/left/up/down.	15.0

Source code in pyxations/pre_processing.py

def saccades_direction(self, tol_deg: float = 15.0) -> None:
    """
    Compute saccade angle (deg) and cardinal direction with tolerance bands.

    Parameters
    ----------
    tol_deg : float
        Half-width of the acceptance band around 0°, ±90°, and ±180°
        for classifying right/left/up/down.
    """
    df = self.saccades.copy()
    self._require_columns(
        df, ["xStart", "xEnd", "yStart", "yEnd"], "saccades_direction"
    )

    x_dif = df["xEnd"].astype(float) - df["xStart"].astype(float)
    y_dif = df["yEnd"].astype(float) - df["yStart"].astype(float)
    deg = np.degrees(np.arctan2(y_dif.to_numpy(), x_dif.to_numpy()))
    df["deg"] = deg.astype(float)

    # Tolerant direction bins
    right = (-tol_deg < df["deg"]) & (df["deg"] < tol_deg)
    left = (df["deg"] > 180 - tol_deg) | (df["deg"] < -180 + tol_deg)
    down = ((90 - tol_deg) < df["deg"]) & (df["deg"] < (90 + tol_deg))
    up = ((-90 - tol_deg) < df["deg"]) & (df["deg"] < (-90 + tol_deg))

    df["dir"] = ""
    df.loc[right, "dir"] = "right"
    df.loc[left, "dir"] = "left"
    df.loc[down, "dir"] = "down"
    df.loc[up, "dir"] = "up"

    self.saccades = df

save_metadata(filename='metadata.json')

Persist metadata next to derivatives for reproducibility.

Source code in pyxations/pre_processing.py

def save_metadata(self, filename: str = "metadata.json") -> None:
    """Persist metadata next to derivatives for reproducibility."""
    self._save_json_sidecar(self.metadata.to_dict(), filename)

set_metadata(coords_unit=None, time_unit=None, pupil_unit=None, screen_width=None, screen_height=None, **extra)

Update session-level metadata used in bounds checks and documentation.

Source code in pyxations/pre_processing.py

def set_metadata(
    self,
    coords_unit: Optional[str] = None,
    time_unit: Optional[str] = None,
    pupil_unit: Optional[str] = None,
    screen_width: Optional[int] = None,
    screen_height: Optional[int] = None,
    **extra,
) -> None:
    """Update session-level metadata used in bounds checks and documentation."""
    if coords_unit is not None:
        self.metadata.coords_unit = coords_unit
    if time_unit is not None:
        self.metadata.time_unit = time_unit
    if pupil_unit is not None:
        self.metadata.pupil_unit = pupil_unit
    if screen_width is not None:
        self.metadata.screen_width = screen_width
    if screen_height is not None:
        self.metadata.screen_height = screen_height
    self.metadata.extra.update(extra)

split_all_into_trials(start_times, end_times, trial_labels=None, *, allow_open_last=True, require_nonoverlap=True)

Segment samples/fixations/saccades/blinks using explicit times.

Source code in pyxations/pre_processing.py

def split_all_into_trials(
    self,
    start_times: Dict[str, List[int]],
    end_times: Dict[str, List[int]],
    trial_labels: Optional[Dict[str, List[str]]] = None,
    *,
    allow_open_last: bool = True,
    require_nonoverlap: bool = True,
) -> None:
    """Segment samples/fixations/saccades/blinks using explicit times."""
    for df in (self.samples, self.fixations, self.saccades, self.blinks):
        self._split_into_trials_df(
            df, start_times, end_times, trial_labels,
            allow_open_last=allow_open_last,
            require_nonoverlap=require_nonoverlap,
        )

split_all_into_trials_by_durations(start_msgs, durations, trial_labels=None, **msg_kwargs)

Segment using start message patterns and per-trial durations (ms).

Source code in pyxations/pre_processing.py

def split_all_into_trials_by_durations(
    self,
    start_msgs: Dict[str, List[str]],
    durations: Dict[str, List[int]],
    trial_labels: Optional[Dict[str, List[str]]] = None,
    **msg_kwargs,
) -> None:
    """Segment using start message patterns and per-trial durations (ms)."""
    starts = self.get_timestamps_from_messages(start_msgs, **msg_kwargs)
    end_times: Dict[str, List[int]] = {}
    for key, durs in durations.items():
        s = starts.get(key, [])
        if len(durs) < len(s):
            raise ValueError(
                f"[{key}] Provided {len(durs)} durations but found {len(s)} start times "
                f"in session: {self.session_path}"
            )
        end_times[key] = [st + du for st, du in zip(s, durs)]
    self.split_all_into_trials(starts, end_times, trial_labels)

split_all_into_trials_by_msgs(start_msgs, end_msgs, trial_labels=None, **msg_kwargs)

Segment tables using start and end message patterns.

Source code in pyxations/pre_processing.py

def split_all_into_trials_by_msgs(
    self,
    start_msgs: Dict[str, List[str]],
    end_msgs: Dict[str, List[str]],
    trial_labels: Optional[Dict[str, List[str]]] = None,
    **msg_kwargs,
) -> None:
    """Segment tables using start and end message patterns."""
    starts = self.get_timestamps_from_messages(start_msgs, **msg_kwargs)
    ends = self.get_timestamps_from_messages(end_msgs, **msg_kwargs)
    self.split_all_into_trials(starts, ends, trial_labels)

SessionMetadata dataclass

Lightweight metadata container saved alongside derivatives.

Source code in pyxations/pre_processing.py

@dataclass
class SessionMetadata:
    """Lightweight metadata container saved alongside derivatives."""
    coords_unit: str = "px"          # 'px' or 'deg'
    time_unit: str = "ms"            # 'ms'
    pupil_unit: str = "arbitrary"
    screen_width: Optional[int] = None
    screen_height: Optional[int] = None
    extra: Dict[str, Union[str, int, float, bool, None]] = field(default_factory=dict)

    def to_dict(self) -> dict:
        return {
            "coords_unit": self.coords_unit,
            "time_unit": self.time_unit,
            "pupil_unit": self.pupil_unit,
            "screen_width": self.screen_width,
            "screen_height": self.screen_height,
            "extra": self.extra,
        }

methods.eyemovement

Fixation and saccade detection algorithms.

eye_movement_detection

REMoDNaV

RemodnavDetection

Bases: EyeMovementDetection

Source code in pyxations/methods/eyemovement/REMoDNaV.py

class RemodnavDetection(EyeMovementDetection):

    def __init__(self, session_folder_path, samples):
        self.session_folder_path = session_folder_path
        self.out_folder = (session_folder_path / 'remodnav_events')
        self.samples = samples

    def detect_eye_movements(self, min_pursuit_dur:float=10., max_pso_dur:float=0.0, min_fix_dur:float=0.05,
                                 sac_max_vel:float=1000., fix_max_amp:float=1.5, sac_time_thresh:float=0.002,
                                 drop_fix_from_blink:bool=True, screen_size:float=38., screen_width:int=1920, screen_distance:float=60,
                                 savgol_length=0.195):

        """
        Detects fixations and saccades from eye-tracking data for both left and right eyes using REMoDNaV, a velocity based eye movement event detection algorithm 
        that is based on, but extends the adaptive Nyström & Holmqvist algorithm (Nyström & Holmqvist, 2010).

        Parameters
        ----------
        min_pursuit_dur : float, optional
            Minimum pursuit duration in seconds for Remodnav detection (default is 10.0).
        max_pso_dur : float, optional
            Maximum post-saccadic oscillation duration in seconds for Remodnav detection (default is 0.0 -No PSO events detection-).
        min_fix_dur : float, optional
            Minimum fixation duration in seconds for Remodnav detection (default is 0.05).
        sac_max_vel : float, optional
            Maximum saccade velocity in deg/s (default is 1000.0).
        fix_max_amp : float, optional
            Maximum fixation amplitude in deg (default is 1.5).
        sac_time_thresh : float, optional
            Time threshold in seconds to consider a saccade as neighboring a fixation (default is 0.002).
        drop_fix_from_blink : bool, optional
            Whether to drop fixations that do not have a previous saccade within the time threshold (default is True).
        screen_size : float, optional
            Size of the screen in cm (default is 38.0).
        screen_width : int, optional
            Horizontal resolution of the screen in pixels (default is 1920).
        screen_distance : float, optional
            Distance from the screen to the participant's eyes in cm (default is 60).

        Returns
        -------
        fixations : dict
            Dictionary containing DataFrames of detected fixations for both left and right eyes with additional columns for mean x, y positions, and pupil size.
        saccades : dict
            Dictionary containing DataFrames of detected saccades for both left and right eyes.

        Raises
        ------
        ValueError
            If Remodnav detection fails.
        """

        # Move eye data, detections file and image to subject results directory
        self.out_folder.mkdir(parents=True, exist_ok=True)

        # Dictionaries to store fixations and saccades DataFrames
        fixations = []
        saccades = []

        # Divide samples by continuous data chunks
        # If the difference in the column "tSample" from row i to row i+1 is greater than 1000 / sample_rate, then it is a new chunk

        # Logic to find chunks
        chunks_start_indexes = np.where(np.diff(self.samples['tSample']) > (1000 / self.samples['Rate_recorded'][0]))[0] + 1

        chunks = np.zeros(len(self.samples))
        chunks[chunks_start_indexes] = 1
        chunks = np.cumsum(chunks)

        # Divide samples by chunks and apply to each chunk the detection algorithm
        for chunk in np.unique(chunks): 
            chunk_samples = self.samples[chunks == chunk].reset_index(drop=True)
            fixations_chunk, saccades_chunk = self.detect_on_chunk(chunk_samples, min_pursuit_dur, max_pso_dur, min_fix_dur, 
                            sac_max_vel, fix_max_amp, sac_time_thresh, drop_fix_from_blink, screen_size, screen_width, screen_distance,
                            savgol_length)
            fixations.append(fixations_chunk)
            saccades.append(saccades_chunk)

        return pd.concat(fixations).sort_values(by='tEnd', ignore_index=True), pd.concat(saccades).sort_values(by='tEnd', ignore_index=True)



    def run_eye_movement_from_samples(self, sample_rate, x_label='X', y_label='Y', config={}, **kwargs):
        '''
        Recieves a pandas dataframe, a sample rate, and optional configuration
        :param dfSamples: Pandas dataframe including x and y columns
        :param sample_rate: an integer representing the data sample rate.
        :param x_label: X column name
        :param y_label: Y column name
        :param config:
        '''
        starting_time = self.samples['tSample'].min()
        eye_data = {
            'x': self.samples[x_label], 
            'y': self.samples[y_label]
        }
        eye_data = np.rec.fromarrays(list(eye_data.values()), names=list(eye_data.keys()))

        if 'pupil_data' not in config.keys():
            config['pupil_data'] = pd.Series([0]* len(eye_data['x']))
        times = np.arange(stop=len(eye_data['x'])) / sample_rate

        return self.run_eye_movement(eye_data['x'], eye_data['y'], sample_rate, 
            times=times, starting_time=starting_time, **config, **kwargs)


    def run_eye_movement(self, gazex_data, gazey_data, sample_rate,
             min_pursuit_dur:float=10., max_pso_dur:float=0.0, min_fix_dur:float=0.05,
             min_saccade_duration=0.04,
             sac_max_vel:float=1000., fix_max_amp:float=1.5, sac_time_thresh:float=0.002,
             drop_fix_from_blink:bool=True, screen_size:float=38., screen_width:int=1920, 
             screen_distance:float=60, calib_index=0, savgol_length=0.19,
             eyes_recorded=None, starting_time=None, times=None, pupil_data=None, eye=None):

        # If not pre run data, run
        print(f'\nRunning eye movements detection for {eye} eye...')

        # Define data to save to excel file needed to run the saccades detection program Remodnav
        eye_data = {'x':gazex_data, 'y':gazey_data}  # eye_data to numpy record array
        eye_data = np.rec.fromarrays(list(eye_data.values()), names=list(eye_data.keys()))

        # Remodnav parameters
        px2deg = math.degrees(math.atan2(.5 * screen_size, screen_distance)) / (.5 * screen_width)  # Pixel to degree conversion

        # Run Remodnav not considering pursuit class and min fixations 50 ms
        clf = EyegazeClassifier(px2deg=px2deg,
            sampling_rate=sample_rate,
            min_pursuit_duration=min_pursuit_dur,
            max_pso_duration=max_pso_dur,
            min_fixation_duration=min_fix_dur,
            min_saccade_duration=min_saccade_duration)
        pp = clf.preproc(
            eye_data,
            savgol_length=savgol_length
        )


        sac_fix = clf(pp, classify_isp=True, sort_events=True)


        # sac_fix to pandas DataFrame
        sac_fix = pd.DataFrame(sac_fix, columns=['start_time', 'end_time', 'label', 'start_x', 'start_y', 'end_x', 'end_y', 'amp', 'peak_vel', 'med_vel', 'avg_vel'])
        sac_fix['duration'] = sac_fix['end_time'] - sac_fix['start_time']
        sac_fix.rename(columns={'start_time':'onset'}, inplace=True)
        # Get saccades and fixations
        saccades_eye_all = sac_fix.loc[(sac_fix['label'] == 'SACC') | (sac_fix['label'] == 'ISAC')]
        fixations_eye_all = sac_fix.loc[sac_fix['label'] == 'FIXA']
        # Drop saccades and fixations based on conditions
        print(f'Dropping saccades with average vel > {sac_max_vel} deg/s, and fixations with amplitude > {fix_max_amp} deg')
        fixations_eye = fixations_eye_all[fixations_eye_all['amp'] <= fix_max_amp].sort_values(by='start_x', ignore_index=True)
        saccades_eye = saccades_eye_all[saccades_eye_all['peak_vel'] <= sac_max_vel].sort_values(by='start_x', ignore_index=True)
        fixations_eye.drop(columns=['label'], inplace=True)
        saccades_eye.drop(columns=['label'], inplace=True)
        print(f'Kept {len(fixations_eye)} out of {len(fixations_eye_all)} fixations')
        print(f'Kept {len(saccades_eye)} out of {len(saccades_eye_all)} saccades')
        # Variables to save fixations features
        mean_x = []
        mean_y = []
        pupil_size = []
        # Drop when no previous saccade detected in sac_time_thresh
        if drop_fix_from_blink:
            prev_sac = []
            # Identify neighbour saccades to each fixation (considering sac_time_thresh)
            print('Finding previous and next saccades')
            for fix_idx, fixation in tqdm(fixations_eye.iterrows(), total=len(fixations_eye)):
                fix_time = fixation['onset']
                fix_dur = fixation['duration']
                # Previous and next saccades
                try:
                    sac0 = saccades_eye.loc[(saccades_eye['onset'] + saccades_eye['duration'] > fix_time - sac_time_thresh) & (saccades_eye['onset'] + saccades_eye['duration'] < fix_time + sac_time_thresh)].index.values[-1]
                except:
                    sac0 = -1
                prev_sac.append(sac0)

            # Add columns
            fixations_eye['prev_sac'] = prev_sac
            fixations_eye.drop(fixations_eye[fixations_eye['prev_sac'] == -1].index, inplace=True)
            print(f'\nKept {len(fixations_eye)} fixations with previous saccade')
            fixations_eye.drop(columns=['prev_sac'], inplace=True)
        # Fixations features
        print('Computing average pupil size, and x and y position')
        for fix_idx, fixation in tqdm(fixations_eye.iterrows(), total=len(fixations_eye)):
            fix_time = fixation['onset']
            fix_dur = fixation['duration']
            # Average pupil size, x and y position
            fix_time_idx = np.where(np.logical_and(times > fix_time, times < fix_time + fix_dur))[0]
            pupil_data_fix = pupil_data[fix_time_idx]
            gazex_data_fix = gazex_data[fix_time_idx]
            gazey_data_fix = gazey_data[fix_time_idx]
            pupil_size.append(np.nanmean(pupil_data_fix))
            mean_x.append(np.nanmean(gazex_data_fix))
            mean_y.append(np.nanmean(gazey_data_fix))

        fixations_eye['xAvg'] = mean_x
        fixations_eye['yAvg'] = mean_y
        fixations_eye['pupilAvg'] = pupil_size
        fixations_eye = fixations_eye.astype({'xAvg':float, 'yAvg':float, 'pupilAvg':float})
        fixations_eye['onset'] = fixations_eye['onset'] * 1000 + starting_time  # Convert to ms
        fixations_eye['end_time'] = fixations_eye['end_time'] * 1000 + starting_time  # Convert to ms
        fixations_eye['duration'] = fixations_eye['duration'] * 1000  # Convert to ms
        saccades_eye['onset'] = saccades_eye['onset'] * 1000 + starting_time  # Convert to ms
        saccades_eye['end_time'] = saccades_eye['end_time'] * 1000 + starting_time  # Convert to ms
        saccades_eye['duration'] = saccades_eye['duration'] * 1000  # Convert to ms
        # Rename columns to match samples columns names
        fixations_eye.rename(columns={'start_x':'xStart', 'start_y':'yStart', 'end_x':'xEnd', 'end_y':'yEnd', 'onset':'tStart', 'end_time':'tEnd', 'amp':'ampDeg', 'peak_vel':'vPeak'}, inplace=True)
        saccades_eye.rename(columns={'start_x':'xStart', 'start_y':'yStart', 'end_x':'xEnd', 'end_y':'yEnd', 'onset':'tStart', 'end_time':'tEnd', 'amp':'ampDeg', 'peak_vel':'vPeak'}, inplace=True)
        fixations_eye['Calib_index'] = calib_index
        saccades_eye['Calib_index'] = calib_index
        fixations_eye['Eyes_recorded'] = eyes_recorded
        saccades_eye['Eyes_recorded'] = eyes_recorded
        fixations_eye['Rate_recorded'] = sample_rate
        saccades_eye['Rate_recorded'] = sample_rate
        fixations_eye['eye'] = eye
        saccades_eye['eye'] = eye  # Save to dictionary

        return fixations_eye, saccades_eye

    def detect_on_chunk(self, chunk, min_pursuit_dur:float=10., max_pso_dur:float=0.0, min_fix_dur:float=0.05,
                                 sac_max_vel:float=1000., fix_max_amp:float=1.5, sac_time_thresh:float=0.002,
                                 drop_fix_from_blink:bool=True, screen_size:float=38., screen_width:int=1920, screen_distance:float=60,
                                 savgol_length:float=0.19):

        sample_rate = chunk["Rate_recorded"].iloc[0]
        calib_index = chunk["Calib_index"].iloc[0]
        eyes_recorded = chunk["Eyes_recorded"].iloc[0]
        starting_time = chunk["tSample"].iloc[0]

        # Check that the columns Rate_recorded, Calib_index and Eyes_recorded are the same for all the samples in the chunk
        assert (chunk["Rate_recorded"] == sample_rate).all()
        assert (chunk["Calib_index"] == calib_index).all()
        assert (chunk["Eyes_recorded"] == eyes_recorded).all()

        times = np.arange(stop=len(chunk)) / sample_rate

        fixations = []
        saccades = []

        min_saccade_duration=0.04

        for gazex_data, gazey_data, pupil_data, eye in zip((chunk['LX'], chunk['RX']),
                                                        (chunk['LY'], chunk['RY']),
                                                        (chunk['LPupil'], chunk['RPupil']),
                                                        ('L', 'R')):

            fixations_eye, saccades_eye = self.run_eye_movement(
                gazex_data, gazey_data, sample_rate,
                min_pursuit_dur, max_pso_dur, min_fix_dur, 
                min_saccade_duration,
                sac_max_vel, fix_max_amp, sac_time_thresh,
                drop_fix_from_blink, screen_size, screen_width, screen_distance, 
                calib_index, savgol_length, eyes_recorded, starting_time, times, pupil_data, eye)


            fixations.append(fixations_eye)
            saccades.append(saccades_eye)

        return pd.concat(fixations), pd.concat(saccades)

detect_eye_movements(min_pursuit_dur=10.0, max_pso_dur=0.0, min_fix_dur=0.05, sac_max_vel=1000.0, fix_max_amp=1.5, sac_time_thresh=0.002, drop_fix_from_blink=True, screen_size=38.0, screen_width=1920, screen_distance=60, savgol_length=0.195)

Detects fixations and saccades from eye-tracking data for both left and right eyes using REMoDNaV, a velocity based eye movement event detection algorithm that is based on, but extends the adaptive Nyström & Holmqvist algorithm (Nyström & Holmqvist, 2010).

Parameters:

Name	Type	Description	Default
min_pursuit_dur	float	Minimum pursuit duration in seconds for Remodnav detection (default is 10.0).	10.0
max_pso_dur	float	Maximum post-saccadic oscillation duration in seconds for Remodnav detection (default is 0.0 -No PSO events detection-).	0.0
min_fix_dur	float	Minimum fixation duration in seconds for Remodnav detection (default is 0.05).	0.05
sac_max_vel	float	Maximum saccade velocity in deg/s (default is 1000.0).	1000.0
fix_max_amp	float	Maximum fixation amplitude in deg (default is 1.5).	1.5
sac_time_thresh	float	Time threshold in seconds to consider a saccade as neighboring a fixation (default is 0.002).	0.002
drop_fix_from_blink	bool	Whether to drop fixations that do not have a previous saccade within the time threshold (default is True).	True
screen_size	float	Size of the screen in cm (default is 38.0).	38.0
screen_width	int	Horizontal resolution of the screen in pixels (default is 1920).	1920
screen_distance	float	Distance from the screen to the participant's eyes in cm (default is 60).	60

Returns:

Name	Type	Description
fixations	dict	Dictionary containing DataFrames of detected fixations for both left and right eyes with additional columns for mean x, y positions, and pupil size.
saccades	dict	Dictionary containing DataFrames of detected saccades for both left and right eyes.

Raises:

Type	Description
ValueError	If Remodnav detection fails.

Source code in pyxations/methods/eyemovement/REMoDNaV.py

def detect_eye_movements(self, min_pursuit_dur:float=10., max_pso_dur:float=0.0, min_fix_dur:float=0.05,
                             sac_max_vel:float=1000., fix_max_amp:float=1.5, sac_time_thresh:float=0.002,
                             drop_fix_from_blink:bool=True, screen_size:float=38., screen_width:int=1920, screen_distance:float=60,
                             savgol_length=0.195):

    """
    Detects fixations and saccades from eye-tracking data for both left and right eyes using REMoDNaV, a velocity based eye movement event detection algorithm 
    that is based on, but extends the adaptive Nyström & Holmqvist algorithm (Nyström & Holmqvist, 2010).

    Parameters
    ----------
    min_pursuit_dur : float, optional
        Minimum pursuit duration in seconds for Remodnav detection (default is 10.0).
    max_pso_dur : float, optional
        Maximum post-saccadic oscillation duration in seconds for Remodnav detection (default is 0.0 -No PSO events detection-).
    min_fix_dur : float, optional
        Minimum fixation duration in seconds for Remodnav detection (default is 0.05).
    sac_max_vel : float, optional
        Maximum saccade velocity in deg/s (default is 1000.0).
    fix_max_amp : float, optional
        Maximum fixation amplitude in deg (default is 1.5).
    sac_time_thresh : float, optional
        Time threshold in seconds to consider a saccade as neighboring a fixation (default is 0.002).
    drop_fix_from_blink : bool, optional
        Whether to drop fixations that do not have a previous saccade within the time threshold (default is True).
    screen_size : float, optional
        Size of the screen in cm (default is 38.0).
    screen_width : int, optional
        Horizontal resolution of the screen in pixels (default is 1920).
    screen_distance : float, optional
        Distance from the screen to the participant's eyes in cm (default is 60).

    Returns
    -------
    fixations : dict
        Dictionary containing DataFrames of detected fixations for both left and right eyes with additional columns for mean x, y positions, and pupil size.
    saccades : dict
        Dictionary containing DataFrames of detected saccades for both left and right eyes.

    Raises
    ------
    ValueError
        If Remodnav detection fails.
    """

    # Move eye data, detections file and image to subject results directory
    self.out_folder.mkdir(parents=True, exist_ok=True)

    # Dictionaries to store fixations and saccades DataFrames
    fixations = []
    saccades = []

    # Divide samples by continuous data chunks
    # If the difference in the column "tSample" from row i to row i+1 is greater than 1000 / sample_rate, then it is a new chunk

    # Logic to find chunks
    chunks_start_indexes = np.where(np.diff(self.samples['tSample']) > (1000 / self.samples['Rate_recorded'][0]))[0] + 1

    chunks = np.zeros(len(self.samples))
    chunks[chunks_start_indexes] = 1
    chunks = np.cumsum(chunks)

    # Divide samples by chunks and apply to each chunk the detection algorithm
    for chunk in np.unique(chunks): 
        chunk_samples = self.samples[chunks == chunk].reset_index(drop=True)
        fixations_chunk, saccades_chunk = self.detect_on_chunk(chunk_samples, min_pursuit_dur, max_pso_dur, min_fix_dur, 
                        sac_max_vel, fix_max_amp, sac_time_thresh, drop_fix_from_blink, screen_size, screen_width, screen_distance,
                        savgol_length)
        fixations.append(fixations_chunk)
        saccades.append(saccades_chunk)

    return pd.concat(fixations).sort_values(by='tEnd', ignore_index=True), pd.concat(saccades).sort_values(by='tEnd', ignore_index=True)

run_eye_movement_from_samples(sample_rate, x_label='X', y_label='Y', config={}, **kwargs)

Recieves a pandas dataframe, a sample rate, and optional configuration :param dfSamples: Pandas dataframe including x and y columns :param sample_rate: an integer representing the data sample rate. :param x_label: X column name :param y_label: Y column name :param config:

Source code in pyxations/methods/eyemovement/REMoDNaV.py

def run_eye_movement_from_samples(self, sample_rate, x_label='X', y_label='Y', config={}, **kwargs):
    '''
    Recieves a pandas dataframe, a sample rate, and optional configuration
    :param dfSamples: Pandas dataframe including x and y columns
    :param sample_rate: an integer representing the data sample rate.
    :param x_label: X column name
    :param y_label: Y column name
    :param config:
    '''
    starting_time = self.samples['tSample'].min()
    eye_data = {
        'x': self.samples[x_label], 
        'y': self.samples[y_label]
    }
    eye_data = np.rec.fromarrays(list(eye_data.values()), names=list(eye_data.keys()))

    if 'pupil_data' not in config.keys():
        config['pupil_data'] = pd.Series([0]* len(eye_data['x']))
    times = np.arange(stop=len(eye_data['x'])) / sample_rate

    return self.run_eye_movement(eye_data['x'], eye_data['y'], sample_rate, 
        times=times, starting_time=starting_time, **config, **kwargs)

Engbert–Kliegl

Created on 5 nov 2024

@author: placiana

EngbertDetection

Bases: EyeMovementDetection

Python implementation for https://github.com/olafdimigen/eye-eeg/blob/master/detecteyemovements.m

Source code in pyxations/methods/eyemovement/engbert.py

class EngbertDetection(EyeMovementDetection):
    '''
    Python implementation for
    https://github.com/olafdimigen/eye-eeg/blob/master/detecteyemovements.m

    '''

    def __init__(self, session_folder_path, samples):
        self.session_folder_path = session_folder_path
        self.out_folder = (session_folder_path / 'engbert_events')
        self.samples = samples

    def detect_eye_movements(
            self,
            vfac: float = 5.0,
            mindur_ms: float = 6.0,
            smoothlevel: int = 1,
            globalthresh: bool = True,
            # deg/px: either give degperpixel OR let it be computed from screen params
            degperpixel: float | None = None,
            screen_size_cm: float = 38.0,
            screen_width_px: int = 1920,
            screen_distance_cm: float = 60.0,
            sample_rate_fallback: float | None = None,
        ):
            """
            Returns (fixations_df, saccades_df) with times in **ms**.
            Columns:
            Saccades: ['tStart','tEnd','duration','xStart','yStart','xEnd','yEnd','ampDeg','vPeak','distDeg','thetaDeg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
            Fixations: ['tStart','tEnd','duration','xAvg','yAvg','pupilAvg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
            """
            df = self.samples.copy()
            if degperpixel is None:
                degperpixel = _compute_px2deg(screen_size_cm, screen_distance_cm, screen_width_px)

            chunk_ids, fs_series = _split_into_chunks(df, fallback_fs=sample_rate_fallback)
            df['_chunk'] = chunk_ids
            df['_fs'] = fs_series.values

            sac_rows = []
            fix_rows = []

            eye_cols = _available_eye_columns(df)

            # Process each chunk & each available eye
            for chunk_id, g in df.groupby('_chunk', sort=True):
                fs = float(g['_fs'].iloc[0])
                t0_ms = float(g['tSample'].iloc[0])
                calib = g['Calib_index'].iloc[0] if 'Calib_index' in g else np.nan
                eyes_rec = g['Eyes_recorded'].iloc[0] if 'Eyes_recorded' in g else np.nan

                # choose which eye streams exist in this chunk
                streams = []
                if eye_cols['has_L'] and g[['LX','LY']].notna().any().any():
                    streams.append(('L', 'LX', 'LY', 'LPupil' if 'LPupil' in g else None))
                if eye_cols['has_R'] and g[['RX','RY']].notna().any().any():
                    streams.append(('R', 'RX', 'RY', 'RPupil' if 'RPupil' in g else None))
                if (not streams) and eye_cols['has_generic'] and g[['X','Y']].notna().any().any():
                    streams.append(('U', 'X', 'Y', 'Pupil' if 'Pupil' in g else None))  # U = unknown/unspecified eye

                if not streams:
                    continue  # nothing to do in this chunk

                # Prepare global thresholds (if requested)
                global_sigmas = {}
                if globalthresh:
                    for eye_label, cx, cy, _ in streams:
                        xy = g[[cx, cy]].to_numpy(dtype=float)
                        vel = vecvel(xy, fs, smoothlevel=smoothlevel)
                        sdx, sdy = velthresh(vel)
                        global_sigmas[eye_label] = (sdx, sdy)

                # Minimum duration in samples
                mindur_samples = max(1, int(round(mindur_ms * fs / 1000.0)))

                for eye_label, cx, cy, cp in streams:
                    xy = g[[cx, cy]].to_numpy(dtype=float)
                    # mark bad rows where both x & y are missing
                    valid = np.isfinite(xy).all(axis=1)
                    # If everything is invalid, skip
                    if not np.any(valid):
                        continue

                    vel = vecvel(xy, fs, smoothlevel=smoothlevel)

                    if globalthresh:
                        sdx, sdy = global_sigmas[eye_label]
                    else:
                        sdx, sdy = velthresh(vel)

                    # handle degenerate (all-NaN) sigmas
                    if not np.isfinite(sdx) or sdx == 0:
                        sdx = np.nanmedian(np.abs(vel[:,0])) or 1.0
                    if not np.isfinite(sdy) or sdy == 0:
                        sdy = np.nanmedian(np.abs(vel[:,1])) or 1.0

                    sac = microsacc_plugin(xy, vel, vfac=vfac, mindur_samples=mindur_samples, sdx=sdx, sdy=sdy)
                    if sac.size == 0:
                        # still may have one big fixation (whole chunk)
                        # synthesize a single fixation over all finite samples
                        idx = np.where(np.isfinite(xy[:,0]) & np.isfinite(xy[:,1]))[0]
                        if idx.size:
                            s_idx = int(idx[0]); e_idx = int(idx[-1])
                            # ms times aligned to t0_ms
                            tStart = t0_ms + (s_idx / fs) * 1000.0
                            tEnd   = t0_ms + (e_idx / fs) * 1000.0
                            dur    = (e_idx - s_idx + 1) / fs * 1000.0
                            xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                            yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                            pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                            fix_rows.append([tStart, tEnd, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])
                        continue

                    # Fill epoch column (not used here) and convert to ms/deg
                    # sac columns: [0 onset, 1 offset, 2 dur(samples), 3 avgV(px/s), 4 vPeak(px/s), 5 dist(px), 6 theta(rad), 7 amp(px), 8 dir(rad), 9 epoch, 10 x0, 11 y0, 12 x1, 13 y1]
                    onset_idx = sac[:,0].astype(int)
                    offset_idx = sac[:,1].astype(int)

                    tStart = t0_ms + (onset_idx / fs) * 1000.0
                    tEnd   = t0_ms + (offset_idx / fs) * 1000.0
                    dur_ms = (sac[:,2] / fs) * 1000.0

                    vPeak_deg = sac[:,4] * degperpixel            # px/s → deg/s
                    dist_deg  = sac[:,5] * degperpixel            # px → deg
                    amp_deg   = sac[:,7] * degperpixel            # px → deg
                    theta_deg = sac[:,6] * (180.0 / np.pi)

                    x0, y0, x1, y1 = sac[:,10], sac[:,11], sac[:,12], sac[:,13]

                    # Build saccade rows
                    for i in range(sac.shape[0]):
                        sac_rows.append([
                            float(tStart[i]), float(tEnd[i]), float(dur_ms[i]),
                            float(x0[i]), float(y0[i]), float(x1[i]), float(y1[i]),
                            float(amp_deg[i]), float(vPeak_deg[i]), float(dist_deg[i]), float(theta_deg[i]),
                            eye_label, calib, eyes_rec, fs, chunk_id
                        ])

                    # Build inter-saccadic fixations inside the chunk
                    # Sort events by onset
                    order = np.argsort(onset_idx)
                    onset_sorted = onset_idx[order]
                    offset_sorted = offset_idx[order]

                    # Leading fixation
                    if onset_sorted[0] > 0:
                        s_idx = 0
                        e_idx = onset_sorted[0] - 1
                        tS = t0_ms + (s_idx / fs) * 1000.0
                        tE = t0_ms + (e_idx / fs) * 1000.0
                        dur = (e_idx - s_idx + 1) / fs * 1000.0
                        xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                        yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                        pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                        fix_rows.append([tS, tE, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])

                    # Middle fixations
                    for k in range(len(onset_sorted) - 1):
                        s_idx = offset_sorted[k] + 1
                        e_idx = onset_sorted[k+1] - 1
                        if e_idx < s_idx:
                            continue
                        tS = t0_ms + (s_idx / fs) * 1000.0
                        tE = t0_ms + (e_idx / fs) * 1000.0
                        dur = (e_idx - s_idx + 1) / fs * 1000.0
                        xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                        yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                        pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                        fix_rows.append([tS, tE, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])

                    # Trailing fixation
                    last_off = int(offset_sorted[-1])
                    if last_off < (len(g)-1):
                        s_idx = last_off + 1
                        e_idx = len(g) - 1
                        tS = t0_ms + (s_idx / fs) * 1000.0
                        tE = t0_ms + (e_idx / fs) * 1000.0
                        dur = (e_idx - s_idx + 1) / fs * 1000.0
                        xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                        yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                        pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                        fix_rows.append([tS, tE, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])

            # Assemble DataFrames
            saccades = pd.DataFrame(
                sac_rows,
                columns=['tStart','tEnd','duration','xStart','yStart','xEnd','yEnd','ampDeg','vPeak','distDeg','thetaDeg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
            ).sort_values('tEnd', ignore_index=True)

            fixations = pd.DataFrame(
                fix_rows,
                columns=['tStart','tEnd','duration','xAvg','yAvg','pupilAvg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
            ).sort_values('tEnd', ignore_index=True)

            return fixations, saccades

detect_eye_movements(vfac=5.0, mindur_ms=6.0, smoothlevel=1, globalthresh=True, degperpixel=None, screen_size_cm=38.0, screen_width_px=1920, screen_distance_cm=60.0, sample_rate_fallback=None)

Returns (fixations_df, saccades_df) with times in ms. Columns: Saccades: ['tStart','tEnd','duration','xStart','yStart','xEnd','yEnd','ampDeg','vPeak','distDeg','thetaDeg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk'] Fixations: ['tStart','tEnd','duration','xAvg','yAvg','pupilAvg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']

Source code in pyxations/methods/eyemovement/engbert.py

def detect_eye_movements(
        self,
        vfac: float = 5.0,
        mindur_ms: float = 6.0,
        smoothlevel: int = 1,
        globalthresh: bool = True,
        # deg/px: either give degperpixel OR let it be computed from screen params
        degperpixel: float | None = None,
        screen_size_cm: float = 38.0,
        screen_width_px: int = 1920,
        screen_distance_cm: float = 60.0,
        sample_rate_fallback: float | None = None,
    ):
        """
        Returns (fixations_df, saccades_df) with times in **ms**.
        Columns:
        Saccades: ['tStart','tEnd','duration','xStart','yStart','xEnd','yEnd','ampDeg','vPeak','distDeg','thetaDeg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
        Fixations: ['tStart','tEnd','duration','xAvg','yAvg','pupilAvg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
        """
        df = self.samples.copy()
        if degperpixel is None:
            degperpixel = _compute_px2deg(screen_size_cm, screen_distance_cm, screen_width_px)

        chunk_ids, fs_series = _split_into_chunks(df, fallback_fs=sample_rate_fallback)
        df['_chunk'] = chunk_ids
        df['_fs'] = fs_series.values

        sac_rows = []
        fix_rows = []

        eye_cols = _available_eye_columns(df)

        # Process each chunk & each available eye
        for chunk_id, g in df.groupby('_chunk', sort=True):
            fs = float(g['_fs'].iloc[0])
            t0_ms = float(g['tSample'].iloc[0])
            calib = g['Calib_index'].iloc[0] if 'Calib_index' in g else np.nan
            eyes_rec = g['Eyes_recorded'].iloc[0] if 'Eyes_recorded' in g else np.nan

            # choose which eye streams exist in this chunk
            streams = []
            if eye_cols['has_L'] and g[['LX','LY']].notna().any().any():
                streams.append(('L', 'LX', 'LY', 'LPupil' if 'LPupil' in g else None))
            if eye_cols['has_R'] and g[['RX','RY']].notna().any().any():
                streams.append(('R', 'RX', 'RY', 'RPupil' if 'RPupil' in g else None))
            if (not streams) and eye_cols['has_generic'] and g[['X','Y']].notna().any().any():
                streams.append(('U', 'X', 'Y', 'Pupil' if 'Pupil' in g else None))  # U = unknown/unspecified eye

            if not streams:
                continue  # nothing to do in this chunk

            # Prepare global thresholds (if requested)
            global_sigmas = {}
            if globalthresh:
                for eye_label, cx, cy, _ in streams:
                    xy = g[[cx, cy]].to_numpy(dtype=float)
                    vel = vecvel(xy, fs, smoothlevel=smoothlevel)
                    sdx, sdy = velthresh(vel)
                    global_sigmas[eye_label] = (sdx, sdy)

            # Minimum duration in samples
            mindur_samples = max(1, int(round(mindur_ms * fs / 1000.0)))

            for eye_label, cx, cy, cp in streams:
                xy = g[[cx, cy]].to_numpy(dtype=float)
                # mark bad rows where both x & y are missing
                valid = np.isfinite(xy).all(axis=1)
                # If everything is invalid, skip
                if not np.any(valid):
                    continue

                vel = vecvel(xy, fs, smoothlevel=smoothlevel)

                if globalthresh:
                    sdx, sdy = global_sigmas[eye_label]
                else:
                    sdx, sdy = velthresh(vel)

                # handle degenerate (all-NaN) sigmas
                if not np.isfinite(sdx) or sdx == 0:
                    sdx = np.nanmedian(np.abs(vel[:,0])) or 1.0
                if not np.isfinite(sdy) or sdy == 0:
                    sdy = np.nanmedian(np.abs(vel[:,1])) or 1.0

                sac = microsacc_plugin(xy, vel, vfac=vfac, mindur_samples=mindur_samples, sdx=sdx, sdy=sdy)
                if sac.size == 0:
                    # still may have one big fixation (whole chunk)
                    # synthesize a single fixation over all finite samples
                    idx = np.where(np.isfinite(xy[:,0]) & np.isfinite(xy[:,1]))[0]
                    if idx.size:
                        s_idx = int(idx[0]); e_idx = int(idx[-1])
                        # ms times aligned to t0_ms
                        tStart = t0_ms + (s_idx / fs) * 1000.0
                        tEnd   = t0_ms + (e_idx / fs) * 1000.0
                        dur    = (e_idx - s_idx + 1) / fs * 1000.0
                        xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                        yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                        pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                        fix_rows.append([tStart, tEnd, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])
                    continue

                # Fill epoch column (not used here) and convert to ms/deg
                # sac columns: [0 onset, 1 offset, 2 dur(samples), 3 avgV(px/s), 4 vPeak(px/s), 5 dist(px), 6 theta(rad), 7 amp(px), 8 dir(rad), 9 epoch, 10 x0, 11 y0, 12 x1, 13 y1]
                onset_idx = sac[:,0].astype(int)
                offset_idx = sac[:,1].astype(int)

                tStart = t0_ms + (onset_idx / fs) * 1000.0
                tEnd   = t0_ms + (offset_idx / fs) * 1000.0
                dur_ms = (sac[:,2] / fs) * 1000.0

                vPeak_deg = sac[:,4] * degperpixel            # px/s → deg/s
                dist_deg  = sac[:,5] * degperpixel            # px → deg
                amp_deg   = sac[:,7] * degperpixel            # px → deg
                theta_deg = sac[:,6] * (180.0 / np.pi)

                x0, y0, x1, y1 = sac[:,10], sac[:,11], sac[:,12], sac[:,13]

                # Build saccade rows
                for i in range(sac.shape[0]):
                    sac_rows.append([
                        float(tStart[i]), float(tEnd[i]), float(dur_ms[i]),
                        float(x0[i]), float(y0[i]), float(x1[i]), float(y1[i]),
                        float(amp_deg[i]), float(vPeak_deg[i]), float(dist_deg[i]), float(theta_deg[i]),
                        eye_label, calib, eyes_rec, fs, chunk_id
                    ])

                # Build inter-saccadic fixations inside the chunk
                # Sort events by onset
                order = np.argsort(onset_idx)
                onset_sorted = onset_idx[order]
                offset_sorted = offset_idx[order]

                # Leading fixation
                if onset_sorted[0] > 0:
                    s_idx = 0
                    e_idx = onset_sorted[0] - 1
                    tS = t0_ms + (s_idx / fs) * 1000.0
                    tE = t0_ms + (e_idx / fs) * 1000.0
                    dur = (e_idx - s_idx + 1) / fs * 1000.0
                    xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                    yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                    pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                    fix_rows.append([tS, tE, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])

                # Middle fixations
                for k in range(len(onset_sorted) - 1):
                    s_idx = offset_sorted[k] + 1
                    e_idx = onset_sorted[k+1] - 1
                    if e_idx < s_idx:
                        continue
                    tS = t0_ms + (s_idx / fs) * 1000.0
                    tE = t0_ms + (e_idx / fs) * 1000.0
                    dur = (e_idx - s_idx + 1) / fs * 1000.0
                    xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                    yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                    pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                    fix_rows.append([tS, tE, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])

                # Trailing fixation
                last_off = int(offset_sorted[-1])
                if last_off < (len(g)-1):
                    s_idx = last_off + 1
                    e_idx = len(g) - 1
                    tS = t0_ms + (s_idx / fs) * 1000.0
                    tE = t0_ms + (e_idx / fs) * 1000.0
                    dur = (e_idx - s_idx + 1) / fs * 1000.0
                    xAvg = float(np.nanmean(xy[s_idx:e_idx+1,0]))
                    yAvg = float(np.nanmean(xy[s_idx:e_idx+1,1]))
                    pupilAvg = float(np.nanmean(g[cp].values[s_idx:e_idx+1])) if cp in g else np.nan
                    fix_rows.append([tS, tE, dur, xAvg, yAvg, pupilAvg, eye_label, calib, eyes_rec, fs, chunk_id])

        # Assemble DataFrames
        saccades = pd.DataFrame(
            sac_rows,
            columns=['tStart','tEnd','duration','xStart','yStart','xEnd','yEnd','ampDeg','vPeak','distDeg','thetaDeg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
        ).sort_values('tEnd', ignore_index=True)

        fixations = pd.DataFrame(
            fix_rows,
            columns=['tStart','tEnd','duration','xAvg','yAvg','pupilAvg','eye','Calib_index','Eyes_recorded','Rate_recorded','chunk']
        ).sort_values('tEnd', ignore_index=True)

        return fixations, saccades

microsacc_plugin(pos_xy, vel_xy, vfac, mindur_samples, sdx, sdy)

Return saccades array with columns: 0 onset, 1 offset, 2 dur, 3 avgvel(px/s), 4 vpeak(px/s), 5 dist(px), 6 theta(rad), 7 amp(px), 8 dir(rad), 9 epoch (filled outside), 10 x0, 11 y0, 12 x1, 13 y1

Source code in pyxations/methods/eyemovement/engbert.py

def microsacc_plugin(pos_xy, vel_xy, vfac, mindur_samples, sdx, sdy):
    """
    Return saccades array with columns:
    0 onset, 1 offset, 2 dur, 3 avgvel(px/s), 4 vpeak(px/s), 5 dist(px),
    6 theta(rad), 7 amp(px), 8 dir(rad), 9 epoch (filled outside),
    10 x0, 11 y0, 12 x1, 13 y1
    """
    vx, vy = vel_xy[:,0], vel_xy[:,1]
    with np.errstate(invalid='ignore'):
        crit = (vx / sdx)**2 + (vy / sdy)**2
    suprath = crit > (vfac**2)
    runs = _find_runs(suprath)

    sac = []
    for s0, s1 in runs:
        if (s1 - s0 + 1) < mindur_samples:
            continue
        seg_v = np.hypot(vx[s0:s1+1], vy[s0:s1+1])
        vpeak = np.nanmax(seg_v)
        avgv  = np.nanmean(seg_v)

        x0, y0 = pos_xy[s0,0], pos_xy[s0,1]
        x1, y1 = pos_xy[s1,0], pos_xy[s1,1]
        amp = np.hypot(x1 - x0, y1 - y0)
        theta = np.arctan2(y1 - y0, x1 - x0)  # main vector angle

        seg = pos_xy[s0:s1+1]
        dist = np.nansum(np.hypot(np.diff(seg[:,0]), np.diff(seg[:,1])))

        sac.append([s0, s1, (s1 - s0 + 1), avgv, vpeak, dist, theta, amp, theta, np.nan, x0, y0, x1, y1])
    return np.array(sac, dtype=float)

vecvel(gaze_xy, fs, smoothlevel=1)

gaze_xy: (N,2) positions in pixels (NaN for missing) returns velocities (N,2) in px/s using central differences (+ optional smoothing of positions)

Source code in pyxations/methods/eyemovement/engbert.py

def vecvel(gaze_xy, fs, smoothlevel=1):
    """
    gaze_xy: (N,2) positions in pixels (NaN for missing)
    returns velocities (N,2) in px/s using central differences (+ optional smoothing of positions)
    """
    x = _smooth_1d(gaze_xy[:,0].astype(float), smoothlevel)
    y = _smooth_1d(gaze_xy[:,1].astype(float), smoothlevel)

    vx = np.empty_like(x); vy = np.empty_like(y)
    vx[1:-1] = (x[2:] - x[:-2]) * (fs / 2.0)
    vy[1:-1] = (y[2:] - y[:-2]) * (fs / 2.0)
    vx[0] = (x[1] - x[0]) * fs
    vy[0] = (y[1] - y[0]) * fs
    vx[-1] = (x[-1] - x[-2]) * fs
    vy[-1] = (y[-1] - y[-2]) * fs
    vx[~np.isfinite(vx)] = np.nan
    vy[~np.isfinite(vy)] = np.nan
    return np.column_stack([vx, vy])

analysis

High-level experiment objects for loading and iterating over derivatives.

generic

Experiment

Source code in pyxations/analysis/generic.py

class Experiment:

    def __init__(self, dataset_path: str, excluded_subjects: list = [], excluded_sessions: dict = {}, excluded_trials: dict = {}, export_format = FEATHER_EXPORT):
        self.dataset_path = Path(dataset_path)
        self.derivatives_path = self.dataset_path.with_name(self.dataset_path.name + "_derivatives")
        self.metadata = pl.read_csv(self.dataset_path / "participants.tsv", separator="\t", 
                                    schema_overrides={"subject_id": pl.Utf8, "old_subject_id": pl.Utf8})
        self.subjects = { subject_id:
            Subject(subject_id, old_subject_id, self, 
                     excluded_sessions.get(subject_id, []), excluded_trials.get(subject_id, {}),export_format)
            for subject_id, old_subject_id in zip(self.metadata.select("subject_id").to_series(),
                                                  self.metadata.select("old_subject_id").to_series())
            if subject_id not in excluded_subjects and old_subject_id not in excluded_subjects
        }
        self.export_format = export_format

    def __iter__(self):
        return iter(self.subjects)

    def __getitem__(self, index):
        return self.subjects[index]

    def __len__(self):
        return len(self.subjects)

    def __repr__(self):
        return f"Experiment = '{self.dataset_path.name}'"

    def __next__(self):
        return next(self.subjects)

    def load_data(self, detection_algorithm: str):
        self.detection_algorithm = detection_algorithm
        for subject in self.subjects.values():
            subject.load_data(detection_algorithm)

    def plot_multipanel(self, display: bool):
        fixations = pl.concat([subject.fixations() for subject in self.subjects.values()])
        saccades = pl.concat([subject.saccades() for subject in self.subjects.values()])

        vis = Visualization(self.derivatives_path, self.detection_algorithm)
        vis.plot_multipanel(fixations, saccades, display)

    def filter_fixations(self, min_fix_dur=50, print_flag=True):
        amount_fix = self.fixations().shape[0]
        for subject in self.subjects.values():
            subject.filter_fixations(min_fix_dur)

        if print_flag:
            print(f"Removed {amount_fix - self.fixations().shape[0]} fixations shorter than {min_fix_dur} ms.")
    def collapse_fixations(self, threshold_px: float, print_flag=True):
        amount_fix = self.fixations().shape[0]
        for subject in self.subjects.values():
            subject.collapse_fixations(threshold_px)
        if print_flag:
            print(f"Removed {amount_fix - self.fixations().shape[0]} fixations that were merged.")

    def drop_trials_with_nan_threshold(self, phase, threshold=0.1,print_flag=True):
        amount_trials_total = self.rts().shape[0]
        for subject in list(self.subjects.values()):
            subject.drop_trials_with_nan_threshold(phase,threshold,False)
        if print_flag:
            print(f"Removed {amount_trials_total - self.rts().shape[0]} trials with NaN values.")

    def drop_trials_longer_than(self, seconds,phase, print_flag=True):
        amount_trials_total = self.rts().shape[0]
        for subject in list(self.subjects.values()):
            subject.drop_trials_longer_than(seconds,phase,False)
        if print_flag:
            print(f"Removed {amount_trials_total - self.rts().shape[0]} trials longer than {seconds} seconds.")

    def plot_scanpaths(self,screen_height,screen_width,display: bool = False):
        with ProcessPoolExecutor(8) as executor:
            futures = [executor.submit(subject.plot_scanpaths,screen_height,screen_width,display) for subject in self.subjects.values()]
            for future in as_completed(futures):
                future.result()

    def drop_poor_or_non_calibrated_trials(self, threshold=1.0, print_flag=True):
        '''
        Drop trials that are not calibrated or have a poor calibration.
        A trial is considered not calibrated if there is no validation data for its calibration index.
        A trial is considered poorly calibrated if the average error is greater than the threshold.
        '''
        amount_trials_total = self.rts().shape[0]
        for subject in list(self.subjects.values()):
            subject.drop_poor_or_non_calibrated_trials(threshold,False)
        if print_flag:
            print(f"Removed {amount_trials_total - self.rts().shape[0]} trials with poor calibration.")

    def rts(self):
        rts = [subject.rts() for subject in self.subjects.values()]
        return pl.concat(rts)

    def get_subject(self, subject_id):
        return self.subjects[subject_id]

    def get_session(self, subject_id, session_id):
        subject = self.get_subject(subject_id)
        return subject.get_session(session_id)

    def get_trial(self, subject_id, session_id, trial_number):
        session = self.get_session(subject_id, session_id)
        return session.get_trial(trial_number)

    def fixations(self):
        return pl.concat([subject.fixations() for subject in self.subjects.values()])

    def saccades(self):
        return pl.concat([subject.saccades() for subject in self.subjects.values()])

    def samples(self):
        return pl.concat([subject.samples() for subject in self.subjects.values()])

    def remove_subject(self, subject_id):
        if subject_id in self.subjects:
            del self.subjects[subject_id]

    def calib_data(self):
        calib_data = [subject.calib_data() for subject in self.subjects.values()]
        calib_indexes = pl.concat([calib_data[1] for calib_data in calib_data])
        calib_data = pl.concat([calib_data[0] for calib_data in calib_data])
        return calib_data, calib_indexes

    def plot_calib_data(self):
        # Step 0: Load and preprocess
        calib_data = self.calib_data()
        trial_numbers = calib_data[1]
        calib_data = calib_data[0].select([
            "subject_id", "session_id", "Calib_index", "eye", "avg_error", "validation_id"
        ])

        # Step 1: Get only rows with max validation_id per group
        max_vals = (
            calib_data
            .group_by(["subject_id", "session_id", "Calib_index", "eye"])
            .agg(pl.col("validation_id").max().alias("max_validation_id"))
        )

        calib_data = (
            calib_data
            .join(max_vals, on=["subject_id", "session_id", "Calib_index", "eye"])
            .filter(pl.col("validation_id") == pl.col("max_validation_id"))
            .drop(["max_validation_id", "validation_id"])
        )

        # Step 2: Choose best eye (lowest avg_error) per calibration
        best_eyes = (
            calib_data
            .group_by(["subject_id", "session_id", "Calib_index"])
            .agg(pl.col("avg_error").min().alias("best_eye_error"))
        )

        calib_data = (
            calib_data
            .join(best_eyes, on=["subject_id", "session_id", "Calib_index"])
            .filter(pl.col("avg_error") == pl.col("best_eye_error"))
            .drop(["eye", "best_eye_error"])
        )

        # Step 3: Add trial number and clean up
        calib_data = (
            calib_data
            .join(trial_numbers, on=["subject_id", "session_id", "Calib_index"],how="right")
            .drop("Calib_index")
        )
        # Replace nans in avg_error with -1
        calib_data = calib_data.with_columns(
            pl.when(pl.col("avg_error").is_null()).then(-1).otherwise(pl.col("avg_error")).alias("avg_error")
        )

        # Step 4: Combine the columns "subject_id" and "session_id" into a single column
        calib_data = (
            calib_data
            .with_columns(
                (pl.col("subject_id").cast(pl.Utf8) + "_" + pl.col("session_id").cast(pl.Utf8)).alias("subject_id")
            )
            .drop("session_id")
        )
        # Create a copy of the colormap and set 'under' color for -1s
        cmap = colormaps["rocket_r"].copy()
        cmap.set_under("yellow")  # or any color you prefer, e.g., "black", "white"

        heatmap_data = (
            calib_data
            .pivot(
                values="avg_error",
                index="subject_id",
                on="trial_number",
                aggregate_function="first"  # safe if unique per cell
            )
            .sort("subject_id")
            .to_pandas()
            .set_index("subject_id")
        )
        heatmap_data = heatmap_data[sorted(heatmap_data.columns, key=lambda x: int(x))]

        # Step 5: Plot with adaptive sizing
        n_subjects = heatmap_data.shape[0]
        n_trials = heatmap_data.shape[1]

        # Define a base size per cell, then scale it
        cell_width = 0.5   # width per trial column
        cell_height = 0.2  # height per subject row

        # Limit extremes so it doesn’t explode with huge data
        fig_width = max(10, min(cell_width * n_trials, 40))
        fig_height = max(8, min(cell_height * n_subjects, 40))

        plt.figure(figsize=(fig_width, fig_height))
        sns.heatmap(
            heatmap_data,
            cmap=cmap,
            center=0.5,
            vmin=0,
            linewidths=0.3,
            linecolor="grey",
            cbar_kws=dict(label="Avg. error (°)")
        )
        plt.xlabel("Trial #", fontsize=14)
        plt.ylabel("Subject", fontsize=14)
        plt.title("Calibration Error per Subject and Trial", fontsize=16)

        # Rotate labels
        plt.xticks(rotation=45, ha="right", fontsize=10)
        plt.yticks(rotation=0, ha="right", va="center", fontsize=10)  # horizontal y labels

        plt.tight_layout()
        plt.show()
        plt.close()

drop_poor_or_non_calibrated_trials(threshold=1.0, print_flag=True)

Drop trials that are not calibrated or have a poor calibration. A trial is considered not calibrated if there is no validation data for its calibration index. A trial is considered poorly calibrated if the average error is greater than the threshold.

Source code in pyxations/analysis/generic.py

def drop_poor_or_non_calibrated_trials(self, threshold=1.0, print_flag=True):
    '''
    Drop trials that are not calibrated or have a poor calibration.
    A trial is considered not calibrated if there is no validation data for its calibration index.
    A trial is considered poorly calibrated if the average error is greater than the threshold.
    '''
    amount_trials_total = self.rts().shape[0]
    for subject in list(self.subjects.values()):
        subject.drop_poor_or_non_calibrated_trials(threshold,False)
    if print_flag:
        print(f"Removed {amount_trials_total - self.rts().shape[0]} trials with poor calibration.")

Session

Source code in pyxations/analysis/generic.py

class Session():

    def __init__(self, session_id: str, subject: Subject, excluded_trials: list = [],export_format = FEATHER_EXPORT):
        self.session_id = session_id
        self.subject = weakref.ref(subject)
        self.excluded_trials = excluded_trials
        self.session_dataset_path = self.subject().subject_dataset_path / f"ses-{self.session_id}"
        self.session_derivatives_path = self.subject().subject_derivatives_path / f"ses-{self.session_id}"
        self._trials = None  # Lazy load trials
        self.export_format = export_format

        if not self.session_derivatives_path.exists():
            raise FileNotFoundError(f"Session path not found: {self.session_derivatives_path}")

    def __getstate__(self):
        state = self.__dict__.copy()
        sess_parent = state.pop("subject", None)
        state["_subject_id"] = id(sess_parent()) if sess_parent else None
        return state

    def __setstate__(self, state):
        self.__dict__.update(state)
        self.subject = lambda: None

    @property
    def trials(self):
        if self._trials is None:
            raise ValueError("Trials not loaded. Please load data first.")
        return self._trials

    def __repr__(self):
        return f"Session = '{self.session_id}', " + self.subject().__repr__()

    def drop_trials_with_nan_threshold(self, phase, threshold=0.1, print_flag=True):
        total_trials = len(self.trials)
        # Filter bad trials

        bad_trials = [trial for trial in self.trials.keys() if self.trials[trial].is_trial_bad(phase, threshold)]
        if len(bad_trials)/total_trials > threshold:
            self.subject().remove_session(self.session_id)


        if print_flag:
            print(f"Removed {len(bad_trials)} trials with NaN values.")

    def drop_poor_or_non_calibrated_trials(self, threshold=1.0, print_flag=True):
        '''
        Drop trials that are not calibrated or have a poor calibration.
        A trial is considered not calibrated if there is no validation data for its calibration index.
        A trial is considered poorly calibrated if the average error is greater than the threshold.
        '''
        trial_numbers = [trial for trial in self.trials.keys()]
        # Step 1: Get only rows with max validation_id per group
        calib_data, trial_numbers = self.calib_data()
        calib_data = calib_data.drop("session_id")
        max_vals = (
            calib_data
            .group_by(["Calib_index", "eye"])
            .agg(pl.col("validation_id").max().alias("max_validation_id"))
        )

        calib_data = (
            calib_data
            .join(max_vals, on=["Calib_index", "eye"])
            .filter(pl.col("validation_id") == pl.col("max_validation_id"))
            .drop(["max_validation_id", "validation_id"])
        )

        # Step 2: Choose best eye (lowest avg_error) per calibration
        best_eyes = (
            calib_data
            .group_by(["Calib_index"])
            .agg(pl.col("avg_error").min().alias("best_eye_error"))
        )

        calib_data = (
            calib_data
            .join(best_eyes, on=["Calib_index"])
            .filter(pl.col("avg_error") == pl.col("best_eye_error"))
            .drop(["eye", "best_eye_error"])
        )

        calib_data = (
            calib_data
            .join(trial_numbers, on=["Calib_index"],how="right")
            .drop("Calib_index")
        )
        # Bad trials are those with avg_error > threshold, or those that have NaN values in avg_error
        bad_trials = calib_data.filter((pl.col("avg_error") > threshold) | (pl.col("avg_error").is_null())).select("trial_number").to_series().unique().to_list()

        for trial in bad_trials:
            self.remove_trial(trial)

        if print_flag:
            print(f"Removed {len(bad_trials)} trials with poor calibration.")

    def drop_trials_longer_than(self, seconds,phase, print_flag=True):

        # Filter bad trials

        bad_trials = [trial for trial in self.trials.keys() if self.trials[trial].is_trial_longer_than(seconds,phase)]
        for trial in bad_trials:
            self.remove_trial(trial)

        if print_flag:
            print(f"Removed {len(bad_trials)} trials longer than {seconds} seconds.")

    def load_behavior_data(self):
        # This should be implemented for each type of experiment
        pass

    def load_data(self, detection_algorithm: str):
        self.detection_algorithm = detection_algorithm
        events_path = self.session_derivatives_path / f"{self.detection_algorithm}_events"


        exporter = get_exporter(self.export_format)
        file_extension = exporter.extension()


        # Check paths and load files efficiently

        samples = exporter.read(self.session_derivatives_path, 'samples')
        fix = exporter.read(events_path, 'fix')
        sacc = exporter.read(events_path, 'sacc')
        blink = exporter.read(events_path, "blink") if (events_path / ("blink" + file_extension)).exists() else None
        self._calib_data = _parse_validations(exporter.read(self.session_derivatives_path, "calib")) if (self.session_derivatives_path / ("calib" + file_extension)).exists() else None

        # Initialize trials
        self._init_trials(samples,fix,sacc,blink,events_path)

    def calib_data(self):
        if self._calib_data is None:
            raise ValueError(f"Calibration data for session {self.session_id} and subject {self.subject().subject_id} not loaded. Please load data first.")

        calib_indexes = [(trial.trial_number,trial.calib_index) for trial in self.trials.values() if trial.calib_index is not None]
        calib_indexes = pl.DataFrame(calib_indexes, schema=["trial_number","Calib_index"],orient="row").with_columns([
            (pl.lit(self.session_id)).alias("session_id")])
        return self._calib_data.with_columns([
            (pl.lit(self.session_id)).alias("session_id")]), calib_indexes

    def _init_trials(self,samples,fix,sacc,blink,events_path):
        cosas = [trial for trial in samples.select("trial_number").to_series().unique() if trial != -1 and trial not in self.excluded_trials]
        self._trials = {trial:
            Trial(trial, self, samples, fix, sacc, blink, events_path)
            for trial in cosas
        } 

    def plot_scanpaths(self,screen_height,screen_width, display: bool = False):
        for trial in self.trials.values():
            trial.plot_scanpath(screen_height,screen_width,display=display)

    def __iter__(self):
        return iter(self.trials)

    def __getitem__(self, index):
        return self.trials[index]

    def __len__(self):
        return len(self.trials)

    def get_trial(self, trial_number):
        return self._trials[trial_number]

    def filter_fixations(self, min_fix_dur=50):
        for trial in self.trials.values():
            trial.filter_fixations(min_fix_dur)

    def collapse_fixations(self, threshold_px: float):
        for trial in self.trials.values():
            trial.collapse_fixations(threshold_px)


    def rts(self):
        rts = [trial.rts() for trial in self.trials.values()]
        rts = pl.concat(rts).with_columns([
            (pl.lit(self.session_id)).alias("session_id"),])
        return rts


    def fixations(self):
        df = pl.concat([trial.fixations() for trial in self.trials.values()]).with_columns([
            (pl.lit(self.session_id)).alias("session_id"),])
        return df

    def saccades(self):
        df = pl.concat([trial.saccades() for trial in self.trials.values()]).with_columns([
            (pl.lit(self.session_id)).alias("session_id"),])
        return df


    def samples(self):
        df = pl.concat([trial.samples() for trial in self.trials.values()]).with_columns([
            (pl.lit(self.session_id)).alias("session_id"),])
        return df

    def remove_trial(self, trial_number):
        if self._trials and trial_number in self._trials:
            del self._trials[trial_number]
            if len(self._trials) == 0:
                subj = self.subject()
                if subj:
                    subj.remove_session(self.session_id)
                self._trials = None
                self.subject = lambda: None

drop_poor_or_non_calibrated_trials(threshold=1.0, print_flag=True)

Drop trials that are not calibrated or have a poor calibration. A trial is considered not calibrated if there is no validation data for its calibration index. A trial is considered poorly calibrated if the average error is greater than the threshold.

Source code in pyxations/analysis/generic.py

def drop_poor_or_non_calibrated_trials(self, threshold=1.0, print_flag=True):
    '''
    Drop trials that are not calibrated or have a poor calibration.
    A trial is considered not calibrated if there is no validation data for its calibration index.
    A trial is considered poorly calibrated if the average error is greater than the threshold.
    '''
    trial_numbers = [trial for trial in self.trials.keys()]
    # Step 1: Get only rows with max validation_id per group
    calib_data, trial_numbers = self.calib_data()
    calib_data = calib_data.drop("session_id")
    max_vals = (
        calib_data
        .group_by(["Calib_index", "eye"])
        .agg(pl.col("validation_id").max().alias("max_validation_id"))
    )

    calib_data = (
        calib_data
        .join(max_vals, on=["Calib_index", "eye"])
        .filter(pl.col("validation_id") == pl.col("max_validation_id"))
        .drop(["max_validation_id", "validation_id"])
    )

    # Step 2: Choose best eye (lowest avg_error) per calibration
    best_eyes = (
        calib_data
        .group_by(["Calib_index"])
        .agg(pl.col("avg_error").min().alias("best_eye_error"))
    )

    calib_data = (
        calib_data
        .join(best_eyes, on=["Calib_index"])
        .filter(pl.col("avg_error") == pl.col("best_eye_error"))
        .drop(["eye", "best_eye_error"])
    )

    calib_data = (
        calib_data
        .join(trial_numbers, on=["Calib_index"],how="right")
        .drop("Calib_index")
    )
    # Bad trials are those with avg_error > threshold, or those that have NaN values in avg_error
    bad_trials = calib_data.filter((pl.col("avg_error") > threshold) | (pl.col("avg_error").is_null())).select("trial_number").to_series().unique().to_list()

    for trial in bad_trials:
        self.remove_trial(trial)

    if print_flag:
        print(f"Removed {len(bad_trials)} trials with poor calibration.")

Subject

Source code in pyxations/analysis/generic.py

class Subject:

    def __init__(self, subject_id: str, old_subject_id: str, experiment: Experiment,
                 excluded_sessions: list = [], excluded_trials: dict = {}, export_format = FEATHER_EXPORT):
        self.subject_id = subject_id
        self.old_subject_id = old_subject_id
        self.experiment = weakref.ref(experiment)
        self._sessions = None  # Lazy load sessions
        self.excluded_sessions = excluded_sessions
        self.excluded_trials = excluded_trials
        self.subject_dataset_path = self.experiment().dataset_path / f"sub-{self.subject_id}"
        self.subject_derivatives_path = self.experiment().derivatives_path / f"sub-{self.subject_id}"
        self.export_format = export_format

    def __getstate__(self):
        state = self.__dict__.copy()
        # Save the parent *id* instead of the weakref itself
        exp = state.pop("experiment", None)
        state["_experiment_id"] = id(exp()) if exp else None
        return state

    def __setstate__(self, state):
        self.__dict__.update(state)
        # In a worker the real Experiment instance isn’t available
        # – keep a placeholder or rebuild if you can look it up.
        self.experiment = lambda: None      # callable that returns None

    @property
    def sessions(self):
        if self._sessions is None:
            self._sessions = { session_folder.name.split("-")[-1] :
                Session(session_folder.name.split("-")[-1], self,
                        self.excluded_trials.get(session_folder.name.split("-")[-1], {}),self.export_format) 
                for session_folder in self.subject_derivatives_path.glob("ses-*") 
                if session_folder.name.split("-")[-1] not in self.excluded_sessions
            }
        return self._sessions

    def __iter__(self):
        return iter(self.sessions)

    def __getitem__(self, index):
        return self.sessions[index]

    def __len__(self):
        return len(self.sessions)

    def __repr__(self):
        return f"Subject = '{self.subject_id}', " + self.experiment().__repr__()

    def __next__(self):
        return next(self.sessions)


    def remove_session(self, session_id):
        if self._sessions and session_id in self._sessions:
            del self._sessions[session_id]
            if len(self._sessions) == 0:
                exp = self.experiment()
                if exp:
                    exp.remove_subject(self.subject_id)
                self._sessions = None
                self.experiment = lambda: None

    def load_data(self, detection_algorithm: str):
        self.detection_algorithm = detection_algorithm
        for session in self.sessions.values():
            session.load_data(detection_algorithm)


    def filter_fixations(self, min_fix_dur=50):
        for session in self.sessions.values():
            session.filter_fixations(min_fix_dur)

    def collapse_fixations(self, threshold_px: float):
        for session in self.sessions.values():
            session.collapse_fixations(threshold_px)

    def drop_trials_with_nan_threshold(self,phase, threshold=0.1, print_flag=True):
        total_sessions = len(self.sessions)
        amount_trials_total = self.rts().shape[0]
        for session in list(self.sessions.values()):
            session.drop_trials_with_nan_threshold(phase,threshold,False)
        bad_sessions_count = total_sessions - len(self.sessions)


        # If the proportion of bad sessions exceeds the threshold, remove all sessions
        if bad_sessions_count / total_sessions > threshold:
            self.experiment().remove_subject(self.subject_id)

        if print_flag:
            print(f"Removed {amount_trials_total - self.rts().shape[0]} trials with NaN values.")

    def drop_poor_or_non_calibrated_trials(self, threshold=1.0, print_flag=True):
        '''
        Drop trials that are not calibrated or have a poor calibration.
        A trial is considered not calibrated if there is no validation data for its calibration index.
        A trial is considered poorly calibrated if the average error is greater than the threshold.
        '''
        amount_trials_total = self.rts().shape[0]
        for session in list(self.sessions.values()):
            session.drop_poor_or_non_calibrated_trials(threshold,False)
        if print_flag:
            print(f"Removed {amount_trials_total - self.rts().shape[0]} trials with poor calibration.")

    def drop_trials_longer_than(self, seconds,phase, print_flag=True):
        amount_trials_total = self.rts().shape[0]
        for session in list(self.sessions.values()):
            session.drop_trials_longer_than(seconds,phase,False)
        if print_flag:
            print(f"Removed {amount_trials_total - self.rts().shape[0]} trials longer than {seconds} seconds.")

    def plot_scanpaths(self,screen_height,screen_width, display: bool = False):
        for session in self.sessions.values():
            session.plot_scanpaths(screen_height,screen_width,display)

    def rts(self):
        rts = [session.rts() for session in self.sessions.values()]
        rts = pl.concat(rts).with_columns([
            (pl.lit(self.subject_id)).alias("subject_id"),])
        return rts

    def get_session(self, session_id):
        return self.sessions[session_id]

    def get_trial(self, session_id, trial_number):
        session = self.get_session(session_id)
        return session.get_trial(trial_number)


    def fixations(self):
        df = pl.concat([session.fixations() for session in self.sessions.values()]).with_columns([
            (pl.lit(self.subject_id)).alias("subject_id"),])
        return df

    def saccades(self):
        df = pl.concat([session.saccades() for session in self.sessions.values()]).with_columns([
            (pl.lit(self.subject_id)).alias("subject_id"),])
        return df

    def samples(self):
        df = pl.concat([session.samples() for session in self.sessions.values()]).with_columns([
            (pl.lit(self.subject_id)).alias("subject_id"),])
        return df

    def calib_data(self):
        calib_data = [session.calib_data() for session in self.sessions.values()]
        calib_indexes = pl.concat([calib_data[1] for calib_data in calib_data]).with_columns([
            (pl.lit(self.subject_id)).alias("subject_id"),])
        calib_data = pl.concat([calib_data[0] for calib_data in calib_data]).with_columns([
            (pl.lit(self.subject_id)).alias("subject_id"),])
        return calib_data, calib_indexes

drop_poor_or_non_calibrated_trials(threshold=1.0, print_flag=True)

Drop trials that are not calibrated or have a poor calibration. A trial is considered not calibrated if there is no validation data for its calibration index. A trial is considered poorly calibrated if the average error is greater than the threshold.

Source code in pyxations/analysis/generic.py

def drop_poor_or_non_calibrated_trials(self, threshold=1.0, print_flag=True):
    '''
    Drop trials that are not calibrated or have a poor calibration.
    A trial is considered not calibrated if there is no validation data for its calibration index.
    A trial is considered poorly calibrated if the average error is greater than the threshold.
    '''
    amount_trials_total = self.rts().shape[0]
    for session in list(self.sessions.values()):
        session.drop_poor_or_non_calibrated_trials(threshold,False)
    if print_flag:
        print(f"Removed {amount_trials_total - self.rts().shape[0]} trials with poor calibration.")

Trial

Source code in pyxations/analysis/generic.py

class Trial:

    def __init__(self, trial_number: int, session: Session, samples: pl.DataFrame, fix: pl.DataFrame, 
                sacc: pl.DataFrame, blink: pl.DataFrame, events_path: Path):
        self.trial_number = trial_number
        self.session = session

        # Filter per trial
        # If "Calib_index" is a column in samples, set self._calib_index to the value of that column

        self._calib_index = samples.filter(pl.col("trial_number") == trial_number).select("Calib_index").to_series()[0] if "Calib_index" in samples.columns else None


        self._samples = samples.filter(pl.col("trial_number") == trial_number).drop("Calib_index",strict=False)
        self._fix = fix.filter(pl.col("trial_number") == trial_number).drop("Calib_index",strict=False)
        self._sacc = sacc.filter(pl.col("trial_number") == trial_number).drop("Calib_index",strict=False)
        self._blink = blink.filter(pl.col("trial_number") == trial_number).drop("Calib_index",strict=False) if blink is not None else None

        # Get the start time
        start_time = self._samples.select("tSample").to_series()[0]

        # Time normalization
        self._samples = self._samples.with_columns([
            (pl.col("tSample") - start_time).alias("tSample")
        ])

        self._fix = self._fix.with_columns([
            (pl.col("tStart") - start_time).alias("tStart"),
            (pl.col("tEnd") - start_time).alias("tEnd")
        ])

        self._sacc = self._sacc.with_columns([
            (pl.col("tStart") - start_time).alias("tStart"),
            (pl.col("tEnd") - start_time).alias("tEnd")
        ])

        if self._blink is not None:
            self._blink = self._blink.with_columns([
                (pl.col("tStart") - start_time).alias("tStart"),
                (pl.col("tEnd") - start_time).alias("tEnd")
            ])

        self.events_path = events_path
        self.detection_algorithm = events_path.name[:-7]


    def fixations(self):
        return self._fix

    @property
    def calib_index(self):
        return self._calib_index


    def saccades(self):
        return self._sacc


    def samples(self):
        return self._samples

    def __repr__(self):
        return f"Trial = '{self.trial_number}', " + self.session.__repr__()

    def plot_scanpath(self,screen_height,screen_width, **kwargs):
        vis = Visualization(self.events_path, self.detection_algorithm)
        (self.events_path / "plots").mkdir(parents=True, exist_ok=True)
        vis.scanpath(fixations=self._fix, saccades=self._sacc, samples=self._samples, screen_height=screen_height, screen_width=screen_width, 
                      folder_path=self.events_path / "plots", **kwargs)

    def plot_animation(self, screen_height, screen_width, video_path=None, background_image_path=None, **kwargs):
        """
        Create an animated visualization of eye-tracking data for this trial.

        When a video is provided, the animation syncs gaze samples with video frames.
        When no video is provided, gaze points are animated on a grey background or
        a provided background image, using the sample timestamps for timing.

        Parameters
        ----------
        screen_height, screen_width
            Stimulus resolution in pixels.
        video_path
            Path to a video file. If provided, gaze is overlaid on video frames.
        background_image_path
            Path to a background image. Only used when video_path is None.
            If both are None, a grey background is used.
        **kwargs
            Additional arguments passed to Visualization.plot_animation():
            - folder_path: Directory to save the animation
            - tmin, tmax: Time window in ms
            - seconds_to_show: Limit animation to first N seconds
            - scale_factor: Resolution scaling (default 0.5)
            - gaze_radius: Gaze point radius in pixels
            - gaze_color: RGB tuple for gaze color
            - fps: Animation frames per second
            - output_format: "html" (default), "mp4", "gif", or "matplotlib"
            - display: If True, return HTML for notebook display

        Returns
        -------
        IPython.display.HTML or None
            Returns HTML animation if display=True and output_format="html".
            For output_format="matplotlib", displays in a GUI window and returns None.
        """
        vis = Visualization(self.events_path, self.detection_algorithm)
        (self.events_path / "plots").mkdir(parents=True, exist_ok=True)

        return vis.plot_animation(
            samples=self._samples,
            screen_height=screen_height,
            screen_width=screen_width,
            video_path=video_path,
            background_image_path=background_image_path,
            **kwargs
        )

    def filter_fixations(self, min_fix_dur: int = 50):
        """
        1.  Delete fixations shorter than `min_fix_dur` (ms).
        2.  Merge the two saccades that flank each deleted fixation
            into one longer saccade, always staying inside a single
            (“phase”, “eye”) stream.

        Returns
        -------
        self   # so you can do:  trial.filter_fixations().is_trial_bad()
        """
        # ─────────────────────── 0 · split keep / drop ──────────────────────
        short_fix = self._fix.filter(pl.col("duration") < min_fix_dur)
        keep_fix  = self._fix.filter(pl.col("duration") >= min_fix_dur)

        if short_fix.is_empty():
            return                                # nothing to do

        # ─────────────────────── 1 · prepare saccades ───────────────────────
        sacc = (self._sacc       # add an integer key that survives every shuffle
                .with_row_count("idx")
                .sort(["phase", "eye", "tStart"]))

        prev_src = sacc.select(["idx", "phase", "eye",
                                pl.col("tEnd").alias("t")])
        next_src = sacc.select(["idx", "phase", "eye",
                                pl.col("tStart").alias("t")])

        # ─────────────────────── 2 · find neighbour IDs ─────────────────────
        short_fix = short_fix.rename({"tStart": "tStart_fix",
                                    "tEnd":   "tEnd_fix"})



        short_fix = short_fix.sort(["phase", "eye", "tStart_fix"])
        prev_src  = prev_src.sort(["phase", "eye", "t"])
        next_src  = next_src.sort(["phase", "eye", "t"])
        with warnings.catch_warnings():
            warnings.filterwarnings(
                "ignore",
                message="Sortedness of columns cannot be checked when 'by' groups provided",
                category=UserWarning,
            )

            short_fix = (
                short_fix
                .join_asof(
                    prev_src,
                    left_on="tStart_fix",
                    right_on="t",
                    by=["phase", "eye"],
                    strategy="backward"
                )
                .rename({"idx": "idx_prev"})
                .drop("t")
                .join_asof(
                    next_src,
                    left_on="tEnd_fix",
                    right_on="t",
                    by=["phase", "eye"],
                    strategy="forward"
                )
                .rename({"idx": "idx_next"})
                .drop("t")
            )

        # only keep rows where we found BOTH neighbours
        short_fix_pairs = short_fix.select(["idx_prev", "idx_next"]).drop_nulls()
        if short_fix_pairs.is_empty():
            # we could not build any (prev,next) pair → only delete fixations
            self._fix = keep_fix.sort(["phase", "tStart"])
            return self

        # ───────────────────── 3 · join the two saccades ────────────────────
        pair_df = (short_fix_pairs.unique()
                .join(sacc, left_on="idx_prev", right_on="idx", how="inner")
                .join(sacc, left_on="idx_next", right_on="idx", suffix="_nxt"))

        # keep **prev** row plus ONLY the four _nxt columns that we still need
        prev_cols = [c for c in pair_df.columns if not c.endswith("_nxt")]
        need_nxt  = ["tEnd_nxt", "xEnd_nxt", "yEnd_nxt", "vPeak_nxt"]
        merged = pair_df.select(prev_cols + need_nxt)

        # ───────── overwrite / derive fields that span both flanks ──────────
        merged = merged.with_columns([
            pl.col("tEnd_nxt").alias("tEnd"),
            (pl.col("tEnd_nxt") - pl.col("tStart")).alias("duration"),
            pl.col("xEnd_nxt").alias("xEnd"),
            pl.col("yEnd_nxt").alias("yEnd"),
            pl.max_horizontal("vPeak", "vPeak_nxt").alias("vPeak"),
            (
                (pl.col("xEnd_nxt") - pl.col("xStart"))**2
            + (pl.col("yEnd_nxt") - pl.col("yStart"))**2
            ).sqrt().alias("ampDeg"),
        ])

        # drop helper columns that end in _nxt (no longer needed)
        merged = merged.drop([c for c in merged.columns if c.endswith("_nxt")])

        # 4 · bring schema in line with original  --------------------------------
        base_cols = sacc.drop("idx").columns

        for col in base_cols:
            if col not in merged.columns:
                if f"{col}_nxt" in pair_df.columns:
                    merged = merged.with_columns(pl.col(f"{col}_nxt").alias(col))
                else:
                    merged = merged.with_columns(
                        pl.lit(None).cast(sacc[col].dtype).alias(col)
                    )

        # --- NEW: make sure every dtype matches the canonical sacc table ----
        for col in base_cols:
            if merged[col].dtype != sacc[col].dtype:
                merged = merged.with_columns(pl.col(col).cast(sacc[col].dtype))

        merged = merged.select(base_cols)

        # ───────────────────── 5 · build the final saccade table ────────────
        to_drop = pl.concat([short_fix_pairs["idx_prev"],
                            short_fix_pairs["idx_next"]]).unique()
        new_sacc = (sacc
                    .filter(~pl.col("idx").is_in(to_drop))
                    .drop("idx")          # helper column gone
                    .vstack(merged)       # add fused rows
                    .sort(["phase", "eye", "tStart"]))

        # ───────────────────── 6 · store back and return ────────────────────
        self._fix  = keep_fix.sort(["phase", "tStart"])
        self._sacc = new_sacc



    def collapse_fixations(self, threshold_px: float) -> None:
        """
        Collapse consecutive fixations that lie ≤ `threshold_px` apart
        *within each phase separately*.  Saccades whose whole time‑span
        falls between the first and last fixation of a pool are discarded.
        The saccade immediately before the pool has its (xEnd, yEnd)
        adjusted to the merged‑fixation centroid; the saccade immediately
        after the pool has its (xStart, yStart) adjusted likewise.

        After running:
            self._fix   → collapsed fixations
            self._sacc  → original saccades minus the discarded ones,
                        plus the updated coordinates for the two
                        bordering saccades.
        """

        # ────────────────── 0 · prepare helpers ──────────────────
        fix = self._fix.sort("tStart").with_row_count("fix_idx")
        sac = self._sacc.sort("tStart").with_row_count("sac_idx")

        new_fix_rows: list[dict] = []
        drop_sac_idx: set[int]   = set()
        mod_sac: dict[int, dict] = {}          # idx → partial‑row updates

        # ────────────────── 1 · loop over phases ─────────────────
        for phase_val in fix["phase"].unique():               # ① per phase
            # Loop over eyes if needed
            for eye in fix["eye"].unique():
                fix_p = fix.filter((pl.col("phase") == phase_val) & (pl.col("eye") == eye))
                sac_p = sac.filter((pl.col("phase") == phase_val) & (pl.col("eye") == eye))

                i, n_fix = 0, len(fix_p)
                while i < n_fix:

                    # ── grow one pool ───────────────────────────────
                    pool = [fix_p.row(i, named=True)]
                    j = i + 1
                    while j < n_fix:
                        dx = fix_p["xAvg"][j] - fix_p["xAvg"][j - 1]
                        dy = fix_p["yAvg"][j] - fix_p["yAvg"][j - 1]
                        if hypot(dx, dy) <= threshold_px:
                            pool.append(fix_p.row(j, named=True))
                            j += 1
                        else:
                            break

                    # ── pool of size 1: keep as‑is ──────────────────
                    if len(pool) == 1:
                        new_fix_rows.append(pool[0].copy())        # unchanged
                        i = j
                        continue

                    # ── merge the pool (>1 fix) ─────────────────────
                    first_fix, last_fix = pool[0], pool[-1]

                    merged_fix = first_fix.copy()
                    merged_fix.update({
                        "tEnd":     last_fix["tEnd"],
                        "duration": sum(f["duration"] for f in pool),
                        "xAvg":     np.mean([f["xAvg"] for f in pool]),
                        "yAvg":     np.mean([f["yAvg"] for f in pool]),
                        "pupilAvg": np.mean([f["pupilAvg"] for f in pool]),
                    })
                    new_fix_rows.append(merged_fix)

                    # ── identify & drop fully‑internal saccades ─────
                    inside = sac_p.filter(
                        (pl.col("tStart") >= first_fix["tEnd"]) &
                        (pl.col("tEnd")   <= last_fix["tStart"])
                    )
                    drop_sac_idx.update(inside["sac_idx"].to_list())

                    # ── adjust bordering saccades ───────────────────
                    merged_x = merged_fix["xAvg"]
                    merged_y = merged_fix["yAvg"]

                    # previous saccade (ends at first_fix.tStart)
                    prev_df = sac_p.filter(pl.col("tEnd") <= first_fix["tStart"]).tail(1)
                    if prev_df.height:
                        prev = prev_df.row(0, named=True)
                        idx  = prev["sac_idx"]
                        upd  = {
                            "xEnd": merged_x,
                            "yEnd": merged_y,
                            "dx":   merged_x - prev["xStart"],
                            "dy":   merged_y - prev["yStart"],
                        }
                        upd["amplitude"] = hypot(upd["dx"], upd["dy"])
                        mod_sac.setdefault(idx, {}).update(upd)

                    # next saccade (starts at last_fix.tEnd)
                    next_df = sac_p.filter(pl.col("tStart") >= last_fix["tEnd"]).head(1)
                    if next_df.height:
                        nxt = next_df.row(0, named=True)
                        idx = nxt["sac_idx"]
                        upd = {
                            "xStart": merged_x,
                            "yStart": merged_y,
                            "dx":     nxt["xEnd"] - merged_x,
                            "dy":     nxt["yEnd"] - merged_y,
                        }
                        upd["amplitude"] = hypot(upd["dx"], upd["dy"])
                        mod_sac.setdefault(idx, {}).update(upd)

                    i = j                                         # advance

        # ────────────────── 2 · rebuild tables ──────────────────
        # 2‑a  fixations
        new_fix = (
            pl.DataFrame(new_fix_rows,
                        schema=fix.drop("fix_idx").schema,
                        orient="row")
            .sort(["phase", "tStart"])
        )

        # 2‑b  saccades: drop + modify in one pass
        new_sac_rows = []
        for row in sac.iter_rows(named=True):
            idx = row["sac_idx"]
            if idx in drop_sac_idx:
                continue                                     # discard
            if idx in mod_sac:                               # apply edits
                row.update(mod_sac[idx])
                # re‑compute amplitude in case only dx/dy were provided
                if "amplitude" not in mod_sac[idx]:
                    row["amplitude"] = hypot(row["dx"], row["dy"])
            new_sac_rows.append({k: v for k, v in row.items() if k != "sac_idx"})

        new_sac = (
            pl.DataFrame(new_sac_rows,
                        schema=sac.drop("sac_idx").schema,
                        orient="row")
            .sort(["phase", "tStart"])
        )

        # ────────────────── 3 · store back ──────────────────────
        self._fix  = new_fix
        self._sacc = new_sac

    def save_rts(self):
        if hasattr(self, "_rts"):
            return

        # Filter out empty phase rows
        filtered = self._samples.filter(pl.col("phase") != "")

        # Calculate RT as the difference between last and first tSample per phase
        rts = (
            filtered
            .group_by("phase")
            .agg([
                (pl.col("tSample").max() - pl.col("tSample").min()).alias("rt")
            ])
            .with_columns([
                pl.lit(self.trial_number).alias("trial_number")
            ])
        )

        self._rts = rts


    def rts(self):
        if not hasattr(self, "_rts"):
            self.save_rts()
        return self._rts

    def is_trial_bad(self, phase, threshold=0.1):
        # Filter samples for the given phase
        samples = self._samples.filter(pl.col("phase") == phase)

        # Remove samples during blinks
        if self._blink is not None and self._blink.height > 0:
            for blink in self._blink.iter_rows(named=True):
                start, end = blink["tStart"], blink["tEnd"]
                samples = samples.filter(~((pl.col("tSample") > start) & (pl.col("tSample") < end)))

        total_samples = samples.height
        if total_samples == 0:
            return True  # If no samples remain, consider it bad

        # Count total NaNs across all columns
        nan_counts = samples.select([pl.col(c).is_null().sum().alias(c) for c in samples.columns])
        nan_total = sum(nan_counts.row(0))

        # Count "bad" values
        bad_values = samples.select(pl.col("bad").sum()).item()

        bad_and_nan_percentage = (nan_total + bad_values) / total_samples

        return bad_and_nan_percentage > threshold


    def is_trial_longer_than(self, seconds, phase):
        rt_row = self.rts().filter(pl.col("phase") == phase)
        if rt_row.is_empty():
            return False  # Or True if no data should be considered long
        return rt_row.select("rt").item() > seconds * 1000.0

    def compute_multimatch(self,other_trial: "Trial",screen_height,screen_width):
        trial_scanpath = self.search_fixations().to_pandas()
        trial_to_compare_scanpath = other_trial.search_fixations().to_pandas()
        # Turn trial scanpath into list of tuples
        trial_scanpath = [tuple(row) for row in trial_scanpath[["xAvg", "yAvg", "duration"]].values]
        trial_to_compare_scanpath = [tuple(row) for row in trial_to_compare_scanpath[["xAvg", "yAvg", "duration"]].values]

        # Convert the list of tuples into a numpy array with the format needed for the multimatch function
        trial_scanpath = np.array(trial_scanpath, dtype=[('start_x', '<f8'), ('start_y', '<f8'), ('duration', '<f8')])
        trial_to_compare_scanpath = np.array(trial_to_compare_scanpath, dtype=[('start_x', '<f8'), ('start_y', '<f8'), ('duration', '<f8')])

        return mm.docomparison(trial_scanpath, trial_to_compare_scanpath, (screen_width, screen_height))

collapse_fixations(threshold_px)

Collapse consecutive fixations that lie ≤ threshold_px apart within each phase separately. Saccades whose whole time‑span falls between the first and last fixation of a pool are discarded. The saccade immediately before the pool has its (xEnd, yEnd) adjusted to the merged‑fixation centroid; the saccade immediately after the pool has its (xStart, yStart) adjusted likewise.

After running: self._fix → collapsed fixations self._sacc → original saccades minus the discarded ones, plus the updated coordinates for the two bordering saccades.

Source code in pyxations/analysis/generic.py

def collapse_fixations(self, threshold_px: float) -> None:
    """
    Collapse consecutive fixations that lie ≤ `threshold_px` apart
    *within each phase separately*.  Saccades whose whole time‑span
    falls between the first and last fixation of a pool are discarded.
    The saccade immediately before the pool has its (xEnd, yEnd)
    adjusted to the merged‑fixation centroid; the saccade immediately
    after the pool has its (xStart, yStart) adjusted likewise.

    After running:
        self._fix   → collapsed fixations
        self._sacc  → original saccades minus the discarded ones,
                    plus the updated coordinates for the two
                    bordering saccades.
    """

    # ────────────────── 0 · prepare helpers ──────────────────
    fix = self._fix.sort("tStart").with_row_count("fix_idx")
    sac = self._sacc.sort("tStart").with_row_count("sac_idx")

    new_fix_rows: list[dict] = []
    drop_sac_idx: set[int]   = set()
    mod_sac: dict[int, dict] = {}          # idx → partial‑row updates

    # ────────────────── 1 · loop over phases ─────────────────
    for phase_val in fix["phase"].unique():               # ① per phase
        # Loop over eyes if needed
        for eye in fix["eye"].unique():
            fix_p = fix.filter((pl.col("phase") == phase_val) & (pl.col("eye") == eye))
            sac_p = sac.filter((pl.col("phase") == phase_val) & (pl.col("eye") == eye))

            i, n_fix = 0, len(fix_p)
            while i < n_fix:

                # ── grow one pool ───────────────────────────────
                pool = [fix_p.row(i, named=True)]
                j = i + 1
                while j < n_fix:
                    dx = fix_p["xAvg"][j] - fix_p["xAvg"][j - 1]
                    dy = fix_p["yAvg"][j] - fix_p["yAvg"][j - 1]
                    if hypot(dx, dy) <= threshold_px:
                        pool.append(fix_p.row(j, named=True))
                        j += 1
                    else:
                        break

                # ── pool of size 1: keep as‑is ──────────────────
                if len(pool) == 1:
                    new_fix_rows.append(pool[0].copy())        # unchanged
                    i = j
                    continue

                # ── merge the pool (>1 fix) ─────────────────────
                first_fix, last_fix = pool[0], pool[-1]

                merged_fix = first_fix.copy()
                merged_fix.update({
                    "tEnd":     last_fix["tEnd"],
                    "duration": sum(f["duration"] for f in pool),
                    "xAvg":     np.mean([f["xAvg"] for f in pool]),
                    "yAvg":     np.mean([f["yAvg"] for f in pool]),
                    "pupilAvg": np.mean([f["pupilAvg"] for f in pool]),
                })
                new_fix_rows.append(merged_fix)

                # ── identify & drop fully‑internal saccades ─────
                inside = sac_p.filter(
                    (pl.col("tStart") >= first_fix["tEnd"]) &
                    (pl.col("tEnd")   <= last_fix["tStart"])
                )
                drop_sac_idx.update(inside["sac_idx"].to_list())

                # ── adjust bordering saccades ───────────────────
                merged_x = merged_fix["xAvg"]
                merged_y = merged_fix["yAvg"]

                # previous saccade (ends at first_fix.tStart)
                prev_df = sac_p.filter(pl.col("tEnd") <= first_fix["tStart"]).tail(1)
                if prev_df.height:
                    prev = prev_df.row(0, named=True)
                    idx  = prev["sac_idx"]
                    upd  = {
                        "xEnd": merged_x,
                        "yEnd": merged_y,
                        "dx":   merged_x - prev["xStart"],
                        "dy":   merged_y - prev["yStart"],
                    }
                    upd["amplitude"] = hypot(upd["dx"], upd["dy"])
                    mod_sac.setdefault(idx, {}).update(upd)

                # next saccade (starts at last_fix.tEnd)
                next_df = sac_p.filter(pl.col("tStart") >= last_fix["tEnd"]).head(1)
                if next_df.height:
                    nxt = next_df.row(0, named=True)
                    idx = nxt["sac_idx"]
                    upd = {
                        "xStart": merged_x,
                        "yStart": merged_y,
                        "dx":     nxt["xEnd"] - merged_x,
                        "dy":     nxt["yEnd"] - merged_y,
                    }
                    upd["amplitude"] = hypot(upd["dx"], upd["dy"])
                    mod_sac.setdefault(idx, {}).update(upd)

                i = j                                         # advance

    # ────────────────── 2 · rebuild tables ──────────────────
    # 2‑a  fixations
    new_fix = (
        pl.DataFrame(new_fix_rows,
                    schema=fix.drop("fix_idx").schema,
                    orient="row")
        .sort(["phase", "tStart"])
    )

    # 2‑b  saccades: drop + modify in one pass
    new_sac_rows = []
    for row in sac.iter_rows(named=True):
        idx = row["sac_idx"]
        if idx in drop_sac_idx:
            continue                                     # discard
        if idx in mod_sac:                               # apply edits
            row.update(mod_sac[idx])
            # re‑compute amplitude in case only dx/dy were provided
            if "amplitude" not in mod_sac[idx]:
                row["amplitude"] = hypot(row["dx"], row["dy"])
        new_sac_rows.append({k: v for k, v in row.items() if k != "sac_idx"})

    new_sac = (
        pl.DataFrame(new_sac_rows,
                    schema=sac.drop("sac_idx").schema,
                    orient="row")
        .sort(["phase", "tStart"])
    )

    # ────────────────── 3 · store back ──────────────────────
    self._fix  = new_fix
    self._sacc = new_sac

filter_fixations(min_fix_dur=50)

1. Delete fixations shorter than min_fix_dur (ms).
2. Merge the two saccades that flank each deleted fixation into one longer saccade, always staying inside a single (“phase”, “eye”) stream.

Returns:

Type	Description
self

Source code in pyxations/analysis/generic.py

def filter_fixations(self, min_fix_dur: int = 50):
    """
    1.  Delete fixations shorter than `min_fix_dur` (ms).
    2.  Merge the two saccades that flank each deleted fixation
        into one longer saccade, always staying inside a single
        (“phase”, “eye”) stream.

    Returns
    -------
    self   # so you can do:  trial.filter_fixations().is_trial_bad()
    """
    # ─────────────────────── 0 · split keep / drop ──────────────────────
    short_fix = self._fix.filter(pl.col("duration") < min_fix_dur)
    keep_fix  = self._fix.filter(pl.col("duration") >= min_fix_dur)

    if short_fix.is_empty():
        return                                # nothing to do

    # ─────────────────────── 1 · prepare saccades ───────────────────────
    sacc = (self._sacc       # add an integer key that survives every shuffle
            .with_row_count("idx")
            .sort(["phase", "eye", "tStart"]))

    prev_src = sacc.select(["idx", "phase", "eye",
                            pl.col("tEnd").alias("t")])
    next_src = sacc.select(["idx", "phase", "eye",
                            pl.col("tStart").alias("t")])

    # ─────────────────────── 2 · find neighbour IDs ─────────────────────
    short_fix = short_fix.rename({"tStart": "tStart_fix",
                                "tEnd":   "tEnd_fix"})



    short_fix = short_fix.sort(["phase", "eye", "tStart_fix"])
    prev_src  = prev_src.sort(["phase", "eye", "t"])
    next_src  = next_src.sort(["phase", "eye", "t"])
    with warnings.catch_warnings():
        warnings.filterwarnings(
            "ignore",
            message="Sortedness of columns cannot be checked when 'by' groups provided",
            category=UserWarning,
        )

        short_fix = (
            short_fix
            .join_asof(
                prev_src,
                left_on="tStart_fix",
                right_on="t",
                by=["phase", "eye"],
                strategy="backward"
            )
            .rename({"idx": "idx_prev"})
            .drop("t")
            .join_asof(
                next_src,
                left_on="tEnd_fix",
                right_on="t",
                by=["phase", "eye"],
                strategy="forward"
            )
            .rename({"idx": "idx_next"})
            .drop("t")
        )

    # only keep rows where we found BOTH neighbours
    short_fix_pairs = short_fix.select(["idx_prev", "idx_next"]).drop_nulls()
    if short_fix_pairs.is_empty():
        # we could not build any (prev,next) pair → only delete fixations
        self._fix = keep_fix.sort(["phase", "tStart"])
        return self

    # ───────────────────── 3 · join the two saccades ────────────────────
    pair_df = (short_fix_pairs.unique()
            .join(sacc, left_on="idx_prev", right_on="idx", how="inner")
            .join(sacc, left_on="idx_next", right_on="idx", suffix="_nxt"))

    # keep **prev** row plus ONLY the four _nxt columns that we still need
    prev_cols = [c for c in pair_df.columns if not c.endswith("_nxt")]
    need_nxt  = ["tEnd_nxt", "xEnd_nxt", "yEnd_nxt", "vPeak_nxt"]
    merged = pair_df.select(prev_cols + need_nxt)

    # ───────── overwrite / derive fields that span both flanks ──────────
    merged = merged.with_columns([
        pl.col("tEnd_nxt").alias("tEnd"),
        (pl.col("tEnd_nxt") - pl.col("tStart")).alias("duration"),
        pl.col("xEnd_nxt").alias("xEnd"),
        pl.col("yEnd_nxt").alias("yEnd"),
        pl.max_horizontal("vPeak", "vPeak_nxt").alias("vPeak"),
        (
            (pl.col("xEnd_nxt") - pl.col("xStart"))**2
        + (pl.col("yEnd_nxt") - pl.col("yStart"))**2
        ).sqrt().alias("ampDeg"),
    ])

    # drop helper columns that end in _nxt (no longer needed)
    merged = merged.drop([c for c in merged.columns if c.endswith("_nxt")])

    # 4 · bring schema in line with original  --------------------------------
    base_cols = sacc.drop("idx").columns

    for col in base_cols:
        if col not in merged.columns:
            if f"{col}_nxt" in pair_df.columns:
                merged = merged.with_columns(pl.col(f"{col}_nxt").alias(col))
            else:
                merged = merged.with_columns(
                    pl.lit(None).cast(sacc[col].dtype).alias(col)
                )

    # --- NEW: make sure every dtype matches the canonical sacc table ----
    for col in base_cols:
        if merged[col].dtype != sacc[col].dtype:
            merged = merged.with_columns(pl.col(col).cast(sacc[col].dtype))

    merged = merged.select(base_cols)

    # ───────────────────── 5 · build the final saccade table ────────────
    to_drop = pl.concat([short_fix_pairs["idx_prev"],
                        short_fix_pairs["idx_next"]]).unique()
    new_sacc = (sacc
                .filter(~pl.col("idx").is_in(to_drop))
                .drop("idx")          # helper column gone
                .vstack(merged)       # add fused rows
                .sort(["phase", "eye", "tStart"]))

    # ───────────────────── 6 · store back and return ────────────────────
    self._fix  = keep_fix.sort(["phase", "tStart"])
    self._sacc = new_sacc

plot_animation(screen_height, screen_width, video_path=None, background_image_path=None, **kwargs)

Create an animated visualization of eye-tracking data for this trial.

When a video is provided, the animation syncs gaze samples with video frames. When no video is provided, gaze points are animated on a grey background or a provided background image, using the sample timestamps for timing.

Parameters:

Name	Type	Description	Default
screen_height		Stimulus resolution in pixels.	required
screen_width		Stimulus resolution in pixels.	required
video_path		Path to a video file. If provided, gaze is overlaid on video frames.	None
background_image_path		Path to a background image. Only used when video_path is None. If both are None, a grey background is used.	None
**kwargs		Additional arguments passed to Visualization.plot_animation(): - folder_path: Directory to save the animation - tmin, tmax: Time window in ms - seconds_to_show: Limit animation to first N seconds - scale_factor: Resolution scaling (default 0.5) - gaze_radius: Gaze point radius in pixels - gaze_color: RGB tuple for gaze color - fps: Animation frames per second - output_format: "html" (default), "mp4", "gif", or "matplotlib" - display: If True, return HTML for notebook display	{}

Returns:

Type	Description
HTML or None	Returns HTML animation if display=True and output_format="html". For output_format="matplotlib", displays in a GUI window and returns None.

Source code in pyxations/analysis/generic.py

def plot_animation(self, screen_height, screen_width, video_path=None, background_image_path=None, **kwargs):
    """
    Create an animated visualization of eye-tracking data for this trial.

    When a video is provided, the animation syncs gaze samples with video frames.
    When no video is provided, gaze points are animated on a grey background or
    a provided background image, using the sample timestamps for timing.

    Parameters
    ----------
    screen_height, screen_width
        Stimulus resolution in pixels.
    video_path
        Path to a video file. If provided, gaze is overlaid on video frames.
    background_image_path
        Path to a background image. Only used when video_path is None.
        If both are None, a grey background is used.
    **kwargs
        Additional arguments passed to Visualization.plot_animation():
        - folder_path: Directory to save the animation
        - tmin, tmax: Time window in ms
        - seconds_to_show: Limit animation to first N seconds
        - scale_factor: Resolution scaling (default 0.5)
        - gaze_radius: Gaze point radius in pixels
        - gaze_color: RGB tuple for gaze color
        - fps: Animation frames per second
        - output_format: "html" (default), "mp4", "gif", or "matplotlib"
        - display: If True, return HTML for notebook display

    Returns
    -------
    IPython.display.HTML or None
        Returns HTML animation if display=True and output_format="html".
        For output_format="matplotlib", displays in a GUI window and returns None.
    """
    vis = Visualization(self.events_path, self.detection_algorithm)
    (self.events_path / "plots").mkdir(parents=True, exist_ok=True)

    return vis.plot_animation(
        samples=self._samples,
        screen_height=screen_height,
        screen_width=screen_width,
        video_path=video_path,
        background_image_path=background_image_path,
        **kwargs
    )

visual_search

VisualSearchExperiment

Bases: Experiment

Source code in pyxations/analysis/visual_search.py

class VisualSearchExperiment(Experiment):
    def __init__(self, dataset_path: str,search_phase_name: str,memorization_phase_name: str, excluded_subjects: list = [], excluded_sessions: dict = {}, excluded_trials: dict = {}, export_format = FEATHER_EXPORT):
        self.dataset_path = Path(dataset_path)
        self.derivatives_path = self.dataset_path.with_name(self.dataset_path.name + "_derivatives")
        self.metadata = pl.read_csv(self.dataset_path / "participants.tsv", separator="\t", 
                                    dtypes={"subject_id": pl.Utf8, "old_subject_id": pl.Utf8})
        self.subjects = { subject_id:
            VisualSearchSubject(subject_id, old_subject_id, self, search_phase_name, memorization_phase_name,
                     excluded_sessions.get(subject_id, []), excluded_trials.get(subject_id, {}),export_format)
            for subject_id, old_subject_id in zip(self.metadata["subject_id"], self.metadata["old_subject_id"])
            if subject_id not in excluded_subjects and old_subject_id not in excluded_subjects
        }
        self.export_format = export_format
        self._search_phase_name = search_phase_name
        self._memorization_phase_name = memorization_phase_name

    def accuracy(self):
        accuracy = pl.concat([subject.accuracy() for subject in self.subjects.values()])

        return accuracy

    def plot_accuracy_by_subject(self):

        correct_responses = self.search_rts()
        # Sort by the sum of correct responses of each subject
        correct_responses_aux = (
            correct_responses
            .group_by(["subject_id", "memory_set_size", "target_present"])
            .agg(pl.col("correct_response").mean().alias("correct_response_mean"))
        ).select(["subject_id", "memory_set_size", "target_present", "correct_response_mean"])
        # Merge the correct_responses with the correct_responses_aux
        correct_responses = correct_responses.join(
            correct_responses_aux,
            on=["subject_id", "memory_set_size", "target_present"],
            how="left"
        ).sort(by=["memory_set_size", "target_present","correct_response_mean"])

        correct_responses = correct_responses.to_pandas()
        # target present to bool
        correct_responses["target_present"] = correct_responses["target_present"].astype(bool)
        # There should be an ax for each memory set size

        mem_set_sizes = correct_responses["memory_set_size"].unique()
        mem_set_sizes.sort()

        width_size = max(0.25 * len(correct_responses["subject_id"].unique()),10)

        n_rows = len(mem_set_sizes)
        fig, axs = plt.subplots(n_rows, 1, figsize=(width_size, 5 * n_rows),sharey=True)

        if n_rows == 1:
            axs = np.array([axs])

        for i, row in enumerate(mem_set_sizes):
            data = correct_responses[(correct_responses["memory_set_size"] == row)]
            sns.lineplot(x='subject_id',y='correct_response',data=data,hue='target_present',errorbar='se',ax=axs[i],estimator='mean')
            axs[i].set_title(f"Memory Set Size {row}")
            axs[i].tick_params(axis='x', rotation=90)
            axs[i].set_xlabel("Subject ID")
            axs[i].set_ylabel("Accuracy")

        plt.tight_layout()
        plt.show()
        plt.close()

    def plot_accuracy_by_stimulus(self):
        # Convert to pandas for Seaborn
        correct_responses = self.search_rts()
        correct_responses_aux = (
            correct_responses
            .group_by(["stimulus", "memory_set_size", "target_present"])
            .agg(pl.col("correct_response").mean().alias("correct_response_mean"))
        ).select(["stimulus", "memory_set_size", "target_present", "correct_response_mean"])

        # Merge the correct_responses with the correct_responses_aux
        correct_responses = correct_responses.join(
            correct_responses_aux,
            on=["stimulus", "memory_set_size", "target_present"],
            how="left"
        ).to_pandas().sort_values(by=["memory_set_size", "target_present", "correct_response_mean"])

        # Convert target_present to bool (in case it's int 0/1)
        correct_responses["target_present"] = correct_responses["target_present"].astype(bool)

        # One subplot per memory set size
        mem_set_sizes = sorted(correct_responses["memory_set_size"].unique())
        n_rows = len(mem_set_sizes)

        width_size = max(0.25 * len(correct_responses["stimulus"].unique()),10)

        fig, axs = plt.subplots(n_rows, 1, figsize=(width_size, 5 * n_rows), sharey=True)

        if n_rows == 1:
            axs = np.array([axs])

        for i, mem_size in enumerate(mem_set_sizes):
            data = correct_responses[correct_responses["memory_set_size"] == mem_size]
            sns.lineplot(x='stimulus',y='correct_response',data=data,hue='target_present',errorbar='se',ax=axs[i],estimator='mean')
            axs[i].set_title(f"Memory Set Size {mem_size}")
            axs[i].tick_params(axis='x', rotation=90)
            axs[i].set_xlabel("Stimulus")
            axs[i].set_ylabel("Accuracy")

        plt.tight_layout()
        plt.show()
        plt.close()

    def search_rts(self):
        rts = self.rts().filter(pl.col("phase") == self._search_phase_name)
        return rts

    def search_saccades(self):
        saccades = self.saccades().filter(pl.col("phase") == self._search_phase_name)
        return saccades

    def search_fixations(self):
        fixations = self.fixations().filter(pl.col("phase") == self._search_phase_name)
        return fixations

    def plot_speed_accuracy_tradeoff_by_subject(self):
        # 1) Aggregate the data
        speed_accuracy = (
            self.search_rts()
            .group_by(["target_present", "memory_set_size", "subject_id"])
            .agg([
                pl.col("rt").mean().alias("rt"),
                pl.col("correct_response").mean().alias("accuracy")
            ])
            .with_columns([
                pl.col("rt") / 1000,  # Convert to seconds
                pl.col("target_present").cast(pl.Boolean)
            ])
            .sort("memory_set_size")
        ).to_pandas()

        # 2) Unique memory set sizes
        mem_set_sizes = np.sort(speed_accuracy["memory_set_size"].unique())
        n_rows = len(mem_set_sizes)

        # 3) Prepare grid layout
        fig = plt.figure(figsize=(6, 1 + 6 * n_rows))
        gs = fig.add_gridspec(
            2 * n_rows, 2,
            width_ratios=(4, 1),
            height_ratios=[1, 4] * n_rows,
            left=0.1, right=0.9, bottom=0.07, top=0.85,
            wspace=0.05, hspace=0.05
        )

        # 4) Loop over memory set sizes
        for i, mem_size in enumerate(mem_set_sizes):
            data = speed_accuracy[speed_accuracy["memory_set_size"] == mem_size]

            row_top = 2 * i
            row_bottom = 2 * i + 1

            ax = fig.add_subplot(gs[row_bottom, 0])
            ax_histx = fig.add_subplot(gs[row_top, 0], sharex=ax)
            ax_histy = fig.add_subplot(gs[row_bottom, 1], sharey=ax)

            # (A) Scatter plot
            sns.scatterplot(
                x="accuracy",
                y="rt",
                data=data,
                hue="target_present",
                ax=ax,
                palette="deep"
            )

            # (B) Connection lines per subject
            for subj_id in data["subject_id"].unique():
                subj_data = data[data["subject_id"] == subj_id]
                if len(subj_data) != 2:
                    continue
                p0 = subj_data[subj_data["target_present"] == False]
                p1 = subj_data[subj_data["target_present"] == True]
                if not p0.empty and not p1.empty:
                    ax.plot(
                        [p0["accuracy"].values[0], p1["accuracy"].values[0]],
                        [p0["rt"].values[0], p1["rt"].values[0]],
                        color="black", alpha=0.3, linewidth=0.5, zorder=0
                    )

            # (C) Marginal histograms
            ax_histx.hist(data["accuracy"], bins=np.linspace(0, 1, 21), color="gray")
            ax_histy.hist(data["rt"], bins=20, orientation='horizontal', color="gray")

            ax_histx.tick_params(axis="x", labelbottom=False)
            ax_histy.tick_params(axis="y", labelleft=False)

            ax_histx.set_title(f"Memory Set Size {mem_size}")
            ax.set_xlim(0, 1)
            ax.set_ylim(0, speed_accuracy["rt"].max() * 1.1)
            ax.set_xlabel("Accuracy")
            ax.set_ylabel("Mean RT (s)")

        plt.suptitle("Speed-Accuracy Tradeoff by Subject", fontsize=14)
        plt.show()
        plt.close()

    def plot_speed_accuracy_tradeoff_by_stimulus(self):
        # 1) Aggregate the data
        speed_accuracy = (
            self.search_rts()
            .group_by(["target_present", "memory_set_size", "stimulus"])
            .agg([
                pl.col("rt").mean().alias("rt"),
                pl.col("correct_response").mean().alias("accuracy")
            ])
            .with_columns([
                (pl.col("rt") / 1000).alias("rt"),  # convert ms → s
                pl.col("target_present").cast(pl.Boolean)
            ])
            .sort("memory_set_size")
            .to_pandas()
        )

        # 2) Unique memory set sizes
        mem_set_sizes = np.sort(speed_accuracy["memory_set_size"].unique())
        n_rows = len(mem_set_sizes)

        # 3) Prepare grid layout
        fig = plt.figure(figsize=(6, 1 + 6 * n_rows))
        gs = fig.add_gridspec(
            2 * n_rows, 2,
            width_ratios=(4, 1),
            height_ratios=[1, 4] * n_rows,
            left=0.1, right=0.9, bottom=0.07, top=0.85,
            wspace=0.05, hspace=0.05
        )

        # 4) Loop over memory set sizes
        for i, mem_size in enumerate(mem_set_sizes):
            data = speed_accuracy[speed_accuracy["memory_set_size"] == mem_size]

            row_top = 2 * i
            row_bottom = 2 * i + 1

            ax = fig.add_subplot(gs[row_bottom, 0])
            ax_histx = fig.add_subplot(gs[row_top, 0], sharex=ax)
            ax_histy = fig.add_subplot(gs[row_bottom, 1], sharey=ax)

            # (A) Main scatter plot
            sns.scatterplot(
                x="accuracy",
                y="rt",
                data=data,
                hue="target_present",
                ax=ax,
                palette="deep"
            )

            # (B) Connect stimulus points (False → True)
            for stim in data["stimulus"].unique():
                stim_data = data[data["stimulus"] == stim]
                if len(stim_data) != 2:
                    continue
                p0 = stim_data[stim_data["target_present"] == False]
                p1 = stim_data[stim_data["target_present"] == True]
                if not p0.empty and not p1.empty:
                    ax.plot(
                        [p0["accuracy"].values[0], p1["accuracy"].values[0]],
                        [p0["rt"].values[0], p1["rt"].values[0]],
                        color="black", alpha=0.3, linewidth=0.5, zorder=0
                    )

            # (C) Marginal histograms
            ax_histx.hist(data["accuracy"], bins=np.linspace(0, 1, 21), color="gray")
            ax_histy.hist(data["rt"], bins=20, orientation='horizontal', color="gray")

            ax_histx.tick_params(axis="x", labelbottom=False)
            ax_histy.tick_params(axis="y", labelleft=False)

            # (D) Titles, limits, labels
            ax_histx.set_title(f"Memory Set Size {mem_size}")
            ax.set_xlim(0, 1)
            ax.set_ylim(0, speed_accuracy["rt"].max() * 1.1)
            ax.set_xlabel("Accuracy")
            ax.set_ylabel("Mean RT (s)")

        # 5) Final touch
        plt.suptitle("Speed-Accuracy Tradeoff by Stimulus", fontsize=14)
        plt.show()
        plt.close()

    def remove_non_answered_trials(self, print_flag=True):
        amount_trials_before_removal = self.search_rts().shape[0]
        for subject in list(self.subjects.values()):
            subject.remove_non_answered_trials(False)

        if print_flag:
            print(f"Removed {amount_trials_before_removal - self.search_rts().shape[0]} non answered trials")

    def remove_poor_accuracy_sessions(self, threshold=0.5, print_flag=True):
        amount_sessions_total = sum([len(subject.sessions) for subject in self.subjects.values()])
        for subject in list(self.subjects.keys()):
            self.subjects[subject].remove_poor_accuracy_sessions(threshold,False)

        if print_flag:
            print(f"Removed {amount_sessions_total - sum([len(subject.sessions) for subject in self.subjects.values()])} sessions with poor accuracy")                



    def scanpaths_by_stimuli(self):
        return pl.concat([subject.scanpaths_by_stimuli() for subject in self.subjects.values()])



    def find_fixation_cutoff(self, percentile=1.0):
        # 1. Gather fixation counts
        fix_counts = [
            {
                "fix_count": trial.search_fixations().height,
                "target_present": trial.target_present,
                "memory_set_size": trial.memory_set_size
            }
            for subject in self.subjects.values()
            for session in subject.sessions.values()
            for trial in session.trials.values()
        ]
        fix_counts = pl.DataFrame(fix_counts)

        # 2. Get all unique group keys
        group_keys = fix_counts.select(["target_present", "memory_set_size"]).unique().to_dicts()

        # 3. Compute cutoff per group
        rows = []
        for group in group_keys:
            tp = group["target_present"]
            mem_size = group["memory_set_size"]

            group_df = fix_counts.filter(
                (pl.col("target_present") == tp) &
                (pl.col("memory_set_size") == mem_size)
            )

            fix_counts_list = group_df["fix_count"].to_list()
            total_fixations = sum(fix_counts_list)
            threshold = total_fixations * percentile
            max_possible = max(fix_counts_list)

            fix_cutoff = _find_fixation_cutoff(
                fix_count_list=fix_counts_list,
                threshold=threshold,
                max_possible=max_possible
            )

            rows.append({
                "target_present": tp,
                "memory_set_size": mem_size,
                "fix_cutoff": fix_cutoff
            })

        return pl.DataFrame(rows)


    def remove_trials_for_stimuli(self,stimuli,print_flag=True):
        '''
        Remove trials for stimuli that are in the list of stimuli.
        Parameters:
            - stimuli: list of stimuli to remove
            - print_flag: if True, print the number of trials removed
        '''
        # Get the trials for the stimuli to remove
        amount_trials_removed = 0
        subj_keys = list(self.subjects.keys())
        for subject_key in subj_keys:
            subject = self.subjects[subject_key]
            session_keys = list(subject.sessions.keys())
            for session_key in session_keys:
                session = subject.sessions[session_key]
                trial_keys = list(session.trials.keys())
                for trial_key in trial_keys:
                    trial = session.trials[trial_key]
                    if trial.stimulus in stimuli:
                        session.remove_trial(trial_key)
                        amount_trials_removed += 1
        if print_flag:
            print(f"Removed {amount_trials_removed} trials for stimuli {stimuli}")



    def remove_trials_for_stimuli_with_poor_accuracy(self, threshold=0.5, print_flag=True):
        '''For now this will be done without grouping by target_present'''
        scanpaths_by_stimuli = self.scanpaths_by_stimuli()
        grouped = scanpaths_by_stimuli.group_by(["stimulus", "memory_set_size"])
        poor_accuracy_stimuli = (
                                grouped.agg(pl.col("correct_response").mean().alias("accuracy"))
                                .filter(pl.col("accuracy") < threshold)
                            )
        # Get the stimulus and memory set size of poor_accuracy_stimuli into a list of tuples
        poor_accuracy_stimuli = poor_accuracy_stimuli.select(pl.col("stimulus"), pl.col("memory_set_size")).to_dicts()
        poor_accuracy_stimuli = [(stimulus["stimulus"], stimulus["memory_set_size"]) for stimulus in poor_accuracy_stimuli]
        amount_trials_removed = 0
        subj_keys = list(self.subjects.keys())
        for subject_key in subj_keys:
            subject = self.subjects[subject_key]
            session_keys = list(subject.sessions.keys())
            for session_key in session_keys:
                session = subject.sessions[session_key]
                trial_keys = list(session.trials.keys())
                for trial_key in trial_keys:
                    trial = session.trials[trial_key]
                    if (trial.stimulus, trial.memory_set_size) in poor_accuracy_stimuli:
                        session.remove_trial(trial_key)
                        amount_trials_removed += 1
        if print_flag:
            print(f"Removed {amount_trials_removed} trials from stimuli with less than {threshold} accuracy.")

    def cumulative_correct_trials_by_fixation(self, group_cutoffs=None):
        if group_cutoffs is None:
            group_cutoffs = self.find_fixation_cutoff()
        cumulative_correct = pl.concat([subject.cumulative_correct_trials_by_fixation(group_cutoffs) for subject in self.subjects.values()])

        return cumulative_correct



    def plot_cumulative_performance(self, group_cutoffs=None):
        if group_cutoffs is None:
            group_cutoffs = self.find_fixation_cutoff()
        cumulative_performance = self.cumulative_correct_trials_by_fixation(group_cutoffs).join(
            group_cutoffs,
            on=["target_present", "memory_set_size"],
            how="left"
        )

        tp_ta = cumulative_performance.select(pl.col("target_present")).unique().to_series()
        tp_ta.sort()
        mem_set_sizes = cumulative_performance.select(pl.col("memory_set_size")).unique().to_series()
        mem_set_sizes.sort()

        # Convert to pandas for Seaborn
        cumulative_performance = cumulative_performance.to_pandas()

        n_cols = len(tp_ta)
        n_rows = len(mem_set_sizes)
        fig, axs = plt.subplots(n_rows, n_cols, figsize=(6 * n_cols, 5 * n_rows),sharey=True)
        fig.suptitle("Cumulative Performance")
        if n_cols == 1:
            axs = np.array([axs])

        if n_rows == 1:
            axs = np.array([axs])

        # For each fixation number (i.e. first "max_fixations" columns), we need the mean and the standard error
        # The X axis will be the fixation number, the Y axis will be the accuracy
        # The area around the mean will be the standard error

        for i, row in enumerate(mem_set_sizes):
            for j, col in enumerate(tp_ta):
                # Get the max fix for the current group, groups_cutoff is in polars

                data = cumulative_performance[(cumulative_performance["memory_set_size"] == row) & (cumulative_performance["target_present"] == col)]
                max_fix = int(data["fix_cutoff"].iloc[0])

                # 1. Trim every array to the same length (optional but handy)
                trimmed = data["cumulative_correct"].apply(lambda arr: arr[:max_fix])

                # 2. Turn the Series-of-lists into long form
                exploded = trimmed.explode().reset_index(drop=True).to_frame("cumulative_correct")

                # 3. Add a 1-based fixation index
                exploded["fixation_number"] = (np.tile(np.arange(1, max_fix + 1), len(data))  # repeat 1..max_fix for every original row
    )
                sns.lineplot(
                    x="fixation_number",
                    y="cumulative_correct",
                    data=exploded,
                    ax=axs[i, j],
                    errorbar='se',
                    estimator='mean',
                    color="black"
                )
                axs[i, j].set_title(f"Memory Set Size {int(row)}, Target Present {bool(col)}")
                # Ticks every 5 fixations
                axs[i, j].set_xticks(range(0, max_fix, 5))
                axs[i, j].set_xticklabels(range(1, max_fix+1, 5))
                axs[i, j].set_xlabel("Fixation Number")
                axs[i, j].set_ylabel("Accuracy")

        plt.ylim(0, 1)
        plt.tight_layout()
        plt.show()
        plt.close()

    def trials_by_rt_bins(self, bin_end, bin_step):
        # 1. Get and filter RTs
        rts = self.rts().filter(pl.col("phase") == self._search_phase_name)
        rts = rts.with_columns([
            (pl.col("rt") / 1000).alias("rt")
        ])

        # 2. Compute bin edges
        bin_edges = np.arange(0, bin_end + bin_step, bin_step)

        # 3. Bin RTs using numpy (returns indices)
        bin_indices = np.digitize(rts["rt"].to_numpy(), bin_edges, right=False)

        # 4. Convert to left edge values
        rt_bin_labels = [bin_edges[i - 1] if i > 0 and i < len(bin_edges) else None for i in bin_indices]

        # 5. Assign back to the DataFrame
        rts = rts.with_columns([
            pl.Series("rt_bin", rt_bin_labels)
        ])

        return rts


    def plot_correct_trials_by_rt_bins(self, bin_end, bin_step):
        # Get relevant trial info with binned RTs
        correct_trials_per_bin = (
            self.trials_by_rt_bins(bin_end, bin_step)
            .select(["rt_bin", "target_present", "memory_set_size", "correct_response"])
            .group_by(["rt_bin", "target_present", "memory_set_size"])
            .agg(pl.col("correct_response").sum().alias("correct_response"))
            .sort(["memory_set_size", "target_present", "rt_bin"])
        ).to_pandas()

        # Ensure sorted and unique values
        tp_ta = sorted(correct_trials_per_bin["target_present"].unique())
        mem_set_sizes = sorted(correct_trials_per_bin["memory_set_size"].unique())

        n_cols = len(tp_ta)
        n_rows = len(mem_set_sizes)

        fig, axs = plt.subplots(n_rows, n_cols, figsize=(6 * n_cols, 5 * n_rows), sharey=True, sharex=True)
        fig.suptitle("Correct Trials by RT Bins")

        # Normalize axis shape for consistent indexing
        if n_cols == 1:
            axs = np.expand_dims(axs, axis=1)
        if n_rows == 1:
            axs = np.expand_dims(axs, axis=0)

        for i, mem_size in enumerate(mem_set_sizes):
            for j, tp in enumerate(tp_ta):
                data = correct_trials_per_bin[
                    (correct_trials_per_bin["memory_set_size"] == mem_size) &
                    (correct_trials_per_bin["target_present"] == tp)
                ]
                sns.barplot(x="rt_bin", y="correct_response", data=data, ax=axs[i, j])
                axs[i, j].set_title(f"Memory Set Size {mem_size}, Target Present {bool(tp)}")
                axs[i, j].set_xlabel("RT Bins (s)")
                axs[i, j].set_ylabel("Correct Trials")
                axs[i, j].set_xticks(range(0, int(bin_end/bin_step)+3, 3))

        plt.tight_layout()
        plt.show()
        plt.close()

    def plot_incorrect_trials_by_rt_bins(self, bin_end, bin_step):
        # Get RT binned trial info
        incorrect_trials_per_bin = (
            self.trials_by_rt_bins(bin_end, bin_step)
            .select(["rt_bin", "target_present", "memory_set_size", "correct_response"])
            .with_columns([
                (1 - pl.col("correct_response")).alias("incorrect_response")
            ])
            .group_by(["rt_bin", "target_present", "memory_set_size"])
            .agg(pl.col("incorrect_response").sum().alias("incorrect_response"))
            .sort(["memory_set_size", "target_present", "rt_bin"])
            .to_pandas()
        )

        # Setup for plotting
        tp_ta = sorted(incorrect_trials_per_bin["target_present"].unique())
        mem_set_sizes = sorted(incorrect_trials_per_bin["memory_set_size"].unique())

        n_cols = len(tp_ta)
        n_rows = len(mem_set_sizes)

        fig, axs = plt.subplots(n_rows, n_cols, figsize=(6 * n_cols, 5 * n_rows), sharey=True, sharex=True)
        fig.suptitle("Incorrect Trials by RT Bins")

        # Normalize shape for subplots
        if n_cols == 1:
            axs = np.expand_dims(axs, axis=1)
        if n_rows == 1:
            axs = np.expand_dims(axs, axis=0)

        for i, mem_size in enumerate(mem_set_sizes):
            for j, tp in enumerate(tp_ta):
                data = incorrect_trials_per_bin[
                    (incorrect_trials_per_bin["memory_set_size"] == mem_size) &
                    (incorrect_trials_per_bin["target_present"] == tp)
                ]
                sns.barplot(x="rt_bin", y="incorrect_response", data=data, ax=axs[i, j])
                axs[i, j].set_title(f"Memory Set Size {mem_size}, Target Present {bool(tp)}")
                axs[i, j].set_xlabel("RT Bins (s)")
                axs[i, j].set_ylabel("Incorrect Trials")
                axs[i, j].set_xticks(range(0, int(bin_end/bin_step)+3, 3))

        plt.tight_layout()
        plt.show()
        plt.close()

    def plot_probability_of_deciding_by_rt_bin(self, bin_end, bin_step):
        # Get RT-binned trials
        trials = (
            self.trials_by_rt_bins(bin_end, bin_step)
            .select(["rt_bin", "target_present", "memory_set_size", "correct_response"])
        )

        # Unique labels
        tp_ta = sorted(trials["target_present"].unique().to_list())
        mem_set_sizes = sorted(trials["memory_set_size"].unique().to_list())

        n_cols = len(tp_ta)
        n_rows = len(mem_set_sizes)

        # Count occurrences per bin
        grouped = (
            trials
            .group_by(["rt_bin", "target_present", "correct_response", "memory_set_size"])
            .agg(pl.count().alias("count"))
            .sort(["correct_response", "target_present", "memory_set_size", "rt_bin"])
        )

        # Compute totals per (correctness, target, memory)
        totals = (
            grouped
            .group_by(["correct_response", "target_present", "memory_set_size"])
            .agg(pl.col("count").sum().alias("total_per_group"))
        )

        # Merge total counts back
        grouped = grouped.join(
            totals,
            on=["correct_response", "target_present", "memory_set_size"],
            how="left"
        )

        # Compute cumulative sums within groups
        grouped = (
            grouped
            .with_columns([
                pl.col("count").cum_sum().over(["correct_response", "target_present", "memory_set_size"]).alias("cumsum"),
            ])
            .with_columns([
                (pl.col("total_per_group") - pl.col("cumsum") + pl.col("count")).alias("total_per_bin"),
                (pl.col("count") / (pl.col("total_per_group") - pl.col("cumsum") + pl.col("count"))).alias("count_normalized"),
                pl.col("correct_response").cast(pl.Boolean)
            ])
        )

        # Convert to pandas for seaborn
        grouped_pd = grouped.to_pandas()

        # Plot setup
        fig, axs = plt.subplots(n_rows, n_cols, figsize=(6 * n_cols, 5 * n_rows), sharey=True, sharex=True)
        fig.suptitle("Probability of Deciding by RT Bins")

        if n_cols == 1:
            axs = np.expand_dims(axs, axis=1)
        if n_rows == 1:
            axs = np.expand_dims(axs, axis=0)

        for i, mem_size in enumerate(mem_set_sizes):
            for j, tp in enumerate(tp_ta):
                data = grouped_pd[
                    (grouped_pd["memory_set_size"] == mem_size) &
                    (grouped_pd["target_present"] == tp)
                ]
                sns.barplot(x="rt_bin", y="count_normalized", hue="correct_response", data=data, ax=axs[i, j])
                axs[i, j].set_title(f"Memory Set Size {mem_size}, Target Present {bool(tp)}")
                axs[i, j].set_xlabel("RT Bins (s)")
                axs[i, j].set_ylabel("Probability of Deciding")
                axs[i, j].set_xticks(range(0, int(bin_end/bin_step)+3, 3))

        plt.tight_layout()
        plt.show()
        plt.close()

remove_trials_for_stimuli(stimuli, print_flag=True)

Remove trials for stimuli that are in the list of stimuli. Parameters: - stimuli: list of stimuli to remove - print_flag: if True, print the number of trials removed

Source code in pyxations/analysis/visual_search.py

def remove_trials_for_stimuli(self,stimuli,print_flag=True):
    '''
    Remove trials for stimuli that are in the list of stimuli.
    Parameters:
        - stimuli: list of stimuli to remove
        - print_flag: if True, print the number of trials removed
    '''
    # Get the trials for the stimuli to remove
    amount_trials_removed = 0
    subj_keys = list(self.subjects.keys())
    for subject_key in subj_keys:
        subject = self.subjects[subject_key]
        session_keys = list(subject.sessions.keys())
        for session_key in session_keys:
            session = subject.sessions[session_key]
            trial_keys = list(session.trials.keys())
            for trial_key in trial_keys:
                trial = session.trials[trial_key]
                if trial.stimulus in stimuli:
                    session.remove_trial(trial_key)
                    amount_trials_removed += 1
    if print_flag:
        print(f"Removed {amount_trials_removed} trials for stimuli {stimuli}")

remove_trials_for_stimuli_with_poor_accuracy(threshold=0.5, print_flag=True)

For now this will be done without grouping by target_present

Source code in pyxations/analysis/visual_search.py

def remove_trials_for_stimuli_with_poor_accuracy(self, threshold=0.5, print_flag=True):
    '''For now this will be done without grouping by target_present'''
    scanpaths_by_stimuli = self.scanpaths_by_stimuli()
    grouped = scanpaths_by_stimuli.group_by(["stimulus", "memory_set_size"])
    poor_accuracy_stimuli = (
                            grouped.agg(pl.col("correct_response").mean().alias("accuracy"))
                            .filter(pl.col("accuracy") < threshold)
                        )
    # Get the stimulus and memory set size of poor_accuracy_stimuli into a list of tuples
    poor_accuracy_stimuli = poor_accuracy_stimuli.select(pl.col("stimulus"), pl.col("memory_set_size")).to_dicts()
    poor_accuracy_stimuli = [(stimulus["stimulus"], stimulus["memory_set_size"]) for stimulus in poor_accuracy_stimuli]
    amount_trials_removed = 0
    subj_keys = list(self.subjects.keys())
    for subject_key in subj_keys:
        subject = self.subjects[subject_key]
        session_keys = list(subject.sessions.keys())
        for session_key in session_keys:
            session = subject.sessions[session_key]
            trial_keys = list(session.trials.keys())
            for trial_key in trial_keys:
                trial = session.trials[trial_key]
                if (trial.stimulus, trial.memory_set_size) in poor_accuracy_stimuli:
                    session.remove_trial(trial_key)
                    amount_trials_removed += 1
    if print_flag:
        print(f"Removed {amount_trials_removed} trials from stimuli with less than {threshold} accuracy.")

VisualSearchSession

Bases: Session

Source code in pyxations/analysis/visual_search.py

class VisualSearchSession(Session):
    BEH_COLUMNS: list[str] = [
        "trial_number", "stimulus", "stimulus_coords", "memory_set", "memory_set_locations",
        "target_present", "target", "target_location", "correct_response", "was_answered"
    ]
    """
    Columns explanation:
    - trial_number: The number of the trial, in the order they were presented. They start from 0.
    - stimulus: The filename of the stimulus presented.
    - stimulus_coords: The coordinates of the stimulus presented. It should be a tuple containing the x, y of the top-left corner of the stimulus and the x, y of the bottom-right corner.
    - memory_set: The set of items memorized by the participant. It should be a list of strings. Each string should be the filename of the stimulus.
    - memory_set_locations: The locations of the items memorized by the participant. It should be a list of tuples. Each tuple should contain bounding
      boxes of the items memorized by the participant. The bounding boxes should be in the format (x1, y1, x2, y2), where (x1, y1) is the top-left corner and
      (x2, y2) is the bottom-right corner.
    - target_present: Whether one of the items is present in the stimulus. It should be a boolean.
    - target: The filename of the target item. It should be a string. If target_present is False, the value for this column will
      not be taken into account.
    - target_location: The location of the target item. It should be a tuple containing the bounding box of the target item. The bounding box should be in
      the format (x1, y1, x2, y2), where (x1, y1) is the top-left corner and (x2, y2) is the bottom-right corner. If target_present is False, the value for this column will
      not be taken into account.
    - correct_response: The correct response for the trial. It should be a boolean.
    - was_answered: Whether the trial was answered by the participant. It should be a boolean.

    Notice that you can get the actual response of the user by using the "correct_response" and "target_present" columns.
    For all of the heights, widths and locations of the items, the values should be in pixels and according to the screen itself.
    """

    COLLECTION_COLUMNS: dict = {
        "stimulus_coords": tuple,           # Parse as a tuple
        "memory_set": list,                 # Parse as a list
        "memory_set_locations": list,       # Parse as a list of tuples
        "target_location": tuple          # Parse as a tuple
    }

    def __init__(
        self, 
        session_id: str, 
        subject: VisualSearchSubject,  
        search_phase_name: str,
        memorization_phase_name: str,
        excluded_trials: list = None,
        export_format = FEATHER_EXPORT
    ):
        excluded_trials = [] if excluded_trials is None else excluded_trials
        super().__init__(session_id, subject, excluded_trials, export_format)
        self._search_phase_name = search_phase_name
        self._memorization_phase_name = memorization_phase_name
        self.behavior_data = None



    def load_behavior_data(self):
        # Get the name of the only csv file in the behavior path
        behavior_path = self.session_dataset_path / "behavioral"

        behavior_files = list(behavior_path.glob("*.csv"))

        if len(behavior_files) != 1:
            raise ValueError(
                f"There should only be one CSV file in the behavior path for session {self.session_id} "
                f"of subject {self.subject.subject_id}. Found files: {[file.name for file in behavior_files]}"
            )

        # Load the CSV file
        name = behavior_files[0].name
        self.behavior_data = pl.read_csv(
            behavior_path / name,
            dtypes={
                "trial_number": pl.Int32,
                "stimulus": pl.Utf8,
                "target_present": pl.Int32,
                "target": pl.Utf8,
                "correct_response": pl.Int32,
                "was_answered": pl.Int32
            }
        )

        # Validate that all required columns are present
        missing_columns = set(self.BEH_COLUMNS) - set(self.behavior_data.columns)
        if missing_columns:
            raise ValueError(f"Missing columns in behavior data: {missing_columns} for session {self.session_id} of subject {self.subject.subject_id}")

    def _init_trials(self,samples,fix,sacc,blink,events_path):
        self._trials = {trial:
            VisualSearchTrial(trial, self, samples, fix, sacc, blink, events_path, self.behavior_data,self._search_phase_name,self._memorization_phase_name)
            for trial in samples["trial_number"].unique() 
            if trial != -1 and trial not in self.excluded_trials and trial in self.behavior_data.select(pl.col("trial_number")).to_series().to_list()
        }

    def load_data(self, detection_algorithm: str):
        self.load_behavior_data()
        super().load_data(detection_algorithm)


    def search_rts(self):
        rts = self.rts().filter(pl.col("phase") == self._search_phase_name)
        return rts

    def search_saccades(self):
        saccades = self.saccades().filter(pl.col("phase") == self._search_phase_name)
        return saccades

    def search_fixations(self):
        fixations = self.fixations().filter(pl.col("phase") == self._search_phase_name)
        return fixations

    def accuracy(self):
        # Accuracy should be grouped by target present and memory set size
        correct_trials = self.search_rts()[["target_present", "correct_response", "memory_set_size"]]
        accuracy = correct_trials.groupby(["target_present", "memory_set_size"]).mean().reset_index()
        # Change the column name to accuracy
        accuracy.rename(columns={"correct_response": "accuracy"}, inplace=True)
        accuracy["session_id"] = self.session_id

        return accuracy

    def remove_non_answered_trials(self,print_flag=True):
        # Remove trials that were not answered
        non_answered_trials = [trial for trial in self.trials if not self.trials[trial].was_answered]
        for trial in non_answered_trials:
            self.remove_trial(trial)
        if print_flag:
            print(f"Removed {len(non_answered_trials)} non answered trials from session {self.session_id}")

    def has_poor_accuracy(self, threshold=0.5):
        correct_trials = self.search_rts()[["target_present", "correct_response", "memory_set_size"]]
        accuracy = correct_trials["correct_response"].sum() / correct_trials["correct_response"].count()
        return accuracy < threshold

    def find_fixation_cutoff(self, percentile=1.0):
        # 1. Gather fixation counts
        fix_counts = [
            {
                "fix_count": trial.search_fixations().height,
                "target_present": trial.target_present,
                "memory_set_size": trial.memory_set_size
            }

            for trial in self.trials.values()
        ]
        fix_counts = pl.DataFrame(fix_counts)

        # 2. Get all unique group keys
        group_keys = fix_counts.select(["target_present", "memory_set_size"]).unique().to_dicts()

        # 3. Compute cutoff per group
        rows = []
        for group in group_keys:
            tp = group["target_present"]
            mem_size = group["memory_set_size"]

            group_df = fix_counts.filter(
                (pl.col("target_present") == tp) &
                (pl.col("memory_set_size") == mem_size)
            )

            fix_counts_list = group_df["fix_count"].to_list()
            total_fixations = sum(fix_counts_list)
            threshold = total_fixations * percentile
            max_possible = max(fix_counts_list)

            fix_cutoff = _find_fixation_cutoff(
                fix_count_list=fix_counts_list,
                threshold=threshold,
                max_possible=max_possible
            )

            rows.append({
                "target_present": tp,
                "memory_set_size": mem_size,
                "fix_cutoff": fix_cutoff
            })

        return pl.DataFrame(rows)


    def cumulative_correct_trials_by_fixation(self, group_cutoffs=None):
        if group_cutoffs is None:
            group_cutoffs = self.find_fixation_cutoff()  # this should return a pl.DataFrame

        records = []

        for trial in self.trials.values():
            scanpath_length = len(trial.search_fixations())

            # ✅ Filter the appropriate fix_cutoff value
            fix_cutoff = (
                group_cutoffs
                .filter(
                    (pl.col("memory_set_size") == trial.memory_set_size) &
                    (pl.col("target_present") == trial.target_present)
                )
                .select("fix_cutoff")
                .item()
            )

            cumulative_correct = np.zeros(fix_cutoff)

            if trial.correct_response and scanpath_length - 1 <= fix_cutoff:
                cumulative_correct[scanpath_length - 1:] = 1

            records.append({
                "cumulative_correct": cumulative_correct,
                "target_present": trial.target_present,
                "memory_set_size": trial.memory_set_size,
            })

        df = pl.DataFrame(records)
        return df

    def scanpaths_by_stimuli(self):
        return pl.DataFrame([trial.scanpath_by_stimuli() for trial in self.trials.values()])

BEH_COLUMNS = ['trial_number', 'stimulus', 'stimulus_coords', 'memory_set', 'memory_set_locations', 'target_present', 'target', 'target_location', 'correct_response', 'was_answered'] class-attribute instance-attribute

Columns explanation: - trial_number: The number of the trial, in the order they were presented. They start from 0. - stimulus: The filename of the stimulus presented. - stimulus_coords: The coordinates of the stimulus presented. It should be a tuple containing the x, y of the top-left corner of the stimulus and the x, y of the bottom-right corner. - memory_set: The set of items memorized by the participant. It should be a list of strings. Each string should be the filename of the stimulus. - memory_set_locations: The locations of the items memorized by the participant. It should be a list of tuples. Each tuple should contain bounding boxes of the items memorized by the participant. The bounding boxes should be in the format (x1, y1, x2, y2), where (x1, y1) is the top-left corner and (x2, y2) is the bottom-right corner. - target_present: Whether one of the items is present in the stimulus. It should be a boolean. - target: The filename of the target item. It should be a string. If target_present is False, the value for this column will not be taken into account. - target_location: The location of the target item. It should be a tuple containing the bounding box of the target item. The bounding box should be in the format (x1, y1, x2, y2), where (x1, y1) is the top-left corner and (x2, y2) is the bottom-right corner. If target_present is False, the value for this column will not be taken into account. - correct_response: The correct response for the trial. It should be a boolean. - was_answered: Whether the trial was answered by the participant. It should be a boolean.

Notice that you can get the actual response of the user by using the "correct_response" and "target_present" columns. For all of the heights, widths and locations of the items, the values should be in pixels and according to the screen itself.

VisualSearchTrial

Bases: Trial

Source code in pyxations/analysis/visual_search.py

class VisualSearchTrial(Trial):

    def __init__(self, trial_number, session, samples, fix, sacc, blink, events_path, behavior_data, search_phase_name, memorization_phase_name):
        super().__init__(trial_number, session, samples, fix, sacc, blink, events_path)

        trial_data = behavior_data.filter(pl.col("trial_number") == trial_number)

        self._target_present = trial_data.select("target_present").item()
        self._target = trial_data.select("target").item()

        if self._target_present:
            self._target_location = _as(trial_data.select("target_location").item(), tuple)

        self._correct_response = trial_data.select("correct_response").item()
        self._stimulus = trial_data.select("stimulus").item()
        self._stimulus_coords = _as(trial_data.select("stimulus_coords").item(), tuple)

        self._memory_set = _as(trial_data.select("memory_set").item(), list)
        self._memory_set_locations = _as(trial_data.select("memory_set_locations").item(), list)
        self._search_phase_name = search_phase_name
        self._memorization_phase_name = memorization_phase_name
        self._was_answered = trial_data.select("was_answered").item()

    @property
    def target(self):
        return self._target

    @property
    def target_location(self):
        return self._target_location

    @property
    def target_present(self):
        return self._target_present

    @property
    def correct_response(self):
        return self._correct_response

    @property
    def memory_set_size(self):
        return len(self._memory_set)

    @property
    def memory_set_locations(self):
        return self._memory_set_locations

    @property
    def memory_set(self):
        return self._memory_set

    @property
    def stimulus(self):
        return self._stimulus

    @property
    def stimulus_coords(self):
        return self._stimulus_coords

    @property
    def was_answered(self):
        return self._was_answered

    def save_rts(self):
        if hasattr(self, "_rts"):
            return

        # Filter out empty phase rows
        filtered = self._samples.filter(pl.col("phase") != "")

        # Calculate RT as the difference between last and first tSample per phase
        self._rts = (
            filtered
            .group_by("phase")
            .agg((pl.col("tSample").max() - pl.col("tSample").min()).alias("rt"))
            .with_columns([
                pl.lit(self.trial_number).alias("trial_number"),
                pl.lit(len(self._memory_set)).alias("memory_set_size"),
                pl.lit(self._target_present).alias("target_present"),
                pl.lit(self._correct_response).alias("correct_response"),
                pl.lit(self._stimulus).alias("stimulus"),
                pl.lit(self._target).alias("target"),
                pl.lit(self._was_answered).alias("was_answered"),
            ])
        )


    def fixations(self):
        fixations = super().fixations().with_columns([
            pl.lit(self._target_present).alias("target_present"),
            pl.lit(self._correct_response).alias("correct_response"),
            pl.lit(self._stimulus).alias("stimulus"),
            pl.lit(self._target).alias("target"),
            pl.lit(self._memory_set).alias("memory_set"),
        ])
        return fixations


    def saccades(self):
        saccades = super().saccades().with_columns([
            pl.lit(self._target_present).alias("target_present"),
            pl.lit(self._correct_response).alias("correct_response"),
            pl.lit(self._stimulus).alias("stimulus"),
            pl.lit(self._target).alias("target"),
            pl.lit(self._memory_set).alias("memory_set"),
        ])
        return saccades



    def search_fixations(self):
        return self.fixations().filter(pl.col("phase") == self._search_phase_name).sort(by="tStart")

    def memorization_fixations(self):
        return self.fixations().filter(pl.col("phase") == self._memorization_phase_name).sort(by="tStart")

    def search_saccades(self):
        return self.saccades().filter(pl.col("phase") == self._search_phase_name).sort(by="tStart")

    def memorization_saccades(self):
        return self.saccades().filter(pl.col("phase") == self._memorization_phase_name).sort(by="tStart")

    def search_samples(self):
        return self.samples().filter(pl.col("phase") == self._search_phase_name).sort(by="tSample")

    def memorization_samples(self):
        return self.samples().filter(pl.col("phase") == self._memorization_phase_name).sort(by="tSample")

    def scanpath_by_stimuli(self):
        return {"fixations": self.search_fixations(), "stimulus": self._stimulus,"correct_response":self._correct_response,"target_present":self._target_present,"memory_set_size":len(self._memory_set)}

    def plot_scanpath(self, screen_height, screen_width, **kwargs):
        '''
        Plots the scanpath of the trial. The scanpath will be plotted in two phases: the search phase and the memorization phase.
        The search phase will be plotted with the stimulus and the memorization phase will be plotted with the items memorized by the participant.
        The search phase will have the fixations and saccades of the trial, while the memorization phase will only have the fixations.
        The names of the phases should be the same ones used in the computation of the derivatives.
        If you don't really care about the memorization phase, you can pass None as an argument.

        '''
        vis = Visualization(self.events_path, self.detection_algorithm)
        (self.events_path / "plots").mkdir(parents=True, exist_ok=True)


        phase_data = {self._search_phase_name:{}, self._memorization_phase_name:{}}
        dataset_parent_folder = self.events_path.parent.parent.parent.parent
        phase_data[self._search_phase_name]["img_paths"] = [dataset_parent_folder / STIMULI_FOLDER / self._stimulus]
        phase_data[self._search_phase_name]["img_plot_coords"] = [self._stimulus_coords]
        if self._memorization_phase_name is not None:
            phase_data[self._memorization_phase_name]["img_paths"] = [dataset_parent_folder / ITEMS_FOLDER / img for img in self._memory_set]
            phase_data[self._memorization_phase_name]["img_plot_coords"] = self._memory_set_locations

        # If the target is present add the "bbox" to the search_phase phase as a key-value pair
        if self._target_present:
            phase_data[self._search_phase_name]["bbox"] = self._target_location
        vis.scanpath(fixations=self._fix,phase_data=phase_data, saccades=self._sacc, samples=self._samples, screen_height=screen_height, screen_width=screen_width, 
                      folder_path=self.events_path / "plots", **kwargs)

    def plot_animation(self, screen_height, screen_width, video_path=None, background_image_path=None, **kwargs):
        """
        Create an animated visualization of eye-tracking data for this trial.

        When a video is provided, the animation syncs gaze samples with video frames.
        When no video is provided, gaze points are animated on a grey background,
        or a provided background image (e.g., the stimulus image), using the sample 
        timestamps for timing.

        Parameters
        ----------
        screen_height, screen_width
            Stimulus resolution in pixels.
        video_path
            Path to a video file. If provided, gaze is overlaid on video frames.
        background_image_path
            Path to a background image. Only used when video_path is None.
            If None and no video, uses the search stimulus as background if available,
            otherwise uses a grey background.
        **kwargs
            Additional arguments passed to Visualization.plot_animation():
            - folder_path: Directory to save the animation
            - tmin, tmax: Time window in ms
            - seconds_to_show: Limit animation to first N seconds
            - scale_factor: Resolution scaling (default 0.5)
            - gaze_radius: Gaze point radius in pixels
            - gaze_color: RGB tuple for gaze color
            - fps: Animation frames per second
            - output_format: "html" (default), "mp4", "gif", or "matplotlib"
            - display: If True, return HTML for notebook display

        Returns
        -------
        IPython.display.HTML or None
            Returns HTML animation if display=True and output_format="html".
            For output_format="matplotlib", displays in a GUI window and returns None.
        """
        vis = Visualization(self.events_path, self.detection_algorithm)
        (self.events_path / "plots").mkdir(parents=True, exist_ok=True)

        # Set default folder_path if not provided
        if 'folder_path' not in kwargs:
            kwargs['folder_path'] = self.events_path / "plots"

        # If no background image provided and no video, try to use the stimulus
        if video_path is None and background_image_path is None:
            dataset_parent_folder = self.events_path.parent.parent.parent.parent
            stimulus_path = dataset_parent_folder / STIMULI_FOLDER / self._stimulus
            if stimulus_path.exists():
                background_image_path = stimulus_path

        return vis.plot_animation(
            samples=self._samples,
            screen_height=screen_height,
            screen_width=screen_width,
            video_path=video_path,
            background_image_path=background_image_path,
            **kwargs
        )

plot_animation(screen_height, screen_width, video_path=None, background_image_path=None, **kwargs)

Create an animated visualization of eye-tracking data for this trial.

When a video is provided, the animation syncs gaze samples with video frames. When no video is provided, gaze points are animated on a grey background, or a provided background image (e.g., the stimulus image), using the sample timestamps for timing.

Parameters:

Name	Type	Description	Default
screen_height		Stimulus resolution in pixels.	required
screen_width		Stimulus resolution in pixels.	required
video_path		Path to a video file. If provided, gaze is overlaid on video frames.	None
background_image_path		Path to a background image. Only used when video_path is None. If None and no video, uses the search stimulus as background if available, otherwise uses a grey background.	None
**kwargs		Additional arguments passed to Visualization.plot_animation(): - folder_path: Directory to save the animation - tmin, tmax: Time window in ms - seconds_to_show: Limit animation to first N seconds - scale_factor: Resolution scaling (default 0.5) - gaze_radius: Gaze point radius in pixels - gaze_color: RGB tuple for gaze color - fps: Animation frames per second - output_format: "html" (default), "mp4", "gif", or "matplotlib" - display: If True, return HTML for notebook display	{}

Returns:

Type	Description
HTML or None	Returns HTML animation if display=True and output_format="html". For output_format="matplotlib", displays in a GUI window and returns None.

Source code in pyxations/analysis/visual_search.py

def plot_animation(self, screen_height, screen_width, video_path=None, background_image_path=None, **kwargs):
    """
    Create an animated visualization of eye-tracking data for this trial.

    When a video is provided, the animation syncs gaze samples with video frames.
    When no video is provided, gaze points are animated on a grey background,
    or a provided background image (e.g., the stimulus image), using the sample 
    timestamps for timing.

    Parameters
    ----------
    screen_height, screen_width
        Stimulus resolution in pixels.
    video_path
        Path to a video file. If provided, gaze is overlaid on video frames.
    background_image_path
        Path to a background image. Only used when video_path is None.
        If None and no video, uses the search stimulus as background if available,
        otherwise uses a grey background.
    **kwargs
        Additional arguments passed to Visualization.plot_animation():
        - folder_path: Directory to save the animation
        - tmin, tmax: Time window in ms
        - seconds_to_show: Limit animation to first N seconds
        - scale_factor: Resolution scaling (default 0.5)
        - gaze_radius: Gaze point radius in pixels
        - gaze_color: RGB tuple for gaze color
        - fps: Animation frames per second
        - output_format: "html" (default), "mp4", "gif", or "matplotlib"
        - display: If True, return HTML for notebook display

    Returns
    -------
    IPython.display.HTML or None
        Returns HTML animation if display=True and output_format="html".
        For output_format="matplotlib", displays in a GUI window and returns None.
    """
    vis = Visualization(self.events_path, self.detection_algorithm)
    (self.events_path / "plots").mkdir(parents=True, exist_ok=True)

    # Set default folder_path if not provided
    if 'folder_path' not in kwargs:
        kwargs['folder_path'] = self.events_path / "plots"

    # If no background image provided and no video, try to use the stimulus
    if video_path is None and background_image_path is None:
        dataset_parent_folder = self.events_path.parent.parent.parent.parent
        stimulus_path = dataset_parent_folder / STIMULI_FOLDER / self._stimulus
        if stimulus_path.exists():
            background_image_path = stimulus_path

    return vis.plot_animation(
        samples=self._samples,
        screen_height=screen_height,
        screen_width=screen_width,
        video_path=video_path,
        background_image_path=background_image_path,
        **kwargs
    )

plot_scanpath(screen_height, screen_width, **kwargs)

Plots the scanpath of the trial. The scanpath will be plotted in two phases: the search phase and the memorization phase. The search phase will be plotted with the stimulus and the memorization phase will be plotted with the items memorized by the participant. The search phase will have the fixations and saccades of the trial, while the memorization phase will only have the fixations. The names of the phases should be the same ones used in the computation of the derivatives. If you don't really care about the memorization phase, you can pass None as an argument.

Source code in pyxations/analysis/visual_search.py

def plot_scanpath(self, screen_height, screen_width, **kwargs):
    '''
    Plots the scanpath of the trial. The scanpath will be plotted in two phases: the search phase and the memorization phase.
    The search phase will be plotted with the stimulus and the memorization phase will be plotted with the items memorized by the participant.
    The search phase will have the fixations and saccades of the trial, while the memorization phase will only have the fixations.
    The names of the phases should be the same ones used in the computation of the derivatives.
    If you don't really care about the memorization phase, you can pass None as an argument.

    '''
    vis = Visualization(self.events_path, self.detection_algorithm)
    (self.events_path / "plots").mkdir(parents=True, exist_ok=True)


    phase_data = {self._search_phase_name:{}, self._memorization_phase_name:{}}
    dataset_parent_folder = self.events_path.parent.parent.parent.parent
    phase_data[self._search_phase_name]["img_paths"] = [dataset_parent_folder / STIMULI_FOLDER / self._stimulus]
    phase_data[self._search_phase_name]["img_plot_coords"] = [self._stimulus_coords]
    if self._memorization_phase_name is not None:
        phase_data[self._memorization_phase_name]["img_paths"] = [dataset_parent_folder / ITEMS_FOLDER / img for img in self._memory_set]
        phase_data[self._memorization_phase_name]["img_plot_coords"] = self._memory_set_locations

    # If the target is present add the "bbox" to the search_phase phase as a key-value pair
    if self._target_present:
        phase_data[self._search_phase_name]["bbox"] = self._target_location
    vis.scanpath(fixations=self._fix,phase_data=phase_data, saccades=self._sacc, samples=self._samples, screen_height=screen_height, screen_width=screen_width, 
                  folder_path=self.events_path / "plots", **kwargs)

visualization

Plotting utilities for scanpaths, fixations, saccades and raw samples.

visualization

Visualization

Source code in pyxations/visualization/visualization.py

class Visualization():
    def __init__(self, derivatives_folder_path,events_detection_algorithm):
        self.derivatives_folder_path = Path(derivatives_folder_path)
        if events_detection_algorithm not in EYE_MOVEMENT_DETECTION_DICT and events_detection_algorithm != 'eyelink':
            raise ValueError(f"Detection algorithm {events_detection_algorithm} not found.")
        self.events_detection_folder = Path(events_detection_algorithm+'_events')

    def scanpath(
        self,
        fixations: pl.DataFrame,
        screen_height: int,
        screen_width: int,
        folder_path: str | Path | None = None,
        tmin: int | None = None,
        tmax: int | None = None,
        saccades: pl.DataFrame | None = None,
        samples: pl.DataFrame | None = None,
        phase_data: dict[str, dict] | None = None,
        display: bool = True,
    ):
        """
        Fast scan‑path visualiser.

        • **Vectorised**: no per‑row Python loops  
        • **Single pass** phase grouping  
        • Uses `BrokenBarHCollection` for fixation spans  
        • Optional asynchronous PNG write via ThreadPoolExecutor (drop‑in‑ready, see comment)

        Parameters
        ----------
        fixations
            Polars DataFrame with at least `tStart`, `duration`, `xAvg`, `yAvg`, `phase`.
        screen_height, screen_width
            Stimulus resolution in pixels.
        folder_path
            Directory where 1 PNG per phase will be stored.  If *None*, nothing is saved.
        tmin, tmax
            Time window in **ms**.  If both `None`, the whole trial is plotted.
        saccades
            Polars DataFrame with `tStart`, `phase`, …  (optional).
        samples
            Polars DataFrame with gaze traces (`tSample`, `LX`, `LY`, `RX`, `RY` or
            `X`, `Y`) (optional).
        phase_data
            Per‑phase extras::

                {
                    "search": {
                        "img_paths": [...],
                        "img_plot_coords": [(x1,y1,x2,y2), ...],
                        "bbox": (x1,y1,x2,y2),
                    },
                    ...
                }

        display
            If *False* the figure canvas is never shown (faster for batch jobs).
        """


        # ------------- small helpers ------------------------------------------------
        def _make_axes(plot_samples: bool):
            if plot_samples:
                fig, (ax_main, ax_gaze) = plt.subplots(
                    2, 1, height_ratios=(4, 1), figsize=(10, 6), sharex=False
                )
            else:
                fig, ax_main = plt.subplots(figsize=(10, 6))
                ax_gaze = None
            ax_main.set_xlim(0, screen_width)
            ax_main.set_ylim(screen_height, 0)
            return fig, ax_main, ax_gaze

        def _maybe_cache_img(path):
            """Load image from disk with a small LRU cache."""

            # Cache hit: move to the end (most recently used)
            if path in _img_cache:
                img = _img_cache.pop(path)
                _img_cache[path] = img
                return img

            # Cache miss: load image
            img = mpimg.imread(path)

            # Optional: reduce memory if image is float64 in [0, 1]
            if isinstance(img, np.ndarray) and img.dtype == np.float64:
                img = (img * 255).clip(0, 255).astype(np.uint8)

            # Insert into cache
            _img_cache[path] = img

            # If cache too big, drop least recently used item
            if len(_img_cache) > _MAX_CACHE_ITEMS:
                _img_cache.popitem(last=False)  # pops the oldest inserted item

            return img

        # ---------------------------------------------------------------------------
        plot_saccades = saccades is not None
        plot_samples = samples is not None
        _img_cache = OrderedDict()
        _MAX_CACHE_ITEMS = 8  # or 5, 10, etc. Tune as you like.

        trial_idx = fixations["trial_number"][0]

        # ---- time filter ----------------------------------------------------------
        if tmin is not None and tmax is not None:
            fixations = fixations.filter(pl.col("tStart").is_between(tmin, tmax))
            if plot_saccades:
                saccades = saccades.filter(pl.col("tStart").is_between(tmin, tmax))
            if plot_samples:
                samples = samples.filter(pl.col("tSample").is_between(tmin, tmax))

        # remove empty phase markings
        fixations = fixations.filter(pl.col("phase") != "")
        if plot_saccades:
            saccades = saccades.filter(pl.col("phase") != "")
        if plot_samples:
            samples = samples.filter(pl.col("phase") != "")

        # ---- split once by phase --------------------------------------------------
        fix_by_phase = fixations.partition_by("phase", as_dict=True)
        sac_by_phase = (
            saccades.partition_by("phase", as_dict=True) if plot_saccades else {}
        )
        samp_by_phase = (
            samples.partition_by("phase", as_dict=True) if plot_samples else {}
        )

        # colour map shared across phases
        cmap = plt.cm.rainbow

        # ---- build & draw ---------------------------------------------------------
        # optional async saver (uncomment if you save hundreds of files)
        from concurrent.futures import ThreadPoolExecutor
        saver = ThreadPoolExecutor(max_workers=4) if folder_path else None

        if not display:
            plt.ioff()

        for phase, phase_fix in fix_by_phase.items():
            if phase_fix.is_empty():
                continue

            # ---------- vectors (zero‑copy) -----------------
            fx, fy, fdur = phase_fix.select(["xAvg", "yAvg", "duration"]).to_numpy().T
            n_fix = fx.size
            fix_idx = np.arange(1, n_fix + 1)

            norm = mplcolors.BoundaryNorm(np.arange(1, n_fix + 2), cmap.N)

            # saccades
            sac_t = (
                sac_by_phase[phase]["tStart"].to_numpy()
                if plot_saccades and phase in sac_by_phase
                else np.empty(0)
            )

            # samples
            if plot_samples and phase in samp_by_phase and samp_by_phase[phase].height:
                samp_phase = samp_by_phase[phase]
                t0 = samp_phase["tSample"][0]
                ts = (samp_phase["tSample"].to_numpy() - t0) 
                get = samp_phase.get_column
                lx = get("LX").to_numpy() if "LX" in samp_phase.columns else None
                ly = get("LY").to_numpy() if "LY" in samp_phase.columns else None
                rx = get("RX").to_numpy() if "RX" in samp_phase.columns else None
                ry = get("RY").to_numpy() if "RY" in samp_phase.columns else None
                gx = get("X").to_numpy() if "X" in samp_phase.columns else None
                gy = get("Y").to_numpy() if "Y" in samp_phase.columns else None
            else:
                t0 = None

            # ---------- figure -----------------------------
            fig, ax_main, ax_gaze = _make_axes(plot_samples and t0 is not None)
            # scatter fixations
            sc = ax_main.scatter(
                fx,
                fy,
                c=fix_idx,
                s=fdur,
                cmap=cmap,
                norm=norm,
                alpha=0.5,
                zorder=2,
            )
            fig.colorbar(
                sc,
                ax=ax_main,
                ticks=[1, n_fix // 2 if n_fix > 2 else n_fix, n_fix],
                fraction=0.046,
                pad=0.04,
            ).set_label("# of fixation")

            # ---------- stimulus imagery / bbox ------------
            if phase_data and phase[0] in phase_data:
                pdict = phase_data[phase[0]]
                coords = pdict.get("img_plot_coords") or []
                bbox = pdict.get('bbox',None) 
                for img_path, box in zip(pdict.get("img_paths", []), coords):

                    ax_main.imshow(_maybe_cache_img(img_path), extent=[box[0], box[2], box[3], box[1]], zorder=0)
                if bbox is not None:
                    x1, y1, x2, y2 = bbox
                    ax_main.plot([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], color='red', linewidth=1.5, zorder=3)

            # ---------- gaze traces ------------------------
            if ax_gaze is not None:
                if lx is not None:
                    ax_main.plot(lx, ly, "--", color="C0", zorder=1)
                    ax_gaze.plot(ts, lx, label="Left X")
                    ax_gaze.plot(ts, ly, label="Left Y")
                if rx is not None:
                    ax_main.plot(rx, ry, "--", color="k", zorder=1)
                    ax_gaze.plot(ts, rx, label="Right X")
                    ax_gaze.plot(ts, ry, label="Right Y")
                if gx is not None:
                    ax_main.plot(gx, gy, "--", color="k", zorder=1, alpha=0.6)
                    ax_gaze.plot(ts, gx, label="X")
                    ax_gaze.plot(ts, gy, label="Y")

                # fixation spans
                bars   = np.c_[phase_fix['tStart'].to_numpy() - t0,
                            phase_fix['duration'].to_numpy()]
                height = ax_gaze.get_ylim()[1] - ax_gaze.get_ylim()[0]
                colors = cmap(norm(fix_idx))

                # Draw all bars in one call; no BrokenBarHCollection import needed
                ax_gaze.broken_barh(bars, (0, height), facecolors=colors, alpha=0.4)
                # saccades
                if sac_t.size:
                    ymin, ymax = ax_gaze.get_ylim()
                    ax_gaze.vlines(
                        sac_t - t0,
                        ymin,
                        ymax,
                        colors="red",
                        linestyles="--",
                        linewidth=0.8,
                    )

                # tidy gaze axis
                h, l = ax_gaze.get_legend_handles_labels()
                by_label = {lab: hdl for hdl, lab in zip(h, l)}
                ax_gaze.legend(
                    by_label.values(),
                    by_label.keys(),
                    loc="center left",
                    bbox_to_anchor=(1, 0.5),
                )
                ax_gaze.set_ylabel("Gaze")
                ax_gaze.set_xlabel("Time [s]")

            fig.tight_layout()

            # ---------- save / show ------------------------
            if folder_path:
                scan_name = f"scanpath_{trial_idx}"
                if tmin is not None and tmax is not None:
                    scan_name += f"_{tmin}_{tmax}"
                out = Path(folder_path) / f"{scan_name}_{phase[0]}.png"
                fig.savefig(out, dpi=150)
                if saver:  saver.submit(fig.savefig, out, dpi=150)

            if display:
                plt.show()
            plt.close(fig)

        if not display:
            plt.ion()


    def fix_duration(self,fixations:pl.DataFrame,axs=None):

        ax = axs
        if ax is None:
            fig, ax = plt.subplots()

        ax.hist(fixations.select(pl.col('duration')).to_numpy().ravel(), bins=100, edgecolor='black', linewidth=1.2, density=True)
        ax.set_title('Fixation duration')
        ax.set_xlabel('Time (ms)')
        ax.set_ylabel('Density')


    def sacc_amplitude(self,saccades:pl.DataFrame,axs=None):

        ax = axs
        if ax is None:
            fig, ax = plt.subplots()

        saccades_amp = saccades.select(pl.col('ampDeg')).to_numpy().ravel()
        ax.hist(saccades_amp, bins=100, range=(0, 20), edgecolor='black', linewidth=1.2, density=True)
        ax.set_title('Saccades amplitude')
        ax.set_xlabel('Amplitude (deg)')
        ax.set_ylabel('Density')


    def sacc_direction(self,saccades:pl.DataFrame,axs=None,figs=None):

        ax = axs
        if ax is None:
            fig = plt.figure()
            ax = plt.subplot(polar=True)
        else:
            ax.set_axis_off()
            ax = figs.add_subplot(2, 2, 3, projection='polar')
        if 'deg' not in saccades.columns or 'dir' not in saccades.columns:
            raise ValueError('Compute saccades direction first by using saccades_direction function from the PreProcessing module.')
        # Convert from deg to rad
        saccades_rad = saccades.select(pl.col('deg')).to_numpy().ravel() * np.pi / 180

        n_bins = 24
        ang_hist, bin_edges = np.histogram(saccades_rad, bins=24, density=True)
        bin_centers = [np.mean((bin_edges[i], bin_edges[i+1])) for i in range(len(bin_edges) - 1)]

        bars = ax.bar(bin_centers, ang_hist, width=2*np.pi/n_bins, bottom=0.0, alpha=0.4, edgecolor='black')
        ax.set_title('Saccades direction')
        ax.set_yticklabels([])

        for r, bar in zip(ang_hist, bars):
            bar.set_facecolor(plt.cm.Blues(r / np.max(ang_hist)))


    def sacc_main_sequence(self,saccades:pl.DataFrame,axs=None, hline=None):

        ax = axs
        if ax is None:
            fig, ax = plt.subplots()
        # Logarithmic bins
        XL = np.log10(25)  # Adjusted to fit the xlim
        YL = np.log10(1000)  # Adjusted to fit the ylim

        saccades_peak_vel = saccades.select(pl.col('vPeak')).to_numpy().ravel()
        saccades_amp = saccades.select(pl.col('ampDeg')).to_numpy().ravel()

        # Create a 2D histogram with logarithmic bins
        ax.hist2d(saccades_amp, saccades_peak_vel, bins=[np.logspace(-1, XL, 50), np.logspace(0, YL, 50)])

        if hline:
            ax.hlines(y=hline, xmin=ax.get_xlim()[0], xmax=ax.get_xlim()[1], colors='grey', linestyles='--', label=hline)
            ax.legend()
        ax.set_yscale('log')
        ax.set_xscale('log')
        ax.set_title('Main sequence')
        ax.set_xlabel('Amplitude (deg)')
        ax.set_ylabel('Peak velocity (deg)')
         # Set the limits of the axes
        ax.set_xlim(0.1, 25)
        ax.set_ylim(10, 1000)
        ax.set_aspect('equal')


    def plot_multipanel(
            self,
            fixations: pl.DataFrame,
            saccades: pl.DataFrame,
            display: bool = True
        ) -> None:
        """
        Create a 2×2 multi‑panel diagnostic plot for every non‑empty
        phase label and save it as PNG in
        <derivatives_folder_path>/<events_detection_folder>/plots/.
        """
        # ── paths & matplotlib style ────────────────────────────────
        folder_path: Path = (
            self.derivatives_folder_path
            / self.events_detection_folder
            / "plots"
        )
        folder_path.mkdir(parents=True, exist_ok=True)
        plt.rcParams.update({"font.size": 12})

        # ── drop practice / invalid trials ─────────────────────────
        fixations = fixations.filter(pl.col("trial_number") != -1)
        saccades  = saccades.filter(pl.col("trial_number") != -1)

        # ── collect valid phase labels (skip empty string) ─────────
        phases = (
            fixations
            .select(pl.col("phase").filter(pl.col("phase") != ""))
            .unique()           # unique values in this Series
            .to_series()
            .to_list()          # plain Python list of strings
        )

        # ── one figure per phase ───────────────────────────────────
        for phase in phases:
            fix_phase   = fixations.filter(pl.col("phase") == phase)
            sacc_phase  = saccades.filter(pl.col("phase") == phase)

            fig, axs = plt.subplots(2, 2, figsize=(12, 7))

            self.fix_duration(fix_phase , axs=axs[0, 0])
            self.sacc_main_sequence(sacc_phase, axs=axs[1, 1])
            self.sacc_direction(sacc_phase, axs=axs[1, 0], figs=fig)
            self.sacc_amplitude(sacc_phase, axs=axs[0, 1])

            fig.tight_layout()
            plt.savefig(folder_path / f"multipanel_{phase}.png")
            if display:
                plt.show()
            plt.close()

    def plot_animation(
        self,
        samples: pl.DataFrame,
        screen_height: int,
        screen_width: int,
        video_path: str | Path | None = None,
        background_image_path: str | Path | None = None,
        folder_path: str | Path | None = None,
        tmin: int | None = None,
        tmax: int | None = None,
        seconds_to_show: float | None = None,
        scale_factor: float = 0.5,
        gaze_radius: int = 10,
        gaze_color: tuple = (255, 0, 0),
        fps: float | None = None,
        output_format: str = "matplotlib",
        display: bool = True,
    ):
        """
        Create an animated visualization of eye-tracking data.

        When a video is provided, the animation syncs gaze samples with video frames.
        When no video is provided, gaze points are animated on a grey background or
        a provided background image, using the sample timestamps for timing.

        Parameters
        ----------
        samples
            Polars DataFrame with gaze samples. Must contain 'tSample' and gaze
            position columns ('X', 'Y' or 'LX', 'LY', 'RX', 'RY').
        screen_height, screen_width
            Stimulus resolution in pixels.
        video_path
            Path to a video file. If provided, gaze is overlaid on video frames.
        background_image_path
            Path to a background image. Only used when video_path is None.
            If both are None, a grey background is used.
        folder_path
            Directory where the animation will be saved. If None, nothing is saved.
            The file format depends on `output_format`.
        tmin, tmax
            Time window in **ms**. If both None, the whole trial is plotted.
        seconds_to_show
            Limit the animation to the first N seconds. If None, shows all available data.
        scale_factor
            Resolution scaling factor (1.0 = original, 0.5 = half resolution).
        gaze_radius
            Radius of the gaze point circle in pixels (before scaling).
        gaze_color
            RGB tuple for gaze point color.
        fps
            Frames per second for the animation. If None:
            - With video: uses the video's native FPS
            - Without video: defaults to 60 FPS
        output_format
            Output format for saved animations:
            - "html": Interactive HTML file (default, works in browsers)
            - "mp4": Video file (requires ffmpeg)
            - "gif": Animated GIF file (requires pillow)
            - "matplotlib": Show in matplotlib GUI window (blocking)
        display
            If True and output_format is "html", returns an HTML object for notebooks.
            If output_format is "matplotlib", this is ignored (always shows window).
            If False, only saves to file (if folder_path is provided).

        Returns
        -------
        IPython.display.HTML or None
            Returns HTML animation if display=True and output_format="html", otherwise None.
        """
        try:
            import cv2
            from matplotlib.animation import FuncAnimation
            import matplotlib as mpl
            mpl.rcParams['animation.embed_limit'] = 100
        except ImportError as e:
            raise ImportError(
                f"Missing required dependency for animation: {e}. "
                "Please install cv2 (opencv-python)."
            )

        # Validate output_format
        valid_formats = ["html", "mp4", "gif", "matplotlib"]
        if output_format not in valid_formats:
            raise ValueError(f"output_format must be one of {valid_formats}, got '{output_format}'")

        # ---- Determine gaze columns ----
        if "X" in samples.columns and "Y" in samples.columns:
            x_col, y_col = "X", "Y"
        elif "LX" in samples.columns and "LY" in samples.columns:
            x_col, y_col = "LX", "LY"
        elif "RX" in samples.columns and "RY" in samples.columns:
            x_col, y_col = "RX", "RY"
        else:
            raise ValueError("Samples DataFrame must contain gaze columns (X, Y) or (LX, LY) or (RX, RY)")

        # ---- Time filter ----
        if tmin is not None and tmax is not None:
            samples = samples.filter(pl.col("tSample").is_between(tmin, tmax))

        if samples.is_empty():
            raise ValueError("No samples available after time filtering")

        # ---- Drop NaN gaze values ----
        samples = samples.filter(pl.col(x_col).is_not_null() & pl.col(y_col).is_not_null())

        # ---- Calculate scaled dimensions ----
        scaled_width = int(screen_width * scale_factor)
        scaled_height = int(screen_height * scale_factor)

        trial_idx = samples["trial_number"][0] if "trial_number" in samples.columns else 0

        # ================= WITH VIDEO =================
        if video_path is not None:
            video_path = Path(video_path)
            if not video_path.exists():
                raise FileNotFoundError(f"Video file not found: {video_path}")

            cap = cv2.VideoCapture(str(video_path))
            video_fps = cap.get(cv2.CAP_PROP_FPS)
            total_frames = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))
            video_width = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
            video_height = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))

            if fps is None:
                fps = video_fps

            # Calculate time to frame mapping
            t_start = samples["tSample"].min()
            t_end = samples["tSample"].max()
            trial_duration = t_end - t_start

            # Create frame-to-time mapping
            frame_edges = np.linspace(t_start, t_end, total_frames + 1)
            frame_times = ((frame_edges[:-1] + frame_edges[1:]) / 2).astype(int)

            # Build a lookup: frame_index -> list of gaze points
            samples_np = samples.select([x_col, y_col, "tSample"]).to_numpy()
            gaze_by_frame = {i: [] for i in range(total_frames)}

            for x, y, t in samples_np:
                # Find the closest frame
                frame_idx = np.searchsorted(frame_times, t, side='right') - 1
                frame_idx = max(0, min(frame_idx, total_frames - 1))
                gaze_by_frame[frame_idx].append((x, y))

            # Limit frames if seconds_to_show is set
            frames_to_show = total_frames
            if seconds_to_show is not None:
                frames_to_show = min(int(fps * seconds_to_show), total_frames)

            # Reset video
            cap.set(cv2.CAP_PROP_POS_FRAMES, 0)

            # Create figure
            fig, ax = plt.subplots(figsize=(10 * scale_factor, 6 * scale_factor))
            ax.axis('off')

            # Initialize with first frame
            ret, frame = cap.read()
            if not ret:
                cap.release()
                raise RuntimeError("Could not read first frame from video")

            frame_resized = cv2.resize(frame, (scaled_width, scaled_height), interpolation=cv2.INTER_AREA)
            frame_rgb = cv2.cvtColor(frame_resized, cv2.COLOR_BGR2RGB)
            im = ax.imshow(frame_rgb)

            def update_frame_video(frame_idx):
                ret, frame = cap.read()
                if not ret:
                    return [im]

                frame_resized = cv2.resize(frame, (scaled_width, scaled_height), interpolation=cv2.INTER_AREA)
                frame_rgb = cv2.cvtColor(frame_resized, cv2.COLOR_BGR2RGB)

                # Draw gaze points for this frame
                for gx, gy in gaze_by_frame.get(frame_idx, []):
                    scaled_x = int(gx * scale_factor)
                    scaled_y = int(gy * scale_factor)
                    if 0 <= scaled_x < scaled_width and 0 <= scaled_y < scaled_height:
                        radius = max(3, int(gaze_radius * scale_factor))
                        cv2.circle(frame_rgb, (scaled_x, scaled_y), radius=radius, color=gaze_color, thickness=-1)

                im.set_array(frame_rgb)
                return [im]

            anim = FuncAnimation(fig, update_frame_video, frames=frames_to_show,
                                 interval=1000/fps, blit=True, repeat=True)

        # ================= WITHOUT VIDEO =================
        else:
            if fps is None:
                fps = 60  # Default FPS for sample-based animation

            # Prepare background
            if background_image_path is not None:
                bg_path = Path(background_image_path)
                if not bg_path.exists():
                    raise FileNotFoundError(f"Background image not found: {bg_path}")
                bg_img = mpimg.imread(str(bg_path))
                if bg_img.dtype == np.float64:
                    bg_img = (bg_img * 255).clip(0, 255).astype(np.uint8)
                # Resize background to match screen dimensions then scale
                bg_img = cv2.resize(bg_img, (scaled_width, scaled_height), interpolation=cv2.INTER_AREA)
            else:
                # Grey background
                bg_img = np.ones((scaled_height, scaled_width, 3), dtype=np.uint8) * 128

            # Get time range
            t_start = samples["tSample"].min()
            t_end = samples["tSample"].max()
            trial_duration = t_end - t_start

            # Limit duration if seconds_to_show is set
            if seconds_to_show is not None:
                t_end = min(t_end, t_start + int(seconds_to_show * 1000))
                samples = samples.filter(pl.col("tSample") <= t_end)
                trial_duration = t_end - t_start

            # Calculate total frames based on duration and fps
            total_frames = int((trial_duration / 1000) * fps)
            if total_frames < 1:
                total_frames = 1

            # Create time bins for each animation frame
            frame_times = np.linspace(t_start, t_end, total_frames + 1)

            # Build gaze lookup by frame
            samples_np = samples.select([x_col, y_col, "tSample"]).to_numpy()
            gaze_by_frame = {i: [] for i in range(total_frames)}

            for x, y, t in samples_np:
                frame_idx = np.searchsorted(frame_times, t, side='right') - 1
                frame_idx = max(0, min(frame_idx, total_frames - 1))
                gaze_by_frame[frame_idx].append((x, y))

            # Create figure
            fig, ax = plt.subplots(figsize=(10 * scale_factor, 6 * scale_factor))
            ax.axis('off')

            # Initialize with background
            im = ax.imshow(bg_img.copy())

            def update_frame_no_video(frame_idx):
                # Start with fresh background copy
                frame_rgb = bg_img.copy()

                # Draw gaze points for this frame
                for gx, gy in gaze_by_frame.get(frame_idx, []):
                    scaled_x = int(gx * scale_factor)
                    scaled_y = int(gy * scale_factor)
                    if 0 <= scaled_x < scaled_width and 0 <= scaled_y < scaled_height:
                        radius = max(3, int(gaze_radius * scale_factor))
                        cv2.circle(frame_rgb, (scaled_x, scaled_y), radius=radius, color=gaze_color, thickness=-1)

                im.set_array(frame_rgb)
                return [im]

            anim = FuncAnimation(fig, update_frame_no_video, frames=total_frames,
                                 interval=1000/fps, blit=True, repeat=True)

        # ================= SAVE / DISPLAY =================
        result = None
        trial_idx_val = trial_idx

        # Build output filename
        anim_name = f"animation_{trial_idx_val}"
        if tmin is not None and tmax is not None:
            anim_name += f"_{tmin}_{tmax}"

        # Handle different output formats
        if output_format == "matplotlib":
            # Show in matplotlib GUI window (blocking)
            plt.show()
            # Cleanup video capture if used
            if video_path is not None:
                cap.release()
            return None

        elif output_format == "mp4":
            if folder_path:
                folder_path = Path(folder_path)
                folder_path.mkdir(parents=True, exist_ok=True)
                out_path = folder_path / f"{anim_name}.mp4"
                try:
                    anim.save(str(out_path), writer='ffmpeg', fps=fps)
                    print(f"Animation saved to: {out_path}")
                except Exception as e:
                    raise RuntimeError(
                        f"Failed to save MP4. Make sure ffmpeg is installed. Error: {e}"
                    )
            plt.close(fig)

        elif output_format == "gif":
            if folder_path:
                folder_path = Path(folder_path)
                folder_path.mkdir(parents=True, exist_ok=True)
                out_path = folder_path / f"{anim_name}.gif"
                try:
                    anim.save(str(out_path), writer='pillow', fps=fps)
                    print(f"Animation saved to: {out_path}")
                except Exception as e:
                    raise RuntimeError(
                        f"Failed to save GIF. Make sure pillow is installed. Error: {e}"
                    )
            plt.close(fig)

        else:  # html (default)
            if folder_path:
                folder_path = Path(folder_path)
                folder_path.mkdir(parents=True, exist_ok=True)
                out_path = folder_path / f"{anim_name}.html"
                with open(out_path, 'w') as f:
                    f.write(anim.to_jshtml())
                print(f"Animation saved to: {out_path}")

            if display:
                try:
                    from IPython.display import HTML
                    plt.close(fig)
                    result = HTML(anim.to_jshtml())
                except ImportError:
                    print("IPython not available. Use output_format='matplotlib' for GUI display.")
                    plt.close(fig)
            else:
                plt.close(fig)

        # Cleanup video capture if used
        if video_path is not None:
            cap.release()

        return result

plot_animation(samples, screen_height, screen_width, video_path=None, background_image_path=None, folder_path=None, tmin=None, tmax=None, seconds_to_show=None, scale_factor=0.5, gaze_radius=10, gaze_color=(255, 0, 0), fps=None, output_format='matplotlib', display=True)

Create an animated visualization of eye-tracking data.

When a video is provided, the animation syncs gaze samples with video frames. When no video is provided, gaze points are animated on a grey background or a provided background image, using the sample timestamps for timing.

Parameters:

Name	Type	Description	Default
samples	DataFrame	Polars DataFrame with gaze samples. Must contain 'tSample' and gaze position columns ('X', 'Y' or 'LX', 'LY', 'RX', 'RY').	required
screen_height	int	Stimulus resolution in pixels.	required
screen_width	int	Stimulus resolution in pixels.	required
video_path	str | Path | None	Path to a video file. If provided, gaze is overlaid on video frames.	None
background_image_path	str | Path | None	Path to a background image. Only used when video_path is None. If both are None, a grey background is used.	None
folder_path	str | Path | None	Directory where the animation will be saved. If None, nothing is saved. The file format depends on output_format.	None
tmin	int | None	Time window in ms. If both None, the whole trial is plotted.	None
tmax	int | None	Time window in ms. If both None, the whole trial is plotted.	None
seconds_to_show	float | None	Limit the animation to the first N seconds. If None, shows all available data.	None
scale_factor	float	Resolution scaling factor (1.0 = original, 0.5 = half resolution).	0.5
gaze_radius	int	Radius of the gaze point circle in pixels (before scaling).	10
gaze_color	tuple	RGB tuple for gaze point color.	(255, 0, 0)
fps	float | None	Frames per second for the animation. If None: - With video: uses the video's native FPS - Without video: defaults to 60 FPS	None
output_format	str	Output format for saved animations: - "html": Interactive HTML file (default, works in browsers) - "mp4": Video file (requires ffmpeg) - "gif": Animated GIF file (requires pillow) - "matplotlib": Show in matplotlib GUI window (blocking)	'matplotlib'
display	bool	If True and output_format is "html", returns an HTML object for notebooks. If output_format is "matplotlib", this is ignored (always shows window). If False, only saves to file (if folder_path is provided).	True

Returns:

Type	Description
HTML or None	Returns HTML animation if display=True and output_format="html", otherwise None.

Source code in pyxations/visualization/visualization.py

def plot_animation(
    self,
    samples: pl.DataFrame,
    screen_height: int,
    screen_width: int,
    video_path: str | Path | None = None,
    background_image_path: str | Path | None = None,
    folder_path: str | Path | None = None,
    tmin: int | None = None,
    tmax: int | None = None,
    seconds_to_show: float | None = None,
    scale_factor: float = 0.5,
    gaze_radius: int = 10,
    gaze_color: tuple = (255, 0, 0),
    fps: float | None = None,
    output_format: str = "matplotlib",
    display: bool = True,
):
    """
    Create an animated visualization of eye-tracking data.

    When a video is provided, the animation syncs gaze samples with video frames.
    When no video is provided, gaze points are animated on a grey background or
    a provided background image, using the sample timestamps for timing.

    Parameters
    ----------
    samples
        Polars DataFrame with gaze samples. Must contain 'tSample' and gaze
        position columns ('X', 'Y' or 'LX', 'LY', 'RX', 'RY').
    screen_height, screen_width
        Stimulus resolution in pixels.
    video_path
        Path to a video file. If provided, gaze is overlaid on video frames.
    background_image_path
        Path to a background image. Only used when video_path is None.
        If both are None, a grey background is used.
    folder_path
        Directory where the animation will be saved. If None, nothing is saved.
        The file format depends on `output_format`.
    tmin, tmax
        Time window in **ms**. If both None, the whole trial is plotted.
    seconds_to_show
        Limit the animation to the first N seconds. If None, shows all available data.
    scale_factor
        Resolution scaling factor (1.0 = original, 0.5 = half resolution).
    gaze_radius
        Radius of the gaze point circle in pixels (before scaling).
    gaze_color
        RGB tuple for gaze point color.
    fps
        Frames per second for the animation. If None:
        - With video: uses the video's native FPS
        - Without video: defaults to 60 FPS
    output_format
        Output format for saved animations:
        - "html": Interactive HTML file (default, works in browsers)
        - "mp4": Video file (requires ffmpeg)
        - "gif": Animated GIF file (requires pillow)
        - "matplotlib": Show in matplotlib GUI window (blocking)
    display
        If True and output_format is "html", returns an HTML object for notebooks.
        If output_format is "matplotlib", this is ignored (always shows window).
        If False, only saves to file (if folder_path is provided).

    Returns
    -------
    IPython.display.HTML or None
        Returns HTML animation if display=True and output_format="html", otherwise None.
    """
    try:
        import cv2
        from matplotlib.animation import FuncAnimation
        import matplotlib as mpl
        mpl.rcParams['animation.embed_limit'] = 100
    except ImportError as e:
        raise ImportError(
            f"Missing required dependency for animation: {e}. "
            "Please install cv2 (opencv-python)."
        )

    # Validate output_format
    valid_formats = ["html", "mp4", "gif", "matplotlib"]
    if output_format not in valid_formats:
        raise ValueError(f"output_format must be one of {valid_formats}, got '{output_format}'")

    # ---- Determine gaze columns ----
    if "X" in samples.columns and "Y" in samples.columns:
        x_col, y_col = "X", "Y"
    elif "LX" in samples.columns and "LY" in samples.columns:
        x_col, y_col = "LX", "LY"
    elif "RX" in samples.columns and "RY" in samples.columns:
        x_col, y_col = "RX", "RY"
    else:
        raise ValueError("Samples DataFrame must contain gaze columns (X, Y) or (LX, LY) or (RX, RY)")

    # ---- Time filter ----
    if tmin is not None and tmax is not None:
        samples = samples.filter(pl.col("tSample").is_between(tmin, tmax))

    if samples.is_empty():
        raise ValueError("No samples available after time filtering")

    # ---- Drop NaN gaze values ----
    samples = samples.filter(pl.col(x_col).is_not_null() & pl.col(y_col).is_not_null())

    # ---- Calculate scaled dimensions ----
    scaled_width = int(screen_width * scale_factor)
    scaled_height = int(screen_height * scale_factor)

    trial_idx = samples["trial_number"][0] if "trial_number" in samples.columns else 0

    # ================= WITH VIDEO =================
    if video_path is not None:
        video_path = Path(video_path)
        if not video_path.exists():
            raise FileNotFoundError(f"Video file not found: {video_path}")

        cap = cv2.VideoCapture(str(video_path))
        video_fps = cap.get(cv2.CAP_PROP_FPS)
        total_frames = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))
        video_width = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
        video_height = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))

        if fps is None:
            fps = video_fps

        # Calculate time to frame mapping
        t_start = samples["tSample"].min()
        t_end = samples["tSample"].max()
        trial_duration = t_end - t_start

        # Create frame-to-time mapping
        frame_edges = np.linspace(t_start, t_end, total_frames + 1)
        frame_times = ((frame_edges[:-1] + frame_edges[1:]) / 2).astype(int)

        # Build a lookup: frame_index -> list of gaze points
        samples_np = samples.select([x_col, y_col, "tSample"]).to_numpy()
        gaze_by_frame = {i: [] for i in range(total_frames)}

        for x, y, t in samples_np:
            # Find the closest frame
            frame_idx = np.searchsorted(frame_times, t, side='right') - 1
            frame_idx = max(0, min(frame_idx, total_frames - 1))
            gaze_by_frame[frame_idx].append((x, y))

        # Limit frames if seconds_to_show is set
        frames_to_show = total_frames
        if seconds_to_show is not None:
            frames_to_show = min(int(fps * seconds_to_show), total_frames)

        # Reset video
        cap.set(cv2.CAP_PROP_POS_FRAMES, 0)

        # Create figure
        fig, ax = plt.subplots(figsize=(10 * scale_factor, 6 * scale_factor))
        ax.axis('off')

        # Initialize with first frame
        ret, frame = cap.read()
        if not ret:
            cap.release()
            raise RuntimeError("Could not read first frame from video")

        frame_resized = cv2.resize(frame, (scaled_width, scaled_height), interpolation=cv2.INTER_AREA)
        frame_rgb = cv2.cvtColor(frame_resized, cv2.COLOR_BGR2RGB)
        im = ax.imshow(frame_rgb)

        def update_frame_video(frame_idx):
            ret, frame = cap.read()
            if not ret:
                return [im]

            frame_resized = cv2.resize(frame, (scaled_width, scaled_height), interpolation=cv2.INTER_AREA)
            frame_rgb = cv2.cvtColor(frame_resized, cv2.COLOR_BGR2RGB)

            # Draw gaze points for this frame
            for gx, gy in gaze_by_frame.get(frame_idx, []):
                scaled_x = int(gx * scale_factor)
                scaled_y = int(gy * scale_factor)
                if 0 <= scaled_x < scaled_width and 0 <= scaled_y < scaled_height:
                    radius = max(3, int(gaze_radius * scale_factor))
                    cv2.circle(frame_rgb, (scaled_x, scaled_y), radius=radius, color=gaze_color, thickness=-1)

            im.set_array(frame_rgb)
            return [im]

        anim = FuncAnimation(fig, update_frame_video, frames=frames_to_show,
                             interval=1000/fps, blit=True, repeat=True)

    # ================= WITHOUT VIDEO =================
    else:
        if fps is None:
            fps = 60  # Default FPS for sample-based animation

        # Prepare background
        if background_image_path is not None:
            bg_path = Path(background_image_path)
            if not bg_path.exists():
                raise FileNotFoundError(f"Background image not found: {bg_path}")
            bg_img = mpimg.imread(str(bg_path))
            if bg_img.dtype == np.float64:
                bg_img = (bg_img * 255).clip(0, 255).astype(np.uint8)
            # Resize background to match screen dimensions then scale
            bg_img = cv2.resize(bg_img, (scaled_width, scaled_height), interpolation=cv2.INTER_AREA)
        else:
            # Grey background
            bg_img = np.ones((scaled_height, scaled_width, 3), dtype=np.uint8) * 128

        # Get time range
        t_start = samples["tSample"].min()
        t_end = samples["tSample"].max()
        trial_duration = t_end - t_start

        # Limit duration if seconds_to_show is set
        if seconds_to_show is not None:
            t_end = min(t_end, t_start + int(seconds_to_show * 1000))
            samples = samples.filter(pl.col("tSample") <= t_end)
            trial_duration = t_end - t_start

        # Calculate total frames based on duration and fps
        total_frames = int((trial_duration / 1000) * fps)
        if total_frames < 1:
            total_frames = 1

        # Create time bins for each animation frame
        frame_times = np.linspace(t_start, t_end, total_frames + 1)

        # Build gaze lookup by frame
        samples_np = samples.select([x_col, y_col, "tSample"]).to_numpy()
        gaze_by_frame = {i: [] for i in range(total_frames)}

        for x, y, t in samples_np:
            frame_idx = np.searchsorted(frame_times, t, side='right') - 1
            frame_idx = max(0, min(frame_idx, total_frames - 1))
            gaze_by_frame[frame_idx].append((x, y))

        # Create figure
        fig, ax = plt.subplots(figsize=(10 * scale_factor, 6 * scale_factor))
        ax.axis('off')

        # Initialize with background
        im = ax.imshow(bg_img.copy())

        def update_frame_no_video(frame_idx):
            # Start with fresh background copy
            frame_rgb = bg_img.copy()

            # Draw gaze points for this frame
            for gx, gy in gaze_by_frame.get(frame_idx, []):
                scaled_x = int(gx * scale_factor)
                scaled_y = int(gy * scale_factor)
                if 0 <= scaled_x < scaled_width and 0 <= scaled_y < scaled_height:
                    radius = max(3, int(gaze_radius * scale_factor))
                    cv2.circle(frame_rgb, (scaled_x, scaled_y), radius=radius, color=gaze_color, thickness=-1)

            im.set_array(frame_rgb)
            return [im]

        anim = FuncAnimation(fig, update_frame_no_video, frames=total_frames,
                             interval=1000/fps, blit=True, repeat=True)

    # ================= SAVE / DISPLAY =================
    result = None
    trial_idx_val = trial_idx

    # Build output filename
    anim_name = f"animation_{trial_idx_val}"
    if tmin is not None and tmax is not None:
        anim_name += f"_{tmin}_{tmax}"

    # Handle different output formats
    if output_format == "matplotlib":
        # Show in matplotlib GUI window (blocking)
        plt.show()
        # Cleanup video capture if used
        if video_path is not None:
            cap.release()
        return None

    elif output_format == "mp4":
        if folder_path:
            folder_path = Path(folder_path)
            folder_path.mkdir(parents=True, exist_ok=True)
            out_path = folder_path / f"{anim_name}.mp4"
            try:
                anim.save(str(out_path), writer='ffmpeg', fps=fps)
                print(f"Animation saved to: {out_path}")
            except Exception as e:
                raise RuntimeError(
                    f"Failed to save MP4. Make sure ffmpeg is installed. Error: {e}"
                )
        plt.close(fig)

    elif output_format == "gif":
        if folder_path:
            folder_path = Path(folder_path)
            folder_path.mkdir(parents=True, exist_ok=True)
            out_path = folder_path / f"{anim_name}.gif"
            try:
                anim.save(str(out_path), writer='pillow', fps=fps)
                print(f"Animation saved to: {out_path}")
            except Exception as e:
                raise RuntimeError(
                    f"Failed to save GIF. Make sure pillow is installed. Error: {e}"
                )
        plt.close(fig)

    else:  # html (default)
        if folder_path:
            folder_path = Path(folder_path)
            folder_path.mkdir(parents=True, exist_ok=True)
            out_path = folder_path / f"{anim_name}.html"
            with open(out_path, 'w') as f:
                f.write(anim.to_jshtml())
            print(f"Animation saved to: {out_path}")

        if display:
            try:
                from IPython.display import HTML
                plt.close(fig)
                result = HTML(anim.to_jshtml())
            except ImportError:
                print("IPython not available. Use output_format='matplotlib' for GUI display.")
                plt.close(fig)
        else:
            plt.close(fig)

    # Cleanup video capture if used
    if video_path is not None:
        cap.release()

    return result

plot_multipanel(fixations, saccades, display=True)

Create a 2×2 multi‑panel diagnostic plot for every non‑empty phase label and save it as PNG in //plots/.

Source code in pyxations/visualization/visualization.py

def plot_multipanel(
        self,
        fixations: pl.DataFrame,
        saccades: pl.DataFrame,
        display: bool = True
    ) -> None:
    """
    Create a 2×2 multi‑panel diagnostic plot for every non‑empty
    phase label and save it as PNG in
    <derivatives_folder_path>/<events_detection_folder>/plots/.
    """
    # ── paths & matplotlib style ────────────────────────────────
    folder_path: Path = (
        self.derivatives_folder_path
        / self.events_detection_folder
        / "plots"
    )
    folder_path.mkdir(parents=True, exist_ok=True)
    plt.rcParams.update({"font.size": 12})

    # ── drop practice / invalid trials ─────────────────────────
    fixations = fixations.filter(pl.col("trial_number") != -1)
    saccades  = saccades.filter(pl.col("trial_number") != -1)

    # ── collect valid phase labels (skip empty string) ─────────
    phases = (
        fixations
        .select(pl.col("phase").filter(pl.col("phase") != ""))
        .unique()           # unique values in this Series
        .to_series()
        .to_list()          # plain Python list of strings
    )

    # ── one figure per phase ───────────────────────────────────
    for phase in phases:
        fix_phase   = fixations.filter(pl.col("phase") == phase)
        sacc_phase  = saccades.filter(pl.col("phase") == phase)

        fig, axs = plt.subplots(2, 2, figsize=(12, 7))

        self.fix_duration(fix_phase , axs=axs[0, 0])
        self.sacc_main_sequence(sacc_phase, axs=axs[1, 1])
        self.sacc_direction(sacc_phase, axs=axs[1, 0], figs=fig)
        self.sacc_amplitude(sacc_phase, axs=axs[0, 1])

        fig.tight_layout()
        plt.savefig(folder_path / f"multipanel_{phase}.png")
        if display:
            plt.show()
        plt.close()

scanpath(fixations, screen_height, screen_width, folder_path=None, tmin=None, tmax=None, saccades=None, samples=None, phase_data=None, display=True)

Fast scan‑path visualiser.

• Vectorised: no per‑row Python loops
• Single pass phase grouping
• Uses BrokenBarHCollection for fixation spans
• Optional asynchronous PNG write via ThreadPoolExecutor (drop‑in‑ready, see comment)

Parameters:

Name	Type	Description	Default
fixations	DataFrame	Polars DataFrame with at least tStart, duration, xAvg, yAvg, phase.	required
screen_height	int	Stimulus resolution in pixels.	required
screen_width	int	Stimulus resolution in pixels.	required
folder_path	str | Path | None	Directory where 1 PNG per phase will be stored. If None, nothing is saved.	None
tmin	int | None	Time window in ms. If both None, the whole trial is plotted.	None
tmax	int | None	Time window in ms. If both None, the whole trial is plotted.	None
saccades	DataFrame | None	Polars DataFrame with tStart, phase, … (optional).	None
samples	DataFrame | None	Polars DataFrame with gaze traces (tSample, LX, LY, RX, RY or X, Y) (optional).	None
phase_data	dict[str, dict] | None	Per‑phase extras:: { "search": { "img_paths": [...], "img_plot_coords": [(x1,y1,x2,y2), ...], "bbox": (x1,y1,x2,y2), }, ... }	None
display	bool	If False the figure canvas is never shown (faster for batch jobs).	True

Source code in pyxations/visualization/visualization.py

def scanpath(
    self,
    fixations: pl.DataFrame,
    screen_height: int,
    screen_width: int,
    folder_path: str | Path | None = None,
    tmin: int | None = None,
    tmax: int | None = None,
    saccades: pl.DataFrame | None = None,
    samples: pl.DataFrame | None = None,
    phase_data: dict[str, dict] | None = None,
    display: bool = True,
):
    """
    Fast scan‑path visualiser.

    • **Vectorised**: no per‑row Python loops  
    • **Single pass** phase grouping  
    • Uses `BrokenBarHCollection` for fixation spans  
    • Optional asynchronous PNG write via ThreadPoolExecutor (drop‑in‑ready, see comment)

    Parameters
    ----------
    fixations
        Polars DataFrame with at least `tStart`, `duration`, `xAvg`, `yAvg`, `phase`.
    screen_height, screen_width
        Stimulus resolution in pixels.
    folder_path
        Directory where 1 PNG per phase will be stored.  If *None*, nothing is saved.
    tmin, tmax
        Time window in **ms**.  If both `None`, the whole trial is plotted.
    saccades
        Polars DataFrame with `tStart`, `phase`, …  (optional).
    samples
        Polars DataFrame with gaze traces (`tSample`, `LX`, `LY`, `RX`, `RY` or
        `X`, `Y`) (optional).
    phase_data
        Per‑phase extras::

            {
                "search": {
                    "img_paths": [...],
                    "img_plot_coords": [(x1,y1,x2,y2), ...],
                    "bbox": (x1,y1,x2,y2),
                },
                ...
            }

    display
        If *False* the figure canvas is never shown (faster for batch jobs).
    """


    # ------------- small helpers ------------------------------------------------
    def _make_axes(plot_samples: bool):
        if plot_samples:
            fig, (ax_main, ax_gaze) = plt.subplots(
                2, 1, height_ratios=(4, 1), figsize=(10, 6), sharex=False
            )
        else:
            fig, ax_main = plt.subplots(figsize=(10, 6))
            ax_gaze = None
        ax_main.set_xlim(0, screen_width)
        ax_main.set_ylim(screen_height, 0)
        return fig, ax_main, ax_gaze

    def _maybe_cache_img(path):
        """Load image from disk with a small LRU cache."""

        # Cache hit: move to the end (most recently used)
        if path in _img_cache:
            img = _img_cache.pop(path)
            _img_cache[path] = img
            return img

        # Cache miss: load image
        img = mpimg.imread(path)

        # Optional: reduce memory if image is float64 in [0, 1]
        if isinstance(img, np.ndarray) and img.dtype == np.float64:
            img = (img * 255).clip(0, 255).astype(np.uint8)

        # Insert into cache
        _img_cache[path] = img

        # If cache too big, drop least recently used item
        if len(_img_cache) > _MAX_CACHE_ITEMS:
            _img_cache.popitem(last=False)  # pops the oldest inserted item

        return img

    # ---------------------------------------------------------------------------
    plot_saccades = saccades is not None
    plot_samples = samples is not None
    _img_cache = OrderedDict()
    _MAX_CACHE_ITEMS = 8  # or 5, 10, etc. Tune as you like.

    trial_idx = fixations["trial_number"][0]

    # ---- time filter ----------------------------------------------------------
    if tmin is not None and tmax is not None:
        fixations = fixations.filter(pl.col("tStart").is_between(tmin, tmax))
        if plot_saccades:
            saccades = saccades.filter(pl.col("tStart").is_between(tmin, tmax))
        if plot_samples:
            samples = samples.filter(pl.col("tSample").is_between(tmin, tmax))

    # remove empty phase markings
    fixations = fixations.filter(pl.col("phase") != "")
    if plot_saccades:
        saccades = saccades.filter(pl.col("phase") != "")
    if plot_samples:
        samples = samples.filter(pl.col("phase") != "")

    # ---- split once by phase --------------------------------------------------
    fix_by_phase = fixations.partition_by("phase", as_dict=True)
    sac_by_phase = (
        saccades.partition_by("phase", as_dict=True) if plot_saccades else {}
    )
    samp_by_phase = (
        samples.partition_by("phase", as_dict=True) if plot_samples else {}
    )

    # colour map shared across phases
    cmap = plt.cm.rainbow

    # ---- build & draw ---------------------------------------------------------
    # optional async saver (uncomment if you save hundreds of files)
    from concurrent.futures import ThreadPoolExecutor
    saver = ThreadPoolExecutor(max_workers=4) if folder_path else None

    if not display:
        plt.ioff()

    for phase, phase_fix in fix_by_phase.items():
        if phase_fix.is_empty():
            continue

        # ---------- vectors (zero‑copy) -----------------
        fx, fy, fdur = phase_fix.select(["xAvg", "yAvg", "duration"]).to_numpy().T
        n_fix = fx.size
        fix_idx = np.arange(1, n_fix + 1)

        norm = mplcolors.BoundaryNorm(np.arange(1, n_fix + 2), cmap.N)

        # saccades
        sac_t = (
            sac_by_phase[phase]["tStart"].to_numpy()
            if plot_saccades and phase in sac_by_phase
            else np.empty(0)
        )

        # samples
        if plot_samples and phase in samp_by_phase and samp_by_phase[phase].height:
            samp_phase = samp_by_phase[phase]
            t0 = samp_phase["tSample"][0]
            ts = (samp_phase["tSample"].to_numpy() - t0) 
            get = samp_phase.get_column
            lx = get("LX").to_numpy() if "LX" in samp_phase.columns else None
            ly = get("LY").to_numpy() if "LY" in samp_phase.columns else None
            rx = get("RX").to_numpy() if "RX" in samp_phase.columns else None
            ry = get("RY").to_numpy() if "RY" in samp_phase.columns else None
            gx = get("X").to_numpy() if "X" in samp_phase.columns else None
            gy = get("Y").to_numpy() if "Y" in samp_phase.columns else None
        else:
            t0 = None

        # ---------- figure -----------------------------
        fig, ax_main, ax_gaze = _make_axes(plot_samples and t0 is not None)
        # scatter fixations
        sc = ax_main.scatter(
            fx,
            fy,
            c=fix_idx,
            s=fdur,
            cmap=cmap,
            norm=norm,
            alpha=0.5,
            zorder=2,
        )
        fig.colorbar(
            sc,
            ax=ax_main,
            ticks=[1, n_fix // 2 if n_fix > 2 else n_fix, n_fix],
            fraction=0.046,
            pad=0.04,
        ).set_label("# of fixation")

        # ---------- stimulus imagery / bbox ------------
        if phase_data and phase[0] in phase_data:
            pdict = phase_data[phase[0]]
            coords = pdict.get("img_plot_coords") or []
            bbox = pdict.get('bbox',None) 
            for img_path, box in zip(pdict.get("img_paths", []), coords):

                ax_main.imshow(_maybe_cache_img(img_path), extent=[box[0], box[2], box[3], box[1]], zorder=0)
            if bbox is not None:
                x1, y1, x2, y2 = bbox
                ax_main.plot([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], color='red', linewidth=1.5, zorder=3)

        # ---------- gaze traces ------------------------
        if ax_gaze is not None:
            if lx is not None:
                ax_main.plot(lx, ly, "--", color="C0", zorder=1)
                ax_gaze.plot(ts, lx, label="Left X")
                ax_gaze.plot(ts, ly, label="Left Y")
            if rx is not None:
                ax_main.plot(rx, ry, "--", color="k", zorder=1)
                ax_gaze.plot(ts, rx, label="Right X")
                ax_gaze.plot(ts, ry, label="Right Y")
            if gx is not None:
                ax_main.plot(gx, gy, "--", color="k", zorder=1, alpha=0.6)
                ax_gaze.plot(ts, gx, label="X")
                ax_gaze.plot(ts, gy, label="Y")

            # fixation spans
            bars   = np.c_[phase_fix['tStart'].to_numpy() - t0,
                        phase_fix['duration'].to_numpy()]
            height = ax_gaze.get_ylim()[1] - ax_gaze.get_ylim()[0]
            colors = cmap(norm(fix_idx))

            # Draw all bars in one call; no BrokenBarHCollection import needed
            ax_gaze.broken_barh(bars, (0, height), facecolors=colors, alpha=0.4)
            # saccades
            if sac_t.size:
                ymin, ymax = ax_gaze.get_ylim()
                ax_gaze.vlines(
                    sac_t - t0,
                    ymin,
                    ymax,
                    colors="red",
                    linestyles="--",
                    linewidth=0.8,
                )

            # tidy gaze axis
            h, l = ax_gaze.get_legend_handles_labels()
            by_label = {lab: hdl for hdl, lab in zip(h, l)}
            ax_gaze.legend(
                by_label.values(),
                by_label.keys(),
                loc="center left",
                bbox_to_anchor=(1, 0.5),
            )
            ax_gaze.set_ylabel("Gaze")
            ax_gaze.set_xlabel("Time [s]")

        fig.tight_layout()

        # ---------- save / show ------------------------
        if folder_path:
            scan_name = f"scanpath_{trial_idx}"
            if tmin is not None and tmax is not None:
                scan_name += f"_{tmin}_{tmax}"
            out = Path(folder_path) / f"{scan_name}_{phase[0]}.png"
            fig.savefig(out, dpi=150)
            if saver:  saver.submit(fig.savefig, out, dpi=150)

        if display:
            plt.show()
        plt.close(fig)

    if not display:
        plt.ion()

samples

Created on Nov 22, 2024

@author: placiana

formats

Vendor-specific input readers. Selected via the dataset_format argument of compute_derivatives_for_dataset.

Created on 2 dic 2024

@author: placiana

export

Writers for persisting derivatives in different on-disk formats.

utils

Helpers for log alignment and miscellaneous utilities.