BIDS_process_physio_ses_2.py

"""
BIDS_process_physio_ses_2.py

Description:
This script processes physiological data, specifically  [ECG, respiration, EDA, PPG], acquired through 
Biopac Acqknowledge during Task fMRI sessions, into BIDS-compliant files. It includes functions for loading data, 
segmenting based on run and event types, writing metadata, and generating quality control plots. The script is designed 
to work with data in .mat format format and expects a BIDS specific directory structure for input and output.

Usage:
Run the script from the command line with the physiological data directory and BIDS dataset directory as arguments.
e.g., 
python BIDS_process_physio_ses_2.py /path/to/physio_data /path/to/BIDS_dataset
Output files will be saved in the specified BIDS dataset directory.

Author: PAMcConnell
Created on: 20231111
Last Modified: 20231111

License: MIT License

Dependencies:
- Python 3.12
- numpy 
- pandas 
- matplotliby)
- scipy 

Environment Setup:
Use the provided environment.yml (datalad.yml) file with Conda for setting up the required environment. 

Change Log:
- 20231111: Initial version

Error Handling and Logging:
The script includes comprehensive error handling and logging. Errors during execution are logged, providing clarity for debugging.
Logs are saved in dataset_root/docs/logs.
"""

import os                                             # Used for interacting with the operating system, like file path operations.
import re                                             # Regular expressions, useful for text matching and manipulation.
import logging                                        # Logging library, for tracking events that happen when running the software.
import argparse                                       # Parser for command-line options, arguments, and sub-commands.
import scipy.io as sio                                # SciPy module for reading and writing MATLAB files.
import numpy as np                                    # NumPy library for numerical operations on arrays and matrices.
import pandas as pd                                   # Pandas library for data manipulation and analysis.
import matplotlib.pyplot as plt                       # Matplotlib's pyplot, a plotting library.
import json                                           # Library for working with JSON data.
import glob                                           # Used for Unix style pathname pattern expansion.
from collections import OrderedDict                   # Dictionary subclass that remembers the insertion order of keys.
import sys                                            # System-specific parameters and functions.
from matplotlib.lines import Line2D                   # Import Line2D from matplotlib to create custom line objects.
from matplotlib.backends.backend_pdf import PdfPages    # for creating multipage PDFs with matplotlib plots.
import shutil                                           # for copying files and directories.

# Sets up archival logging for the script, directing log output to both a file and the console.
def setup_logging(subject_id, session_id, bids_root_dir):
    """

    The function configures logging to capture informational, warning, and error messages. It creates a unique log file for each 
    subject-session combination, stored in a 'logs' directory within the 'doc' folder adjacent to the BIDS root directory.

  
    Parameters:
    - subject_id (str): The identifier for the subject.
    - session_id (str): The identifier for the session.
    - bids_root_dir (str): The root directory of the BIDS dataset.

    Returns:
    - log_file_path (str): The path to the log file.

    This function sets up a logging system that writes logs to both a file and the console. 
    The log file is named based on the subject ID, session ID, and the script name. 
    It's stored in a 'logs' directory within the 'doc' folder by subject ID, which is located at the same 
    level as the BIDS root directory.

    The logging level is set to INFO, meaning it captures all informational, warning, and error messages.

    Usage Example:
    setup_logging('sub-01', 'ses-1', '/path/to/bids_root_dir')
    """

    try: 
        # Extract the base name of the script without the .py extension.
        script_name = os.path.basename(__file__).replace('.py', '')

        # Construct the log directory path within 'doc/logs'
        log_dir = os.path.join(os.path.dirname(bids_root_dir), 'doc', 'logs', script_name, subject_id)

        # Create the log directory if it doesn't exist.
        if not os.path.exists(log_dir):
            os.makedirs(log_dir)

        # Construct the log file name using subject ID, session ID, and script name.
        log_file_name = f"{subject_id}_{session_id}_{script_name}.log"
        log_file_path = os.path.join(log_dir, log_file_name)

        # Configure file logging.
        logging.basicConfig(
            level=logging.INFO,
            format='%(asctime)s - %(levelname)s - %(message)s',
            filename=log_file_path,
            filemode='w' # 'w' mode overwrites existing log file
        )

        # If you also want to log to console.
        console_handler = logging.StreamHandler(sys.stdout)
        console_handler.setLevel(logging.INFO)
        formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
        console_handler.setFormatter(formatter)
        logging.getLogger().addHandler(console_handler)

        logging.info(f"Logging setup complete. Log file: {log_file_path}")

        return log_file_path
    
    except Exception as e:
            print(f"Error setting up logging: {e}")
            sys.exit(1) # Exiting the script due to logging setup failure.

# Extract the subject and session IDs from the provided physio root directory path.
def extract_subject_session(physio_root_dir):
    """
    Parameters:
    - physio_root_dir (str): The directory path that includes subject and session information. 
                             This path should follow the BIDS convention, containing 'sub-' and 'ses-' prefixes.

    Returns:
    - subject_id (str): The extracted subject ID.
    - session_id (str): The extracted session ID.

    Raises:
    - ValueError: If the subject_id and session_id cannot be extracted from the path.

    This function assumes that the directory path follows the Brain Imaging Data Structure (BIDS) naming convention. 
    It uses regular expressions to find and extract the subject and session IDs from the path.

    Usage Example:
    subject_id, session_id = extract_subject_session('/path/to/data/sub-01/ses-1/physio')

    Note: This function will raise an error if it cannot find a pattern matching the BIDS convention in the path.
    """

    # Normalize the path to remove any trailing slashes for consistency.
    physio_root_dir = os.path.normpath(physio_root_dir)

    # The pattern looks for 'sub-' followed by any characters until a slash, and similar for 'ses-'.
    match = re.search(r'(sub-[^/]+)/(ses-[^/]+)', physio_root_dir)
    if not match:
        raise ValueError("Unable to extract subject_id and session_id from path: %s", physio_root_dir)
    
    subject_id, session_id = match.groups()

    return subject_id, session_id

# Load a MATLAB (.mat) file containing physiological data and extracts labels, data, and units for physiological data.
def load_mat_file(mat_file_path):
    """
    Loads a MATLAB (.mat) file and extracts physiological data labels, data, and units.

    Parameters:
    - mat_file_path (str): Path to the .mat file.

    Returns:
    - labels (ndarray): Array of data channel names.
    - data (ndarray): Array containing physiological data.
    - units (ndarray): Array of units corresponding to each data channel.

    Raises:
    - FileNotFoundError: If the .mat file is not found at the specified path.
    - KeyError: If the .mat file lacks required keys ('labels', 'data', 'units').
    - Exception: For any other issues encountered during loading.

    The function verifies the existence of the .mat file, loads its contents, and checks for the presence 
    of required keys. It provides detailed logging for each step and potential errors for troubleshooting.

    Usage Example:
    labels, data, units = load_mat_file('/path/to/physio_data.mat')

    Note:
    - The .mat file must contain 'labels', 'data', and 'units' keys.
    - Compatibility with MATLAB file formats should be verified for different versions.
    """
    
    # Verify if the specified .mat file exists.
    if not os.path.isfile(mat_file_path):
        error_msg = f"MAT file does not exist at {mat_file_path}"
        logging.error(error_msg)
        raise FileNotFoundError(error_msg)
    
    try:
        # Attempt to load the .mat file.
        logging.info(f"Loading MAT file from: {mat_file_path}")
        mat_contents = sio.loadmat(mat_file_path)
        
        # Verify that required keys are in the loaded .mat file.
        required_keys = ['labels', 'data', 'units']
        if not all(key in mat_contents for key in required_keys):
            error_msg = f"MAT file at {mat_file_path} is missing required keys: {required_keys}"
            logging.error(error_msg)
            raise KeyError(error_msg)
        
        # Extract labels, data, and units.
        labels = mat_contents['labels'].flatten()  # Flatten in case it's a 2D array.
        data = mat_contents['data']
        units = mat_contents['units'].flatten()  # Flatten in case it's a 2D array
        logging.info(f"Successfully loaded MAT file from: {mat_file_path}")
              
    except Exception as e:
        # Log the exception and re-raise to handle it upstream.
        logging.error("Failed to load MAT file from %s: %s", mat_file_path, e)
        raise
    
    return labels, data, units

#   Renames physiological data channels according to BIDS (Brain Imaging Data Structure) conventions.
def rename_channels(labels):
    """
    Parameters:
    - labels (array): Original names of the physiological data channels.

    Returns:
    - bids_labels_dictionary (dict): Mapping from original labels to BIDS-compliant labels.
    - bids_labels_list (list): A list of the renamed, BIDS-compliant labels.

    The function iterates through the provided original labels and renames them based on predefined 
    mapping to meet BIDS standards. Channels not recognized are omitted with a warning. This function
    is essential for ensuring that physiological data aligns with the BIDS format for further processing.

    Usage Example:
    bids_labels_dict, bids_labels_list = rename_channels(['ECG', 'RSP', 'EDA', 'PPG', 'Digital input'])

    Note: The function expects a specific set of channel names. Make sure to update the mapping 
    if different channels or naming conventions are used.
    """

    logging.info("Renaming channels according to BIDS conventions")

    # Define the mapping from original labels to BIDS labels.
    original_label_mapping = {
        'ECG': 'cardiac',
        'RSP': 'respiratory',
        'EDA': 'eda',
        'Trigger': 'trigger',
        'PPG': 'ppg',  # Only if exists
        # Add other mappings as required
    }

    # Initialize an empty dictionary and list to store the renamed labels.
    bids_labels_dictionary = {}
    bids_labels_list = []

    # Initialize dictionary and list for storing BIDS-compliant labels.
    for label in labels:
        # Skip any labels for digital inputs
        if 'Digital input' in label:
            continue
        
        # Check and rename the label if it matches one of the keys in original_label_mapping.
        for original, bids in original_label_mapping.items():
            if original in label:
                bids_labels_dictionary[label] = bids
                bids_labels_list.append(bids)
                break
        else:
            logging.warning("Label '%s' does not match any BIDS convention and will be omitted.", label)
    
    logging.info("BIDS labels list after renaming: %s", bids_labels_list)

    # Debug log to print the renamed labels in the list.
    return bids_labels_dictionary, bids_labels_list
   
#  Extracts metadata from a JSON file and identifies the associated run.
def extract_metadata_from_json(json_file_path, processed_jsons):
    """
    Extracts specific metadata from a JSON file and returns the run identifier and metadata.

    Parameters:
    - json_file_path (str): Path to the JSON file containing metadata for a specific fMRI run.
    - processed_jsons (set): A set that tracks already processed JSON files to avoid duplication.

    Returns:
    - tuple: (run_id, run_metadata), where:
        - run_id (str): The identifier of the fMRI run, extracted from the file name.
        - run_metadata (dict): A dictionary containing extracted metadata.

    Raises:
    - FileNotFoundError: If the specified JSON file does not exist.
    - ValueError: If the run_id cannot be determined or required metadata fields are missing.
    - json.JSONDecodeError: If the JSON file contains invalid JSON.

    This function is critical for parsing and organizing metadata necessary for subsequent data processing steps.
    It verifies the existence of essential fields and logs detailed information for debugging and audit purposes.
    
    Usage Example:
    run_id, metadata = extract_metadata_from_json('/path/to/run-01_bold.json', processed_jsons_set)
    """
   
    # Log the start of metadata extraction.
    logging.info(f"Extracting metadata from JSON file: {json_file_path}")

    # Avoid reprocessing the same file.
    if json_file_path in processed_jsons:
        logging.info("JSON file %s has already been processed.", json_file_path)
        return None, None  # No new metadata to return.
    
    # Check if the JSON file exists.
    if not os.path.isfile(json_file_path):
        logging.error("JSON file does not exist at %s", json_file_path)
        raise FileNotFoundError(f"No JSON file found at the specified path: {json_file_path}")
    
    # Load the JSON file content.
    try:
        with open(json_file_path, 'r') as file:
            metadata = json.load(file)
    except json.JSONDecodeError as e:
        logging.error(f"Error decoding JSON: {e}")
        raise
    
    try:
        # Extract run_id from the file name.
        run_id_match = re.search(r'run-\d+', json_file_path)
        if not run_id_match:
            raise ValueError(f"Run identifier not found in JSON file name: {json_file_path}")
        run_id = run_id_match.group()

        # Verify essential metadata fields are present.
        required_fields = ['TaskName', 'RepetitionTime', 'NumVolumes', 'SeriesDescription']
        run_metadata = {field: metadata.get(field) for field in required_fields}
        if not all(run_metadata.values()):
            missing_fields = [key for key, value in run_metadata.items() if value is None]
            raise ValueError(f"JSON file {json_file_path} is missing required fields: {missing_fields}")

        # Add the JSON file to the set of processed files.
        processed_jsons.add(json_file_path)
        logging.info(f"Successfully extracted metadata for {run_id}: {run_metadata}")

        # Return the run ID and extracted metadata.
        return run_id, run_metadata
    
    except Exception as e:
        logging.error(f"Unexpected error occurred while extracting metadata from {json_file_path}: {e}")
        raise

# Extracts the starting points of triggers in MRI data based on a specified threshold and minimum consecutive points.
def extract_trigger_points(data, threshold, min_consecutive):
    """
    Parameters:
    - data (numpy.ndarray): The data array from the MRI trigger channel.
    - threshold (float): The threshold value above which a data point is considered part of a trigger.
    - min_consecutive (int): The minimum number of consecutive data points above the threshold required to confirm a trigger start.

    Returns:
    - list: A list of indices in the data array where valid trigger starts are detected.

    This function is integral to preprocessing fMRI data, where accurate identification of trigger points is essential for aligning physiological data with MRI volumes. It uses a simple thresholding approach combined with a criterion for a minimum number of consecutive points to reduce false positives.

    The function first converts the data to a binary sequence based on the threshold, then identifies changes 
    from 0 to 1 as potential trigger starts. It further checks these starts to ensure they have the specified 
    minimum number of consecutive points above the threshold.

    Usage Example:
    trigger_starts = extract_trigger_points(mri_trigger_data, threshold=5, min_consecutive=5)

    Note:
    - This function assumes the trigger channel data is a 1D numpy array.
    - The choice of threshold and min_consecutive parameters may vary based on the specifics of the MRI acquisition and should be validated for each dataset.
    """
    logging.info("Extracting trigger points with threshold: %s and minimum consecutive points: %s", threshold, min_consecutive)

    # Convert the data points to binary values based on the threshold.
    triggers = (data > threshold).astype(int)
    logging.info(f"Trigger data converted to binary based on threshold {threshold}")

    # Calculate the difference to identify rising edges of the trigger signal.
    diff_triggers = np.diff(triggers, prepend=0)
    potential_trigger_starts = np.where(diff_triggers == 1)[0]

    # Validate trigger starts based on the minimum consecutive points criterion.
    valid_trigger_starts = []
    for start in potential_trigger_starts:
        # Check if the start has the required minimum number of consecutive points above the threshold
        if np.sum(triggers[start:start+min_consecutive]) == min_consecutive:
            valid_trigger_starts.append(start)

    logging.info("Identified %d valid trigger starts", len(valid_trigger_starts))
    return valid_trigger_starts

# Finds the index of the next MRI trigger start after a given index in the sequence of trigger starts.
def find_trigger_start(trigger_starts, current_index):
    """
    Parameters:
    - trigger_starts (numpy.ndarray): An array of indices where MRI trigger starts are detected.
    - current_index (int): The current index in the physiological data, used as a reference to find the next trigger start.

    Returns:
    - int or None: The index of the next trigger start if found; otherwise, None.

    This function is essential for processing fMRI physiological data, where identifying the start of each scan run is crucial for accurate data analysis. It works by searching for the nearest index in the trigger_starts array that comes after the current_index.

    Usage Example:
    next_trigger = find_trigger_start(trigger_starts, current_index)

    Note:
    - Ensure trigger_starts is a sorted numpy array.
    - If no trigger start is found after the current_index, the function returns None, indicating the end of the sequence or missing data.
    """
    
    logging.info(f"Finding trigger start after index {current_index}.")

    # Finding the appropriate index in the trigger_starts array.
    next_trigger_index = np.searchsorted(trigger_starts, current_index, side='right')
    
    # Verify if the next trigger index is within the bounds of the trigger_starts array.
    if next_trigger_index < len(trigger_starts):
        next_trigger_start = trigger_starts[next_trigger_index]
        logging.info(f"Next trigger start found at index {next_trigger_start}.")
        return next_trigger_start
    else:
        logging.error(f"No trigger start found after index {current_index}.")
        return None  # Indicate that there are no more triggers to process.

# Segments physiological data into individual fMRI runs based on metadata and trigger starts.
def find_runs(data, all_runs_metadata, trigger_starts, sampling_rate, invalid_runs):
    """
    Parameters:
    - data (numpy.ndarray): The complete physiological data set.
    - all_runs_metadata (dict): Metadata for each run, keyed by run identifiers.
    - trigger_starts (list): Indices marking the start of each potential fMRI run.
    - sampling_rate (int): The frequency at which data points were sampled.

    Returns:
    - list of dicts: Each dictionary contains the start and end indices of a run and its metadata.

    This function is critical for neuroimaging research where analyses are often conducted on data aligned with individual fMRI runs. It ensures that each run is accurately captured based on the timing information provided in the metadata.

    Usage Example:
    runs_info = find_runs(data, all_runs_metadata, trigger_starts, sampling_rate)

    Note:
    - Ensure the data array and all_runs_metadata are correctly formatted and complete.
    - The function logs detailed information about the process, facilitating troubleshooting and verification.
    """
    runs_info = []
    start_from_index = 0
    skipped_invalid_runs = set()

    for run_id, metadata in sorted(all_runs_metadata.items()):
        try:
            repetition_time = metadata['RepetitionTime']
            num_volumes = metadata['NumVolumes']

            logging.info(f"Searching for {run_id} with {num_volumes} volumes and TR={repetition_time}s")
            samples_per_volume = int(sampling_rate * repetition_time)

            # Continue searching for the next valid run, or skip if it's marked as invalid
            while start_from_index < len(trigger_starts):
                valid_run_start_index = None

                # Find the start index of the next valid run
                for i, start in enumerate(trigger_starts[start_from_index:], start_from_index):
                    if i + num_volumes > len(trigger_starts):
                        break  # Not enough triggers left for the full run

                    expected_end = trigger_starts[i] + (num_volumes - 1) * samples_per_volume
                    if trigger_starts[i + num_volumes - 1] <= expected_end:
                        valid_run_start_index = i
                        break

                if valid_run_start_index is not None:
                    # Check if the run is marked as invalid and should be skipped
                    if run_id in invalid_runs and run_id not in skipped_invalid_runs:
                        logging.info(f"Skipping invalid run {run_id} starting at index {valid_run_start_index}.")
                        skipped_invalid_runs.add(run_id)
                        start_from_index = valid_run_start_index + num_volumes
                        continue

                    # If it's not invalid, or has already been skipped once, process it
                    start_idx = trigger_starts[valid_run_start_index]
                    end_idx = start_idx + num_volumes * samples_per_volume

                    runs_info.append({
                        'run_id': run_id,
                        'start_index': start_idx,
                        'end_index': end_idx,
                        'metadata': metadata
                    })
                    logging.info(f"Found valid run {run_id} from index {start_idx} to {end_idx}.")
                    start_from_index = valid_run_start_index + num_volumes
                    break

                else:
                    logging.warning(f"No more valid trigger sequences found for {run_id} after index {start_from_index}.")
                    break

        except KeyError as e:
            logging.error(f"Key error for {run_id}: {e}")
        except IndexError as e:
            logging.error(f"Index error while processing {run_id}: {e}")
        except Exception as e:
            logging.error(f"Unexpected error while processing {run_id}: {e}")

    logging.info(f"Total valid runs identified: {len(runs_info)}")
    if not runs_info:
        logging.warning("No valid runs found. Check the trigger starts and metadata.")

    return runs_info

# Segments the runs based on the information provided by find_runs() and writes them to output files. 
def segment_runs(runs_info, output_dir, metadata_dict, labels, subject_id, session_id):
    """
    Parameters:
    - runs_info (list): A list of dictionaries, each containing information for a segmented run.
    - output_dir (str): Directory where output files will be saved.
    - metadata_dict (dict): Additional metadata for all runs.
    - labels (list of str): Labels of the data channels.
    - subject_id (str): Identifier for the subject.
    - session_id (str): Identifier for the session.

    Returns:
    - list: A list of file paths to the written output files.

    This function plays a crucial role in organizing and saving segmented physiological data into individual files, 
    each corresponding to a specific fMRI run. It's a vital step for further data analysis in neuroimaging studies.

    Usage Example:
    output_files = segment_runs(runs_info, output_dir, metadata_dict, labels, subject_id, session_id)

    Note:
    - Ensure that the output directory exists and is writable.
    - The function assumes that the data processing method (e.g., filtering, normalization) is defined elsewhere.
    """
    
    output_files = []
    for run in runs_info:
        run_id = run['run_id']
        data = run['data']
        start_index = run['start_index']
        end_index = run['end_index']
        run_metadata = run['run_metadata']
        task_name = run_metadata['TaskName']  # Extract the TaskName from the metadata.
        
        logging.info("Segmenting full run %s from index %d to %d", run_id, start_index, end_index)

        # Write the processed data to an output file.
        try:
            # Call the existing write_output_files function with the appropriate parameters.
            output_file_path = os.path.join(output_dir, f"{subject_id}_{session_id}_task-learn_{run_id}_physio.tsv.gz")
            write_output_files(data, run_metadata, metadata_dict, labels, output_dir, subject_id, session_id, run_id)
            
            output_files.append(output_file_path)
            logging.info("Output file for run %s written to %s", run_id, output_file_path)

        except IOError as e:
            logging.error("Failed to write output file for full run %s: %s", run_id, e, exc_info=True)
            raise # Propagate the exception for further handling.
    
    logging.info(f"Completed writing output files for all segmented full runs.")
    return output_files

# Creates a metadata dictionary for a run with channel-specific information.
def create_metadata_dict(run_info, sampling_rate, bids_labels_list, units_dict):
    """
    Parameters:
    - run_info (dict): Information about the run, including start index and run ID.
    - sampling_rate (int): The sampling rate of the physiological data.
    - bids_labels_list (list): List of BIDS-compliant labels for data channels.
    - units_dict (dict): Mapping from BIDS labels to their units.

    Returns:
    - dict: A dictionary containing metadata for the run.

    This function is vital for associating physiological data with its metadata, ensuring the data's integrity 
    and facilitating its use in research and analysis.

    Note:
    - The function assumes that the run_info dictionary and units_dict are correctly formatted and provided.
    - Channel-specific metadata is predefined and may need updates based on new information or different setups.
    
    Usage Example:
    run_metadata = create_metadata_dict(run_info, sampling_rate, bids_labels_list, units_dict)
    """

    # Initialize the metadata dictionary with common information.
    metadata_dict = {
        "RunID": run_info['run_id'],
        "SamplingFrequency": sampling_rate, # BIDS label
        "SamplingRate": {
            "Value": sampling_rate,  
            "Units": "Hz" # provide units
        },
        "StartTime": run_info['start_index'] / sampling_rate, # BIDS label
        "StartTimeSec": {
            "Value": run_info['start_index'] / sampling_rate, 
            "Description": "Start time of the current run relative to recording onset",
            "Units": "seconds" # provide units
        },
        "StartTimeMin": {
            "Value": (run_info['start_index'] / sampling_rate)/60, 
            "Description": "Start time of the current run relative to recording onset",
            "Units": "minutes"
        },
        "Columns": bids_labels_list,
        "Manufacturer": "Biopac",
        "Acquisition Software": "Acqknowledge",
    }

    # Channel-specific metadata.
    channel_metadata = {
        "cardiac": {
            "Description": "Continuous ECG measurement",
            "Placement": "Lead 1",
            "Gain": 500,
            "35HzLPN": "off / 150HzLP",
            "HPF": "0.05 Hz",
        },
        "respiratory": {
            "Description": "Continuous measurements by respiration belt",
            "Gain": 10,
            "LPF": "10 Hz",
            "HPF1": "DC",
            "HPF2": "0.05 Hz",
        },
        "eda": {
            "Description": "Continuous EDA measurement",
            "Placement": "Right plantar instep",
            "Gain": 5,
            "LPF": "1.0 Hz",
            "HPF1": "DC",
            "HPF2": "DC",
        },
        "trigger": {
            "Description": "fMRI Volume Marker",
        },
        "ppg": {
            "Description": "Continuous PPG measurement",
            "Placement": "Left index toe",
            "Gain": 10,
            "LPF": "3.0 Hz",
            "HPF1": "0.5 Hz",
            "HPF2": "0.05 Hz",
        }
    }

    # Add channel-specific metadata to the dictionary if the channel is present.
    for channel in channel_metadata:

    # Assuming 'channel' is a variable in this function that holds the current channel being processed.
        if channel in bids_labels_list:
            channel_specific_metadata = channel_metadata[channel]
            
            # Set the 'Units' dynamically based on the units_dict
            channel_specific_metadata['Units'] = units_dict.get(channel, "Unknown")
            metadata_dict[channel] = channel_specific_metadata
    
    logging.info(f"Full run metadata dictionary created for {run_info['run_id']}")
    return metadata_dict

# Creates metadata for an event segment, including run details and specific event information.
def create_event_metadata_dict(run_info, segment_length, sampling_rate, repetition_time, 
                               bids_labels_list, units_dict, trial_type, 
                               segment_start_time, event_onset, segment_duration):
    """
    Parameters:
    - run_info (dict): Information about the run, including the start index and run ID.
    - segment_length (int): Length of the segment in data points.
    - sampling_rate (int): Sampling rate of the data in Hz.
    - repetition_time (float): Repetition time for the fMRI volumes in seconds.
    - bids_labels_list (list): List of BIDS-compliant labels for data channels.
    - units_dict (dict): Mapping from BIDS labels to their units.
    - trial_type (str): Type of trial, e.g., 'sequence', 'random'.
    - segment_start_time (float): Start time of the segment in seconds.
    - event_onset (float): Onset time of the event relative to the start of the run.
    - segment_duration (float): Duration of the segment in seconds.

    Returns:
    - dict: A dictionary containing metadata for the event segment.

    The function meticulously details each aspect of an event segment, crucial for accurate data analysis 
    and subsequent research applications. It assumes the input parameters are correctly formatted and provided.

    Note:
    - This function can be adapted to different experimental setups and data structures.
    - The channel-specific metadata may require updates to reflect specific measurement techniques or equipment used.
    
    Usage Example:
    event_metadata = create_event_metadata_dict(run_info, segment_length, sampling_rate, repetition_time, 
                                                bids_labels_list, units_dict, trial_type, 
                                                segment_start_time, event_onset, segment_duration)
    """

    # Calculate the number of volumes and the duration of the segment.
    num_volumes_events = int(segment_length / (sampling_rate * repetition_time))
    segment_duration_min = (segment_length / sampling_rate) / 60
    
    # Logging for debugging and verification.
    logging.info(f"Number of volumes covered by '{trial_type}' block: {num_volumes_events}")
    #logging.info(f"Segment length: {segment_length} data points")
    logging.info(f"The total duration of the '{trial_type}' segment in {run_info['run_id']} is {segment_duration_min} minutes")
    #logging.info((run_info['start_index'] / sampling_rate) + event_onset)

    # Initialize the metadata dictionary with segment-specific information.
    metadata_dict_events = {
        "RunID": run_info['run_id'],
        "NumVolumes": num_volumes_events,
        "RepetitionTime": repetition_time,
        "SamplingFrequency": sampling_rate,
        "StartTime": (run_info['start_index'] / sampling_rate) + event_onset,
        "StartTimeSec": {
            "Value": (run_info['start_index'] / sampling_rate) + event_onset,
            "Description": "Start time of the segment relative to recording onset",
            "Units": "seconds"
        },
        "StartTimeMin": {
            "Value": ((run_info['start_index'] / sampling_rate) + event_onset) / 60,
            "Description": "Start time of the segment relative to recording onset",
            "Units": "minutes"
        },
        "DurationMin": {
            "Value": segment_duration_min,
            "Description": "Duration of the segment",
            "Units": "minutes"
        },
        "TrialType": trial_type,  # Include trial type in the metadata.
        "Columns": bids_labels_list
    }
  
    # Channel-specific metadata.
    channel_event_metadata = {
        "cardiac": {
            "Description": "Continuous ECG measurement",
            "Placement": "Lead 1",
            "Gain": 500,
            "35HzLPN": "off / 150HzLP",
            "HPF": "0.05 Hz",
        },
        "respiratory": {
            "Description": "Continuous measurements by respiration belt",
            "Gain": 10,
            "LPF": "10 Hz",
            "HPF1": "DC",
            "HPF2": "0.05 Hz",
        },
        "eda": {
            "Description": "Continuous EDA measurement",
            "Placement": "Right plantar instep",
            "Gain": 5,
            "LPF": "1.0 Hz",
            "HPF1": "DC",
            "HPF2": "DC",
        },
        "trigger": {
            "Description": "fMRI Volume Marker",
        },
        "ppg": {
            "Description": "Continuous PPG measurement",
            "Placement": "Left index toe",
            "Gain": 10,
            "LPF": "3.0 Hz",
            "HPF1": "0.5 Hz",
            "HPF2": "0.05 Hz",
        }
    }

    # Add channel-specific metadata to the dictionary if the channel is present.
    for channel in channel_event_metadata:
    # Assuming 'channel' is a variable in this function that holds the current channel being processed.
        if channel in bids_labels_list:
            channel_specific_event_metadata = channel_event_metadata[channel]
            # Set the 'Units' dynamically based on the units_dict.
            channel_specific_event_metadata['Units'] = units_dict.get(channel, "Unknown")
            metadata_dict_events[channel] = channel_specific_event_metadata

    logging.info(f"Event metadata dictionary created for {run_info['run_id']} {trial_type} segment.")
    return metadata_dict_events

# Writes segmented data to TSV and JSON files in the Brain Imaging Data Structure (BIDS) format.
def write_output_files(segmented_data, run_metadata, metadata_dict, bids_labels_list, output_dir, subject_id, session_id, run_id):
    """
    Parameters:
    - segmented_data (numpy.ndarray): The data segmented for a run.
    - run_metadata (dict): The metadata for the run.
    - metadata_dict (dict): The additional metadata for the run.
    - bids_labels_list (list of str): The BIDS-compliant labels of the data channels.
    - output_dir (str): The directory to which the files should be written.
    - subject_id (str): The identifier for the subject.
    - session_id (str): The identifier for the session.
    - run_id (str): The identifier for the run.

    Raises:
    - Exception: If any error occurs during the file writing process, it logs the error and re-raises the exception.

    The function creates TSV and JSON files named according to BIDS naming conventions.
    It writes the segmented data to a compressed TSV file and the combined metadata to a JSON file.

    Usage Example:
    write_output_files(segmented_data, run_meta, meta_dict, labels, '/output/path', 'sub-01', 'ses-1', 'run-01')

    Example Usage:
    write_output_files(segmented_data_array, run_metadata_dict, additional_metadata_dict, channel_labels, '/output/path', 'sub-01', 'sess-01', 'run-01')
    """
    try:
        
        # Move up four levels to get to dataset_root_dir
        dataset_root_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(output_dir))))

        #logging.info(f"Dataset root directory: {dataset_root_dir}")

        output_derivatives_dir = os.path.join(dataset_root_dir, 'derivatives', 'physio', 'learn', 'runs')

        # Ensure the output directories exist
        os.makedirs(output_dir, exist_ok=True)
        os.makedirs(output_derivatives_dir, exist_ok=True)

        # Define filenames based on the BIDS format.
        tsv_filename = f"{subject_id}_{session_id}_task-learn_{run_id}_physio.tsv.gz"
        json_filename = f"{subject_id}_{session_id}_task-learn_{run_id}_physio.json"
        
        # Prepare the full file paths.
        tsv_file_path = os.path.join(output_dir, tsv_filename)
        json_file_path = os.path.join(output_dir, json_filename)

        # Create a DataFrame with the segmented data and correct labels.
        df = pd.DataFrame(segmented_data, columns=bids_labels_list)

        # Save the DataFrame to a TSV file with GZip compression.
        df.to_csv(tsv_file_path, sep='\t', index=False, compression='gzip')
        logging.info(f"TSV file written: {tsv_file_path}")

        # Copy the file
        shutil.copy2(tsv_file_path, output_derivatives_dir)

        # Log the action
        logging.info(f"TSV file copied to: {output_derivatives_dir}")

        # Merge the run-specific metadata with the additional metadata.
        combined_metadata = {**run_metadata, **metadata_dict}

        # Write the combined metadata to a JSON file.
        with open(json_file_path, 'w') as json_file:
            json.dump(combined_metadata, json_file, indent=4)
        logging.info(f"JSON file written: {json_file_path}")

        # Copy the file
        shutil.copy2(json_file_path, output_derivatives_dir)

        # Log the action
        logging.info(f"JSON file copied to: {output_derivatives_dir}")


    except Exception as e:
        # Log any exceptions that occur during the file writing process.
        logging.error(f"Failed to write output files for {run_id}: {e}", exc_info=True)
        raise
    
    # Log the successful writing of files.
    logging.info(f"Full length run output files for {run_id} also written successfully to {output_dir}")

# Writes output files (.tsv and .json)for each event segment of physiological data in BIDS format.
def write_event_output_files(event_segments, run_metadata, metadata_dict_events, bids_labels_list, output_dir, subject_id, session_id, run_id):
    """
    Parameters:
    - event_segments (list of tuples): Each tuple contains a segment of data (numpy array) and its corresponding trial type (str), start time (float), and duration (float).
    - run_metadata (dict): Metadata for the run.
    - metadata_dict_events (dict): Additional metadata for the events.
    - bids_labels_list (list of str): BIDS-compliant labels of the data channels.
    - output_dir (str): Directory to write the files.
    - subject_id (str): Identifier for the subject.
    - session_id (str): Identifier for the session.
    - run_id (str): Identifier for the run.

    Raises:
    - Exception: If an error occurs during file writing.

    Notes:
    - This function requires 'pandas' for DataFrame operations and 'os', 'json' for file handling.
    - It assumes that event_segments contains tuples with exactly four elements: data, trial type, start time, and duration.
    - The function creates TSV and JSON files for each event segment, naming them according to the BIDS format and including relevant metadata.
    - Error handling is included to manage unexpected segment formats or file writing issues.

    Example Usage:
    write_event_output_files(event_segments, run_metadata, metadata_dict_events, channel_labels, '/output/path', 'sub-01', 'sess-01', 'run-01')
    """
    try:
        
        # Move up four levels to get to dataset_root_dir
        dataset_root_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(output_dir))))

        logging.info(f"Dataset root directory: {dataset_root_dir}")

        output_derivatives_dir = os.path.join(dataset_root_dir, 'derivatives', 'physio', 'learn', 'events')

        # Ensure the output directories exist
        os.makedirs(output_dir, exist_ok=True)
        os.makedirs(output_derivatives_dir, exist_ok=True)
      
        # Loop through each segment and write to individual files.
        for i, segment_info in enumerate(event_segments):

            # Adjust the segment index for logging and naming (starting from 1 instead of 0).
            segment_index = i + 1

            # Unpack the segment_info tuple.
            if len(segment_info) == 4:  # Assuming there are four elements in the tuple.
                segment_data, trial_type, segment_start_time, segment_duration = segment_info
            else:
                logging.error(f"Unexpected event segment format for {run_id}: {segment_info}")
                continue  # Skip this segment.

            # Define filenames based on the BIDS format.
            tsv_event_filename = f"{subject_id}_{session_id}_task-learn_{run_id}_recording-{trial_type}_physio.tsv.gz"
            json_event_filename = f"{subject_id}_{session_id}_task-learn_{run_id}_recording-{trial_type}_physio.json"

            # Prepare file paths.
            tsv_event_file_path = os.path.join(output_dir, tsv_event_filename)
            json_event_file_path = os.path.join(output_dir, json_event_filename)
       
            # Check if segment data is not empty.
            if segment_data.size == 0:
                logging.warning(f"No data found for event segment {segment_index} of trial type '{trial_type}' in run '{run_id}'. Skipping file writing.")
                continue

            # Create a DataFrame and save to a TSV file.
            df = pd.DataFrame(segment_data, columns=bids_labels_list)
            df.to_csv(tsv_event_file_path, sep='\t', index=False, compression='gzip')
            logging.info(f"TSV file written: {tsv_event_file_path}")

            # Copy the file
            shutil.copy2(tsv_event_file_path, output_derivatives_dir)

            # Log the action
            logging.info(f"TSV file copied to: {output_derivatives_dir}")

            # Write the event metadata to a JSON file.
            with open(json_event_file_path, 'w') as json_file:
                json.dump(metadata_dict_events, json_file, indent=4)
               
            logging.info(f"JSON file written: {json_event_file_path}")

            # Copy the file
            shutil.copy2(json_event_file_path, output_derivatives_dir)

            # Log the action
            logging.info(f"JSON file copied to: {output_derivatives_dir}")

            # Log the successful writing of files.
            logging.info(f"Event segment output files for {run_id}, segment {segment_index}, trial type '{trial_type}' written successfully to {output_dir}")

    except Exception as e:
        # Log any exceptions during the file writing process.
        logging.error(f"Failed to write event output files for {run_id}: {e}", exc_info=True)
        raise

# Plots the segmented data for each run and saves the plots to a PNG file.
def plot_runs(original_data, segmented_data_list, runs_info, bids_labels_list, sampling_rate, plot_file_path, units_dict):
    """
    Parameters:
    - original_data (np.ndarray): The entire original data to be used as a background.
    - segmented_data_list (list of np.ndarray): Each element is an array of data for a run.
    - runs_info (list): Information about each run, including start and end indices and metadata.
    - bids_labels_list (list of str): BIDS-compliant channel labels.
    - sampling_rate (int): The rate at which data was sampled.
    - plot_file_path (str): The file path to save the plot.
    - units_dict (dict): A dictionary mapping labels to their units.

    Raises:
    - Exception: If an error occurs during plotting.

    Notes:
    - This function requires 'numpy' for numerical operations and 'matplotlib.pyplot' for plotting.
    - It creates a multi-panel plot, each panel representing a channel of the data.
    - The original data is plotted as a grey background, and segmented data are overlaid with different colors.
    - The function handles the generation of a dynamic legend and ensures proper layout adjustment.
    - Error handling includes logging the error and re-raising the exception for further handling.

    Example Usage:
    plot_runs(original_data_array, segmented_data_list, runs_info_list, channel_labels, 250, '/path/to/save/plot.png', units_dictionary)
    """

    try:
        # Define a list of colors for different runs.
        colors = [
            'red', 'green', 'blue', 'cyan', 'magenta', 'yellow', 'black', 'orange',
            'pink', 'purple', 'lime', 'indigo', 'violet', 'gold', 'grey', 'brown'
        ]

        # Create a figure and a set of subplots.
        fig, axes = plt.subplots(nrows=len(bids_labels_list), ncols=1, figsize=(20, 10))

        # Time axis for the original data.
        time_axis_original = np.arange(original_data.shape[0]) / sampling_rate / 60

        # Plot the entire original data as background.
        for i, label in enumerate(bids_labels_list):
            unit = units_dict.get(label, 'Unknown unit')
            axes[i].plot(time_axis_original, original_data[:, i], color='grey', alpha=0.5, label='Background' if i == 0 else "")
            axes[i].set_ylabel(f'{label}\n({unit})')

        # Overlay each segmented run on the background.
        for segment_index, (segment_data, run_info) in enumerate(zip(segmented_data_list, runs_info)):

            # Define time_axis_segment for each segmented run.
            time_axis_segment = np.arange(run_info['start_index'], run_info['end_index']) / sampling_rate / 60
            color = colors[segment_index % len(colors)]  # Cycle through colors.
            for i, label in enumerate(bids_labels_list):
                axes[i].plot(time_axis_segment, segment_data[:, i], color=color, label=f'{run_info["run_id"]}' if i == 0 else "")

        # Set the x-axis label for the bottom subplot.
        axes[-1].set_xlabel('Time (min)')

        # Collect handles and labels for the legend from all axes.
        handles, labels = [], []
        for ax in axes.flat:
            h, l = ax.get_legend_handles_labels()

            # Add the handle/label if it's not already in the list.
            for hi, li in zip(h, l):
                if li not in labels:
                    handles.append(hi)
                    labels.append(li)

        # Only create a legend if there are items to display.
        if handles and labels:
            ncol = min(len(handles), len(labels))
            fig.legend(handles, labels, loc='lower center', ncol=ncol)
        else:
            logging.info("No legend items to display.")

        # Apply tight layout with padding to make room for the legend and axis labels.
        fig.tight_layout(rect=[0.05, 0.1, 0.95, 0.97])  # Adjust the left and bottom values as needed.

        # Save and show the figure.
        plt.savefig(plot_file_path, dpi=600)
        #plt.show()  # Comment this line if you want to bypass plot display.
        logging.info(f"Plot saved to {plot_file_path}")

    except Exception as e:
        logging.error("Failed to plot runs: %s", e, exc_info=True)
        raise

# Segments the physiological data based on events described in the events DataFrame, considering 'sequence' and 'random' blocks.
def segment_data_by_events(data, events_df, sampling_rate, run_info, run_id):
    """
    Segments physiological data based on event descriptions, specifically for 'sequence' and 'random' events.
    
    This function handles special runs with only 'random' events (like 'run-00', 'run-07') differently from other runs.
    For regular runs, it segments the data into 'sequence' and 'random' blocks based on their occurrence order.

    Parameters:
    - data (np.ndarray): The original data array.
    - events_df (pd.DataFrame): Contains 'onset', 'duration', and 'trial_type' for each event.
    - sampling_rate (int): The rate at which data was sampled.
    - run_info (dict): Information about the current run, including the start index.
    - run_id (str): The identifier of the current run.

    Returns:
    - List of tuples: Each tuple contains a segment of data, its trial type, start time, and duration.
    
    Raises:
    - Exception: If any error occurs during the segmentation process.

    Dependencies:
    - Requires 'numpy' for data handling and 'pandas' for DataFrame operations.

    Example Usage:
    segments = segment_data_by_events(data_array, events_dataframe, 250, run_info_dict, 'run-01')
    """
    
    segments = []
    logging.info(f"Starting data segmentation by events for {run_id}.")

    # Special handling for 'random' block only runs (e.g., 'run-00', 'run-07').
    if run_id in ["run-00", "run-07"]:
        
        # Determine end time of the last 'random' event, configure logging outputs.
        last_random_event = events_df[events_df['trial_type'] == 'random'].iloc[-1]
        segment_end_time = int(last_random_event['onset'] + last_random_event['duration']) # seconds relative to the start of the run
        end_index = int((run_info['start_index'] + segment_end_time) * sampling_rate) # samples
        block_end_time_fixed = int(segment_end_time + (run_info['start_index'] / sampling_rate)) # seconds relative to the start of the session
        block_start_time_fixed = int(run_info['start_index'] / sampling_rate) # seconds relative to the start of the session
        block_duration_fixed = int(block_end_time_fixed - block_start_time_fixed) # seconds relative to the start of the session
        fixed_start_index = int(run_info['start_index']) # samples
        fixed_end_index = int(run_info['start_index'] + (segment_end_time * sampling_rate)) # samples
        fixed_index_length = fixed_end_index - fixed_start_index # samples

        # Special handling for 'random' block only runs (e.g., 'run-00', 'run-07') - process logging. 
        logging.info(f"Event fixed random segment Start Time: {block_start_time_fixed} seconds")
        logging.info(f"Event fixed random segment End Time: {block_end_time_fixed} seconds")
        logging.info(f"Event fixed random segment Duration: {block_duration_fixed} seconds")
        logging.info(f"Event fixed random segment Start index: {run_info['start_index']} samples")
        logging.info(f"Event fixed random segment End index: {fixed_end_index} samples")
        logging.info(f"Event fixed random segment Length: {fixed_index_length} samples")
        logging.info(
            f"Segment for this fixed random block starts at {block_start_time_fixed} seconds relative to the session start, " 
            f"ends at {block_end_time_fixed} seconds, encompassing a total duration of {block_duration_fixed} seconds " 
            f"and covering {fixed_index_length} data points, from start index {run_info['start_index']} to end index {fixed_end_index} in the data array.")

        # Extract segment for 'random' block.
        segment_data = data[fixed_start_index:fixed_end_index]
        segments.append((segment_data, 'random', run_info['start_index'] / sampling_rate, segment_end_time))
        logging.info(f"Segmented 'fixed random' block for {run_id} with random segment ending after {segment_end_time} seconds relative to the start of the run.")
        return segments

    trial_index = 0 # initialize trial index counter

    # Determine the first onset time for both 'sequence' and 'random' events
    first_sequence_event = events_df[events_df['trial_type'] == 'sequence'].iloc[0]['onset'] if not events_df[events_df['trial_type'] == 'sequence'].empty else float('inf')
    first_random_event = events_df[events_df['trial_type'] == 'random'].iloc[0]['onset'] if not events_df[events_df['trial_type'] == 'random'].empty else float('inf')

    # Define the order in which to process the segments
    segment_order = ['sequence', 'random'] if first_sequence_event < first_random_event else ['random', 'sequence']

    # for trial_type in ['sequence', 'random']:
    for trial_type in segment_order:
        block_events = events_df[events_df['trial_type'] == trial_type]
        if block_events.empty:
            logging.warning(f"No events found for trial type '{trial_type}' in {run_id}.")
            continue

        trial_index += 1  # Increment trial index counter

        # Calculate start and end times of the block.
        block_start_onset = block_events.iloc[0]['onset']
        last_event = block_events.iloc[-1]
        block_end_time = last_event['onset'] + last_event['duration']
        block_end_time_int = int(block_end_time)

        # Log start and end times of the block.
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s first onset begins at {block_start_onset} seconds relative to start of the run.")
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s last onset begins at {last_event['onset']} seconds relative to start of the run.")
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s end time is {block_end_time_int} seconds relative to start of the session.")
        
        # Convert times to indices
        start_index = int(block_start_onset * sampling_rate) + run_info['start_index']
        end_index = int(block_end_time * sampling_rate) + run_info['start_index']
        segment_data = data[start_index:end_index]
        index_length = end_index - start_index
        
        # Log start and end indices of the block.
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s start index begins at sample #{start_index}")
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s end index ends at sample #{end_index}")
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s index length spans {index_length} samples")
        
        # Log start and end times of the segment.
        segment_start_time = start_index / sampling_rate  # Convert to seconds
        segment_end_time = end_index / sampling_rate  # Convert to seconds
        segment_duration = segment_end_time - segment_start_time

        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s start time begins at {segment_start_time} seconds relative to start of the session.")
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s ends at {segment_end_time} seconds relative to start of the session.")
        logging.info(f"The {run_id} '{trial_type}' segment {trial_index}'s duration spans {segment_duration} seconds")
        
        # Extract data for the segment.
        segment_data = data[start_index:end_index]

        # Append segment info.
        segments.append((segment_data, trial_type, segment_start_time, segment_duration))
        #logging.info(f"The {run_id} '{trial_type}' segment {trial_index} for Start time {segment_start_time} s, Duration {segment_duration} s")
        logging.info(
                    f"The {run_id} '{trial_type}' segment {trial_index} starts at {segment_start_time} seconds relative to the session start," 
                    f"ends at {segment_end_time} seconds, encompassing a total duration of {segment_duration} seconds " 
                    f"and covering {index_length} samples, from start index {start_index} to end index {end_index} in the data array.")

    logging.info(f"Data segmentation completed for {run_id}. Total number of segments: {len(segments)}")
    return segments

# Plots segmented data over the original data for visual comparison.
def plot_runs_with_events(original_data, event_segments_by_run, events_df, sampling_rate, plot_events_file_path, units_dict, bids_labels_list, run_info_dict):
    """
    Plots segmented data for each run, overlaying it on the background of the original data. 
    This helps in visual comparison between the original and segmented data.

    Parameters:
    - original_data (np.ndarray): The original full dataset.
    - event_segments_by_run (dict): Segments of data for each run, keyed by run ID.
    - events_df (pd.DataFrame): DataFrame containing event information.
    - sampling_rate (int): Data sampling rate.
    - plot_events_file_path (str): Path to save the plot.
    - units_dict (dict): Dictionary mapping data channels to their units.
    - bids_labels_list (list of str): List of BIDS-compliant labels for data channels.
    - run_info_dict (dict): Dictionary containing run-specific information.

    Notes:
    - The function creates a plot with multiple subplots, each representing a data channel.
    - Original data is plotted as a grey background for reference.
    - Segmented data from different runs are overlaid with distinct colors for 'sequence' and 'random' trial types.
    - This function uses 'matplotlib.pyplot' for plotting and 'numpy' for numerical operations.
    - Logging is used to track the progress and status of the plotting process.

    Raises:
    - Exception: If any error occurs during the plotting process.

    Example Usage:
    plot_runs_with_events(original_data, event_segments_by_run, events_df, 250, '/path/to/save/plot.png', units_dict, channel_labels, run_info)
    """

    # Define colors for different trial types.
    colors = {'sequence': 'red', 'random': 'blue'}

    # Create custom legend handles.
    custom_handles = [Line2D([0], [0], color=colors['sequence'], lw=2, label='Sequence'),
                    Line2D([0], [0], color=colors['random'], lw=2, label='Random')]
    
    # Create figure and subplots.
    fig, axes = plt.subplots(nrows=len(bids_labels_list), ncols=1, figsize=(20, 10))
  
    # Plot the original data as a background reference.
    time_axis_original = np.arange(original_data.shape[0]) / sampling_rate / 60  # Convert to minutes.
    for i, label in enumerate(bids_labels_list):
        axes[i].plot(time_axis_original, original_data[:, i], color='grey', alpha=0.5, label='Background')
        axes[i].set_ylabel(f'{label} ({units_dict.get(label, "Unknown unit")})')

    # Overlay segmented data on top of the original data.
    for run_id, segments in event_segments_by_run.items():
        for segment_data, trial_type, segment_start_time, segment_duration in segments:
            # Calculate the time axis for the segment.
            start_time = segment_start_time / 60  # Convert to minutes.
            end_time = start_time + segment_duration / 60
            time_axis_segment = np.linspace(start_time, end_time, len(segment_data))

            # Plot each segment.
            for i, label in enumerate(bids_labels_list):
                axes[i].plot(time_axis_segment, segment_data[:, i], color=colors[trial_type])
            
    # Set x-axis label and legend.
    axes[-1].set_xlabel('Time (min)')
    handles, labels = axes[-1].get_legend_handles_labels()
    fig.legend(handles=custom_handles, loc='lower center', ncol=len(custom_handles))

    # Adjust layout to fit all elements.
    plt.tight_layout()
    plt.subplots_adjust(bottom=0.2)  # Adjust to accommodate the legend.

    # Save and display the plot.
    plt.savefig(plot_events_file_path, dpi=600)
    #plt.show() # Uncomment to display figure. 
    logging.info(f"Plot saved to {plot_events_file_path}")

# SeriesDecription to run ID mapping for invalid run use cases. 
def determine_expected_series_description(run_id):
    series_descriptions = {
    "run-00": "sms3_TASK_seqexec_PRE",
    "run-01": "sms3_TASK_1",
    "run-02": "sms3_TASK_2",
    "run-03": "sms3_TASK_3",
    "run-04": "sms3_TASK_4",
    "run-05": "sms3_TASK_5",
    "run-06": "sms3_TASK_6",
    "run-07": "sms3_TASK_seqexec_POST"
    }

    return series_descriptions.get(run_id, '')
    
# Main function to orchestrate the conversion of physiological data to BIDS format.
def main(physio_root_dir, bids_root_dir, force_process_flag=False, invalid_runs_str=''):
    """
    Main function to orchestrate the conversion of physiological data to BIDS format.

    Parameters:
    - physio_root_dir (str): Directory containing the physiological .mat files.
    - bids_root_dir (str): Root directory of the BIDS dataset where output files will be saved.

    The function performs the following steps:
    1. Load physiological data from .mat files.
    2. Rename channels according to BIDS conventions.
    3. Extract metadata from associated JSON files.
    4. Segment the data into runs based on triggers and metadata.
    5. Write segmented data to output files in BIDS format.
    6. Plot the physiological data for all runs.

    Usage Example:
    main('/path/to/physio_data', '/path/to/bids_dataset')

    Dependencies:
    - Requires several helper functions defined in the same script for loading data, renaming channels,
      extracting metadata, finding runs, segmenting data, writing output files, and plotting.
    """

    # Define the known sampling rate.
    sampling_rate = 5000  # Replace with the actual sampling rate if different.
    
    # Parse the invalid runs string into a list
    invalid_runs = invalid_runs_str.split(',') if invalid_runs_str else []

    try:
        # Extract subject and session IDs from the path.
        subject_id, session_id = extract_subject_session(physio_root_dir)

        # Define output directory for the BIDS dataset.
        output_dir = os.path.join(bids_root_dir, subject_id, session_id, 'func')
        
        # Search for *_physio.tsv files in the output directory
        file_pattern = os.path.join(output_dir, f"{subject_id}_{session_id}_task-learn_run-01_physio.tsv.gz")
        existing_files = glob.glob(file_pattern)
        print(f"File pattern is: {file_pattern}")

        if existing_files:
            # Print the list of existing files to the console
            print("The following *_physio.tsv.gz files already exist:", existing_files)
            return
            # raise ValueError("Files already exist, processing may already be complete")
            # Alternatively, use sys.exit() or return, depending on your requirements

        # Setup logging after extracting subject_id and session_id.
        log_file_path = setup_logging(subject_id, session_id, bids_root_dir)
        logging.info("Processing subject: %s, session: %s", subject_id, session_id)
        logging.info(f"Force processing flag: {force_process_flag}")
        logging.info(f"Invalid runs: {invalid_runs}")
        
        # Load physiological data from the .mat file. 
        mat_file_path = os.path.join(physio_root_dir, f"{subject_id}_{session_id}_task-learn_physio.mat")
        labels, data, units = load_mat_file(mat_file_path)
        if data is None or not data.size:
            raise ValueError("Data is empty after loading.")

        # Rename channels according to BIDS conventions.
        bids_labels_dictionary, bids_labels_list = rename_channels(labels)

        # Create a dictionary mapping from original labels to units.
        original_labels_to_units = dict(zip(labels, units))

        # Create a mapping of original labels to their indices in the data array.
        original_label_indices = {label: idx for idx, label in enumerate(labels)}

        # Now create the units_dict by using the bids_labels_dictionary to look up the original labels.
        units_dict = {
            bids_labels_dictionary[original_label]: unit
            for original_label, unit in original_labels_to_units.items()
            if original_label in bids_labels_dictionary
        }
    
        # Filter the data array to retain only the columns with BIDS labels.
        segmented_data_bids_only = data[:, [original_label_indices[original] for original in bids_labels_dictionary.keys()]]

        # Process JSON files to extract metadata for each run.
        json_file_paths = glob.glob(os.path.join(bids_root_dir, subject_id, session_id, 'func', '*_bold.json'))
 
        # Assume json_file_paths is a list of paths to your JSON files.
        all_runs_metadata = {}
        processed_jsons = set()

        # Sort json_file_paths based on the run ID number extracted from the file name.
        sorted_json_file_paths = sorted(
            json_file_paths,
            key=lambda x: int(re.search(r"run-(\d+)_bold\.json$", x).group(1)))
     
        # Now use the sorted_json_file_paths instead of the unsorted json_file_paths.
        for json_file_path in sorted_json_file_paths:
            run_id, run_metadata = extract_metadata_from_json(json_file_path, processed_jsons)
            if run_id and run_metadata:  # Checks that neither are None.

                # Retrieve the actual SeriesDescription from the run_metadata.
                actual_series_description = run_metadata['SeriesDescription']

                # Logic to determine expected series description based on run_id and run mapping.
                expected_series_description = determine_expected_series_description(run_id)
                
                # Check for mismatch between actual and expected SeriesDescription.
                if actual_series_description != expected_series_description:
                    if run_id not in invalid_runs and invalid_runs:
                        logging.warning(f"Warning: SeriesDescription mismatch for {run_id} but it's marked as invalid. Expected: {expected_series_description}, Found: {actual_series_description}.")
                    elif run_id not in invalid_runs and not invalid_runs:
                        logging.error(f"Error: SeriesDescription mismatch for {run_id}. Expected: {expected_series_description}, Found: {actual_series_description}.")
                        sys.exit(1)  # Terminate the script with an error if a mismatch is found and it's not marked as invalid.
                    
                # Get metadata for the current run.
                all_runs_metadata[run_id] = run_metadata

        # Sort the run IDs based on their numeric part.
        sorted_run_ids = sorted(all_runs_metadata.keys(), key=lambda x: int(x.split('-')[1]))

        # Create an ordered dictionary that respects the sorted run IDs.
        sorted_all_runs_metadata = OrderedDict((run_id, all_runs_metadata[run_id]) for run_id in sorted_run_ids)

        # Sort the expected runs to ensure they are processed in the correct order.
        expected_runs = sorted(all_runs_metadata.keys(), key=lambda x: int(x.split('-')[1]))

        # Find the index for the trigger channel using the BIDS labels list.
        if 'trigger' in bids_labels_list:
            trigger_channel_index = bids_labels_list.index('trigger')
            mri_trigger_data = data[:, trigger_channel_index]
        else:
            raise ValueError("Trigger channel not found in BIDS labels.")

        # Extract trigger points with the appropriate threshold and min_consecutive.
        # Note: Adjust the threshold and min_consecutive based on your specific data
        trigger_starts = extract_trigger_points(mri_trigger_data, threshold=5, min_consecutive=5) # Note: Adjust the threshold and min_consecutive based on your specific data.
        if len(trigger_starts) == 0:
            raise ValueError("No trigger points found, please check the threshold and min_consecutive parameters.")
        logging.info("Trigger starts: %s", len(trigger_starts)) # trigger_starts.

        # Find runs using the extracted trigger points and the list of invalid runs
        runs_info = find_runs(data, all_runs_metadata, trigger_starts, sampling_rate, invalid_runs)
        if len(runs_info) == 0:
            raise ValueError("No runs were found, please check the triggers and metadata.")

        # Catch error if no runs were found.
        if not runs_info:
            raise ValueError("No runs were found. Please check the triggers and metadata.")
 
        # expected_runs = set(run_info['run_id'] for run_info in runs_info)
        # if expected_runs != set(all_runs_metadata.keys()):
        #     raise ValueError("Mismatch between found runs and expected runs based on JSON metadata.")
        found_runs = sorted(run_info['run_id'] for run_info in runs_info)
        # Verify that the found runs match the expected runs from the JSON metadata.
        logging.info(f"Expected runs: {expected_runs}")
        logging.info(f"Found runs: {found_runs}")
        if expected_runs != found_runs:
            if not args.force:
                raise ValueError("Mismatch between found runs and expected runs based on JSON metadata.")
            else:
                logging.warning("Warning: Mismatch between found runs and expected runs. Proceeding due to --force flag.")
        
        # Create a mapping from run_id to run_info.
        run_info_dict = {info['run_id']: info for info in runs_info}

        # # Verify that the found runs match the expected runs from the JSON metadata.
        # if not set(sorted_run_ids) == set(run_info_dict.keys()):
        #     raise ValueError("Mismatch between found runs and expected runs based on JSON metadata.")

        event_segments_by_run = {}  # Dictionary to store event segments for each run.

        # Create a set of valid run IDs from runs_info for quick lookup
        valid_run_ids = {info['run_id'] for info in runs_info}

        # Segment runs and write output files for each run, using sorted_run_ids to maintain order.
        output_files = []
        for run_id in sorted_run_ids:
            if run_id not in valid_run_ids:
                logging.info(f"Skipping event segmentation for invalid or missing run {run_id}.")
                continue  # Skip this run
            try:
                run_info = run_info_dict[run_id]
                repetition_time = all_runs_metadata[run_id]['RepetitionTime']  # Retrieve repetition time from metadata.
                
                # Construct the events file path.
                events_file_path = os.path.join(
                    bids_root_dir, 
                    subject_id, 
                    session_id, 
                    'func', 
                    f"{subject_id}_{session_id}_task-learn_{run_id}_events.tsv"
                )
                logging.info(f"Events file path for {run_id} sourced from: {events_file_path}")
                
                # Read the events file into events_df.
                try:
                    events_df = pd.read_csv(events_file_path, sep='\t')
                except Exception as e:
                    logging.error(f"Error reading events file for {run_id}: {e}", exc_info=True)
                    continue  # Skip this run if events file cannot be read.
                
                # Process each segment within the run.
                event_segments = segment_data_by_events(data, events_df, sampling_rate, run_info, run_id)
                event_segments_by_run[run_id] = event_segments

                start_index, end_index = run_info['start_index'], run_info['end_index']
                segmented_data = data[start_index:end_index]

                output_dir = os.path.join(bids_root_dir, subject_id, session_id, 'func')
                
                # Create the metadata dictionary for the current run.
                metadata_dict = create_metadata_dict(run_info, sampling_rate, bids_labels_list, units_dict)
        
                for segment in event_segments:
                    if len(segment) != 4:
                        logging.error(f"Invalid segment format for {run_id}: {segment}")
                        continue  # Skip malformed segments.

                    segment_data, trial_type, segment_start_time, segment_duration = segment
            
                    segment_length = len(segment_data)
                
                    # Retrieve the onset time for the segment from events_df.
                    event_onset = events_df[events_df['trial_type'] == trial_type]['onset'].iloc[0]
                    logging.info(f"First '{trial_type}' event onset beings at {event_onset} seconds relative to the start of the run")

                    # Create event metadata for each segment.
                    metadata_dict_events = create_event_metadata_dict(
                        run_info, segment_length, sampling_rate, repetition_time, 
                        bids_labels_list, units_dict, trial_type, 
                        segment_start_time, event_onset, segment_duration
                    )

                    # logging.info(f"Metadata event dictionary for {run_id}: {metadata_dict_events}")
                    logging.info(f"Writing event output files for {run_id} '{trial_type}' segment")
                    write_event_output_files(
                        [segment], sorted_all_runs_metadata[run_id], metadata_dict_events, 
                        bids_labels_list, output_dir, subject_id, session_id, run_id
                    )
        
                # Call the write_output_files function with the correct parameters.
                output_files.append(write_output_files(
                    segmented_data_bids_only[start_index:end_index],
                    sorted_all_runs_metadata[run_id],
                    metadata_dict,
                    bids_labels_list,
                    output_dir,
                    subject_id,
                    session_id,
                    run_id
                ))  

            except KeyError as e:
                if not args.force:
                    raise  # Re-raise the exception if --force flag is not set
                else:
                    print(f"Warning: Run {e.args[0]} not found. Skipping due to --force flag.")
                    continue  # Skip this run and continue with the next
        
        # Create a list of segmented data for plotting.
        segmented_data_list = [segmented_data_bids_only[run_info['start_index']:run_info['end_index']] for run_info in runs_info]
        logging.info("Number of runs segmented: %s", len(segmented_data_list))

        # Filter the original data array to retain only the columns with BIDS labels.
        data_bids_only = data[:, [original_label_indices[original] for original in bids_labels_dictionary.keys()]]

        # Plot physiological data for all runs with the filtered background data.
        if segmented_data_list:
            
            # Extract the base name of the script without the .py extension.
            script_name = os.path.basename(__file__).replace('.py', '')

            # Construct the log directory path within 'doc/logs'
            log_dir = os.path.join(os.path.dirname(bids_root_dir), 'doc', 'logs', script_name, subject_id)

            # Create the log directory if it doesn't exist.
            if not os.path.exists(log_dir):
                os.makedirs(log_dir)
            
            # Move up four levels to get to dataset_root_dir
            dataset_root_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(output_dir))))
            output_derivatives_dir_plot_runs = os.path.join(dataset_root_dir, 'derivatives', 'physio', 'learn', 'runs', 'plots')
            output_derivatives_dir_plot_events = os.path.join(dataset_root_dir, 'derivatives', 'physio', 'learn', 'events', 'plots')
            os.makedirs(output_derivatives_dir_plot_runs, exist_ok=True)
            os.makedirs(output_derivatives_dir_plot_events, exist_ok=True)
            

            # Call the plot_runs function with the correct parameters and save plots to the log directory.
            logging.info("Preparing to plot runs.")
            log_file_path_plot_runs = os.path.join(log_dir, f"{subject_id}_{session_id}_task-learn_all_runs_physio.png")
            plot_file_path = log_file_path_plot_runs
            plot_runs(data_bids_only, segmented_data_list, runs_info, bids_labels_list, sampling_rate, plot_file_path, units_dict)

            # Log the action
            logging.info(f"Runs plot png file copied to: {output_derivatives_dir_plot_runs}")

            # Copy the file
            shutil.copy2(plot_file_path, output_derivatives_dir_plot_runs)

            logging.info("Preparing to plot event blocks.")
            log_file_path_plot_events = os.path.join(log_dir, f"{subject_id}_{session_id}_task-learn_all_blocks_physio.png")
            plot_events_file_path = log_file_path_plot_events
            plot_runs_with_events(data_bids_only, event_segments_by_run, events_df, sampling_rate, plot_events_file_path, units_dict, bids_labels_list, run_info_dict)

            # Log the action
            logging.info(f"Events plot png file copied to: {output_derivatives_dir_plot_events}")

            # Copy the file
            shutil.copy2(plot_events_file_path, output_derivatives_dir_plot_events)
        
        else:
            logging.error("No data available to plot.")

    except Exception as e:
        logging.error("An error occurred in the main function: %s", e, exc_info=True)
        raise

# Main function to run the script from the command line.
if __name__ == '__main__':
    """
    Command-line execution entry point for the script.

    This block allows the script to be run directly from the command line. It uses argparse to handle
    command-line arguments, specifically the paths to the directories containing the physiological data
    and the BIDS dataset. These arguments are then passed to the main function of the script.

   Usage:

    python BIDS_process_physio_ses_2.py <physio_root_dir> <bids_root_dir>

    """
    # Set up an argument parser to handle command-line arguments.
    parser = argparse.ArgumentParser(description="Convert physiological data to BIDS format.")

    # The first argument is the root directory containing the matlab physio data.
    parser.add_argument("physio_root_dir", help="Directory containing the physiological .mat file.")

    # The second argument is the root directory of the BIDS dataset.
    parser.add_argument("bids_root_dir", help="Path to the root of the BIDS dataset.")
    
    # The third argument is an optional flag to force processing even if run mismatch occurs.
    parser.add_argument('--force', action='store_true', help='Force processing even if run mismatch occurs')

    # Another optional argument is a comma-separated list of invalid runs (e.g., run-04,run-05).
    parser.add_argument('--invalid', type=str, help='Comma-separated list of invalid runs.', default='')

    # Parse the arguments provided by the user.
    args = parser.parse_args()
    
    # Starting script messages
    print(f"Starting script with provided arguments.")
    print(f"Physiological data directory: {args.physio_root_dir}")
    print(f"BIDS root directory: {args.bids_root_dir}")
    print(f"Force processing: {args.force}")
    print(f"Invalid runs: {args.invalid}")
 
    # Call the main function with the parsed arguments.
    try:
        main(args.physio_root_dir, args.bids_root_dir, args.force, args.invalid)
    except Exception as e:
        logging.error("An error occurred during script execution: %s", e, exc_info=True)
        logging.info("Script execution completed with errors.")
    else:
        logging.info("Script executed successfully.")