Repository: physiopy/phys2bids

Issue:
Chopping files semi-automatically 
A lot of researchers will have one physiological file for the whole scan session, which will need to be chopped into separate files corresponding to each scan to be analyzed. Offering a way of semi-automatically chopping these files, and allowing phys2bids to output multiple physiological files from the same input file, would be useful. 

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index 6044203..48a6f50 100644
--- a/.gitignore
+++ b/.gitignore
@@ -126,3 +126,4 @@ dmypy.json
 
 docs/generated/
 .vscode/
+phys2bids/tests/data/*
diff --git a/docs/howto.rst b/docs/howto.rst
index a08ee7c..526dbc9 100644
--- a/docs/howto.rst
+++ b/docs/howto.rst
@@ -279,6 +279,74 @@ By looking at this figure, we can work out that we need a smaller threshold in o
     Tip: Time 0 is the time of first trigger
     ------------------------------------------------
 
+Splitting your input file into multiple run output files
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If your file contains more than one (f)MRI acquisition (or runs), you can provide multiple values to ``-ntp`` in order to get multiple ``.tsv.gz`` outputs. If the TR of the entire session is consistent (i.e. all your acquisitions have the same TR), then you can specify one value after ``-tr``.
+
+By specifying the number of timepoints in each acquisition, ``phys2bids`` will recursively cut the input file by detecting the first trigger of the entire session and the ones after the number of timepoints you specified.
+
+.. code-block:: shell
+
+    phys2bids -in two_scans_samefreq_all.txt -chtrig 2 -ntp 536 398 -tr 1.2 -thr 2
+
+Now, instead of counting the trigger timepoints once, ``physbids`` will check the trigger channel recursively with all the values listed in ``-ntp``. The logger will inform you about the number of timepoints left at each iteration.
+
+.. code-block:: shell
+
+    INFO:phys2bids.physio_obj:Counting trigger points
+    INFO:phys2bids.physio_obj:The trigger is in channel 2
+    INFO:phys2bids.physio_obj:The number of timepoints found with the manual threshold of 2.0000 is 934
+    INFO:phys2bids.physio_obj:Checking number of timepoints
+    INFO:phys2bids.physio_obj:Found just the right amount of timepoints!
+    WARNING:phys2bids.slice4phys:
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    phys2bids will split the input file according to the given -tr and -ntp arguments
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    INFO:phys2bids.physio_obj:Counting trigger points
+    INFO:phys2bids.physio_obj:The trigger is in channel 2
+    INFO:phys2bids.physio_obj:The number of timepoints found with the manual threshold of 2.0000 is 934
+    INFO:phys2bids.physio_obj:Checking number of timepoints
+    WARNING:phys2bids.physio_obj:Found 398 timepoints more than expected!
+    Assuming extra timepoints are at the end (try again with a more liberal thr)
+    INFO:phys2bids.slice4phys:
+    --------------------------------------------------------------
+    Slicing between 0.0 seconds and 961.381 seconds
+    --------------------------------------------------------------
+    INFO:phys2bids.physio_obj:Counting trigger points
+    INFO:phys2bids.physio_obj:The trigger is in channel 2
+    INFO:phys2bids.physio_obj:The number of timepoints found with the manual threshold of 2.0000 is 400
+    INFO:phys2bids.physio_obj:Checking number of timepoints
+    WARNING:phys2bids.physio_obj:Found 2 timepoints more than expected!
+    Assuming extra timepoints are at the end (try again with a more liberal thr)
+    INFO:phys2bids.slice4phys:
+    --------------------------------------------------------------
+    Slicing between 952.381 seconds and 1817.96 seconds
+    --------------------------------------------------------------
+    INFO:phys2bids.viz:Plot trigger
+    INFO:phys2bids.viz:Plot trigger
+    INFO:phys2bids.phys2bids:Preparing 2 output files.
+    INFO:phys2bids.phys2bids:Exporting files for run 1 freq 1000.0
+
+The logger also notifies you about the slicing points used (the first always being from the beginning of session, until the specified number of timepoints after the first trigger). The user can also check the resulting slice by looking at the plot of the trigger channel for each run. Each slice is adjusted with a padding after the last trigger. Such padding can be specified while calling ``phys2bids`` with ``-pad``. If nothing is specified, the default value of 9 seconds will be used. This padding is also applied at the beginning (before the first trigger of the run) of the 2nd to last run.
+
+What if I have multiple acquisition types ?
+*******************************************
+
+The user can also benefit from this utility when dealing with multiple **acquisition types** such as different functional scans with different TRs. Like ``-ntp``, ``-tr`` can take multiple values. **Though, if more than one value is specified, they require the same amount of values**. The idea is simple: if you only have one acquisition type, the one ``-tr`` input you gave will be broadcast through all runs, but if you have different acquisition types, you have to list all of them in order.
+
+.. warning::
+    There are currently no ``multi-run tutorial files`` available along with the package (under ``phys2bids/tests/data``). Although, you can visit `phys2bids OSF <https://osf.io/3txqr/files/>`_ storage to access a LabChart physiological recording with multiple fMRI acquisitions. Find it under ``labchart/chicago``.
+
+.. note::
+    **Why would I have more than one fMRI acquisition in the physiological recording?**
+
+    The idea is to reduce human error and have a good padding around your fMRI scan!
+
+    Synchronization between start of both fMRI and physiological acquisitions can be difficult, so it is safer to have as few physiological recordings as possible, with multiple imaging sequences.
+
+    Moreover, if you want to correct for recording/physiological delays, you will love *that bit of recorded information* around your fMRI scan!
+
 Generating outputs in BIDs format
 #################################
 
diff --git a/phys2bids/bids.py b/phys2bids/bids.py
index f529693..2cfc80f 100644
--- a/phys2bids/bids.py
+++ b/phys2bids/bids.py
@@ -1,5 +1,5 @@
-import os
 import logging
+import os
 from csv import reader
 from pathlib import Path
 
@@ -25,6 +25,9 @@ UNIT_ALIASES = {
                 # ampere: electric current
                 'ampere': 'A', 'amp': 'A', 'amps': 'A',
                 # second: time and hertzs: frequency
+                # siemens: electric conductance (e.g. EDA)
+                'siemens': 'S',
+                # second: time and hertzs
                 '1/hz': 's', '1/hertz': 's', 'hz': 'Hz',
                 '1/s': 'Hz', '1/second': 'Hz', '1/seconds': 'Hz',
                 '1/sec': 'Hz', '1/secs': 'Hz', 'hertz': 'Hz',
@@ -96,7 +99,7 @@ def bidsify_units(orig_unit):
     return orig_unit
 
 
-def use_heuristic(heur_file, sub, ses, filename, outdir, record_label=''):
+def use_heuristic(heur_file, sub, ses, filename, outdir, run='', record_label=''):
     """
     Import and use the heuristic specified by the user to rename the file.
 
@@ -130,7 +133,7 @@ def use_heuristic(heur_file, sub, ses, filename, outdir, record_label=''):
 
     # Initialise a dictionary of bids_keys that has already "recording"
     bids_keys = {'sub': '', 'ses': '', 'task': '', 'acq': '', 'ce': '',
-                 'dir': '', 'rec': '', 'run': '', 'recording': record_label}
+                 'dir': '', 'rec': '', 'run': run, 'recording': record_label}
 
     # Start filling bids_keys dictionary and path with subject and session
     if sub.startswith('sub-'):
@@ -150,7 +153,7 @@ def use_heuristic(heur_file, sub, ses, filename, outdir, record_label=''):
 
     # Load heuristic and use it to fill dictionary
     heur = utils.load_heuristic(heur_file)
-    bids_keys.update(heur.heur(Path(filename).stem))
+    bids_keys.update(heur.heur(Path(filename).stem, run))
 
     # If bids_keys['task'] is still empty, stop the program
     if not bids_keys['task']:
diff --git a/phys2bids/cli/run.py b/phys2bids/cli/run.py
index adf181f..497eb48 100644
--- a/phys2bids/cli/run.py
+++ b/phys2bids/cli/run.py
@@ -1,7 +1,6 @@
 # -*- coding: utf-8 -*-
-"""
-Parser for phys2bids
-"""
+"""Parser for phys2bids."""
+
 
 import argparse
 
@@ -10,7 +9,7 @@ from phys2bids import __version__
 
 def _get_parser():
     """
-    Parses command line inputs for this function
+    Parse command line inputs for this function.
 
     Returns
     -------
@@ -27,8 +26,8 @@ def _get_parser():
     required.add_argument('-in', '--input-file',
                           dest='filename',
                           type=str,
-                          help='The name of the file containing physiological data, with or '
-                               'without extension.',
+                          help='The name of the file containing physiological '
+                               'data, with or without extension.',
                           required=True)
     optional.add_argument('-info', '--info',
                           dest='info',
@@ -61,11 +60,6 @@ def _get_parser():
                                'in the current folder. Edit the heur_ex.py file in '
                                'heuristics folder.',
                           default=None)
-    # optional.add_argument('-hdir', '--heur-dir',
-    #                       dest='heurdir',
-    #                       type=str,
-    #                       help='Folder containing heuristic file.',
-    #                       default='.')
     optional.add_argument('-sub', '--subject',
                           dest='sub',
                           type=str,
@@ -94,26 +88,39 @@ def _get_parser():
                           default=None)
     optional.add_argument('-ntp', '--numtps',
                           dest='num_timepoints_expected',
+                          nargs='*',
                           type=int,
-                          help='Number of expected timepoints (TRs). '
-                               'Default is 0. Note: the estimation of when the '
-                               'neuroimaging acquisition started cannot take place '
-                               'with this default.',
-                          default=0)
+                          help='Number of expected trigger timepoints (TRs). '
+                               'Default is None. Note: the estimation of beggining of '
+                               'neuroimaging acquisition cannot take place with this default. '
+                               'If you\'re running phys2bids on a multi-run recording, '
+                               'give a list of each expected ntp for each run.',
+                          default=None)
     optional.add_argument('-tr', '--tr',
                           dest='tr',
+                          nargs='*',
                           type=float,
                           help='TR of sequence in seconds. '
-                               'Default is 0 second.',
-                          default=0)
+                               'If you\'re running phys2bids on a multi-run recording, '
+                               'you can give a list of each expected ntp for each run, '
+                               'or just one TR if it is consistent throughout the session.',
+                          default=None)
     optional.add_argument('-thr', '--threshold',
                           dest='thr',
                           type=float,
                           help='Threshold to use for trigger detection. '
-                               'If "ntp" and "TR" are specified, phys2bids automatically computes '
-                               'a threshold to detect the triggers. Use this parameter to set it '
-                               'manually',
-                               default=None)
+                               'If "ntp" and "TR" are specified, phys2bids '
+                               'automatically computes a threshold to detect '
+                               'the triggers. Use this parameter to set it manually. '
+                               'This parameter is necessary for multi-run recordings. ',
+                          default=None)
+    optional.add_argument('-pad', '--padding',
+                          dest='pad',
+                          type=float,
+                          help='Padding in seconds used around a single run '
+                               'when separating multi-run session files. '
+                               'Default is 9 seconds.',
+                          default=9)
     optional.add_argument('-chnames', '--channel-names',
                           dest='ch_name',
                           nargs='*',
diff --git a/phys2bids/heuristics/heur_euskalibur.py b/phys2bids/heuristics/heur_euskalibur.py
index 8c0e1d5..75ae8c6 100644
--- a/phys2bids/heuristics/heur_euskalibur.py
+++ b/phys2bids/heuristics/heur_euskalibur.py
@@ -1,7 +1,7 @@
 import fnmatch
 
 
-def heur(physinfo):
+def heur(physinfo, run=''):
     """
     Set of if .. elif statements to fill BIDS names.
 
diff --git a/phys2bids/heuristics/heur_tutorial.py b/phys2bids/heuristics/heur_tutorial.py
index eebc082..694c633 100644
--- a/phys2bids/heuristics/heur_tutorial.py
+++ b/phys2bids/heuristics/heur_tutorial.py
@@ -1,7 +1,7 @@
 import fnmatch
 
 
-def heur(physinfo):
+def heur(physinfo, run=''):
     """
     Set of if .. elif statements to fill BIDS names.
 
diff --git a/phys2bids/interfaces/acq.py b/phys2bids/interfaces/acq.py
index 7ea3e4c..fd61398 100644
--- a/phys2bids/interfaces/acq.py
+++ b/phys2bids/interfaces/acq.py
@@ -1,9 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
-"""
-phys2bids interface for acqknowledge files.
-"""
+"""phys2bids interface for acqknowledge files."""
 import logging
 import warnings
 
diff --git a/phys2bids/interfaces/txt.py b/phys2bids/interfaces/txt.py
index f0ed4f1..31cc644 100644
--- a/phys2bids/interfaces/txt.py
+++ b/phys2bids/interfaces/txt.py
@@ -1,9 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
-"""
-phys2bids interface for txt files.
-"""
+"""phys2bids interface for txt files."""
 
 import logging
 from collections import Counter
diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py
index e4781dc..1a4e316 100644
--- a/phys2bids/phys2bids.py
+++ b/phys2bids/phys2bids.py
@@ -33,11 +33,12 @@ import sys
 from copy import deepcopy
 from shutil import copy as cp
 
-from numpy import savetxt
+import numpy as np
 
 from phys2bids import utils, viz, _version, bids
 from phys2bids.cli.run import _get_parser
 from phys2bids.physio_obj import BlueprintOutput
+from phys2bids.slice4phys import slice4phys
 
 from . import __version__
 from .due import due, Doi
@@ -126,10 +127,10 @@ def print_json(outfile, samp_freq, time_offset, ch_name):
     description='The BIDS specification',
     cite_module=True)
 def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
-              sub=None, ses=None, chtrig=0, chsel=None, num_timepoints_expected=0,
-              tr=1, thr=None, ch_name=[], yml='', debug=False, quiet=False):
+              sub=None, ses=None, chtrig=0, chsel=None, num_timepoints_expected=None,
+              tr=None, thr=None, pad=9, ch_name=[], yml='', debug=False, quiet=False):
     """
-    Main workflow of phys2bids.
+    Run main workflow of phys2bids.
 
     Runs the parser, does some checks on input, then imports
     the right interface file to read the input. If only info is required,
@@ -141,18 +142,15 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
     ------
     NotImplementedError
         If the file extension is not supported yet.
-
     """
     # Check options to make them internally coherent pt. I
     # #!# This can probably be done while parsing?
     outdir = utils.check_input_dir(outdir)
     utils.path_exists_or_make_it(outdir)
     utils.path_exists_or_make_it(os.path.join(outdir, 'code'))
-    conversion_path = os.path.join(outdir, 'code/conversion')
+    conversion_path = os.path.join(outdir, 'code', 'conversion')
     utils.path_exists_or_make_it(conversion_path)
-    # generate extra path
-    extra_dir = os.path.join(outdir, 'bids_ignore')
-    utils.path_exists_or_make_it(extra_dir)
+
     # Create logfile name
     basename = 'phys2bids_'
     extension = 'tsv'
@@ -203,6 +201,18 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
     infile = os.path.join(indir, filename)
     utils.check_file_exists(infile)
 
+    if isinstance(num_timepoints_expected, int):
+        num_timepoints_expected = [num_timepoints_expected]
+    if isinstance(tr, (int, float)):
+        tr = [tr]
+
+    if tr is not None and num_timepoints_expected is not None:
+        # If tr and ntp were specified, check that tr is either length one or ntp.
+        if len(num_timepoints_expected) != len(tr) and len(tr) != 1:
+            raise Exception('Number of sequence types listed with TR '
+                            'doesn\'t match expected number of runs in '
+                            'the session')
+
     # Read file!
     if ftype == 'acq':
         from phys2bids.interfaces.acq import populate_phys_input
@@ -214,8 +224,11 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
 
     LGR.info(f'Reading the file {infile}')
     phys_in = populate_phys_input(infile, chtrig)
+
+    LGR.info('Checking that units of measure are BIDS compatible')
     for index, unit in enumerate(phys_in.units):
         phys_in.units[index] = bids.bidsify_units(unit)
+
     LGR.info('Reading infos')
     phys_in.print_info(filename)
     # #!# Here the function viz.plot_channel should be called
@@ -225,25 +238,6 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
     if info:
         return
 
-    # Create trigger plot. If possible, to have multiple outputs in the same
-    # place, adds sub and ses label.
-    if tr != 0 and num_timepoints_expected != 0:
-        # Run analysis on trigger channel to get first timepoint and the time offset.
-        # #!# Get option of no trigger! (which is wrong practice or Respiract)
-        phys_in.check_trigger_amount(thr, num_timepoints_expected, tr)
-        LGR.info('Plot trigger')
-        plot_path = os.path.join(conversion_path,
-                                 os.path.splitext(os.path.basename(filename))[0])
-        if sub:
-            plot_path += f'_sub-{sub}'
-        if ses:
-            plot_path += f'_ses-{ses}'
-        viz.plot_trigger(phys_in.timeseries[0], phys_in.timeseries[chtrig],
-                         plot_path, tr, phys_in.thr, num_timepoints_expected, filename)
-    else:
-        LGR.warning('Skipping trigger pulse count. If you want to run it, '
-                    'call phys2bids using "-ntp" and "-tr" arguments')
-
     # The next few lines remove the undesired channels from phys_in.
     if chsel:
         LGR.info('Dropping unselected channels')
@@ -256,40 +250,88 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
         LGR.info('Renaming channels with given names')
         phys_in.rename_channels(ch_name)
 
+    # Checking acquisition type via user's input
+    if tr is not None and num_timepoints_expected is not None:
+
+        #  Multi-run acquisition type section
+        #  Check list length, more than 1 means multi-run
+        if len(num_timepoints_expected) > 1:
+            # if multi-run of same sequence type, pad list with ones
+            # and multiply array with user's input
+            if len(tr) == 1:
+                tr = np.ones(len(num_timepoints_expected)) * tr[0]
+
+            # Sum of values in ntp_list should be equivalent to num_timepoints_found
+            phys_in.check_trigger_amount(thr=thr,
+                                         num_timepoints_expected=sum(num_timepoints_expected),
+                                         tr=1)
+
+            # Check that sum of tp expected is equivalent to num_timepoints_found,
+            # if it passes call slice4phys
+            if phys_in.num_timepoints_found != sum(num_timepoints_expected):
+                raise Exception('The number of triggers found is different '
+                                'than expected. Better stop now than break '
+                                'something.')
+
+            # slice the recording based on user's entries
+            # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY
+            phys_in = slice4phys(phys_in, num_timepoints_expected, tr,
+                                 phys_in.thr, pad)
+            # returns a dictionary in the form {run_idx: phys_in[startpoint, endpoint]}
+
+            # save a figure for each run | give the right acquisition parameters for runs
+            fileprefix = os.path.join(conversion_path,
+                                      os.path.splitext(os.path.basename(filename))[0])
+            for i, run in enumerate(phys_in.keys()):
+                plot_fileprefix = f'{fileprefix}_{run}'
+                viz.export_trigger_plot(phys_in[run], chtrig, plot_fileprefix, tr[i],
+                                        num_timepoints_expected[i], filename,
+                                        sub, ses)
+
+        # Single run acquisition type, or : nothing to split workflow
+        else:
+            # Run analysis on trigger channel to get first timepoint
+            # and the time offset.
+            phys_in.check_trigger_amount(thr, num_timepoints_expected[0], tr[0])
+            # save a figure of the trigger
+            fileprefix = os.path.join(conversion_path,
+                                      os.path.splitext(os.path.basename(filename))[0])
+            viz.export_trigger_plot(phys_in, chtrig, fileprefix, tr[0],
+                                    num_timepoints_expected[0], filename,
+                                    sub, ses)
+
+            # Reassign phys_in as dictionary
+            # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY
+            phys_in = {1: phys_in}
+
+    else:
+        LGR.warning('Skipping trigger pulse count. If you want to run it, '
+                    'call phys2bids using both "-ntp" and "-tr" arguments')
+        # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY
+        phys_in = {1: phys_in}
+
     # The next few lines create a dictionary of different BlueprintInput
-    # objects, one for each unique frequency in phys_in
-    uniq_freq_list = set(phys_in.freq)
-    output_amount = len(uniq_freq_list)
-    if output_amount > 1:
-        LGR.warning(f'Found {output_amount} different frequencies in input!')
-
-    LGR.info(f'Preparing {output_amount} output files.')
-    phys_out = {}  # create phys_out dict that will have a
-    # blueprint object per frequency
-    # for each different frequency
-    for uniq_freq in uniq_freq_list:
-        # copy the phys_in object to the new dict entry
-        phys_out[uniq_freq] = deepcopy(phys_in)
-        # this counter will take into account how many channels are eliminated
-        count = 0
-        # for each channel in the original phys_in object
-        # take the frequency
-        for idx, i in enumerate(phys_in.freq):
-            # if that frequency is different than the frequency of the phys_obj entry
-            if i != uniq_freq:
-                # eliminate that channel from the dict since we only want channels
-                # with the same frequency
-                phys_out[uniq_freq].delete_at_index(idx - count)
-                # take into acount the elimination so in the next eliminated channel we
-                # eliminate correctly
-                count += 1
-        # Also create a BlueprintOutput object for each unique frequency found.
-        # Populate it with the corresponding blueprint input and replace it
-        # in the dictionary.
-        phys_out[uniq_freq] = BlueprintOutput.init_from_blueprint(phys_out[uniq_freq])
-
-    if heur_file and sub:
+    # objects, one for each unique frequency for each run in phys_in
+    # they also save the amount of runs and unique frequencies
+    run_amount = len(phys_in)
+    uniq_freq_list = set(phys_in[1].freq)
+    freq_amount = len(uniq_freq_list)
+    if freq_amount > 1:
+        LGR.info(f'Found {freq_amount} different frequencies in input!')
+
+    if run_amount > 1:
+        LGR.info(f'Found {run_amount} different scans in input!')
+
+    LGR.info(f'Preparing {freq_amount*run_amount} output files.')
+    # Create phys_out dict that will have a blueprint object for each different frequency
+    phys_out = {}
+
+    if heur_file is not None and sub is not None:
         LGR.info(f'Preparing BIDS output using {heur_file}')
+        # If heuristics are used, init a dict of arguments to pass to use_heuristic
+        heur_args = {'heur_file': heur_file, 'sub': sub, 'ses': ses,
+                     'filename': filename, 'outdir': outdir, 'run': '',
+                     'record_label': ''}
         # Generate participants.tsv file if it doesn't exist already.
         # Update the file if the subject is not in the file.
         # Do not update if the subject is already in the file.
@@ -300,39 +342,85 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
         bids.readme_file(outdir)
         cp(heur_file, os.path.join(conversion_path,
            os.path.splitext(os.path.basename(heur_file))[0] + '.py'))
-    elif heur_file and not sub:
+    elif heur_file is not None and sub is None:
         LGR.warning('While "-heur" was specified, option "-sub" was not.\n'
                     'Skipping BIDS formatting.')
 
-    # Preparing output parameters: name and folder.
-    for uniq_freq in uniq_freq_list:
-        # If possible, prepare bids renaming.
-        if heur_file and sub:
-            if output_amount > 1:
-                # Add "recording-freq" to filename if more than one freq
-                outfile = bids.use_heuristic(heur_file, sub, ses, filename,
-                                             outdir, uniq_freq)
-            else:
-                outfile = bids.use_heuristic(heur_file, sub, ses, filename, outdir)
+    # Export a (set of) phys_out for each element in phys_in
+    # run keys start from 1 (human friendly)
+    for run in phys_in.keys():
+        for uniq_freq in uniq_freq_list:
+            # Initialise the key for the (possibly huge amount of) dictionary entries
+            key = f'{run}_{uniq_freq}'
+            # copy the phys_in object to the new dict entry
+            phys_out[key] = deepcopy(phys_in[run])
+            # this counter will take into account how many channels are eliminated
+            count = 0
+            # for each channel in the original phys_in object
+            # take the frequency
+            for idx, i in enumerate(phys_in[run].freq):
+                # if that frequency is different than the frequency of the phys_obj entry
+                if i != uniq_freq:
+                    # eliminate that channel from the dict since we only want channels
+                    # with the same frequency
+                    phys_out[key].delete_at_index(idx - count)
+                    # take into acount the elimination so in the next eliminated channel we
+                    # eliminate correctly
+                    count += 1
+            # Also create a BlueprintOutput object for each unique frequency found.
+            # Populate it with the corresponding blueprint input and replace it
+            # in the dictionary.
+            phys_out[key] = BlueprintOutput.init_from_blueprint(phys_out[key])
+
+        # Preparing output parameters: name and folder.
+        for uniq_freq in uniq_freq_list:
+            key = f'{run}_{uniq_freq}'
+            # If possible, prepare bids renaming.
+            if heur_file is not None and sub is not None:
+                # Add run info to heur_args if more than one run is present
+                if run_amount > 1:
+                    heur_args['run'] = f'{run:02d}'
+
+                # Append "recording-freq" to filename if more than one freq
+                if freq_amount > 1:
+                    heur_args['record_label'] = f'{uniq_freq:.0f}Hz'
+
+                phys_out[key].filename = bids.use_heuristic(**heur_args)
+
+                # If any filename exists already because of multirun, append labels
+                # But warn about the non-validity of this BIDS-like name.
+                if run_amount > 1:
+                    if any([phys.filename == phys_out[key].filename
+                           for phys in phys_out.values()]):
+                        phys_out[key].filename = (f'{phys_out[key].filename}'
+                                                  f'_take-{run}')
+                        LGR.warning('Identified multiple outputs with the same name.\n'
+                                    'Appending fake label to avoid overwriting.\n'
+                                    '!!! ATTENTION !!! the output is not BIDS compliant.\n'
+                                    'Please check heuristics to solve the problem.')
 
-        else:
-            outfile = os.path.join(outdir,
-                                   os.path.splitext(os.path.basename(filename))[0])
-            if output_amount > 1:
+            else:
+                phys_out[key].filename = os.path.join(outdir,
+                                                      os.path.splitext(os.path.basename(filename)
+                                                                       )[0])
+                # Append "run" to filename if more than one run
+                if run_amount > 1:
+                    phys_out[key].filename = f'{phys_out[key].filename}_{run:02d}'
                 # Append "freq" to filename if more than one freq
-                outfile = f'{outfile}_{uniq_freq}'
-
-        LGR.info(f'Exporting files for freq {uniq_freq}')
-        savetxt(outfile + '.tsv.gz', phys_out[uniq_freq].timeseries,
-                fmt='%.8e', delimiter='\t')
-        print_json(outfile, phys_out[uniq_freq].freq,
-                   phys_out[uniq_freq].start_time,
-                   phys_out[uniq_freq].ch_name)
-        print_summary(filename, num_timepoints_expected,
-                      phys_in.num_timepoints_found, uniq_freq,
-                      phys_out[uniq_freq].start_time,
-                      os.path.join(conversion_path,
-                                   os.path.splitext(os.path.basename(outfile))[0]))
+                if freq_amount > 1:
+                    phys_out[key].filename = f'{phys_out[key].filename}_{uniq_freq:.0f}Hz'
+
+            LGR.info(f'Exporting files for run {run} freq {uniq_freq}')
+            np.savetxt(phys_out[key].filename + '.tsv.gz',
+                       phys_out[key].timeseries, fmt='%.8e', delimiter='\t')
+            print_json(phys_out[key].filename, phys_out[key].freq,
+                       phys_out[key].start_time, phys_out[key].ch_name)
+            print_summary(filename, num_timepoints_expected,
+                          phys_in[run].num_timepoints_found, uniq_freq,
+                          phys_out[key].start_time,
+                          os.path.join(conversion_path,
+                                       os.path.splitext(os.path.basename(phys_out[key].filename)
+                                                        )[0]))
 
 
 def _main(argv=None):
diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py
index 3eb1ce1..24fb923 100644
--- a/phys2bids/physio_obj.py
+++ b/phys2bids/physio_obj.py
@@ -298,8 +298,15 @@ class BlueprintInput():
 
             idx = slice(idx, idx + 1)
 
+        # If idx.start or stop are None, make them 0 or trigger length.
+        if idx.start is None:
+            idx = slice(0, idx.stop)
+        if idx.stop is None:
+            idx = slice(idx.start, trigger_length)
+
+        # Check that the indexes are not out of bounds
         if idx.start >= trigger_length or idx.stop > trigger_length:
-            raise IndexError(f'slice ({idx.start}, {idx.stop}) is out of '
+            raise IndexError(f'Slice ({idx.start}, {idx.stop}) is out of '
                              f'bounds for channel {self.trigger_idx} '
                              f'with size {trigger_length}')
 
@@ -549,6 +556,9 @@ class BlueprintOutput():
     start_time : float
         Starting time of acquisition (equivalent to first TR,
         or to the opposite sign of the time offset).
+    filename : string
+        Filename the object will be saved with. Init as empty string
+
 
     Methods
     -------
@@ -564,13 +574,14 @@ class BlueprintOutput():
         method to populate from input blueprint instead of init
     """
 
-    def __init__(self, timeseries, freq, ch_name, units, start_time):
+    def __init__(self, timeseries, freq, ch_name, units, start_time, filename=''):
         """Initialise BlueprintOutput (see class docstring)."""
         self.timeseries = is_valid(timeseries, np.ndarray)
         self.freq = is_valid(freq, (int, float))
         self.ch_name = has_size(ch_name, self.ch_amount, 'unknown')
         self.units = has_size(units, self.ch_amount, '[]')
         self.start_time = start_time
+        self.filename = is_valid(filename, str)
 
     @property
     def ch_amount(self):
diff --git a/phys2bids/slice4phys.py b/phys2bids/slice4phys.py
new file mode 100644
index 0000000..573a7d1
--- /dev/null
+++ b/phys2bids/slice4phys.py
@@ -0,0 +1,155 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import logging
+from copy import deepcopy
+
+import numpy as np
+
+LGR = logging.getLogger(__name__)
+
+
+def find_runs(phys_in, ntp_list, tr_list, thr=None, padding=9):
+    """
+    Find runs slicing index.
+
+    Returns dictionary key for each run in BlueprintInput object based on
+    user's entries. Each key has a tuple of 4 elements, 2 expressing the
+    timestamps of run in nb of samples. Timestamps are the index of first and
+    last triggers of a run, adjusted with padding. run_start and run_end
+    indexes refer to the samples contained in the whole session.
+    First trigger time offset and nb of triggers contained in the run are also indicated.
+
+    Parameters
+    ---------
+    phys_in: BlueprintInput object
+        Object returned by BlueprintInput class
+    ntp_list: list
+        a list of integers given by the user as `ntp` input
+    tr_list: list
+        a list of float given by the user as `tr` input
+    thr: int
+        inherit threshold for detection of trigger given by user
+    padding: int
+        extra time at beginning and end of timeseries, expressed in seconds (s)
+        Default: 9
+
+    Returns
+    --------
+    run_timestamps: dictionary
+        Containing tuples of run start and end indexes for each run, based on
+        trigger channels. It also contains run attributes: time offset from
+        session beggining, and nb of triggers in the form of
+        run_timestamps{1:(start, end, time offset, nb of triggers),
+                       2:(...), ... }
+    """
+    # Initialize dictionaries to save  run timestamps and phys_in attributes
+    run_timestamps = {}
+
+    # Express the padding in samples equivalent
+    padding = padding * phys_in.freq[0]
+
+    # enumerate user input  num_timepoints_expected
+    for run_idx, run_tps in enumerate(ntp_list):
+
+        # correct time offset for this iteration's object
+        phys_in.check_trigger_amount(thr=thr, num_timepoints_expected=run_tps,
+                                     tr=tr_list[run_idx])
+        # If it's the very first run, start the run at sample 0,
+        # otherwise start is first trigger (adjust with padding later)
+        if run_idx == 0:
+            run_start = 0
+        else:
+            run_start = int(np.where(np.isclose(phys_in.timeseries[0], 0))[0])
+
+        # Defining end of acquisition
+        # run length in seconds
+        end_sec = run_tps * tr_list[run_idx]
+
+        # define index of the run's last trigger + padding (HAS TO BE INT type)
+        # pick first value of time array that is over specified run length
+        # where returns list of values over end_sec and its dtype, choose [list][first value]
+        run_end = int(np.where(phys_in.timeseries[0] > end_sec)[0][0] + padding)
+        update = int(run_end - padding + 1)
+
+        # if the padding is too much for the remaining timeseries length
+        # then the padding stops at the end of recording
+        if phys_in.timeseries[0].shape[0] < run_end:
+            run_end = phys_in.timeseries[0].shape[0]
+
+        # Adjust timestamps with previous end_index
+        # Except if it's the first run
+        if run_idx > 0:
+            previous_end_index = run_timestamps[run_idx][1]
+            # adjust time_offset to keep original timing information
+            phys_in.time_offset = phys_in.time_offset + run_timestamps[run_idx][2]
+            # update run_start, removing 2 paddings (one for this run, one for the previous)
+            run_start = int(run_start + previous_end_index - 2 * padding)
+            # update run_end, removing the padding of the previous end
+            run_end = int(run_end + previous_end_index - padding)
+
+        # Save *start* and *end_index* in dictionary along with *time_offset* and *ntp found*
+        # dict key must be readable by human
+        # LGRinfo
+        LGR.info('\n--------------------------------------------------------------\n'
+                 f'Slicing between {(run_start/phys_in.freq[phys_in.trigger_idx])} seconds and '
+                 f'{run_end/phys_in.freq[phys_in.trigger_idx]} seconds\n'
+                 '--------------------------------------------------------------')
+
+        run_timestamps[run_idx + 1] = (run_start, run_end,
+                                       phys_in.time_offset,
+                                       phys_in.num_timepoints_found)
+
+        # update the object so that next iteration will look for the first trigger
+        # after previous run's last trigger. maybe padding extends to next run
+        phys_in = deepcopy(phys_in[update:-1])
+
+    return run_timestamps
+
+
+def slice4phys(phys_in, ntp_list, tr_list, thr, padding=9):
+    """
+    Slice runs for phys2bids.
+
+    Parameters
+    ---------
+    phys_in: BlueprintInput object
+        Object returned by BlueprintInput class
+    ntp_list: list
+        a list of integers given by the user as `ntp` input
+        Default: [0, ]
+    tr_list: list
+        a list of float given by the user as `tr` input
+        Default: [1,]
+    padding: int
+        extra time at beginning and end of timeseries, expressed in seconds (s)
+        Default: 9
+
+    Returns
+    --------
+    phys_in_slices: dict
+        keys start by `run 1` until last (`run n`).
+        items are slices of BlueprintInput objects based on run attributes returned by
+        internal function (`slice4phys` takes the same arguments as `find_runs`)
+    """
+    phys_in_slices = {}
+    # inform the user
+    LGR.warning('\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+                '\nphys2bids will split the input file according to the given -tr and -ntp'
+                ' arguments'
+                '\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
+    # Find the timestamps
+    run_timestamps = find_runs(phys_in, ntp_list, tr_list, thr, padding)
+    for n, run in enumerate(run_timestamps.keys()):
+
+        # tmp variable to collect run's info
+        run_attributes = run_timestamps[run]
+
+        phys_in_slices[run] = deepcopy(phys_in[run_attributes[0]:run_attributes[1]])
+
+        # Run check_trigger amount
+        phys_in_slices[run].check_trigger_amount(thr=thr,
+                                                 num_timepoints_expected=ntp_list[n],
+                                                 tr=tr_list[n])
+
+    return phys_in_slices
diff --git a/phys2bids/viz.py b/phys2bids/viz.py
index 190a61a..0cb0938 100644
--- a/phys2bids/viz.py
+++ b/phys2bids/viz.py
@@ -80,7 +80,7 @@ def plot_trigger(time, trigger, fileprefix, tr, thr, num_timepoints_expected,
     subplot.set_ylabel('Volts')
     subplot.plot(time, trigger, '-', time, thrline, 'r-.', time, block, '-')
     subplot.fill_between(time, block, where=block >= d, interpolate=True, color='#ffbb6e')
-    subplot.legend(["trigger", "Trigger detection threshold", "time block"], loc='upper right')
+    subplot.legend(['trigger', 'Trigger detection threshold', 'time block'], loc='upper right')
     # plot the first spike according to the user threshold
     subplot = fig.add_subplot(223)
     subplot.set_xlim([-tr * 4, tr * 4])
@@ -128,6 +128,57 @@ def plot_trigger(time, trigger, fileprefix, tr, thr, num_timepoints_expected,
     plt.close()
 
 
+def export_trigger_plot(phys_in, chtrig, fileprefix, tr, num_timepoints_expected,
+                        filename, sub=None, ses=None, figsize=FIGSIZE,
+                        dpi=SET_DPI):
+    """
+    Save a trigger plot.
+
+    Used in main workflow (`phys2bids`), this function minimizes repetition in code for parallel
+    workflow (multi-run workflow and default workflow) and maintains readability of code
+
+    Parameters
+    ---------
+    phys_in : object
+        Object returned by BlueprintInput class
+        For multi-run acquisitions, phys_in is a slice of the whole object
+    chtrig : int
+        trigger channel
+        integer representing the index of the trigger on phys_in.timeseries
+    fileprefix: str or path
+        A string representing a file name or a fullpath
+        to a file, WITHOUT extension
+    tr : list
+        a list of float given by the user as `tr` input
+    num_timepoints_expected : list
+        a list of integers given by the user as `ntp` input
+    filename : str
+        name of the input file given by user's entry
+    sub: str or int
+        Name of subject. Default is None
+    ses: str or int or None
+        Name of session. Default is None
+    figsize: tuple or list of floats
+        Size of the figure expressed as (size_x, size_y),
+        Default is {FIGSIZE}
+    dpi: int
+        Desired DPI of the figure,
+        Default is {SET_DPI}
+    """
+    LGR.info('Plot trigger')
+    # Create trigger plot. If possible, to have multiple outputs in the same
+    # place, adds sub and ses label.
+    if sub is not None:
+        fileprefix += f'_sub-{sub}'
+    if ses is not None:
+        fileprefix += f'_ses-{ses}'
+
+    # adjust for multi run arguments, iterate through acquisition attributes
+    plot_trigger(phys_in.timeseries[0], phys_in.timeseries[chtrig],
+                 fileprefix, tr, phys_in.thr, num_timepoints_expected,
+                 filename, figsize, dpi)
+
+
 def plot_all(ch_name, timeseries, units, freq, infile, outfile, dpi=SET_DPI, size=FIGSIZE):
     """
     Plot all the channels for visualizations and saves them in outfile.
@@ -178,7 +229,7 @@ def plot_all(ch_name, timeseries, units, freq, infile, outfile, dpi=SET_DPI, siz
         ax[row].set_ylabel(units[row + 1])
         ax[row].xlim = 30 * 60 * freq[0]  # maximum display of half an hour
         ax[row].grid()
-    ax[row].set_xlabel("seconds")
+    ax[row].set_xlabel('seconds')
     outfile = os.path.join(outfile, os.path.splitext(os.path.basename(infile))[0] + '.png')
     LGR.info(f'saving channel plot to {outfile}')
     fig.savefig(outfile, dpi=dpi, bbox_inches='tight')
diff --git a/requirements.txt b/requirements.txt
index 3118cec..c8772ca 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,3 @@
-numpy>=1.9.3
-matplotlib>=3.1.1
-PyYAML
\ No newline at end of file
+numpy>=1.9.3, !=*rc*
+matplotlib>=3.1.1 ,!=3.3.0rc1
+PyYAML!=*rc*
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
index a941818..1c921e4 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -22,9 +22,9 @@ provides =
 [options]
 python_requires = >=3.6.1
 install_requires =
-    numpy >=1.9.3
-    matplotlib >=3.1.1
-    PyYAML
+    numpy >=1.9.3, !=*rc*
+    matplotlib >=3.1.1, !=3.3.0rc1
+    PyYAML !=*rc*
 tests_require =
     pytest >=3.6
 test_suite = pytest


Test patch:
diff --git a/phys2bids/heuristics/heur_test_acq.py b/phys2bids/heuristics/heur_test_acq.py
index 9f553dc..8836119 100644
--- a/phys2bids/heuristics/heur_test_acq.py
+++ b/phys2bids/heuristics/heur_test_acq.py
@@ -1,7 +1,7 @@
 import fnmatch
 
 
-def heur(physinfo):
+def heur(physinfo, run=''):
     """
     Set of if .. elif statements to fill BIDS names.
 
diff --git a/phys2bids/heuristics/heur_test_multifreq.py b/phys2bids/heuristics/heur_test_multifreq.py
index c6bdbbd..bd91a69 100644
--- a/phys2bids/heuristics/heur_test_multifreq.py
+++ b/phys2bids/heuristics/heur_test_multifreq.py
@@ -1,7 +1,7 @@
 import fnmatch
 
 
-def heur(physinfo):
+def heur(physinfo, run=''):
     """
     Set of if .. elif statements to fill BIDS names.
 
diff --git a/phys2bids/tests/test_bids.py b/phys2bids/tests/test_bids.py
index 8bb6598..6c062ee 100644
--- a/phys2bids/tests/test_bids.py
+++ b/phys2bids/tests/test_bids.py
@@ -1,5 +1,4 @@
 import os
-
 from pkg_resources import resource_filename
 
 import pytest
@@ -8,7 +7,6 @@ from csv import reader
 
 from phys2bids import bids
 from phys2bids.bids import UNIT_ALIASES
-from phys2bids.utils import append_list_as_row
 
 
 def test_bidsify_units():
@@ -42,7 +40,8 @@ def test_use_heuristic(tmpdir, test_sub, test_ses):
     test_record_label = 'test'
 
     heur_path = bids.use_heuristic(test_full_heur_path, test_sub, test_ses,
-                                   test_full_input_path, test_outdir, test_record_label)
+                                   test_full_input_path, test_outdir,
+                                   record_label=test_record_label)
 
     if test_sub[:4] == 'sub-':
         test_sub = test_sub[4:]
diff --git a/phys2bids/tests/test_integration.py b/phys2bids/tests/test_integration.py
index 3aa0ce1..d8f8d79 100644
--- a/phys2bids/tests/test_integration.py
+++ b/phys2bids/tests/test_integration.py
@@ -33,7 +33,7 @@ def test_integration_acq(samefreq_full_acq_file):
     conversion_path = join(test_path, 'code', 'conversion')
 
     phys2bids(filename=test_filename, indir=test_path, outdir=test_path,
-              chtrig=test_chtrig, num_timepoints_expected=1)
+              chtrig=test_chtrig, num_timepoints_expected=1, tr=1)
 
     # Check that files are generated
     for suffix in ['.json', '.tsv.gz']:
@@ -120,8 +120,8 @@ def test_integration_heuristic(multifreq_lab_file):
     # Check that files are generated in conversion path
     for freq in ['40', '100', '500', '1000']:
         assert isfile(join(conversion_path,
-                           'sub-006_ses-01_task-test_rec-biopac_run-01_recording-'
-                           + freq + '.log'))
+                           'sub-006_ses-01_task-test_rec-biopac_run-01_'
+                           f'recording-{freq}Hz_physio.log'))
     assert isfile(join(conversion_path,
                        'Test1_multifreq_onescan_sub-006_ses-01_trigger_time.png'))
     assert isfile(join(conversion_path, 'Test1_multifreq_onescan.png'))
@@ -131,12 +131,13 @@ def test_integration_heuristic(multifreq_lab_file):
     # Check that files are generated
     base_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-'
     for suffix in ['.json', '.tsv.gz']:
-        for freq in ['40.0', '100.0', '500.0', '1000.0']:
-            assert isfile(join(test_path_output, base_filename + freq + '_physio' + suffix))
+        for freq in ['40', '100', '500', '1000']:
+            assert isfile(join(test_path_output,
+                               f'{base_filename}{freq}Hz_physio{suffix}'))
 
     # ##### Checks for 40 Hz files
     # Read log file (note that this file is not the logger file)
-    log_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-40.log'
+    log_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-40Hz_physio.log'
     with open(join(conversion_path, log_filename)) as log_info:
         log_info = log_info.readlines()
 
@@ -152,7 +153,7 @@ def test_integration_heuristic(multifreq_lab_file):
     assert check_string(log_info, 'first trigger', 'Time 0', is_num=False)
 
     # Checks json file
-    json_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-40.0_physio.json'
+    json_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-40Hz_physio.json'
     with open(join(test_path_output, json_filename)) as json_file:
         json_data = json.load(json_file)
 
@@ -163,7 +164,7 @@ def test_integration_heuristic(multifreq_lab_file):
 
     # ##### Checks for 100 Hz files
     # Read log file (note that this file is not the logger file)
-    log_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-100.log'
+    log_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-100Hz_physio.log'
     with open(join(conversion_path, log_filename)) as log_info:
         log_info = log_info.readlines()
 
@@ -179,7 +180,7 @@ def test_integration_heuristic(multifreq_lab_file):
     assert check_string(log_info, 'first trigger', 'Time 0', is_num=False)
 
     # Checks json file
-    json_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-100.0_physio.json'
+    json_filename = 'sub-006_ses-01_task-test_rec-biopac_run-01_recording-100Hz_physio.json'
     with open(join(test_path_output, json_filename)) as json_file:
         json_data = json.load(json_file)
 
diff --git a/phys2bids/tests/test_phys2bids.py b/phys2bids/tests/test_phys2bids.py
index aea1f59..73dd56b 100644
--- a/phys2bids/tests/test_phys2bids.py
+++ b/phys2bids/tests/test_phys2bids.py
@@ -10,7 +10,7 @@ from phys2bids import phys2bids
 
 def test_print_summary(tmpdir):
     test_filename = 'input.txt'
-    test_ntp_expected = 10
+    test_ntp_expected = [10, ]
     test_ntp_found = 5
     test_samp_freq = 0.2
     test_time_offset = 0.8