Repository: SpikeInterface/spikeinterface
Pull request #131 (2021-05-03 09:23:42)

concatenate_recordings / append_recordings
@alejoe91 : tell me what do you think.


Patch:
diff --git a/TODO_big_refactoring.txt b/TODO_big_refactoring.txt
index 5a12024b9..a16a3a759 100644
--- a/TODO_big_refactoring.txt
+++ b/TODO_big_refactoring.txt
@@ -16,7 +16,6 @@ extractors:
     * all event reader (get_ttl()) (alessio)
     * MdaRecordingExtractor  : need to be done in neo ??
     * BiocamRecordingExtractor  : need to be done in neo
-    * MaxOneRecordingExtractor  : need to be done in neo
     * MCSH5RecordingExtractor  : need to be done in neo ???
     * SHYBRIDRecordingExtractor:  need to be done in neo ???
     * NIXIORecordingExtractor: wrappe neo
diff --git a/environment_rtd.yml b/environment_rtd.yml
index 795f5f7b2..a3161880c 100644
--- a/environment_rtd.yml
+++ b/environment_rtd.yml
@@ -23,7 +23,7 @@ dependencies:
   - pip:
     - git+https://github.com/NeuralEnsemble/python-neo.git
     - probeinterface
-    - MEArec
+    - MEArec>=1.7.1
     - isosplit5
     - mountainsort4>=1.0.0
     - hdbscan
diff --git a/examples/modules/core/plot_5_append_concatenate_segments.py b/examples/modules/core/plot_5_append_concatenate_segments.py
new file mode 100644
index 000000000..b1afece76
--- /dev/null
+++ b/examples/modules/core/plot_5_append_concatenate_segments.py
@@ -0,0 +1,65 @@
+"""
+Append and/or concatenate segments
+===================================
+
+Sometimes a recording can be splitted in several sequence, for instance because baseline/intervention/pauses.
+
+The terminology in spikeniterface is "segment" (like in neo).
+
+spikeinterface have tools to manipulate theses segments. There are two ways:
+
+  1. append_recordings()/append_sortings()/append_events()
+  2. concatenate_recordings()/ concatenate_sortings()/ concatenate_events()
+
+Imagine we have 2 recordings with respectively 2 and 3 segments:
+
+  1. In case 1 we will have one recording with 5 segments
+  2. In case 2 we will have one recording with 1 "big" segment that virtually concatenate all segments.
+
+Here a short example.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+from spikeinterface import NumpyRecording, NumpySorting
+from spikeinterface import append_recordings, concatenate_recordings
+
+##############################################################################
+# First let's generate 2 recordings with 2 and 3 segments respectively.
+# 
+
+sampling_frequency = 1000.
+
+trace0 = np.zeros((150, 5), dtype='float32')
+trace1 = np.zeros((100, 5), dtype='float32')
+rec0 = NumpyRecording([trace0, trace1], sampling_frequency)
+print(rec0)
+
+trace2 = np.zeros((50, 5), dtype='float32')
+trace3 = np.zeros((200, 5), dtype='float32')
+trace4 = np.zeros((120, 5), dtype='float32')
+rec1 = NumpyRecording([trace2, trace3, trace4], sampling_frequency)
+print(rec1)
+
+
+##############################################################################
+# lets use `append_recordings()`
+
+recording_list = [rec0, rec1]
+rec = append_recordings(recording_list)
+print(rec)
+for i in range(rec.get_num_segments()):
+    s = rec.get_num_samples(segment_index=i)
+    print(f'segment {i} num_samples {s}')
+
+##############################################################################
+# lets use `concatenate_recordings()`
+
+recording_list = [rec0, rec1]
+rec = concatenate_recordings(recording_list)
+print(rec)
+s = rec.get_num_samples(segment_index=0)
+print(f'segment {0} num_samples {s}')
+
diff --git a/requirements_rtd.txt b/requirements_rtd.txt
deleted file mode 100644
index 790594200..000000000
--- a/requirements_rtd.txt
+++ /dev/null
@@ -1,49 +0,0 @@
-# main dependancy
-numpy
-neo
-joblib
-tqdm
-probeinterface
-
-# sphinx-gallery
-pillow
-sphinx-gallery
-numpydoc
-matplotlib
-
-
-# optional dependancy for sub modules
-h5py
-pandas
-scipy
-sklearn
-matplotlib
-networkx
-datalad
-MEArec
-
-## install sorters ##
-
-## mountainsort4
-mountainsort4>=1.0.0
-
-# klusta
-# Do not work anymore on python3.8
-# Cython
-# tqdm
-# click
-# klusta
-# klustakwik2==0.2.7
-
-
-## tridesclous
-pyqt5
-hdbscan
-loky
-tridesclous
-
-## spyking-circus impossible to run RTD because need mpich
-# spyking-circus
-
-## herdingspikes
-herdingspikes
diff --git a/spikeinterface/core/__init__.py b/spikeinterface/core/__init__.py
index ecb75dc4f..2302c9a57 100644
--- a/spikeinterface/core/__init__.py
+++ b/spikeinterface/core/__init__.py
@@ -22,6 +22,13 @@ from .numpyextractors import NumpyRecording , NumpySorting
 from .channelslicerecording import ChannelSliceRecording
 from .unitsselectionsorting import UnitsSelectionSorting
 
+# utils to append and concatenate segment (equivalent to OLD MultiRecordingTimeExtractor)
+from .segmentutils import (
+    append_recordings, AppendSegmentRecording,
+    concatenate_recordings, ConcatenateSegmentRecording,
+    append_sortings, AppendSegmentSorting)
+
+
 # default folder
 from .default_folders import (set_global_tmp_folder, get_global_tmp_folder,
     is_set_global_tmp_folder, reset_global_tmp_folder, 
diff --git a/spikeinterface/core/segmentutils.py b/spikeinterface/core/segmentutils.py
new file mode 100644
index 000000000..04a881272
--- /dev/null
+++ b/spikeinterface/core/segmentutils.py
@@ -0,0 +1,208 @@
+"""
+Implementation if utils class to manipulate segment with 2 diffrents concept:
+  * append_recordings/append_sortings/append_events
+  * concatenate_recordings/ concatenate_sortings/ concatenate_events
+
+
+For instance:
+  * append_recording: given one recording with 2 segments and one recording with 3 segments give one recording with 5 segments
+  * conatenate_recording: given a list of several recording with several segment give one recording with one segment that sum the total duration
+
+"""
+import numpy as np
+
+from .baserecording import BaseRecording, BaseRecordingSegment
+from .basesorting import BaseSorting, BaseSortingSegment
+from .baseevent import BaseEvent, BaseEventSegment
+
+
+class AppendSegmentRecording(BaseRecording):
+    """
+    Return a recording that "appends" all segments from all recording
+    into one recording multi segment.
+    
+    For instance, given one recording with 2 segments and one recording with
+    3 segments, this class will give one recording with 5 segments
+
+    Parameters
+    ----------
+    recording_list : list of BaseRecording
+        A list of recordings
+    """
+    def __init__(self, recording_list):
+        
+        rec0 = recording_list[0]
+        sampling_frequency = rec0.get_sampling_frequency()
+        dtype = rec0.get_dtype()
+        channel_ids = rec0.channel_ids
+        
+        # check same carracteristics
+        ok1 = all(sampling_frequency == rec.get_sampling_frequency() for rec in recording_list)
+        ok2 = all(dtype == rec.get_dtype() for rec in recording_list)
+        ok3 = all(np.array_equal(channel_ids, rec.channel_ids) for rec in recording_list)
+        if not (ok1 and ok2 and ok3):
+            raise ValueError("Recording don't have the same sampling_frequency/dtype/channel_ids")
+        
+        BaseRecording.__init__(self, sampling_frequency, channel_ids, dtype)
+        self.copy_metadata(rec0)
+        
+        for rec in recording_list:
+            for parent_segment in rec._recording_segments:
+                rec_seg = ProxyAppendRecordingSegment(parent_segment)
+                self.add_recording_segment(rec_seg)
+
+        self._kwargs = {'recording_list': [rec.to_dict() for rec in recording_list]}
+
+
+class ProxyAppendRecordingSegment(BaseRecordingSegment):
+    def __init__(self, parent_segment):
+        BaseRecordingSegment.__init__(self)
+        self.parent_segment = parent_segment
+
+    def get_num_samples(self):
+        return self.parent_segment.get_num_samples()
+
+    def get_traces(self, *args, **kwargs):
+        return self.parent_segment.get_traces(*args, **kwargs)
+
+def append_recordings(*args, **kargs):
+    return AppendSegmentRecording(*args, **kargs)
+append_recordings.__doc__ == AppendSegmentRecording.__doc__
+
+
+class ConcatenateSegmentRecording(BaseRecording):
+    """
+    Return a recording that "concatenates" all segments from all recording
+    into one recording mono segment. The operation is lazy.
+    
+    For instance, given one recording with 2 segments and one recording with
+    3 segments, this class will give one recording with 1 segment
+
+    Parameters
+    ----------
+    recording_list : list of BaseRecording
+        A list of recordings
+    """
+    def __init__(self, recording_list):
+        
+        one_rec = append_recordings(recording_list)
+        
+        BaseRecording.__init__(self, one_rec.get_sampling_frequency(), one_rec.channel_ids, one_rec.get_dtype())
+        self.copy_metadata(one_rec)
+        
+        parent_segments = []
+        for rec in recording_list:
+            for parent_segment in rec._recording_segments:
+                parent_segments.append(parent_segment)
+        rec_seg = ProxyConcatenateRecordingSegment(parent_segments)
+        self.add_recording_segment(rec_seg)
+
+        self._kwargs = {'recording_list': [rec.to_dict() for rec in recording_list]}
+
+
+class ProxyConcatenateRecordingSegment(BaseRecordingSegment):
+    def __init__(self, parent_segments):
+        BaseRecordingSegment.__init__(self)
+        self.parent_segments = parent_segments
+        self.all_length = [rec_seg.get_num_samples() for rec_seg in self.parent_segments]
+        self.cumsum_length = np.cumsum([0] + self.all_length)
+        self.total_length = np.sum(self.all_length)
+
+    def get_num_samples(self):
+        return self.total_length
+
+    def get_traces(self, start_frame, end_frame, channel_indices):
+        if start_frame is None:
+            start_frame = 0
+        if end_frame is None:
+            end_frame = self.get_num_samples()
+            
+        i0 = np.searchsorted(self.cumsum_length, start_frame, side='right') - 1
+        i1 = np.searchsorted(self.cumsum_length, end_frame,side='right') - 1
+        
+        # several case:
+        #  * come from one segment (i0 == i1)
+        #  * come from several segment (i0 < i1)
+        
+        if i0 == i1:
+            # one segment
+            rec_seg = self.parent_segments[i0]
+            seg_start = self.cumsum_length[i0]
+            traces = rec_seg.get_traces(start_frame - seg_start, end_frame - seg_start, channel_indices)
+        else:
+            # sveral segments
+            all_traces = []
+            for i in range(i0, i1+1):
+                if i == len(self.parent_segments):
+                    # limit case
+                    continue
+
+                rec_seg = self.parent_segments[i]
+                seg_start = self.cumsum_length[i]
+                if i == i0:
+                    # first
+                    traces_chunk = rec_seg.get_traces(start_frame - seg_start, None, channel_indices)
+                    all_traces.append(traces_chunk)
+                elif i == i1:
+                    #last
+                    if (end_frame - seg_start) >0:
+                        traces_chunk = rec_seg.get_traces(None, end_frame - seg_start, channel_indices)
+                        all_traces.append(traces_chunk)
+                else:
+                    # in between
+                    traces_chunk = rec_seg.get_traces(None, None, channel_indices)
+                    all_traces.append(traces_chunk)
+            traces = np.concatenate(all_traces, axis=0)
+        
+        return traces
+
+
+def concatenate_recordings(*args, **kargs):
+    return ConcatenateSegmentRecording(*args, **kargs)
+concatenate_recordings.__doc__ == ConcatenateSegmentRecording.__doc__
+
+
+class AppendSegmentSorting(BaseSorting):
+    """
+    Return a sorting that "append" all segments from all sorting
+    into one sorting multi segment.
+    
+    Parameters
+    ----------
+    sorting_list : list of BaseSorting
+        A list of sortings
+    """
+    def __init__(self, sorting_list):
+        
+        sorting0 = sorting_list[0]
+        sampling_frequency = sorting0.get_sampling_frequency()
+        unit_ids = sorting0.unit_ids
+        
+        # check same carracteristics
+        ok1 = all(sampling_frequency == sorting.get_sampling_frequency() for sorting in sorting_list)
+        ok2 = all(np.array_equal(unit_ids, sorting.unit_ids) for sorting in sorting_list)
+        if not (ok1 and  ok2):
+            raise ValueError("Sorting don't have the same sampling_frequency/unit_ids")
+        
+        BaseSorting.__init__(self, sampling_frequency, unit_ids)
+        self.copy_metadata(sorting0)
+        
+        for sorting in sorting_list:
+            for parent_segment in sorting._sorting_segments:
+                sorting_seg = ProxyAppendSortingSegment(parent_segment)
+                self.add_sorting_segment(sorting_seg)
+
+        self._kwargs = {'sorting_list': [sorting.to_dict() for sorting in sorting_list]}
+
+
+class ProxyAppendSortingSegment(BaseSortingSegment):
+    def __init__(self, parent_segment):
+        BaseSortingSegment.__init__(self)
+        self.parent_segment = parent_segment
+
+    def get_unit_spike_train(self, *args, **kwargs):
+        return self.parent_segment.get_unit_spike_train(*args, **kwargs)
+
+def append_sortings(*args, **kargs):
+    return AppendSegmentSorting(*args, **kargs)
+append_sortings.__doc__ == AppendSegmentSorting.__doc__
diff --git a/spikeinterface/extractors/extractorlist.py b/spikeinterface/extractors/extractorlist.py
index 7f7777cca..1ecbd4490 100644
--- a/spikeinterface/extractors/extractorlist.py
+++ b/spikeinterface/extractors/extractorlist.py
@@ -24,6 +24,9 @@ from .neoextractors import (
     BlackrockRecordingExtractor,
     MCSRawRecordingExtractor,
     KiloSortSortingExtractor,
+    
+    # not in python-neo master yet
+    # MaxwellRecordingExtractor,  
 )
 
 recording_extractor_full_list = [
@@ -43,6 +46,9 @@ recording_extractor_full_list = [
     NeuralynxRecordingExtractor,
     BlackrockRecordingExtractor,
     MCSRawRecordingExtractor,
+    
+    # not in python-neo master yet
+    # MaxwellRecordingExtractor,
 
     ##OLD
     # ~ MdaRecordingExtractor,
diff --git a/spikeinterface/extractors/neoextractors/__init__.py b/spikeinterface/extractors/neoextractors/__init__.py
index 5b9d200a7..57d1eb218 100644
--- a/spikeinterface/extractors/neoextractors/__init__.py
+++ b/spikeinterface/extractors/neoextractors/__init__.py
@@ -9,3 +9,6 @@ from .neuralynx import NeuralynxRecordingExtractor
 from .blackrock import BlackrockRecordingExtractor
 from .mscraw import MCSRawRecordingExtractor
 from .kilosort import KiloSortSortingExtractor
+
+# not in python-neo master yet
+# from .maxwell import MaxwellRecordingExtractor
diff --git a/spikeinterface/extractors/neoextractors/maxwell.py b/spikeinterface/extractors/neoextractors/maxwell.py
new file mode 100644
index 000000000..65aa604d9
--- /dev/null
+++ b/spikeinterface/extractors/neoextractors/maxwell.py
@@ -0,0 +1,31 @@
+import neo
+
+from .neobaseextractor import NeoBaseRecordingExtractor
+
+
+class MaxwellRecordingExtractor(NeoBaseRecordingExtractor):
+    """
+    Class for reading data from Maxwell device.
+    It handle MaxOne (old and new format) and MaxTwo.
+    
+    
+    
+    Based on neo.rawio.IntanRawIO
+    
+    Parameters
+    ----------
+    file_path: str
+        
+    stream_id: str or None
+        For Maxtwo when there are several well at the same time you 
+        need to speficify. stream_id='well000' or 'well0001' or ...
+    rec_name: str or None
+        When th file contain several block (aka rec) you need to specify the one
+        you want.  rec_name='rec0000'
+    """ 
+    mode = 'file'
+    NeoRawIOClass = 'MaxwellRawIO'
+
+    def __init__(self, file_path, stream_id=None, rec_name=None):
+        neo_kwargs = {'filename' : str(file_path)}
+        NeoBaseRecordingExtractor.__init__(self, stream_id=stream_id, rec_name=rec_name)
diff --git a/spikeinterface/toolkit/preprocessing/common_reference.py b/spikeinterface/toolkit/preprocessing/common_reference.py
index a69a16fdb..dfb8d461c 100644
--- a/spikeinterface/toolkit/preprocessing/common_reference.py
+++ b/spikeinterface/toolkit/preprocessing/common_reference.py
@@ -134,17 +134,15 @@ class CommonReferenceRecordingSegment(BasePreprocessorSegment):
                 for chan_ind in _channel_indices])
         
         return out_traces
-
+    
     def _groups(self, channel_indices):
         selected_groups = []
         selected_channels = []
         if self.groups:
             for chan_inds in self.groups:
+                selected_groups.append(chan_inds)
                 sel_inds = [ind for ind in channel_indices if ind in chan_inds]
-                # if no channels are in a group, do not return the group
-                if len(sel_inds) > 0:
-                    selected_groups.append(chan_inds)
-                    selected_channels.append(sel_inds)
+                selected_channels.append(sel_inds)
         else:
             # no groups = all channels
             selected_groups = [slice(None)]


Test patch:
diff --git a/requirements_test.txt b/requirements_test.txt
index 6eb6c602e..6141bdb0f 100644
--- a/requirements_test.txt
+++ b/requirements_test.txt
@@ -8,7 +8,7 @@ scikit-learn
 matplotlib
 networkx
 datalad
-MEArec
+MEArec>=1.7.1
 
 
 ## install sorters ##
diff --git a/spikeinterface/core/tests/test_segmentutils.py b/spikeinterface/core/tests/test_segmentutils.py
new file mode 100644
index 000000000..823cc7ecc
--- /dev/null
+++ b/spikeinterface/core/tests/test_segmentutils.py
@@ -0,0 +1,78 @@
+import pytest
+import numpy as np
+
+from spikeinterface.core import(
+    append_recordings, AppendSegmentRecording,
+    concatenate_recordings, ConcatenateSegmentRecording,
+    append_sortings, AppendSegmentSorting)
+
+from spikeinterface.core import NumpyRecording, NumpySorting
+
+
+def test_append_concatenate_recordings():
+    traces = np.zeros((1000, 5), dtype='float64')
+    traces[:] = np.arange(1000)[:, None]
+    sampling_frequency = 30000
+    rec0 = NumpyRecording([traces] * 3, sampling_frequency)
+    rec1 = NumpyRecording([traces] * 2, sampling_frequency)
+    
+    # append
+    rec = append_recordings([rec0, rec1])
+    # print(rec)
+    assert rec.get_num_segments() == 5
+    for segment_index in range(5):
+        traces = rec.get_traces(segment_index=segment_index)
+        assert rec.get_num_samples(segment_index) == 1000
+    
+    # concatenate
+    rec = concatenate_recordings([rec0, rec1])
+    # print(rec)
+    assert rec.get_num_samples(0) == 5* 1000
+    assert rec.get_num_segments() == 1
+    
+    # case one segment
+    traces = rec.get_traces(start_frame=0, end_frame=15)
+    assert np.array_equal(traces[:, 0], np.arange(0, 15))
+    traces = rec.get_traces(start_frame=500, end_frame=750)
+    assert np.array_equal(traces[:, 0], np.arange(500, 750))
+    
+    # case on limit
+    traces = rec.get_traces(start_frame=1000, end_frame=2000)
+    assert traces.shape == (1000, 5)
+    assert np.array_equal(traces[:, 0], np.arange(0, 1000, dtype='float64'))
+    
+    
+    # case total
+    traces = rec.get_traces(start_frame=None, end_frame=None)
+    assert traces.shape == (5000, 5)
+    assert traces[1000, 0] == 0
+    
+    # several segments
+    traces = rec.get_traces(start_frame=50, end_frame=4500)
+    assert traces.shape == (4450, 5)
+    assert traces[0, 0] == 50
+    assert traces[-1, 0] == 499
+    
+
+def test_append_sortings():
+    sampling_frequency = 30000.
+    times = np.arange(0, 1000, 10)
+    labels = np.zeros(times.size, dtype='int64')
+    labels[0::3] = 0
+    labels[1::3] = 1
+    labels[2::3] = 2
+    sorting0 = NumpySorting.from_times_labels([times]*3, [labels]*3, sampling_frequency)
+    sorting1 = NumpySorting.from_times_labels([times]*2, [labels]*2, sampling_frequency)
+    
+    sorting = append_sortings([sorting0, sorting1])
+    # print(sorting)
+    assert sorting.get_num_segments() == 5
+
+    
+    
+if __name__ == '__main__':
+    #~ test_append_concatenate_recordings()
+    test_append_sortings()
+
+
+    
\ No newline at end of file
diff --git a/spikeinterface/extractors/tests/test_neoextractors.py b/spikeinterface/extractors/tests/test_neoextractors.py
index d833e0e70..686a725e5 100644
--- a/spikeinterface/extractors/tests/test_neoextractors.py
+++ b/spikeinterface/extractors/tests/test_neoextractors.py
@@ -107,11 +107,20 @@ class KiloSortSortingTest(SortingCommonTestSuite, unittest.TestCase):
         'phy/phy_example_0',
     ]
 
+#~ @pytest.mark.skip(reason="Maxwell HDF5 compression need a manual installable plugin!!!")
+#~ class MaxwellRecordingTest(RecordingCommonTestSuite, unittest.TestCase):
+    #~ ExtractorClass = MaxwellRecordingExtractor
+    #~ downloads = ['maxwell']
+    #~ entities = [
+        #~ 'maxwell/MaxOne_data/Record/000011/data.raw.h5',
+        #~ ('maxwell/MaxTwo_data/Network/000028/data.raw.h5', {'stream_id':'well0000', 'rec_name': 'rec0000'})
+    #~ ]
+
 
 if __name__ == '__main__':
     #~ test = MearecRecordingTest()
     #~ test = MearecSortingTest()
-    test = SpikeGLXRecordingTest()
+    #~ test = SpikeGLXRecordingTest()
     #~ test = OpenEphysBinaryRecordingTest()
     # test = OpenEphysLegacyRecordingTest()
     # test = ItanRecordingTest()
@@ -121,6 +130,7 @@ if __name__ == '__main__':
     # test = BlackrockRecordingTest()
     # test = MCSRawRecordingTest()
     # test = KiloSortSortingTest()
+    #test = MaxwellRecordingTest()
     
     test.setUp()
     test.test_open()