From 056a523b9458b77f13f0a3d54238bfd8868d6db2 Mon Sep 17 00:00:00 2001
From: Thijs <t.r.ruikes@uva.nl>
Date: Thu, 28 May 2020 10:34:29 +0200
Subject: [PATCH 1/6] Add info about inversion Commented not-used lines in
 `read_txt_header` In `read_txt_header` writes None to 'recording_opened' or
 'recording_closed' if either breaks

---
 neo/rawio/neuralynxrawio.py | 30 +++++++++++++++++++++---------
 1 file changed, 21 insertions(+), 9 deletions(-)

diff --git a/neo/rawio/neuralynxrawio.py b/neo/rawio/neuralynxrawio.py
index 6f6f12217..30c980c36 100644
--- a/neo/rawio/neuralynxrawio.py
+++ b/neo/rawio/neuralynxrawio.py
@@ -37,6 +37,8 @@ class NeuralynxRawIO(BaseRawIO):
     """"
     Class for reading dataset recorded by Neuralynx.
 
+    Data recorded with 'Invert' enabled is autmoatically stratified in 'read_ncs' <--- something like this
+
     Examples:
         >>> reader = NeuralynxRawIO(dirname='Cheetah_v5.5.1/original_data')
         >>> reader.parse_header()
@@ -672,23 +674,33 @@ def read_txt_header(filename):
     if old_date_format:
         datetime1_regex = r'## Time Opened \(m/d/y\): (?P<date>\S+)  \(h:m:s\.ms\) (?P<time>\S+)'
         datetime2_regex = r'## Time Closed \(m/d/y\): (?P<date>\S+)  \(h:m:s\.ms\) (?P<time>\S+)'
-        filename_regex = r'## File Name (?P<filename>\S+)'
+        # filename_regex = r'## File Name (?P<filename>\S+)'  # not used
         datetimeformat = '%m/%d/%Y %H:%M:%S.%f'
     else:
         datetime1_regex = r'-TimeCreated (?P<date>\S+) (?P<time>\S+)'
         datetime2_regex = r'-TimeClosed (?P<date>\S+) (?P<time>\S+)'
-        filename_regex = r'-OriginalFileName "?(?P<filename>\S+)"?'
+        # filename_regex = r'-OriginalFileName "?(?P<filename>\S+)"?'  # not used
         datetimeformat = '%Y/%m/%d %H:%M:%S'
 
-    original_filename = re.search(filename_regex, txt_header).groupdict()['filename']
+    # original_filename = re.search(filename_regex, txt_header).groupdict()['filename']  Not used
+
+    dt1 = re.search(datetime1_regex, txt_header)  # cant call groupdict if re.search fails
+    dt2 = re.search(datetime2_regex, txt_header)
 
-    dt1 = re.search(datetime1_regex, txt_header).groupdict()
-    dt2 = re.search(datetime2_regex, txt_header).groupdict()
+    # Using this at least the code won't break if the regexp fails
+    if dt1 is not None:
+        dt1 = dt1.groupdict()
+        info['recording_opened'] = datetime.datetime.strptime(
+            dt1['date'] + ' ' + dt1['time'], datetimeformat)
+    else:
+        info['recording_opened'] = None
 
-    info['recording_opened'] = datetime.datetime.strptime(
-        dt1['date'] + ' ' + dt1['time'], datetimeformat)
-    info['recording_closed'] = datetime.datetime.strptime(
-        dt2['date'] + ' ' + dt2['time'], datetimeformat)
+    if dt2 is not None:
+        dt2 = dt2.groupdict()
+        info['recording_closed'] = datetime.datetime.strptime(
+            dt2['date'] + ' ' + dt2['time'], datetimeformat)
+    else:
+        info['recording_closed'] = None
 
     return info
 

From 808056310f364e9ba90ba509d1f1fb914d85e041 Mon Sep 17 00:00:00 2001
From: Thijs <t.r.ruikes@uva.nl>
Date: Thu, 28 May 2020 19:18:17 +0200
Subject: [PATCH 2/6] Started adding patching functionality. Setting
 `align_and_patch_ncs_channels` to True show fix my problems, leaving it to
 false will return old behaviour.

Patching will make 1 segment out of the data, not recommended for data with large gaps.
---
 neo/rawio/neuralynxrawio.py | 363 +++++++++++++++++++-----------------
 1 file changed, 189 insertions(+), 174 deletions(-)

diff --git a/neo/rawio/neuralynxrawio.py b/neo/rawio/neuralynxrawio.py
index 30c980c36..238e4d9fc 100644
--- a/neo/rawio/neuralynxrawio.py
+++ b/neo/rawio/neuralynxrawio.py
@@ -67,6 +67,7 @@ def _parse_header(self):
         event_channels = []
 
         self.ncs_filenames = OrderedDict()  # (chan_name, chan_id): filename
+        self._ncs_info = OrderedDict()  # (chan_name, chan_id): info
         self.nse_ntt_filenames = OrderedDict()  # (chan_name, chan_id): filename
         self.nev_filenames = OrderedDict()  # chan_id: filename
 
@@ -115,30 +116,8 @@ def _parse_header(self):
                     sig_channels.append((chan_name, chan_id, info['sampling_rate'],
                                          'int16', units, gain, offset, group_id))
                     self.ncs_filenames[chan_uid] = filename
-                    keys = [
-                        'DspFilterDelay_µs',
-                        'recording_opened',
-                        'FileType',
-                        'DspDelayCompensation',
-                        'recording_closed',
-                        'DspLowCutFilterType',
-                        'HardwareSubSystemName',
-                        'DspLowCutNumTaps',
-                        'DSPLowCutFilterEnabled',
-                        'HardwareSubSystemType',
-                        'DspHighCutNumTaps',
-                        'ADMaxValue',
-                        'DspLowCutFrequency',
-                        'DSPHighCutFilterEnabled',
-                        'RecordSize',
-                        'InputRange',
-                        'DspHighCutFrequency',
-                        'input_inverted',
-                        'NumADChannels',
-                        'DspHighCutFilterType',
-                    ]
-                    d = {k: info[k] for k in keys if k in info}
-                    signal_annotations.append(d)
+                    signal_annotations.append(info)  # is there a reason to not append all information
+                    self._ncs_info[chan_uid] = info  # make info accessible to rest of class (during .read_ncs_files call)
 
                 elif ext in ('nse', 'ntt'):
                     # nse and ntt are pretty similar except for the wavform shape
@@ -204,83 +183,83 @@ def _parse_header(self):
             self._sigs_sampling_rate = sampling_rate[0]
 
         # read ncs files for gaps detection and nb_segment computation
-        self.read_ncs_files(self.ncs_filenames)
-
-        # timestamp limit in nev, nse
-        # so need to scan all spike and event to
-        ts0, ts1 = None, None
-        for _data_memmap in (self._spike_memmap, self._nev_memmap):
-            for _, data in _data_memmap.items():
-                ts = data['timestamp']
-                if ts.size == 0:
-                    continue
-                if ts0 is None:
-                    ts0 = ts[0]
-                    ts1 = ts[-1]
-                ts0 = min(ts0, ts[0])
-                ts1 = max(ts0, ts[-1])
-
-        if self._timestamp_limits is None:
-            # case  NO ncs but HAVE nev or nse
-            self._timestamp_limits = [(ts0, ts1)]
-            self._seg_t_starts = [ts0 / 1e6]
-            self._seg_t_stops = [ts1 / 1e6]
-            self.global_t_start = ts0 / 1e6
-            self.global_t_stop = ts1 / 1e6
-        elif ts0 is not None:
-            # case  HAVE ncs AND HAVE nev or nse
-            self.global_t_start = min(ts0 / 1e6, self._sigs_t_start[0])
-            self.global_t_stop = max(ts1 / 1e6, self._sigs_t_stop[-1])
-            self._seg_t_starts = list(self._sigs_t_start)
-            self._seg_t_starts[0] = self.global_t_start
-            self._seg_t_stops = list(self._sigs_t_stop)
-            self._seg_t_stops[-1] = self.global_t_stop
-        else:
-            # case HAVE ncs but  NO nev or nse
-            self._seg_t_starts = self._sigs_t_start
-            self._seg_t_stops = self._sigs_t_stop
-            self.global_t_start = self._sigs_t_start[0]
-            self.global_t_stop = self._sigs_t_stop[-1]
-
-        if self.keep_original_times:
-            self.global_t_stop = self.global_t_stop - self.global_t_start
-            self.global_t_start = 0
-
-        # fill into header dict
-        self.header = {}
-        self.header['nb_block'] = 1
-        self.header['nb_segment'] = [self._nb_segment]
-        self.header['signal_channels'] = sig_channels
-        self.header['unit_channels'] = unit_channels
-        self.header['event_channels'] = event_channels
-
-        # Annotations
-        self._generate_minimal_annotations()
-        bl_annotations = self.raw_annotations['blocks'][0]
-
-        for seg_index in range(self._nb_segment):
-            seg_annotations = bl_annotations['segments'][seg_index]
-
-            for c in range(sig_channels.size):
-                sig_ann = seg_annotations['signals'][c]
-                sig_ann.update(signal_annotations[c])
-
-            for c in range(unit_channels.size):
-                unit_ann = seg_annotations['units'][c]
-                unit_ann.update(unit_annotations[c])
-
-            for c in range(event_channels.size):
-                # annotations for channel events
-                event_id, ttl_input = self.internal_event_ids[c]
-                chan_id = event_channels[c]['id']
-
-                ev_ann = seg_annotations['events'][c]
-                ev_ann['file_origin'] = self.nev_filenames[chan_id]
-
-                # ~ ev_ann['marker_id'] =
-                # ~ ev_ann['nttl'] =
-                # ~ ev_ann['digital_marker'] =
-                # ~ ev_ann['analog_marker'] =
+        self.read_ncs_files()  # no need to pass class attributes to function
+        #
+        # # timestamp limit in nev, nse
+        # # so need to scan all spike and event to
+        # ts0, ts1 = None, None
+        # for _data_memmap in (self._spike_memmap, self._nev_memmap):
+        #     for _, data in _data_memmap.items():
+        #         ts = data['timestamp']
+        #         if ts.size == 0:
+        #             continue
+        #         if ts0 is None:
+        #             ts0 = ts[0]
+        #             ts1 = ts[-1]
+        #         ts0 = min(ts0, ts[0])
+        #         ts1 = max(ts1, ts[-1])  # TODO: shouldn't this be ts1
+        #
+        # if self._timestamp_limits is None:
+        #     # case  NO ncs but HAVE nev or nse
+        #     self._timestamp_limits = [(ts0, ts1)]
+        #     self._seg_t_starts = [ts0 / 1e6]
+        #     self._seg_t_stops = [ts1 / 1e6]
+        #     self.global_t_start = ts0 / 1e6
+        #     self.global_t_stop = ts1 / 1e6
+        # elif ts0 is not None:
+        #     # case  HAVE ncs AND HAVE nev or nse
+        #     self.global_t_start = min(ts0 / 1e6, self._sigs_t_start[0])
+        #     self.global_t_stop = max(ts1 / 1e6, self._sigs_t_stop[-1])
+        #     self._seg_t_starts = list(self._sigs_t_start)
+        #     self._seg_t_starts[0] = self.global_t_start
+        #     self._seg_t_stops = list(self._sigs_t_stop)
+        #     self._seg_t_stops[-1] = self.global_t_stop
+        # else:
+        #     # case HAVE ncs but  NO nev or nse
+        #     self._seg_t_starts = self._sigs_t_start
+        #     self._seg_t_stops = self._sigs_t_stop
+        #     self.global_t_start = self._sigs_t_start[0]
+        #     self.global_t_stop = self._sigs_t_stop[-1]
+        #
+        # if self.keep_original_times:
+        #     self.global_t_stop = self.global_t_stop - self.global_t_start
+        #     self.global_t_start = 0
+        #
+        # # fill into header dict
+        # self.header = {}
+        # self.header['nb_block'] = 1
+        # self.header['nb_segment'] = [self._nb_segment]
+        # self.header['signal_channels'] = sig_channels
+        # self.header['unit_channels'] = unit_channels
+        # self.header['event_channels'] = event_channels
+        #
+        # # Annotations
+        # self._generate_minimal_annotations()
+        # bl_annotations = self.raw_annotations['blocks'][0]
+        #
+        # for seg_index in range(self._nb_segment):
+        #     seg_annotations = bl_annotations['segments'][seg_index]
+        #
+        #     for c in range(sig_channels.size):
+        #         sig_ann = seg_annotations['signals'][c]
+        #         sig_ann.update(signal_annotations[c])
+        #
+        #     for c in range(unit_channels.size):
+        #         unit_ann = seg_annotations['units'][c]
+        #         unit_ann.update(unit_annotations[c])
+        #
+        #     for c in range(event_channels.size):
+        #         # annotations for channel events
+        #         event_id, ttl_input = self.internal_event_ids[c]
+        #         chan_id = event_channels[c]['id']
+        #
+        #         ev_ann = seg_annotations['events'][c]
+        #         ev_ann['file_origin'] = self.nev_filenames[chan_id]
+        #
+        #         # ~ ev_ann['marker_id'] =
+        #         # ~ ev_ann['nttl'] =
+        #         # ~ ev_ann['digital_marker'] =
+        #         # ~ ev_ann['analog_marker'] =
 
     def _segment_t_start(self, block_index, seg_index):
         return self._seg_t_starts[seg_index] - self.global_t_start
@@ -413,7 +392,7 @@ def _rescale_event_timestamp(self, event_timestamps, dtype):
         event_times -= self.global_t_start
         return event_times
 
-    def read_ncs_files(self, ncs_filenames):
+    def read_ncs_files(self, ncs_filenames=None):
         """
         Given a list of ncs files contrsuct:
             * self._sigs_memmap = [ {} for seg_index in range(self._nb_segment) ]
@@ -433,67 +412,96 @@ def read_ncs_files(self, ncs_filenames):
         gap_indexes can be given (when cached) to avoid full read.
 
         """
+        align_and_patch_ncs_channels = True
+
+        # Default behaviour: use self.ncs_filenames
+        if ncs_filenames is None:
+            ncs_filenames = self.ncs_filenames
+        else:
+            raise ValueError('Havent implemented custom ncs_filenames')  # TODO: remove ncs_filenames keyword?
+
+        # If no ncs_filenames are available
         if len(ncs_filenames) == 0:
             self._nb_segment = 1
             self._timestamp_limits = None
             return
 
+        ncs_info = OrderedDict(
+            DspFilterDelay_us=[],
+            DspDelayCompensation=[],
+            t0=[],
+            t1=[],
+            gap_pairs=[],
+            gap_pair_times=[],
+            data_size=[],
+            timestamp=[],  # only filled if not align_and_path
+        )
         good_delta = int(BLOCK_SIZE * 1e6 / self._sigs_sampling_rate)
-        chan_uid0 = list(ncs_filenames.keys())[0]
-        filename0 = ncs_filenames[chan_uid0]
-
-        data0 = np.memmap(filename0, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
-
-        gap_indexes = None
-        lost_indexes = None
-
-        if self.use_cache:
-            gap_indexes = self._cache.get('gap_indexes')
-            lost_indexes = self._cache.get('lost_indexes')
+        max_tolerance = 2
 
-        # detect gaps on first file
-        if (gap_indexes is None) or (lost_indexes is None):
+        # summary of features per channel
+        for chan_uid, filename in ncs_filenames.items():
+            if 'DspFilterDelay_µs' in self._ncs_info[chan_uid].keys():
+                ncs_info['DspFilterDelay_us'].append(self._ncs_info[chan_uid]['DspFilterDelay_µs'])
+            else:
+                ncs_info['DspFilterDelay_us'].append(None)
 
-            # this can be long!!!!
-            timestamps0 = data0['timestamp']
-            deltas0 = np.diff(timestamps0)
+            if 'DspDelayCompensation' in self._ncs_info[chan_uid].keys():
+                ncs_info['DspDelayCompensation'].append(self._ncs_info[chan_uid]['DspDelayCompensation'])
+            else:
+                ncs_info['DspDelayCompensation'].append(None)
 
-            # It should be that:
-            # gap_indexes, = np.nonzero(deltas0!=good_delta)
-            # but for a file I have found many deltas0==15999, 16000, 16001 (for sampling at 32000)
-            # I guess this is a round problem
-            # So this is the same with a tolerance of 1 or 2 ticks
-            max_tolerance = 2
-            mask = np.abs((deltas0 - good_delta).astype('int64')) > max_tolerance
+            data = np.memmap(filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
+            ncs_info['t0'].append(data['timestamp'][0])
+            ncs_info['t1'].append(data['timestamp'][-1])
+            ncs_info['data_size'].append(data.size)
 
+            timestamps = data['timestamp']
+            deltas = np.diff(timestamps)
+            mask = np.abs((deltas - good_delta).astype('int64')) > max_tolerance
             gap_indexes, = np.nonzero(mask)
+            lost_indexes, = np.nonzero(data['nb_valid'] < BLOCK_SIZE)
+            gap_candidates = np.unique([0]
+                                       + [data.size-1]
+                                       + (gap_indexes + 1).tolist()
+                                       + lost_indexes.tolist())  # linear
+
+            gap_pairs = np.vstack([gap_candidates[:-1], gap_candidates[1:]]).T  # 2D (n_segments, 2)
+
+            # construct proper gap ranges free of lost samples artifacts
+            minimal_segment_length = 1  # in blocks
+            goodpairs = np.diff(gap_pairs, 1).reshape(-1) > minimal_segment_length
+            gap_pairs = gap_pairs[goodpairs]  # ensures a segment is at least a block wide
+            ncs_info['gap_pairs'].append(gap_pairs)
+            gap_times = [[data['timestamp'][i0], data['timestamp'][i1]] for i0, i1 in gap_pairs]
+            ncs_info['gap_pair_times'].append(gap_times)
+
+            if not align_and_patch_ncs_channels:
+                ncs_info['timestamp'].append(data['timestamp'])
+                chan_uid0 = chan_uid
+                break
+
+        if np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]) or \
+                len(np.unique(ncs_info['DspFilterDelay_us'])) == 1:
+            # compensation enabled on all channels, all channels have same lag; the channels are aligned
+            ncs_t0 = min(ncs_info['t0'])
+            ncs_t1 = max(ncs_info['t1'])
+
+        elif not np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]) and \
+                len(np.unique(ncs_info['DspFilterDelay_us'])) > 1:
+            # no compesation enabled on some channels, different delays between channels; channels are not aligned
+            raise ValueError('to be implemented')
+
+        if align_and_patch_ncs_channels:
+            self._nb_segment = 1
+            self._sigs_memmap = [{}]
+            # TODO: make time axis (# expected datapoints, t0 and t1)
+            # TODO: check what channels need patching (len(gap_pairs) > 1)
+        else:
+            self._nb_segment = len(ncs_info['gap_pairs'][0])
+            self._sigs_memmap = [{} for seg_index in range(self._nb_segment)]
+            gap_pairs = ncs_info['gap_pairs'][0]
 
-            if self.use_cache:
-                self.add_in_cache(gap_indexes=gap_indexes)
-
-            # update for lost_indexes
-            # Sometimes NLX writes a faulty block, but it then validates how much samples it wrote
-            # the validation field is in delta0['nb_valid'], it should be equal to BLOCK_SIZE
-
-            lost_indexes, = np.nonzero(data0['nb_valid'] < BLOCK_SIZE)
-
-            if self.use_cache:
-                self.add_in_cache(lost_indexes=lost_indexes)
-
-        gap_candidates = np.unique([0]
-                                   + [data0.size]
-                                   + (gap_indexes + 1).tolist()
-                                   + lost_indexes.tolist())  # linear
-
-        gap_pairs = np.vstack([gap_candidates[:-1], gap_candidates[1:]]).T  # 2D (n_segments, 2)
-
-        # construct proper gap ranges free of lost samples artifacts
-        minimal_segment_length = 1  # in blocks
-        goodpairs = np.diff(gap_pairs, 1).reshape(-1) > minimal_segment_length
-        gap_pairs = gap_pairs[goodpairs]  # ensures a segment is at least a block wide
-
-        self._nb_segment = len(gap_pairs)
-        self._sigs_memmap = [{} for seg_index in range(self._nb_segment)]
         self._sigs_t_start = []
         self._sigs_t_stop = []
         self._sigs_length = []
@@ -503,29 +511,36 @@ def read_ncs_files(self, ncs_filenames):
         for chan_uid, ncs_filename in self.ncs_filenames.items():
 
             data = np.memmap(ncs_filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
-            assert data.size == data0.size, 'ncs files do not have the same data length'
-
-            for seg_index, (i0, i1) in enumerate(gap_pairs):
-
-                assert data[i0]['timestamp'] == data0[i0][
-                    'timestamp'], 'ncs files do not have the same gaps'
-                assert data[i1 - 1]['timestamp'] == data0[i1 - 1][
-                    'timestamp'], 'ncs files do not have the same gaps'
-
-                subdata = data[i0:i1]
-                self._sigs_memmap[seg_index][chan_uid] = subdata
-
-                if chan_uid == chan_uid0:
-                    ts0 = subdata[0]['timestamp']
-                    ts1 = subdata[-1]['timestamp'] \
-                            + np.uint64(BLOCK_SIZE / self._sigs_sampling_rate * 1e6)
-                    self._timestamp_limits.append((ts0, ts1))
-                    t_start = ts0 / 1e6
-                    self._sigs_t_start.append(t_start)
-                    t_stop = ts1 / 1e6
-                    self._sigs_t_stop.append(t_stop)
-                    length = subdata.size * BLOCK_SIZE
-                    self._sigs_length.append(length)
+
+            if not align_and_patch_ncs_channels:
+                assert data.size == ncs_info['data_size'][0], 'ncs files do not have the same data length'
+
+                for seg_index, (i0, i1) in enumerate(gap_pairs):
+
+                    assert data[i0]['timestamp'] == ncs_info['timestamp'][0][i0][
+                        'timestamp'], 'ncs files do not have the same gaps'
+                    assert data[i1 - 1]['timestamp'] == ncs_info['timestamp'][0][i1 - 1][
+                        'timestamp'], 'ncs files do not have the same gaps'
+
+                    subdata = data[i0:i1]
+                    self._sigs_memmap[seg_index][chan_uid] = subdata
+
+                    if chan_uid == chan_uid0:
+                        ts0 = subdata[0]['timestamp']
+                        ts1 = subdata[-1]['timestamp'] \
+                                + np.uint64(BLOCK_SIZE / self._sigs_sampling_rate * 1e6)
+                        self._timestamp_limits.append((ts0, ts1))
+                        t_start = ts0 / 1e6
+                        self._sigs_t_start.append(t_start)
+                        t_stop = ts1 / 1e6
+                        self._sigs_t_stop.append(t_stop)
+                        length = subdata.size * BLOCK_SIZE
+                        self._sigs_length.append(length)
+
+            else:
+                # TODO: find missing data
+                # TODO: fill missing data
+                # TODO: write to new file
 
 
 # Keys funcitons

From 9694a8bd3f5e998954817eec7d971ab7d38e0099 Mon Sep 17 00:00:00 2001
From: Thijs <t.r.ruikes@uva.nl>
Date: Fri, 29 May 2020 12:00:56 +0200
Subject: [PATCH 3/6] Moved all patching funtionality to `patch_ncs_files`.
 This function will scan a directory for ncs files, detect misalignements/gaps
 in data and patches accordingly. Output is written in a new directory.

`read_ncs_files` is now unchanged
---
 neo/rawio/neuralynxrawio.py | 285 ++++++++++++++++++++++--------------
 1 file changed, 174 insertions(+), 111 deletions(-)

diff --git a/neo/rawio/neuralynxrawio.py b/neo/rawio/neuralynxrawio.py
index 238e4d9fc..73e4bf630 100644
--- a/neo/rawio/neuralynxrawio.py
+++ b/neo/rawio/neuralynxrawio.py
@@ -117,7 +117,6 @@ def _parse_header(self):
                                          'int16', units, gain, offset, group_id))
                     self.ncs_filenames[chan_uid] = filename
                     signal_annotations.append(info)  # is there a reason to not append all information
-                    self._ncs_info[chan_uid] = info  # make info accessible to rest of class (during .read_ncs_files call)
 
                 elif ext in ('nse', 'ntt'):
                     # nse and ntt are pretty similar except for the wavform shape
@@ -392,7 +391,104 @@ def _rescale_event_timestamp(self, event_timestamps, dtype):
         event_times -= self.global_t_start
         return event_times
 
-    def read_ncs_files(self, ncs_filenames=None):
+    @staticmethod
+    def patch_ncs_files(dirname):
+        """
+        Scan all ncs files in dir for broken/misaligned data and write the patched
+        data to patchdir.This directory can then be used by neuralynx io.
+
+        For now assumes the data should be 1 continous segment
+        """
+
+        # This should be callable without _parse_header (because it breaks), so
+        # I'll scan for ncs files here aswell
+        sampling_rate = 32000
+        verbose = True
+        ncs_info = OrderedDict(
+            filenames=[],
+            chan_uids=[],
+            chan_names=[],
+            DspFilterDelay_us=[],
+            DspDelayCompensation=[],
+            t0=[],
+            t1=[],
+            gap_pairs=[],
+            gap_pair_times=[],
+            data_size=[],
+            timestamp=[],  # only filled if not align_and_path
+            sampling_rate=[],
+        )
+        good_delta = int(BLOCK_SIZE * 1e6 / sampling_rate)
+        max_tolerance = 2
+
+        if verbose:
+            print('Scanning for ncs files in :\n ', dirname)
+
+        # Scan directory for ncs files, read headers to find inconsistencies
+        # in data
+        for filename in sorted(os.listdir(dirname)):
+
+            fullname = os.path.join(dirname, filename)
+            ncs_info['filenames'].append(fullname)
+            if (os.path.getsize(fullname) <= HEADER_SIZE) or not filename.endswith('ncs'):
+                continue
+
+            if verbose:
+                print(f'\t{filename}')
+
+            info = read_txt_header(fullname)
+            chan_names = info['channel_names']
+            chan_ids = info['channel_ids']
+            info['filename'] = fullname
+            for idx, chan_id in enumerate(chan_ids):
+                chan_name = chan_names[idx]
+                chan_uid = (chan_name, chan_id)
+
+                ncs_info['chan_uids'].append(chan_uid)
+                ncs_info['chan_names'].append(chan_name)
+                ncs_info['sampling_rate'].append(info['sampling_rate'])
+
+                if 'DspFilterDelay_µs' in info.keys():
+                    ncs_info['DspFilterDelay_us'].append(info['DspFilterDelay_µs'])
+                else:
+                    ncs_info['DspFilterDelay_us'].append(None)
+
+                if 'DspDelayCompensation' in ncs_info.keys():
+                    ncs_info['DspDelayCompensation'].append(info['DspDelayCompensation'])
+                else:
+                    ncs_info['DspDelayCompensation'].append(None)
+
+                data = get_ncs_memmap(fullname)
+                ncs_info['t0'].append(data['timestamp'][0])
+                ncs_info['t1'].append(data['timestamp'][-1])
+                ncs_info['data_size'].append(data.size)
+
+                # check good_delta is indeed good delta for this channel
+                assert good_delta in np.unique(np.diff(data['timestamp'])), 'good delta does not match for channel'
+
+        # assert len(np.unique(ncs_info['DspFilterDelay_us'])) == 1 or np.all([ddc == 'Enabled' for ])
+        if len(np.unique(ncs_info['DspFilterDelay_us'])) > 0 and not np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]):
+            raise ValueError('some channels have dsp filtering, but lag compensation is disabled')
+        elif len(np.unique(ncs_info['DspFilterDelay_us'])) > 0 and np.all(
+                [ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]):
+            print('Channels have different lengths due to Dsp filters, padding shorter channels')
+
+        assert len(np.unique(ncs_info['sampling_rate'])) == 1, f'some chans have different samplingrate: {ncs_info["sampling_rate"]}'
+
+        t0_all = np.min(ncs_info['t0'])
+        t1_all = np.max(ncs_info['t1'])
+        rng = t1_all - t0_all
+        tsinterp = np.arange(0, rng, good_delta, dtype='uint64')
+
+        for i, chan_uid in enumerate(ncs_info['chan_uids']):
+            data = get_ncs_memmap(ncs_info['filenames'][i])
+            gaps = np.where(np.diff(data['timestamp']) > good_delta)[0]
+
+            # TODO: open memmap to write patched data
+            # get offset
+            # fill gaps
+
+    def read_ncs_files(self, ncs_filenames):
         """
         Given a list of ncs files contrsuct:
             * self._sigs_memmap = [ {} for seg_index in range(self._nb_segment) ]
@@ -412,136 +508,99 @@ def read_ncs_files(self, ncs_filenames=None):
         gap_indexes can be given (when cached) to avoid full read.
 
         """
-        align_and_patch_ncs_channels = True
-
-        # Default behaviour: use self.ncs_filenames
-        if ncs_filenames is None:
-            ncs_filenames = self.ncs_filenames
-        else:
-            raise ValueError('Havent implemented custom ncs_filenames')  # TODO: remove ncs_filenames keyword?
-
-        # If no ncs_filenames are available
         if len(ncs_filenames) == 0:
             self._nb_segment = 1
             self._timestamp_limits = None
             return
 
-        ncs_info = OrderedDict(
-            DspFilterDelay_us=[],
-            DspDelayCompensation=[],
-            t0=[],
-            t1=[],
-            gap_pairs=[],
-            gap_pair_times=[],
-            data_size=[],
-            timestamp=[],  # only filled if not align_and_path
-        )
         good_delta = int(BLOCK_SIZE * 1e6 / self._sigs_sampling_rate)
-        max_tolerance = 2
+        chan_uid0 = list(ncs_filenames.keys())[0]
+        filename0 = ncs_filenames[chan_uid0]
 
-        # summary of features per channel
-        for chan_uid, filename in ncs_filenames.items():
-            if 'DspFilterDelay_µs' in self._ncs_info[chan_uid].keys():
-                ncs_info['DspFilterDelay_us'].append(self._ncs_info[chan_uid]['DspFilterDelay_µs'])
-            else:
-                ncs_info['DspFilterDelay_us'].append(None)
+        data0 = np.memmap(filename0, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
 
-            if 'DspDelayCompensation' in self._ncs_info[chan_uid].keys():
-                ncs_info['DspDelayCompensation'].append(self._ncs_info[chan_uid]['DspDelayCompensation'])
-            else:
-                ncs_info['DspDelayCompensation'].append(None)
+        gap_indexes = None
+        lost_indexes = None
 
-            data = np.memmap(filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
-            ncs_info['t0'].append(data['timestamp'][0])
-            ncs_info['t1'].append(data['timestamp'][-1])
-            ncs_info['data_size'].append(data.size)
+        if self.use_cache:
+            gap_indexes = self._cache.get('gap_indexes')
+            lost_indexes = self._cache.get('lost_indexes')
+
+        # detect gaps on first file
+        if (gap_indexes is None) or (lost_indexes is None):
+
+            # this can be long!!!!
+            timestamps0 = data0['timestamp']
+            deltas0 = np.diff(timestamps0)
+
+            # It should be that:
+            # gap_indexes, = np.nonzero(deltas0!=good_delta)
+            # but for a file I have found many deltas0==15999, 16000, 16001 (for sampling at 32000)
+            # I guess this is a round problem
+            # So this is the same with a tolerance of 1 or 2 ticks
+            max_tolerance = 2
+            mask = np.abs((deltas0 - good_delta).astype('int64')) > max_tolerance
 
-            timestamps = data['timestamp']
-            deltas = np.diff(timestamps)
-            mask = np.abs((deltas - good_delta).astype('int64')) > max_tolerance
             gap_indexes, = np.nonzero(mask)
-            lost_indexes, = np.nonzero(data['nb_valid'] < BLOCK_SIZE)
-            gap_candidates = np.unique([0]
-                                       + [data.size-1]
-                                       + (gap_indexes + 1).tolist()
-                                       + lost_indexes.tolist())  # linear
-
-            gap_pairs = np.vstack([gap_candidates[:-1], gap_candidates[1:]]).T  # 2D (n_segments, 2)
-
-            # construct proper gap ranges free of lost samples artifacts
-            minimal_segment_length = 1  # in blocks
-            goodpairs = np.diff(gap_pairs, 1).reshape(-1) > minimal_segment_length
-            gap_pairs = gap_pairs[goodpairs]  # ensures a segment is at least a block wide
-            ncs_info['gap_pairs'].append(gap_pairs)
-            gap_times = [[data['timestamp'][i0], data['timestamp'][i1]] for i0, i1 in gap_pairs]
-            ncs_info['gap_pair_times'].append(gap_times)
-
-            if not align_and_patch_ncs_channels:
-                ncs_info['timestamp'].append(data['timestamp'])
-                chan_uid0 = chan_uid
-                break
-
-        if np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]) or \
-                len(np.unique(ncs_info['DspFilterDelay_us'])) == 1:
-            # compensation enabled on all channels, all channels have same lag; the channels are aligned
-            ncs_t0 = min(ncs_info['t0'])
-            ncs_t1 = max(ncs_info['t1'])
-
-        elif not np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]) and \
-                len(np.unique(ncs_info['DspFilterDelay_us'])) > 1:
-            # no compesation enabled on some channels, different delays between channels; channels are not aligned
-            raise ValueError('to be implemented')
-
-        if align_and_patch_ncs_channels:
-            self._nb_segment = 1
-            self._sigs_memmap = [{}]
-            # TODO: make time axis (# expected datapoints, t0 and t1)
-            # TODO: check what channels need patching (len(gap_pairs) > 1)
-        else:
-            self._nb_segment = len(ncs_info['gap_pairs'][0])
-            self._sigs_memmap = [{} for seg_index in range(self._nb_segment)]
-            gap_pairs = ncs_info['gap_pairs'][0]
 
-        self._sigs_t_start = []
-        self._sigs_t_stop = []
-        self._sigs_length = []
-        self._timestamp_limits = []
+            if self.use_cache:
+                self.add_in_cache(gap_indexes=gap_indexes)
 
-        # create segment with subdata block/t_start/t_stop/length
-        for chan_uid, ncs_filename in self.ncs_filenames.items():
+            # update for lost_indexes
+            # Sometimes NLX writes a faulty block, but it then validates how much samples it wrote
+            # the validation field is in delta0['nb_valid'], it should be equal to BLOCK_SIZE
 
-            data = np.memmap(ncs_filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
+            lost_indexes, = np.nonzero(data0['nb_valid'] < BLOCK_SIZE)
 
-            if not align_and_patch_ncs_channels:
-                assert data.size == ncs_info['data_size'][0], 'ncs files do not have the same data length'
+            if self.use_cache:
+                self.add_in_cache(lost_indexes=lost_indexes)
 
-                for seg_index, (i0, i1) in enumerate(gap_pairs):
+        gap_candidates = np.unique([0]
+                                   + [data0.size]
+                                   + (gap_indexes + 1).tolist()
+                                   + lost_indexes.tolist())  # linear
 
-                    assert data[i0]['timestamp'] == ncs_info['timestamp'][0][i0][
-                        'timestamp'], 'ncs files do not have the same gaps'
-                    assert data[i1 - 1]['timestamp'] == ncs_info['timestamp'][0][i1 - 1][
-                        'timestamp'], 'ncs files do not have the same gaps'
+        gap_pairs = np.vstack([gap_candidates[:-1], gap_candidates[1:]]).T  # 2D (n_segments, 2)
 
-                    subdata = data[i0:i1]
-                    self._sigs_memmap[seg_index][chan_uid] = subdata
+        # construct proper gap ranges free of lost samples artifacts
+        minimal_segment_length = 1  # in blocks
+        goodpairs = np.diff(gap_pairs, 1).reshape(-1) > minimal_segment_length
+        gap_pairs = gap_pairs[goodpairs]  # ensures a segment is at least a block wide
 
-                    if chan_uid == chan_uid0:
-                        ts0 = subdata[0]['timestamp']
-                        ts1 = subdata[-1]['timestamp'] \
-                                + np.uint64(BLOCK_SIZE / self._sigs_sampling_rate * 1e6)
-                        self._timestamp_limits.append((ts0, ts1))
-                        t_start = ts0 / 1e6
-                        self._sigs_t_start.append(t_start)
-                        t_stop = ts1 / 1e6
-                        self._sigs_t_stop.append(t_stop)
-                        length = subdata.size * BLOCK_SIZE
-                        self._sigs_length.append(length)
+        self._nb_segment = len(gap_pairs)
+        self._sigs_memmap = [{} for seg_index in range(self._nb_segment)]
+        self._sigs_t_start = []
+        self._sigs_t_stop = []
+        self._sigs_length = []
+        self._timestamp_limits = []
 
-            else:
-                # TODO: find missing data
-                # TODO: fill missing data
-                # TODO: write to new file
+        # create segment with subdata block/t_start/t_stop/length
+        for chan_uid, ncs_filename in self.ncs_filenames.items():
 
+            data = np.memmap(ncs_filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
+            assert data.size == data0.size, 'ncs files do not have the same data length'
+
+            for seg_index, (i0, i1) in enumerate(gap_pairs):
+
+                assert data[i0]['timestamp'] == data0[i0][
+                    'timestamp'], 'ncs files do not have the same gaps'
+                assert data[i1 - 1]['timestamp'] == data0[i1 - 1][
+                    'timestamp'], 'ncs files do not have the same gaps'
+
+                subdata = data[i0:i1]
+                self._sigs_memmap[seg_index][chan_uid] = subdata
+
+                if chan_uid == chan_uid0:
+                    ts0 = subdata[0]['timestamp']
+                    ts1 = subdata[-1]['timestamp'] \
+                          + np.uint64(BLOCK_SIZE / self._sigs_sampling_rate * 1e6)
+                    self._timestamp_limits.append((ts0, ts1))
+                    t_start = ts0 / 1e6
+                    self._sigs_t_start.append(t_start)
+                    t_stop = ts1 / 1e6
+                    self._sigs_t_stop.append(t_stop)
+                    length = subdata.size * BLOCK_SIZE
+                    self._sigs_length.append(length)
 
 # Keys funcitons
 def _to_bool(txt):
@@ -762,3 +821,7 @@ def get_nse_or_ntt_dtype(info, ext):
         dtype += [('samples', 'int16', (nb_sample, nb_chan))]
 
     return dtype
+
+
+def get_ncs_memmap(filename):
+    return np.memmap(filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
\ No newline at end of file

From 34dd9c10ce58131378da2c2f60de3f32d0fcbecc Mon Sep 17 00:00:00 2001
From: Thijs <t.r.ruikes@uva.nl>
Date: Sat, 30 May 2020 20:54:42 +0200
Subject: [PATCH 4/6] Finished first draft of path_ncs_files

---
 neo/rawio/baserawio.py      |   2 +-
 neo/rawio/neuralynxrawio.py | 373 +++++++++++++++++++++++++-----------
 2 files changed, 263 insertions(+), 112 deletions(-)

diff --git a/neo/rawio/baserawio.py b/neo/rawio/baserawio.py
index fa0e98750..291256c09 100644
--- a/neo/rawio/baserawio.py
+++ b/neo/rawio/baserawio.py
@@ -147,7 +147,7 @@ def parse_header(self):
 
         """
         self._parse_header()
-        self._group_signal_channel_characteristics()
+        # self._group_signal_channel_characteristics()
 
     def source_name(self):
         """Return fancy name of file source"""
diff --git a/neo/rawio/neuralynxrawio.py b/neo/rawio/neuralynxrawio.py
index 73e4bf630..df961ccd5 100644
--- a/neo/rawio/neuralynxrawio.py
+++ b/neo/rawio/neuralynxrawio.py
@@ -32,7 +32,6 @@
 BLOCK_SIZE = 512  # nb sample per signal block
 HEADER_SIZE = 2 ** 14  # file have a txt header of 16kB
 
-
 class NeuralynxRawIO(BaseRawIO):
     """"
     Class for reading dataset recorded by Neuralynx.
@@ -182,83 +181,83 @@ def _parse_header(self):
             self._sigs_sampling_rate = sampling_rate[0]
 
         # read ncs files for gaps detection and nb_segment computation
-        self.read_ncs_files()  # no need to pass class attributes to function
-        #
-        # # timestamp limit in nev, nse
-        # # so need to scan all spike and event to
-        # ts0, ts1 = None, None
-        # for _data_memmap in (self._spike_memmap, self._nev_memmap):
-        #     for _, data in _data_memmap.items():
-        #         ts = data['timestamp']
-        #         if ts.size == 0:
-        #             continue
-        #         if ts0 is None:
-        #             ts0 = ts[0]
-        #             ts1 = ts[-1]
-        #         ts0 = min(ts0, ts[0])
-        #         ts1 = max(ts1, ts[-1])  # TODO: shouldn't this be ts1
-        #
-        # if self._timestamp_limits is None:
-        #     # case  NO ncs but HAVE nev or nse
-        #     self._timestamp_limits = [(ts0, ts1)]
-        #     self._seg_t_starts = [ts0 / 1e6]
-        #     self._seg_t_stops = [ts1 / 1e6]
-        #     self.global_t_start = ts0 / 1e6
-        #     self.global_t_stop = ts1 / 1e6
-        # elif ts0 is not None:
-        #     # case  HAVE ncs AND HAVE nev or nse
-        #     self.global_t_start = min(ts0 / 1e6, self._sigs_t_start[0])
-        #     self.global_t_stop = max(ts1 / 1e6, self._sigs_t_stop[-1])
-        #     self._seg_t_starts = list(self._sigs_t_start)
-        #     self._seg_t_starts[0] = self.global_t_start
-        #     self._seg_t_stops = list(self._sigs_t_stop)
-        #     self._seg_t_stops[-1] = self.global_t_stop
-        # else:
-        #     # case HAVE ncs but  NO nev or nse
-        #     self._seg_t_starts = self._sigs_t_start
-        #     self._seg_t_stops = self._sigs_t_stop
-        #     self.global_t_start = self._sigs_t_start[0]
-        #     self.global_t_stop = self._sigs_t_stop[-1]
-        #
-        # if self.keep_original_times:
-        #     self.global_t_stop = self.global_t_stop - self.global_t_start
-        #     self.global_t_start = 0
-        #
-        # # fill into header dict
-        # self.header = {}
-        # self.header['nb_block'] = 1
-        # self.header['nb_segment'] = [self._nb_segment]
-        # self.header['signal_channels'] = sig_channels
-        # self.header['unit_channels'] = unit_channels
-        # self.header['event_channels'] = event_channels
-        #
-        # # Annotations
-        # self._generate_minimal_annotations()
-        # bl_annotations = self.raw_annotations['blocks'][0]
-        #
-        # for seg_index in range(self._nb_segment):
-        #     seg_annotations = bl_annotations['segments'][seg_index]
-        #
-        #     for c in range(sig_channels.size):
-        #         sig_ann = seg_annotations['signals'][c]
-        #         sig_ann.update(signal_annotations[c])
-        #
-        #     for c in range(unit_channels.size):
-        #         unit_ann = seg_annotations['units'][c]
-        #         unit_ann.update(unit_annotations[c])
-        #
-        #     for c in range(event_channels.size):
-        #         # annotations for channel events
-        #         event_id, ttl_input = self.internal_event_ids[c]
-        #         chan_id = event_channels[c]['id']
-        #
-        #         ev_ann = seg_annotations['events'][c]
-        #         ev_ann['file_origin'] = self.nev_filenames[chan_id]
-        #
-        #         # ~ ev_ann['marker_id'] =
-        #         # ~ ev_ann['nttl'] =
-        #         # ~ ev_ann['digital_marker'] =
-        #         # ~ ev_ann['analog_marker'] =
+        self.read_ncs_files(ncs_filenames=self.ncs_filenames)  # no need to pass class attributes to function
+
+        # timestamp limit in nev, nse
+        # so need to scan all spike and event to
+        ts0, ts1 = None, None
+        for _data_memmap in (self._spike_memmap, self._nev_memmap):
+            for _, data in _data_memmap.items():
+                ts = data['timestamp']
+                if ts.size == 0:
+                    continue
+                if ts0 is None:
+                    ts0 = ts[0]
+                    ts1 = ts[-1]
+                ts0 = min(ts0, ts[0])
+                ts1 = max(ts1, ts[-1])  # TODO: shouldn't this be ts1
+
+        if self._timestamp_limits is None:
+            # case  NO ncs but HAVE nev or nse
+            self._timestamp_limits = [(ts0, ts1)]
+            self._seg_t_starts = [ts0 / 1e6]
+            self._seg_t_stops = [ts1 / 1e6]
+            self.global_t_start = ts0 / 1e6
+            self.global_t_stop = ts1 / 1e6
+        elif ts0 is not None:
+            # case  HAVE ncs AND HAVE nev or nse
+            self.global_t_start = min(ts0 / 1e6, self._sigs_t_start[0])
+            self.global_t_stop = max(ts1 / 1e6, self._sigs_t_stop[-1])
+            self._seg_t_starts = list(self._sigs_t_start)
+            self._seg_t_starts[0] = self.global_t_start
+            self._seg_t_stops = list(self._sigs_t_stop)
+            self._seg_t_stops[-1] = self.global_t_stop
+        else:
+            # case HAVE ncs but  NO nev or nse
+            self._seg_t_starts = self._sigs_t_start
+            self._seg_t_stops = self._sigs_t_stop
+            self.global_t_start = self._sigs_t_start[0]
+            self.global_t_stop = self._sigs_t_stop[-1]
+
+        if self.keep_original_times:
+            self.global_t_stop = self.global_t_stop - self.global_t_start
+            self.global_t_start = 0
+
+        # fill into header dict
+        self.header = {}
+        self.header['nb_block'] = 1
+        self.header['nb_segment'] = [self._nb_segment]
+        self.header['signal_channels'] = sig_channels
+        self.header['unit_channels'] = unit_channels
+        self.header['event_channels'] = event_channels
+
+        # Annotations
+        self._generate_minimal_annotations()
+        bl_annotations = self.raw_annotations['blocks'][0]
+
+        for seg_index in range(self._nb_segment):
+            seg_annotations = bl_annotations['segments'][seg_index]
+
+            for c in range(sig_channels.size):
+                sig_ann = seg_annotations['signals'][c]
+                sig_ann.update(signal_annotations[c])
+
+            for c in range(unit_channels.size):
+                unit_ann = seg_annotations['units'][c]
+                unit_ann.update(unit_annotations[c])
+
+            for c in range(event_channels.size):
+                # annotations for channel events
+                event_id, ttl_input = self.internal_event_ids[c]
+                chan_id = event_channels[c]['id']
+
+                ev_ann = seg_annotations['events'][c]
+                ev_ann['file_origin'] = self.nev_filenames[chan_id]
+
+                # ~ ev_ann['marker_id'] =
+                # ~ ev_ann['nttl'] =
+                # ~ ev_ann['digital_marker'] =
+                # ~ ev_ann['analog_marker'] =
 
     def _segment_t_start(self, block_index, seg_index):
         return self._seg_t_starts[seg_index] - self.global_t_start
@@ -392,20 +391,41 @@ def _rescale_event_timestamp(self, event_timestamps, dtype):
         return event_times
 
     @staticmethod
-    def patch_ncs_files(dirname):
+    def patch_ncs_files(dirname, patchdir=None, interpolate=False, verbose=True, scan_only=False):
         """
-        Scan all ncs files in dir for broken/misaligned data and write the patched
-        data to patchdir.This directory can then be used by neuralynx io.
+        Scan a directory containing .ncs files, and detect misaligned data. Then patches the data
+        and saves the fixed .ncs files to another directory.
+
+
+        Parameters
+        ----------
+        dirname: [str]
+            full path containing ncs files scan
+        patchdir: [str] (default: None)
+            when patching data, patched ncs files are written to patchdir
+        interpolate: [bool] (default: False)
+            doesn't do anything
+        verbose: [bool] (default=True)
+            if true, will print out some extra info
+        scan_only: [bool] (default=False)
+            if true, will only scan the data, and list faulty files (will not patch),
+            does not require 'patchdir'.
 
-        For now assumes the data should be 1 continous segment
         """
 
-        # This should be callable without _parse_header (because it breaks), so
-        # I'll scan for ncs files here aswell
-        sampling_rate = 32000
-        verbose = True
+        time_margin = 1  # [ms] allow std of time_margin for time starts and stops across ncs files
+        max_gap_duration_to_fill = 60  # [s] maximum gap size in data to fill
+
+        # --------------------------- SCAN NCS FOR MISALIGNMENTS ------------------------------------
+        # Scan all files in directory and list:
+        #   -gaps due to ringbuffer errors
+        #   -lag by dsp filters
+        #   -dsp comensation
+
+        # Summarize in 'ncs_info' (dict with lists, one entry per file
         ncs_info = OrderedDict(
             filenames=[],
+            fullnames=[],
             chan_uids=[],
             chan_names=[],
             DspFilterDelay_us=[],
@@ -418,28 +438,28 @@ def patch_ncs_files(dirname):
             timestamp=[],  # only filled if not align_and_path
             sampling_rate=[],
         )
-        good_delta = int(BLOCK_SIZE * 1e6 / sampling_rate)
-        max_tolerance = 2
 
         if verbose:
             print('Scanning for ncs files in :\n ', dirname)
+            print('Detected files:')
 
-        # Scan directory for ncs files, read headers to find inconsistencies
-        # in data
+        found_missing_dsp_compensation = False
         for filename in sorted(os.listdir(dirname)):
 
             fullname = os.path.join(dirname, filename)
-            ncs_info['filenames'].append(fullname)
             if (os.path.getsize(fullname) <= HEADER_SIZE) or not filename.endswith('ncs'):
                 continue
 
+            ncs_info['filenames'].append(filename)
+            ncs_info['fullnames'].append(fullname)
+
             if verbose:
-                print(f'\t{filename}')
+                print(f'\t{filename}:')
 
             info = read_txt_header(fullname)
             chan_names = info['channel_names']
             chan_ids = info['channel_ids']
-            info['filename'] = fullname
+
             for idx, chan_id in enumerate(chan_ids):
                 chan_name = chan_names[idx]
                 chan_uid = (chan_name, chan_id)
@@ -458,35 +478,153 @@ def patch_ncs_files(dirname):
                 else:
                     ncs_info['DspDelayCompensation'].append(None)
 
-                data = get_ncs_memmap(fullname)
+                data = read_ncs_memmap(fullname)
                 ncs_info['t0'].append(data['timestamp'][0])
                 ncs_info['t1'].append(data['timestamp'][-1])
                 ncs_info['data_size'].append(data.size)
 
                 # check good_delta is indeed good delta for this channel
+
+                good_delta = int(BLOCK_SIZE * 1e6 / info['sampling_rate'])
                 assert good_delta in np.unique(np.diff(data['timestamp'])), 'good delta does not match for channel'
 
-        # assert len(np.unique(ncs_info['DspFilterDelay_us'])) == 1 or np.all([ddc == 'Enabled' for ])
-        if len(np.unique(ncs_info['DspFilterDelay_us'])) > 0 and not np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]):
-            raise ValueError('some channels have dsp filtering, but lag compensation is disabled')
-        elif len(np.unique(ncs_info['DspFilterDelay_us'])) > 0 and np.all(
-                [ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]):
-            print('Channels have different lengths due to Dsp filters, padding shorter channels')
+                if ncs_info['DspFilterDelay_us'][-1] is not None and int(ncs_info['DspFilterDelay_us'][-1]) > 0 and \
+                    ncs_info['DspDelayCompensation'] != 'Enabled':
+                    print(f'WARNING: {filename} has dsplag ({ncs_info["DspFilterDelay_us"][-1]}) and no delay compensation!')
+                    found_missing_dsp_compensation = True
+
+                if verbose and scan_only:
+                    print(f'\t\tchan id: {chan_uid[1]}')
+                    print(f'\t\tdsp delay: {ncs_info["DspFilterDelay_us"][-1]}, dsp compensation: {ncs_info["DspDelayCompensation"][-1]}')
+                    print(f'\t\tt start: {ncs_info["t0"][-1]}, t stop: {ncs_info["t1"][-1]}')
 
-        assert len(np.unique(ncs_info['sampling_rate'])) == 1, f'some chans have different samplingrate: {ncs_info["sampling_rate"]}'
+        # --------------------------- CHECK IF DATA NEEDS PATCHING ------------------------------------
 
         t0_all = np.min(ncs_info['t0'])
         t1_all = np.max(ncs_info['t1'])
-        rng = t1_all - t0_all
-        tsinterp = np.arange(0, rng, good_delta, dtype='uint64')
 
+        print('\nDetected errors:')
+
+        needs_patching = False
+        if np.any((ncs_info['t0'] - t0_all) > time_margin) or \
+                np.any((t1_all - ncs_info['t1'] > time_margin)):
+            print('\t-misalignment of timestamps')
+            needs_patching = True
+
+        if len(np.unique(ncs_info['data_size'])) > 1:
+            print('\t-Records do not have the same nr of samples')
+            needs_patching = True
+
+        if scan_only:
+            return
+
+        if not needs_patching and not found_missing_dsp_compensation:
+            print('\nNo errors found, no need to patch this data!')
+            return
+        elif not needs_patching and found_missing_dsp_compensation:
+            print('\nNo broken data, but some channels have no Dsp delay compensation')
+            return
+
+        # --------------------------- PATCH DATA ------------------------------------
+        assert patchdir is not None, 'To patch data, provide patchdir keyword'
+        if not os.path.isdir(patchdir):
+            print(f'Patchdir did not exist, making directory:\n{patchdir}')
+            os.makedirs(patchdir)
+
+        # Check if there is no misalginment due to dsp filters and no dsp compensation
+        assert len(np.unique(ncs_info['DspFilterDelay_us'])) == 1 or \
+               np.all([ddc == 'Enabled' for ddc in ncs_info['DspDelayCompensation']]), \
+            'some channels have dsp filtering, but lag compensation is disabled'
+
+        time_range = t1_all - t0_all
+
+        print(f'Saving patched data in: {patchdir}')
+        print('Patching:')
         for i, chan_uid in enumerate(ncs_info['chan_uids']):
-            data = get_ncs_memmap(ncs_info['filenames'][i])
-            gaps = np.where(np.diff(data['timestamp']) > good_delta)[0]
+            print(f'\t-{ncs_info["filenames"][i]}')
+
+            # Infile and outfile
+            ncs_in_name = ncs_info['fullnames'][i]
+            ncs_out_name = os.path.join(patchdir, ncs_info['filenames'][i])
+
+            # Load broken data
+            data = read_ncs_memmap(ncs_info['fullnames'][i])
+
+            good_delta = int(BLOCK_SIZE * 1e6 / ncs_info['sampling_rate'][i])
+            ts_interpolated = np.arange(0, time_range, good_delta, dtype='uint64')
+            chan_id = np.ones(ts_interpolated.shape, dtype='uint32')  # read per ncs file
+            sample_rate = np.ones(ts_interpolated.shape, dtype='uint32')  # read per ncs file
+            nb_valid = np.ones(ts_interpolated.shape, dtype='uint32') * 512  # hardcode as we patch here
+
+            # Fil constants
+            chan_ids_in_ncs = np.unique(data['channel_id'])
+            assert len(chan_ids_in_ncs) == 1, f'found multiple channel ids in ncs: {chan_ids_in_ncs}'
+            chan_id[:] = chan_ids_in_ncs[0]
+
+            sample_rates_in_ncs = np.unique(data['sample_rate'])
+            assert len(sample_rates_in_ncs == 1), f'found multiple sample rates in ncs: {sample_rate}'
+            sample_rate[:] = sample_rates_in_ncs[0]
 
-            # TODO: open memmap to write patched data
-            # get offset
-            # fill gaps
+            sample_dt = 1e6/sample_rates_in_ncs[0]  # [us] time between two sample points
+
+            # Patch data
+            # Make flat array of size (nsamples,) to fill data in (more convienent due to shifting of blocks)
+            # Have fill value default to signal mean
+            records = np.zeros((len(ts_interpolated) * 512), dtype='int16') * np.mean(data['samples'])
+
+            # Get times of ncs file, normalized to common t0
+            ts_norm = data['timestamp'] - t0_all
+
+            assert len(np.where(np.diff(data['timestamp']) < good_delta)[0]) == 0, 'this shouldnt happen'
+
+            # gap indices are detected by dt > good_delta (thus pointing to rows)
+            gap_indices = np.where(np.diff(data['timestamp']) > good_delta)[0]  # row (block) indices
+            for g in gap_indices:
+                assert (data['timestamp'][g+1] - data['timestamp'][g]) <= max_gap_duration_to_fill * 1e6, 'found a big gap, what should we do???'
+
+            # Get the intervals of continuous data (pointing to 'ts_norm' or the orignal (broken) data
+            data_starts = np.concatenate([np.zeros(1), gap_indices+1]).astype('int')
+            data_stops = np.concatenate([gap_indices, [data['timestamp'].size-1]]).astype('int')
+
+            # Construct bins for interpolated time axis
+            ts_bins = np.vstack([ts_interpolated, np.roll(ts_interpolated, -1)])
+            ts_bins[1, -1] = ts_bins[0, -1] + good_delta
+
+            for i0, i1 in zip(data_starts, data_stops):
+
+                # Find in which row of the interpolated time_axis the timepoint of the gap fits
+                i0_itp = np.where((ts_bins[0, :] <= ts_norm[i0]) & (ts_bins[1, :] > ts_norm[i0]))[0]
+                i1_itp = np.where((ts_bins[0, :] <= ts_norm[i1]) & (ts_bins[1, :] > ts_norm[i1]))[0]
+
+                # cant be sure enough
+                assert len(i0_itp) == 1 and len(i1_itp) == 1
+                i0_itp = i0_itp[0]
+                i1_itp = i1_itp[0]
+                assert ts_interpolated[i0_itp] <= ts_norm[i0] < ts_interpolated[i0_itp] + good_delta
+                assert ts_interpolated[i1_itp] <= ts_norm[i1] < ts_interpolated[i1_itp] + good_delta
+
+                # account for disalignment between interpolated time axis and OG time points
+                n_idx_in_block = int((ts_norm[i0] - ts_interpolated[i0_itp]) / sample_dt)
+                lin_i0_itp = sub2lin(i0_itp) + n_idx_in_block
+                lin_i1_itp = sub2lin(i1_itp) + n_idx_in_block
+
+                # write the data to the correct points!
+                records[lin_i0_itp:lin_i1_itp] = data['samples'].flatten()[sub2lin(i0): sub2lin(i1)]
+
+                # remark: I could mark the gap intervals, and interpolate the data in the gap
+
+            # Prepare the data for writing to ncs format
+            records = np.resize(records, (int(records.size/512), 512))
+            patched_data = np.core.records.fromarrays(
+                [ts_interpolated + t0_all, chan_id, sample_rate, nb_valid, records], dtype=np.dtype(ncs_dtype)
+            )
+
+            header = read_ncs_header(ncs_in_name)
+            with open(ncs_out_name, 'wb') as f:
+                f.write(header)
+                f.write(patched_data.tobytes())
+
+        print('Done!')
 
     def read_ncs_files(self, ncs_filenames):
         """
@@ -602,6 +740,7 @@ def read_ncs_files(self, ncs_filenames):
                     length = subdata.size * BLOCK_SIZE
                     self._sigs_length.append(length)
 
+
 # Keys funcitons
 def _to_bool(txt):
     if txt == 'True':
@@ -781,7 +920,8 @@ def read_txt_header(filename):
 
 ncs_dtype = [('timestamp', 'uint64'), ('channel_id', 'uint32'), ('sample_rate', 'uint32'),
              ('nb_valid', 'uint32'), ('samples', 'int16', (BLOCK_SIZE,))]
-
+ncs_dtype2 = [('timestamp', 'uint64', (1,1)), ('channel_id', 'uint32', (1,1)), ('sample_rate', 'uint32', (1,1)),
+             ('nb_valid', 'uint32', (1,1)), ('samples', 'int16', (BLOCK_SIZE, 512))]
 nev_dtype = [
     ('reserved', '<i2'),
     ('system_id', '<i2'),
@@ -823,5 +963,16 @@ def get_nse_or_ntt_dtype(info, ext):
     return dtype
 
 
-def get_ncs_memmap(filename):
-    return np.memmap(filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
\ No newline at end of file
+def read_ncs_memmap(filename):
+    return np.memmap(filename, dtype=ncs_dtype, mode='r', offset=HEADER_SIZE)
+
+
+def read_ncs_header(filename):
+    with open(filename, 'rb') as f:
+        header = f.read(HEADER_SIZE)
+    return header
+
+
+def sub2lin(idx):
+    return int(float(idx) * BLOCK_SIZE)
+

From 066c204c71ea74f0ebd0245f3b061c062e8f5e37 Mon Sep 17 00:00:00 2001
From: Thijs <t.r.ruikes@uva.nl>
Date: Sat, 30 May 2020 21:04:17 +0200
Subject: [PATCH 5/6] Restored baserawio.py

---
 neo/rawio/baserawio.py      |  2 +-
 neo/rawio/neuralynxrawio.py | 32 +++++++++++++++++++++++++++-----
 2 files changed, 28 insertions(+), 6 deletions(-)

diff --git a/neo/rawio/baserawio.py b/neo/rawio/baserawio.py
index 291256c09..fa0e98750 100644
--- a/neo/rawio/baserawio.py
+++ b/neo/rawio/baserawio.py
@@ -147,7 +147,7 @@ def parse_header(self):
 
         """
         self._parse_header()
-        # self._group_signal_channel_characteristics()
+        self._group_signal_channel_characteristics()
 
     def source_name(self):
         """Return fancy name of file source"""
diff --git a/neo/rawio/neuralynxrawio.py b/neo/rawio/neuralynxrawio.py
index df961ccd5..b7055f863 100644
--- a/neo/rawio/neuralynxrawio.py
+++ b/neo/rawio/neuralynxrawio.py
@@ -32,6 +32,7 @@
 BLOCK_SIZE = 512  # nb sample per signal block
 HEADER_SIZE = 2 ** 14  # file have a txt header of 16kB
 
+
 class NeuralynxRawIO(BaseRawIO):
     """"
     Class for reading dataset recorded by Neuralynx.
@@ -66,7 +67,6 @@ def _parse_header(self):
         event_channels = []
 
         self.ncs_filenames = OrderedDict()  # (chan_name, chan_id): filename
-        self._ncs_info = OrderedDict()  # (chan_name, chan_id): info
         self.nse_ntt_filenames = OrderedDict()  # (chan_name, chan_id): filename
         self.nev_filenames = OrderedDict()  # chan_id: filename
 
@@ -115,7 +115,30 @@ def _parse_header(self):
                     sig_channels.append((chan_name, chan_id, info['sampling_rate'],
                                          'int16', units, gain, offset, group_id))
                     self.ncs_filenames[chan_uid] = filename
-                    signal_annotations.append(info)  # is there a reason to not append all information
+                    keys = [
+                        'DspFilterDelay_µs',
+                        'recording_opened',
+                        'FileType',
+                        'DspDelayCompensation',
+                        'recording_closed',
+                        'DspLowCutFilterType',
+                        'HardwareSubSystemName',
+                        'DspLowCutNumTaps',
+                        'DSPLowCutFilterEnabled',
+                        'HardwareSubSystemType',
+                        'DspHighCutNumTaps',
+                        'ADMaxValue',
+                        'DspLowCutFrequency',
+                        'DSPHighCutFilterEnabled',
+                        'RecordSize',
+                        'InputRange',
+                        'DspHighCutFrequency',
+                        'input_inverted',
+                        'NumADChannels',
+                        'DspHighCutFilterType',
+                    ]
+                    d = {k: info[k] for k in keys if k in info}
+                    signal_annotations.append(d)
 
                 elif ext in ('nse', 'ntt'):
                     # nse and ntt are pretty similar except for the wavform shape
@@ -181,7 +204,7 @@ def _parse_header(self):
             self._sigs_sampling_rate = sampling_rate[0]
 
         # read ncs files for gaps detection and nb_segment computation
-        self.read_ncs_files(ncs_filenames=self.ncs_filenames)  # no need to pass class attributes to function
+        self.read_ncs_files(self.ncs_filenames)
 
         # timestamp limit in nev, nse
         # so need to scan all spike and event to
@@ -920,8 +943,7 @@ def read_txt_header(filename):
 
 ncs_dtype = [('timestamp', 'uint64'), ('channel_id', 'uint32'), ('sample_rate', 'uint32'),
              ('nb_valid', 'uint32'), ('samples', 'int16', (BLOCK_SIZE,))]
-ncs_dtype2 = [('timestamp', 'uint64', (1,1)), ('channel_id', 'uint32', (1,1)), ('sample_rate', 'uint32', (1,1)),
-             ('nb_valid', 'uint32', (1,1)), ('samples', 'int16', (BLOCK_SIZE, 512))]
+
 nev_dtype = [
     ('reserved', '<i2'),
     ('system_id', '<i2'),

From 9ce810d7a190e0e0f399a5a81e65e53f276a5c66 Mon Sep 17 00:00:00 2001
From: Thijs <t.r.ruikes@uva.nl>
Date: Fri, 5 Jun 2020 09:49:22 +0200
Subject: [PATCH 6/6] changed big gap assert to warning

---
 neo/rawio/neuralynxrawio.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/neo/rawio/neuralynxrawio.py b/neo/rawio/neuralynxrawio.py
index b7055f863..8eff4ace6 100644
--- a/neo/rawio/neuralynxrawio.py
+++ b/neo/rawio/neuralynxrawio.py
@@ -603,8 +603,8 @@ def patch_ncs_files(dirname, patchdir=None, interpolate=False, verbose=True, sca
             # gap indices are detected by dt > good_delta (thus pointing to rows)
             gap_indices = np.where(np.diff(data['timestamp']) > good_delta)[0]  # row (block) indices
             for g in gap_indices:
-                assert (data['timestamp'][g+1] - data['timestamp'][g]) <= max_gap_duration_to_fill * 1e6, 'found a big gap, what should we do???'
-
+                if data['timestamp'][g+1] - data['timestamp'][g] <= max_gap_duration_to_fill * 1e6:  #, 'found a big gap, what should we do???'
+                    print(f'WARNING - big gap filled {data["timestamp"][g+1] - data["timestamp"][g]} - {ncs_info["filenames"][i]}')
             # Get the intervals of continuous data (pointing to 'ts_norm' or the orignal (broken) data
             data_starts = np.concatenate([np.zeros(1), gap_indices+1]).astype('int')
             data_stops = np.concatenate([gap_indices, [data['timestamp'].size-1]]).astype('int')