Add an automatic algorithm to find group of triggers, i.e. "takes" or "runs", and attempt splitting recording automagically.

.github/workflows/auto-release.yml

@@ -10,7 +10,7 @@ on:
 
 jobs:
   auto-release:
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-22.04
     # Set skip ci to avoid loops
     if: "!contains(github.event.head_commit.message, 'ci skip') && !contains(github.event.head_commit.message, 'skip ci')"
     # Set bash as default shell for jobs
@@ -19,7 +19,7 @@ jobs:
         shell: bash
     steps:
     - name: Checkout source
-      uses: actions/checkout@v3
+      uses: actions/checkout@v4
       with:
         # Fetch all history for all branches and tags
         fetch-depth: 0

.github/workflows/python-publish.yml

@@ -10,15 +10,15 @@ on:
 jobs:
   deploy:
 
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-22.04
 
     steps:
     - name: Checkout source
-      uses: actions/checkout@v3
+      uses: actions/checkout@v4
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
-        python-version: '3.10'
+        python-version: 3.11
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
@@ -29,4 +29,5 @@ jobs:
         TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
       run: |
         python setup.py sdist bdist_wheel
+        twine check dist/*
         twine upload dist/*

phys2bids/__init__.py

@@ -4,11 +4,3 @@
 
 __version__ = get_versions()["version"]
 del get_versions
-
-from . import _version
-
-__version__ = _version.get_versions()["version"]
-
-from . import _version
-
-__version__ = _version.get_versions()["version"]

phys2bids/cli/run.py

@@ -133,6 +133,34 @@ def _get_parser():
         "or just one TR if it is consistent throughout the session.",
         default=None,
     )
+    optional.add_argument(
+        "-esttakes",
+        "--estimate_takes",
+        dest="estimate_takes",
+        action="store_true",
+        help="Run automatic algorithm to estimate clusters of triggers, e.g. the "
+        "'takes' or 'runs' of fMRI. Useful when sequences were stopped and restarted, "
+        "or when you don't know how many triggers or trs you have in each take. "
+        "This should work 95%% of the time, but might fail if the padding "
+        "between takes is too similar to the takes TRs or if the trigger are heavily "
+        "unevenly registered. Default is False.",
+        default=False,
+    )
+    optional.add_argument(
+        "-ci",
+        "--confidence-interval",
+        dest="ci",
+        # Here always as float, later it will check if the float is an integer instead.
+        type=float,
+        help="The Confidence Interval (CI) to use in the estimation of the trigger clusters. "
+        "The cluster algorithm considers triggers with duration (in samples) within this "
+        "CI as part of the same group, thus the same. If CI is an integer, it will consider "
+        "that amount of samples around the distance. If CI is a float and < 1, it will "
+        "consider that percentage of the trigger duration. CI cannot be a float > 1. "
+        "Default is 1. Change to .25 if there is a CMRR DWI sequence or when recording "
+        "sub-triggers.",
+        default=1,
+    )
     optional.add_argument(
         "-thr",
         "--threshold",

phys2bids/phys2bids.py

@@ -38,7 +38,7 @@
 from phys2bids import _version, bids, utils, viz
 from phys2bids.cli.run import _get_parser
 from phys2bids.physio_obj import BlueprintOutput
-from phys2bids.slice4phys import slice4phys
+from phys2bids.slice4phys import estimate_ntp_and_tr, slice4phys
 
 from . import __version__
 from .due import Doi, due
@@ -141,6 +141,8 @@ def phys2bids(
     chsel=None,
     num_timepoints_expected=None,
     tr=None,
+    estimate_takes=False,
+    ci=1,
     thr=None,
     pad=9,
     ch_name=[],
@@ -260,6 +262,13 @@ def phys2bids(
         from phys2bids.io import load_gep
 
         phys_in = load_gep(infile)
+    # We can't check smr due to specific version required
+    elif ftype == "smr":  # pragma: no cover
+        from phys2bids.io import load_smr
+
+        phys_in = load_smr(infile, chtrig)
+    else:
+        raise RuntimeError(f'Unsupported extension "{ftype}"')  # pragma: no cover
 
     LGR.info("Checking that units of measure are BIDS compatible")
     for index, unit in enumerate(phys_in.units):
@@ -298,6 +307,11 @@ def phys2bids(
         LGR.info("Renaming channels with given names")
         phys_in.rename_channels(ch_name)
 
+    # If requested, run the automatic detection of timepoints and groups
+
+    if estimate_takes:
+        num_timepoints_expected, tr = estimate_ntp_and_tr(phys_in, thr=None, ci=1)
+
     # Checking acquisition type via user's input
     if tr is not None and num_timepoints_expected is not None:
         #  Multi-run acquisition type section

phys2bids/physio_obj.py

@@ -465,10 +465,10 @@ def delete_at_index(self, idx):
         del self.units[idx]
 
         if self.trigger_idx == idx:
-            LGR.warning("Removing trigger channel - are you sure you are doing" "the right thing?")
+            LGR.warning("Removing trigger channel - are you sure you are doing the right thing?")
             self.trigger_idx = 0
 
-    def check_trigger_amount(self, thr=None, num_timepoints_expected=0, tr=0):
+    def check_trigger_amount(self, thr=None, num_timepoints_expected=None, tr=None):
         """
         Count trigger points and correct time offset in channel "time".
 
@@ -527,7 +527,7 @@ def check_trigger_amount(self, thr=None, num_timepoints_expected=0, tr=0):
         num_timepoints_found = np.count_nonzero(np.ediff1d(timepoints.astype(np.int8)) > 0)
         if flag == 1:
             LGR.info(
-                f"The number of timepoints according to the std_thr method "
+                f"The number of timepoints according to the average threshold method "
                 f"is {num_timepoints_found}. The computed threshold is {thr:.4f}"
             )
         else:
@@ -550,7 +550,8 @@ def check_trigger_amount(self, thr=None, num_timepoints_expected=0, tr=0):
 
             elif num_timepoints_found < num_timepoints_expected:
                 timepoints_missing = num_timepoints_expected - num_timepoints_found
-                LGR.warning(f"Found {timepoints_missing} timepoints" " less than expected!")
+                LGR.warning(f"Found {timepoints_missing} timepoints less than expected!")
+                # if tr checks for any truthy value, including not 0 nor empty things.
                 if tr:
                     LGR.warning(
                         "Correcting time offset, assuming missing "
@@ -559,14 +560,14 @@ def check_trigger_amount(self, thr=None, num_timepoints_expected=0, tr=0):
                     )
                     time_offset -= timepoints_missing * tr
                 else:
-                    LGR.warning("Can't correct time offset - you should " "specify the TR")
+                    LGR.warning("Can't correct time offset - you should specify the TR")
 
             else:
                 LGR.info("Found just the right amount of timepoints!")
 
         else:
             LGR.warning(
-                "The necessary options to find the amount of timepoints " "were not provided."
+                "The necessary options to find the amount of timepoints were not provided."
             )
         self.thr = thr
         self.time_offset = time_offset

phys2bids/slice4phys.py

@@ -8,6 +8,127 @@
 LGR = logging.getLogger(__name__)
 
 
+def estimate_ntp_and_tr(phys_in, thr=None, ci=1):
+    """
+    Find groups of trigger in a spiky signal like the trigger channel signal.
+
+    Parameters
+    ----------
+    phys_in : BlueprintInput object
+        A BlueprintInput object containing a physiological acquisition
+    thr : None, optional
+        The threshold for automatic spike detection. Default is to use the average of
+        the signal.
+    ci : int or float, optional
+        Confidence Interval (CI) to use in the estimation of the trigger clusters. The
+        cluster algorithm considers triggers with duration (in samples) within this CI
+        as part of the same group, thus the same. If CI is an integer, it will consider
+        that amount of samples. If CI is a float and < 1, it will consider that
+        percentage of the trigger duration. CI cannot be a float > 1. Default is 1.
+        Change to .25 if there is a CMRR DWI sequence or when recording sub-triggers.
+
+    Returns
+    -------
+    ntp
+        The list of number of timepoints found for each take.
+    tr
+        The list of corresponding TR, computed as average samples per group / frequency.
+
+    Raises
+    ------
+    Exception
+        If it doesn't find at least one group.
+    ValueError
+        If CI is a float above 1 or if it is not an int or a float.
+    """
+    LGR.info('Running automatic clustering of triggers to find timepoints and tr of each "take"')
+    trigger = phys_in.timeseries[phys_in.trigger_idx]
+
+    thr = np.mean(trigger) if thr is None else thr
+    timepoints = trigger > thr
+    spikes = np.flatnonzero(np.ediff1d(timepoints.astype(np.int8)) > 0)
+    interspike_interval = np.diff(spikes)
+    unique_isi, counts = np.unique(interspike_interval, return_counts=True)
+
+    # The following line is for python < 3.12. From 3.12, ci.is_integer() is enough.
+    if isinstance(ci, int) or isinstance(ci, float) and ci.is_integer():
+        upper_ci_isi = unique_isi + ci
+    elif isinstance(ci, float) and ci < 1:
+        upper_ci_isi = unique_isi * (1 + ci)
+    elif isinstance(ci, float) and ci > 1:
+        raise ValueError("Confidence intervals percentages above 1 are not supported.")
+    else:
+        raise ValueError("Confidence intervals must be either integers or floats.")
+
+    # Loop through the uniques ISI and group them within the specified CI.
+    # Also compute the average TR of the group.
+    isi_groups = {}
+    average_tr = {}
+    k = 0
+    current_group = [unique_isi[0]]
+
+    # np.unique returns sorted elements → unique_isi[0] == min(unique_isi), so THIS WORKS.
+    for n, i in enumerate(range(1, len(unique_isi))):
+        if unique_isi[i] <= upper_ci_isi[n]:
+            current_group.append(unique_isi[i])
+        else:
+            isi_groups[k] = current_group
+            average_tr[k] = np.mean(current_group) / phys_in.freq[0]
+            k += 1
+            current_group = [unique_isi[i]]
+
+    isi_groups[k] = current_group
+    average_tr[k] = np.mean(current_group) / phys_in.freq[0]
+
+    # Invert the isi_group into value per group
+    group_by_isi = {isi: group for group, isis in isi_groups.items() for isi in isis}
+
+    # Use the found groups to find the number of timepoints and assign the right TR
+    estimated_ntp = []
+    estimated_tr = []
+
+    i = 0
+    while i < interspike_interval.size - 1:
+        current_group = group_by_isi.get(interspike_interval[i])
+        for n in range(i + 1, interspike_interval.size):
+            if current_group != group_by_isi.get(interspike_interval[n]):
+                break
+        # Repeat one last time outside of for loop
+        estimated_ntp += [n - i]
+        estimated_tr += [average_tr[current_group]]
+        i = n
+
+    if len(estimated_ntp) < 1:
+        raise Exception("This should not happen. Something went very wrong.")
+    # The algorithm found n groups, the last of which has two timepoints less due to
+    # diff computations. Each real group of n>1 triggers counts one trigger less but is
+    # followed by a "fake" group of 1 trigger that is actually the interval to the next
+    # group. That does not hold if there is a real group of 1 trigger.
+    # Loop through the estiamtions to fix all that.
+    ntp = []
+    tr = []
+    i = 0
+
+    while i < len(estimated_ntp):
+        if estimated_ntp[i] == 1:
+            ntp.append(estimated_ntp[i])
+            tr.append(estimated_tr[i])
+            i += 1
+        elif i + 1 < len(estimated_ntp):
+            ntp.append(estimated_ntp[i] + estimated_ntp[i + 1])
+            tr.append(estimated_tr[i])
+            i += 2
+        else:
+            ntp.append(estimated_ntp[i] + 2)
+            tr.append(estimated_tr[i])
+            i += 1
+
+    LGR.info(
+        f"The automatic clustering found {len(ntp)} groups of triggers long: {ntp} with respective TR: {tr}"
+    )
+    return ntp, tr
+
+
 def find_takes(phys_in, ntp_list, tr_list, thr=None, padding=9):
     """
     Find takes slicing index.

phys2bids/tests/test_integration.py

@@ -280,6 +280,26 @@ def test_integration_multirun(skip_integration, multi_run_file):
     #     assert isfile(join(conversion_path, f'Test2_samefreq_TWOscans_{run}_trigger_time.png'))
     assert isfile(join(conversion_path, "Test2_samefreq_TWOscans.png"))
 
+    test_outpath = join(test_path, "out")
+    # Run for estimation
+    phys2bids(
+        filename=test_filename,
+        indir=test_path,
+        outdir=test_outpath,
+        chtrig=test_chtrig,
+        ch_name=["An", "Advantage"],
+    )
+
+    # Check that files are generated in outdir
+    base_filename = "Test2_samefreq_TWOscans_"
+    for suffix in [".json", ".tsv.gz"]:
+        for run in ["01", "02"]:
+            assert isfile(join(test_path, f"{base_filename}{run}{suffix}"))
+
+    # Check that files are generated in conversion_path
+    for run in ["01", "02"]:
+        assert isfile(join(conversion_path, f"Test2_samefreq_TWOscans_{run}.log"))
+
 
 def test_integration_gep_onefile(skip_integration, ge_one_gep_file):
     """

phys2bids/tests/test_phys2bids.py

@@ -4,6 +4,7 @@
 
 import json
 import os
+from shutil import copyfile
 
 from pytest import raises
 
@@ -78,3 +79,6 @@ def test_raise_exception(samefreq_full_acq_file):
             outdir=test_path,
         )
     assert "stop now" in str(errorinfo.value)
+
+    test_newfilename = test_filename + ".owowow"
+    copyfile(samefreq_full_acq_file, test_newfilename)

phys2bids/tests/test_physio_obj.py

@@ -118,12 +118,15 @@ def test_BlueprintInput():
     assert blueprint_in.ch_amount == num_channnels - 1
 
     # Tests check_trigger_amount
-    blueprint_in.check_trigger_amount(thr=0.9, num_timepoints_expected=1)
+    blueprint_in.check_trigger_amount(thr=0.9, num_timepoints_expected=1, tr=1)
     assert blueprint_in.num_timepoints_found == 1
     assert blueprint_in.time_offset == 1
     test_offset_time = test_time - 1
     assert np.array_equal(blueprint_in.timeseries[0], test_offset_time)
 
+    blueprint_in.check_trigger_amount(thr=0.9, num_timepoints_expected=1)
+    blueprint_in.check_trigger_amount(thr=0.9)
+
     # Tests delete_at_index with trigger channel
     blueprint_in.delete_at_index(test_chtrig)
     assert blueprint_in.ch_amount == num_channnels - 2