Merge pull request #11 from celprov/fix/debugging

Implementation of the workflow

physioqc/cli/run.py

@@ -0,0 +1,67 @@
+# -*- coding: utf-8 -*-
+"""Parser for physioqc."""
+
+
+import argparse
+
+from physioqc import __version__
+
+
+def _get_parser():
+    """
+    Parse command line inputs for this function.
+
+    Returns
+    -------
+    parser.parse_args() : argparse dict
+
+    Notes
+    -----
+    Default values must be updated in __call__ method from MetricsArgDict class.
+    # Argument parser follow template provided by RalphyZ.
+    # https://stackoverflow.com/a/43456577
+    """
+
+    parser = argparse.ArgumentParser(
+        description=(
+            "%(prog)s, a tool perform quality control on physiological signal "
+            "It generates a visual report and some metrics \n"
+            f"Version {__version__}"
+        ),
+        add_help=False,
+    )
+    # Required arguments
+    required = parser.add_argument_group("Required Argument")
+    required.add_argument(
+        "-in",
+        "--input-file",
+        dest="filename",
+        type=str,
+        help="Full path and name of the file containing "
+        "physiological data, with or without extension.",
+        required=True,
+    )
+
+    # Important optional arguments
+    optional = parser.add_argument_group("Optional arguments")
+    optional.add_argument(
+        "-outdir",
+        "--output-dir",
+        dest="outdir",
+        type=str,
+        help="Folder where output should be placed. Default is current folder.",
+        default=".",
+    )
+    optional.add_argument(
+        "-h", "--help", action="help", help="Show this help message and exit"
+    )
+
+    return parser
+
+
+if __name__ == "__main__":
+    raise RuntimeError(
+        "physioqc/cli/run.py should not be run directly;\n"
+        "Please `pip install` physioqc and use the "
+        "`physioqc` command"
+    )

physioqc/interfaces/interfaces.py

@@ -0,0 +1,116 @@
+import os
+
+import matplotlib.pyplot as plt
+import pandas as pd
+
+from physioqc.metrics.multimodal import peak_amplitude, peak_detection, peak_distance
+
+
+def generate_figures(figures, data, outdir):
+    """
+    Generates all the figure needed to populate the visual report.
+    Save the figures in the 'figures' folder
+
+    Parameters
+    ----------
+    figures : dictionary
+        A list of functions to run to generate all the figures.
+    data : np.array or peakdet Physio object
+        Physiological data
+    """
+    # Create the output directory if not existing
+    os.makedirs(os.path.join(outdir, "figures"), exist_ok=True)
+
+    for figure in figures:
+        # Get the plot name from the name of the function that was ran
+        plot_name = "_".join(figure.__name__.split("_")[1:])
+
+        # Run the function to generate the figure
+        if figure.__name__ == "plot_histogram":
+            # plot histogram takes as input the peak distance or peak height,
+            # while the other function take a peakdet Physio object
+
+            data = peak_detection(data)
+            # Plot histogram of peak distance
+            peak_dist = peak_distance(data)
+            fig, _ = figure(peak_dist)
+
+            plot_name = "histogram_peak_distance"
+            # TO IMPLEMENT the subject name should be automatically read when the data are loaded
+            fig.savefig(os.path.join(outdir, "figures", f"sub-01_desc-{plot_name}.svg"))
+
+            # Plot histogram of peak amplitude
+            peak_ampl = peak_amplitude(data)
+            fig, _ = figure(peak_ampl)
+
+            plot_name = "histogram_peak_distance"
+            fig.savefig(os.path.join(outdir, "figures", f"sub-01_desc-{plot_name}.svg"))
+
+        else:
+            fig, _ = figure(data)
+            # Save the figure
+            fig.savefig(os.path.join(outdir, "figures", f"sub-01_desc-{plot_name}.svg"))
+
+
+def run_metrics(metrics_dict, data):
+    """
+    Goes through the list of metrics and calls them
+
+    Parameters
+    ----------
+    metrics : dictionary
+        A dictionary that maps the physiological signal instance
+        to the function that needs to be called on it.
+    data : np.array or peakdet Physio object
+        Physiological data
+
+    Returns
+    -------
+    Dataframe :obj:`panda.dataframe`
+        A dataframe containing the value of each metric
+    """
+    metrics_df = pd.DataFrame()
+    name_list = []
+    value_list = []
+    for m in metrics_dict.items():
+        # Extract the physiological signal instance (signal, peak amplitude ...)
+        # TO IMPLEMENT : if the function m is peak_amplitude or peak_distance, you should first
+        # run peak_detection and then pass the output of peak_detection to the latter.
+        result = m[0](data)
+        for fct in m[1]:
+            # Run each metric on the physiological signal instance
+            value_list.append(fct(result))
+            # Save metric name and value
+            name_list.append(m[0].__name__ + "_" + fct.__name__)
+
+        metrics_df = pd.DataFrame(
+            list(zip(name_list, value_list)), columns=["Metric", "Value"]
+        )
+
+    return metrics_df
+
+
+def save_metrics(metrics_df, outdir, to_csv=False):
+    """
+    Save the metrics in the defined output path
+
+    Parameters
+    ----------
+    metrics_df : pd.Dataframe
+        A dataframe gathering the value of the metrics
+    outdir : string
+        Physiological data
+    to.csv : boolean
+        If true save in CSV format otherwise save to JSON
+
+    Returns
+    -------
+    Dataframe :obj:`panda.dataframe`
+        A dataframe containing the value of each metric
+    """
+    # TO IMPLEMENT : there may be a bug associated to the next line IsDirectoryError
+    os.makedirs(outdir, exist_ok=True)
+    if to_csv:
+        metrics_df.to_csv(os.path.join(outdir, "metrics.csv"), index=False)
+    else:
+        metrics_df.to_json(os.path.join(outdir, "metrics.json"))

physioqc/interfaces/visualizations.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from typing import List, Literal, Sequence, Tuple
 
 import matplotlib.pyplot as plt

physioqc/metrics/multimodal.py

@@ -1,11 +1,30 @@
 """These functions compute various non-modality dependent signal processing metrics."""
 
+import warnings
+
 import numpy as np
 import peakdet as pk
 from scipy import signal
 from scipy.misc import derivative
 
-from .utils import physio_or_numpy
+from .utils import has_peakfind_physio, physio_or_numpy
+
+
+def signal_fct(signal):
+    """
+    Wrapper that turns the object into a function for loop
+
+    Parameters
+    ----------
+    signal : np.array or peakdet Physio object
+        Physiological data
+
+    Returns
+    -------
+    signal : np.array or peakdet Physio object
+        Physiological data
+    """
+    return signal
 
 
 def std(signal):
@@ -14,7 +33,7 @@ def std(signal):
 
     Parameters
     ----------
-    signal : np.array
+    signal : np.array or peakdet Physio object
         Physiological data
 
     Returns
@@ -33,7 +52,7 @@ def mean(signal: np.array):
 
     Parameters
     ----------
-    signal : np.array
+    signal : np.array or peakdet Physio object
         Physiological data
 
     Returns
@@ -52,7 +71,7 @@ def tSNR(signal):
 
     Parameters
     ----------
-    signal : np.array
+    signal : np.array or peakdet Physio object
         Physiological data
 
     Returns
@@ -71,7 +90,7 @@ def CoV(signal):
 
     Parameters
     ----------
-    signal : np.array
+    signal : np.array or peakdet Physio object
         Physiological data
 
     Returns
@@ -80,8 +99,8 @@ def CoV(signal):
         Temporal signal to noise ratio of signal.
     """
     signal = physio_or_numpy(signal)
-    tSNR_val = np.std(signal, axis=0) / np.mean(signal, axis=0)
-    return tSNR_val
+    cov_val = np.std(signal, axis=0) / np.mean(signal, axis=0)
+    return cov_val
 
 
 def min(signal: np.array):
@@ -90,7 +109,7 @@ def min(signal: np.array):
 
     Parameters
     ----------
-    signal : np.array
+    signal : np.array or peakdet Physio object
         Physiological data
 
     Returns
@@ -109,7 +128,7 @@ def max(signal: np.array):
 
     Parameters
     ----------
-    signal : np.array
+    signal : np.array or peakdet Physio object
         Physiological data
 
     Returns
@@ -207,7 +226,10 @@ def peak_distance(ph: pk.Physio):
     np.array
         np.array of shape [npeaks, ]
     """
-    # TODO Check if peaks have been estimated.
+    if not has_peakfind_physio(ph):
+        warnings.warn("Peaks not estimated, estimating")
+        ph = peak_detection(ph)
+
     diff_peak = np.diff(ph.peaks, axis=0)
 
     return diff_peak
@@ -226,7 +248,9 @@ def peak_amplitude(ph: pk.Physio):
     np.array
         np.array of shape [npeaks - 1, ]
     """
-    # TODO Check if peaks have been estimated.
+    if not has_peakfind_physio(ph):
+        warnings.warn("Peaks not estimated, estimating")
+        ph = peak_detection(ph)
     # Assuming that peaks and troughs are aligned. Last peak has no trough.
     peak_amp = ph.data[ph.peaks[:-1]]
     trough_amp = ph.data[ph.troughs]
@@ -279,8 +303,8 @@ def energy(data, lowf=None, highf=None):
     energy_density = np.square(psd)
 
     if lowf is None or highf is None:
-        # If frequencies are not precised, compute the total power
-        idx_band = np.ones(psd.shape)
+        # If frequencies are not defined, compute the total power
+        idx_band = np.ones(psd.shape).astype(bool)
     else:
         # Define frequency band
         idx_band = np.logical_and(freqs >= lowf, freqs <= highf)
@@ -290,7 +314,7 @@ def energy(data, lowf=None, highf=None):
     return energy
 
 
-def fALFF(data, lowf, highf):
+def fALFF(data, lowf=0, highf=0.5):
     """
     Calculate the fractional amplitude of low-frequency fluctuations (fALFF).
 
@@ -310,20 +334,24 @@ def fALFF(data, lowf, highf):
     -------
     Float :obj:`numpy.ndarray`
         fALFF
+
+    Note
+    -------
+    The default value of lowf and highf were set randomly. Please update them with more meaningful value
     """
     # Extract energy in the frequency band
-    band_energy = energy(data.data, lowf=lowf, highf=highf)
+    band_energy = energy(data, lowf=lowf, highf=highf)
 
     # Extract total energy
-    total_energy = energy(data.data)
+    total_energy = energy(data)
 
     # Compute the relative energy
     rel_energy = band_energy / total_energy
 
     return rel_energy
 
 
-def freq_energy(data, thr):
+def freq_energy(data, thr=0.5):
     """
     Compute the minimum frequency with energy higher than the threshold.
 
@@ -338,30 +366,12 @@ def freq_energy(data, thr):
     -------
     Float :obj:`numpy.ndarray`
         Minimum frequency with power higher than the threshold
+
+    Note
+    ----
+    The value of the threshold has been selected randomly for now. Please update it with a more meaningful value.
     """
-    energy_nd = energy(data.data)
+    energy_nd = energy(data)
     freq = np.argmax(energy_nd > thr)
 
     return freq
-
-
-def smoothness(data):
-    """
-    Compute smoothness as the second derivative of the signal.
-
-    Parameters
-    ----------
-    args : data
-        a peakdet Physio object
-
-    Returns
-    -------
-    Float :obj:`numpy.ndarray`
-        Smoothness
-    """
-    time = np.arange(0, len(data.data) / data.fs, 1 / data.fs)
-    dx2 = np.empty(len(time))
-    for t in time:
-        dx2[t] = derivative(data.data, t, n=2)
-
-    return smoothness