Update example script to new io and requirements

code/example.py

@@ -48,7 +48,6 @@ import quantities as pq
 
 from neo import Block, Segment
 from elephant.signal_processing import butter
-
 from reachgraspio import reachgraspio
 from neo.utils import cut_segment_by_epoch, add_epoch, get_events
 # from neo.utils import add_epoch, cut_segment_by_epoch, get_events
@@ -69,15 +68,11 @@ odml_dir = os.path.join('..', 'datasets')
 # Open the session for reading
 session = reachgraspio.ReachGraspIO(session_name, odml_directory=odml_dir)
 
-# Read the first 300s of data (time series at 1000Hz (ns2) and 30kHz (ns6)
-# scaled to units of voltage, sorted spike trains, spike waveforms and events)
-# from electrode 62 of the recording session and return it as a Neo Block. The
-# time shift of the ns2 signal (LFP) induced by the online filter is
-# automatically corrected for by a heuristic factor stored in the metadata
-# (correct_filter_shifts=True).
-data_block = session.read_block(load_waveforms=True, correct_filter_shifts=True, lazy=True)
+# Read the complete dataset in lazy mode. Neo object will be created, but
+# data are not loaded in to memory.
+data_block = session.read_block(correct_filter_shifts=True, lazy=True)
 
-# Access the single Segment of the data block, reaching up to 300s.
+# Access the single Segment of the data block
 assert len(data_block.segments) == 1
 data_segment = data_block.segments[0]
 
@@ -92,32 +87,49 @@ data_segment = data_block.segments[0]
 # Neo AnalogSignal, which is used for plotting later on in this script.
 # =============================================================================
 
-filtered_anasig = []
-# Loop through all AnalogSignal objects in the loaded data
-# Create LFP representation of the data if not preset already
-for anasig in data_block.segments[0].analogsignals:
-    if all(anasig.array_annotations['nsx'] == 2):
-        # AnalogSignal is LFP from ns2
-        anasig.name = 'LFP (online filtered, ns2)'
-    elif all(np.isin(anasig.array_annotations['nsx'], [5, 6])):
-        # AnalogSignal is raw signal from ns5 or ns6
-        assert len(np.unique(anasig.array_annotations['nsx'])) == 1
-        anasig.name = 'raw (ns%i)' % np.unique(anasig.array_annotations['nsx'])[0]
-
-        # temporarily loading the data into memory
-        anasig_loaded = anasig.load(time_slice=(None, 300*pq.s), channel_indexes=[63])
-
-        # Use the Elephant library to filter the analog signal
-        f_anasig = butter(
-                anasig_loaded,
-                highpass_freq=None,
-                lowpass_freq=250 * pq.Hz,
-                order=4)
-        print('filtering done.')
-        f_anasig.name = 'LFP (offline filtered ns%i)' % np.unique(anasig.array_annotations['nsx'])[0]
-        filtered_anasig.append(f_anasig)
-# Attach all offline filtered LFPs to the segment of data
-data_block.segments[0].analogsignals.extend(filtered_anasig)
+# Iterate through all analog signals and replace these lazy object by new
+# analog signals containing only data of channel 62 (target_channel_id) and
+# provide human readable name for the analog signal (LFP / raw signal type)
+
+target_channel_id = 62
+nsx_to_anasig_name = {2: 'LFP signal (online filtered)',
+                      5: 'raw signal',
+                      6: 'raw signal'}
+
+idx = 0
+while idx < len(data_segment.analogsignals):
+    # remove analog signals, that don't contain target channel
+    channel_ids = data_segment.analogsignals[idx].array_annotations['channel_ids']
+    if target_channel_id not in channel_ids:
+        data_segment.analogsignals.pop(idx)
+        continue
+
+    # replace analog signal with analog signal containing data
+    target_channel_index = np.where(channel_ids == target_channel_id)[0][0]
+    anasig = data_segment.analogsignals[idx].load(
+        channel_indexes=[target_channel_index])
+    data_segment.analogsignals[idx] = anasig
+    idx += 1
+
+    # replace name by label of contained signal type
+    anasig.name = nsx_to_anasig_name[anasig.array_annotations['nsx'][0]]
+
+# load spiketrains of same channel
+channel_spiketrains = data_segment.filter({'channel_id': target_channel_id})
+data_segment.spiketrains = [st.load(load_waveforms=True) for st in channel_spiketrains]
+
+# The LFP is not present in the data fils of both recording. Here, we
+# generate the LFP signal from the raw signal if it's not present already.
+if not data_segment.filter({'name': 'LFP signal (online filtered)'}):
+    raw_signal = data_segment.filter({'name': 'raw signal'})[0]
+
+    # Use the Elephant library to filter the raw analog signal
+    f_anasig = butter(raw_signal,  highpass_freq=None, lowpass_freq=250 * pq.Hz, order=4)
+    print('filtering done.')
+
+    f_anasig.name = 'LFP signal (offline filtered)'
+    # Attach offline filtered LFP to the segment of data
+    data_segment.analogsignals.extend(f_anasig)
 
 
 # =============================================================================
@@ -174,10 +186,11 @@ print('added epoch.')
 # such that each segment starts at time 0 (parameter reset_times); annotations
 # describing the analysis epoch are carried over to the segments. A new Neo
 # Block named "data_cut_to_analysis_epochs" is created to capture all cut
-# analysis epochs.
+# analysis epochs. For execution time reason, we are only considering the
+# first 10 epochs here.
 cut_trial_block = Block(name="data_cut_to_analysis_epochs")
 cut_trial_block.segments = cut_segment_by_epoch(
-    data_segment, epoch, reset_time=True)
+    data_segment, epoch[:10], reset_time=True)
 
 # =============================================================================
 # Plot data
@@ -238,7 +251,7 @@ for event in trial_segment.events:
     if event.name == 'TrialEvents':
         for ev_id, ev in enumerate(event):
                 plt.axvline(
-                    ev,
+                    ev.rescale(time_unit),
                     alpha=0.2,
                     linewidth=3,
                     linestyle='dashed',

code/requirements.txt

@@ -4,6 +4,13 @@ scipy>=0.14.0
 six
 enum34; python_version < '3.4'
 
+# following package is required to use elephant without installation
+tqdm
+# following packages are required to use odML without installation
+rdflib
+pyyaml
+lxml
+
 # Plotting module required for example scripts
 matplotlib>=3.0.3
 
@@ -18,3 +25,5 @@ matplotlib>=3.0.3
 # python-neo=0.9.0
 # elephant=0.9.0
 # odml=1.5.1
+
+