diff --git a/neo/io/__init__.py b/neo/io/__init__.py index e3280da6f..7f5df3078 100644 --- a/neo/io/__init__.py +++ b/neo/io/__init__.py @@ -35,6 +35,7 @@ * :attr:`ElanIO` * :attr:`IgorIO` * :attr:`IntanIO` +* :attr:`IntanBinaryIO` * :attr:`MEArecIO` * :attr:`KlustaKwikIO` * :attr:`KwikIO` @@ -151,6 +152,10 @@ .. autoattribute:: extensions +.. autoclass:: neo.io.IntanBinaryIO + + .. autoattribute:: extensions + .. autoclass:: neo.io.KlustaKwikIO .. autoattribute:: extensions @@ -314,6 +319,7 @@ from neo.io.exampleio import ExampleIO from neo.io.igorproio import IgorIO from neo.io.intanio import IntanIO +from neo.io.intanbinaryio import IntanBinaryIO from neo.io.klustakwikio import KlustaKwikIO from neo.io.kwikio import KwikIO from neo.io.mearecio import MEArecIO @@ -368,6 +374,7 @@ ExampleIO, IgorIO, IntanIO, + IntanBinaryIO, KlustaKwikIO, KwikIO, MEArecIO, diff --git a/neo/io/intanbinaryio.py b/neo/io/intanbinaryio.py new file mode 100644 index 000000000..16a4073b3 --- /dev/null +++ b/neo/io/intanbinaryio.py @@ -0,0 +1,11 @@ +from neo.io.basefromrawio import BaseFromRaw +from neo.rawio.intanbinaryrawio import IntanBinaryRawIO + + +class IntanBinaryIO(IntanBinaryRawIO, BaseFromRaw): + __doc__ = IntanBinaryRawIO.__doc__ + _prefered_signal_group_mode = 'group-by-same-units' + + def __init__(self, dirname): + IntanBinaryRawIO.__init__(self, dirname=dirname) + BaseFromRaw.__init__(self, dirname) diff --git a/neo/rawio/__init__.py b/neo/rawio/__init__.py index 896a491fd..803e5940d 100644 --- a/neo/rawio/__init__.py +++ b/neo/rawio/__init__.py @@ -90,6 +90,9 @@ .. autoattribute:: extensions +.. autoclass:: neo.rawio.IntanBinaryRawIO + .. autoattribute:: extensions + .. autoclass:: neo.rawio.MaxwellRawIO .. autoattribute:: extensions @@ -189,6 +192,7 @@ from neo.rawio.elanrawio import ElanRawIO from neo.rawio.examplerawio import ExampleRawIO from neo.rawio.intanrawio import IntanRawIO +from neo.rawio.intanbinaryrawio import IntanBinaryRawIO from neo.rawio.maxwellrawio import MaxwellRawIO from neo.rawio.mearecrawio import MEArecRawIO from neo.rawio.medrawio import MedRawIO @@ -223,6 +227,7 @@ EDFRawIO, ElanRawIO, IntanRawIO, + IntanBinaryRawIO, MicromedRawIO, MaxwellRawIO, MEArecRawIO, diff --git a/neo/rawio/intanbinaryrawio.py b/neo/rawio/intanbinaryrawio.py new file mode 100644 index 000000000..ee8c7a6ed --- /dev/null +++ b/neo/rawio/intanbinaryrawio.py @@ -0,0 +1,442 @@ +""" +Support for Intan's binary file mode which is turned on by setting the save +option in Intan to 'One File Per Signal Type' or 'One File Per Channel'. +The reader for this format is for rhd only since it uses the rhd header. +Need to confirm if rhs can also use this save format. + +RHD supported 1.x, 2.x, and 3.0, 3.1, 3.2 +See: + * http://intantech.com/files/Intan_RHD2000_data_file_formats.pdf + * http://intantech.com/files/Intan_RHS2000_data_file_formats.pdf + +Author: Zach McKenzie, based on Samuel Garcia's IntanRawIO +""" +from pathlib import Path + +import numpy as np +from collections import OrderedDict +from packaging.version import Version as V + +from .baserawio import ( + BaseRawIO, + _signal_channel_dtype, + _signal_stream_dtype, + _spike_channel_dtype, + _event_channel_dtype, + _common_sig_characteristics, +) + +from .intanrawio import ( + read_variable_header, + rhd_global_header_base, + rhd_global_header_part1, + rhd_global_header_v11, + rhd_global_header_v13, + rhd_global_header_v20, + rhd_global_header_final, + rhd_signal_group_header, + rhd_signal_channel_header, + stream_type_to_name, +) + + +class IntanBinaryRawIO(BaseRawIO): + """ + Class for processing Intan Data when saved in binary format. Requires an + `info.rhd` as well one file per signal stream or one file per channel + + Parameters + ---------- + dirname: str, Path + The root directory containing the info.rhd file + + """ + + extensions = ["rhd", "dat"] + rawmode = "one-dir" + + def __init__(self, dirname=""): + BaseRawIO.__init__(self) + self.dirname = dirname + + def _source_name(self): + return self.dirname + + def _parse_header(self): + + dir_path = Path(self.dirname) + assert (dir_path / "info.rhd").exists(), ( + "IntanBinaryRawIO requires the root directory containing the `info.rhd`" + ) + + header_file = dir_path / "info.rhd" + + one_file_per_signal = any((dir_path / file).exists() for file in one_file_per_signal_filenames) + + self.one_file_per_signal = one_file_per_signal + + if one_file_per_signal: + raw_file_paths_dict = create_one_file_per_signal_dict(dir_path) + else: + raw_file_paths_dict = create_one_file_per_channel_dict(dir_path) + + ( + self._global_info, + self._ordered_channels, + data_dtype, + self._block_size, + ) = read_rhd(header_file) + + self._raw_data = {} + for stream_index, sub_datatype in data_dtype.items(): + if one_file_per_signal: + self._raw_data[stream_index] = np.memmap( + raw_file_paths_dict[stream_index], dtype=sub_datatype, mode="r" + ) + else: + self._raw_data[stream_index] = [] + for channel_index, datatype in enumerate(sub_datatype): + self._raw_data[stream_index].append( + np.memmap( + raw_file_paths_dict[stream_index][channel_index], + dtype=[datatype], + mode="r", + ) + ) + + # check timestamp continuity + if one_file_per_signal: + timestamp = self._raw_data[6]["timestamp"].flatten() + else: + timestamp = self._raw_data[6][0]["timestamp"].flatten() + assert np.all(np.diff(timestamp) == 1), "timestamp have gaps" + + # signals + signal_channels = [] + for c, chan_info in enumerate(self._ordered_channels): + name = chan_info["native_channel_name"] + chan_id = str(c) # the chan_id have no meaning in intan + if chan_info["signal_type"] == 20: + # exception for temperature + sig_dtype = "int16" + else: + sig_dtype = "uint16" + stream_id = str(chan_info["signal_type"]) + signal_channels.append( + ( + name, + chan_id, + chan_info["sampling_rate"], + sig_dtype, + chan_info["units"], + chan_info["gain"], + chan_info["offset"], + stream_id, + ) + ) + signal_channels = np.array(signal_channels, dtype=_signal_channel_dtype) + + stream_ids = np.unique(signal_channels["stream_id"]) + signal_streams = np.zeros(stream_ids.size, dtype=_signal_stream_dtype) + signal_streams["id"] = stream_ids + for stream_index, stream_id in enumerate(stream_ids): + signal_streams["name"][stream_index] = stream_type_to_name.get( + int(stream_id), "" + ) + + self._max_sampling_rate = np.max(signal_channels["sampling_rate"]) + + if one_file_per_signal: + self._max_sigs_length = max( + [ + raw_data.size * self._block_size + for raw_data in self._raw_data.values() + ] + ) + else: + self._max_sigs_length = max( + [ + len(raw_data) * raw_data[0].size * self._block_size + for raw_data in self._raw_data.values() + ] + ) + + # No events + event_channels = [] + event_channels = np.array(event_channels, dtype=_event_channel_dtype) + + # No spikes + spike_channels = [] + spike_channels = np.array(spike_channels, dtype=_spike_channel_dtype) + + # fill into header dict + self.header = {} + self.header["nb_block"] = 1 + self.header["nb_segment"] = [1] + self.header["signal_streams"] = signal_streams + self.header["signal_channels"] = signal_channels + self.header["spike_channels"] = spike_channels + self.header["event_channels"] = event_channels + + self._generate_minimal_annotations() + + def _segment_t_start(self, block_index, seg_index): + return 0.0 + + def _segment_t_stop(self, block_index, seg_index): + t_stop = self._max_sigs_length / self._max_sampling_rate + return t_stop + + def _get_signal_size(self, block_index, seg_index, stream_index): + stream_id = self.header["signal_streams"][stream_index]["id"] + mask = self.header["signal_channels"]["stream_id"] == stream_id + signal_channels = self.header["signal_channels"][mask] + channel_names = signal_channels["name"] + chan_name0 = channel_names[0] + if self.one_file_per_signal: + size = self._raw_data[stream_index][chan_name0].size + else: + size = self._raw_data[stream_index][0][chan_name0].size + return size + + def _get_signal_t_start(self, block_index, seg_index, stream_index): + return 0.0 + + def _get_analogsignal_chunk( + self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes + ): + + if i_start is None: + i_start = 0 + if i_stop is None: + i_stop = self._get_signal_size(block_index, seg_index, stream_index) + + stream_id = self.header["signal_streams"][stream_index]["id"] + mask = self.header["signal_channels"]["stream_id"] == stream_id + signal_channels = self.header["signal_channels"][mask] + if channel_indexes is None: + channel_indexes = slice(None) + channel_names = signal_channels["name"][channel_indexes] + + if self.one_file_per_signal: + shape = self._raw_data[stream_index][channel_names[0]].shape + else: + shape = self._raw_data[stream_index][0][channel_names[0]].shape + + # some channel (temperature) have 1D field so shape 1D + # because 1 sample per block + if len(shape) == 2: + # this is the general case with 2D + block_size = shape[1] + block_start = i_start // block_size + block_stop = i_stop // block_size + 1 + + sl0 = i_start % block_size + sl1 = sl0 + (i_stop - i_start) + + sigs_chunk = np.zeros((i_stop - i_start, len(channel_names)), dtype="uint16") + for i, chan_name in enumerate(channel_names): + if self.one_file_per_signal: + data_chan = self._raw_data[stream_index][chan_name] + else: + data_chan = self._raw_data[stream_index][i][chan_name] + if len(shape) == 1: + sigs_chunk[:, i] = data_chan[i_start:i_stop] + else: + sigs_chunk[:, i] = data_chan[block_start:block_stop].flatten()[sl0:sl1] + + return sigs_chunk + + +############ +# RHD Zone for Binary Files + +# For One File Per Signal +one_file_per_signal_filenames = [ + "amplifier.dat", + "auxiliary.dat", + "supply.dat", + "analogin.dat", + "digitalin.dat", + "digitalout.dat", +] + +# For One File Per Channel +possible_raw_file_prefixes = [ + "amp", + "aux", + "vdd", + "board-ANALOG", + "board-DIGITAL-IN", + "board-DIGITAL-OUT", +] + + +def create_one_file_per_signal_dict(dirname): + """Function for One File Per Signal Type""" + raw_file_dict = {} + for raw_index, raw_file in enumerate(one_file_per_signal_filenames): + if Path(dirname / raw_file).is_file(): + raw_file_dict[raw_index] = Path(dirname / raw_file) + raw_file_dict[6] = Path(dirname / "time.dat") + + return raw_file_dict + + +def create_one_file_per_channel_dict(dirname): + """Utility function for One File Per Channel""" + file_names = dirname.glob("**/*.dat") + files = [file for file in file_names if file.is_file()] + raw_file_dict = {} + for raw_index, prefix in enumerate(possible_raw_file_prefixes): + raw_file_dict[raw_index] = [file for file in files if prefix in file.name] + raw_file_dict[6] = [Path(dirname / "time.dat")] + + return raw_file_dict + + +def read_rhd(filename): + with open(filename, mode="rb") as f: + + global_info = read_variable_header(f, rhd_global_header_base) + version = V("{major_version}.{minor_version}".format(**global_info)) + + # the header size depends on the version :-( + header = list(rhd_global_header_part1) # make a copy + + if version >= V("1.1"): + header = header + rhd_global_header_v11 + else: + global_info["num_temp_sensor_channels"] = 0 + + if version >= V("1.3"): + header = header + rhd_global_header_v13 + else: + global_info["eval_board_mode"] = 0 + + if version >= V("2.0"): + header = header + rhd_global_header_v20 + else: + global_info["reference_channel"] = "" + + header = header + rhd_global_header_final + + global_info.update(read_variable_header(f, header)) + + # read channel group and channel header + channels_by_type = {k: [] for k in [0, 1, 2, 3, 4, 5]} + data_dtype = {k: [] for k in range(7)} # 5 channels + 6 is for time stamps + for g in range(global_info["nb_signal_group"]): + group_info = read_variable_header(f, rhd_signal_group_header) + + if bool(group_info["signal_group_enabled"]): + for c in range(group_info["channel_num"]): + chan_info = read_variable_header(f, rhd_signal_channel_header) + if bool(chan_info["channel_enabled"]): + channels_by_type[chan_info["signal_type"]].append(chan_info) + + sr = global_info["sampling_rate"] + + # construct the data block dtype and reorder channels + if version >= V("2.0"): + BLOCK_SIZE = 128 + else: + BLOCK_SIZE = 60 # 256 channels + + ordered_channels = [] + + # 6: Timestamp stored in time.dat + if version >= V("1.2"): + data_dtype[6] = [("timestamp", "int32", BLOCK_SIZE)] + else: + data_dtype[6] = [("timestamp", "uint32", BLOCK_SIZE)] + + # 0: RHD2000 amplifier channel stored in amplifier.dat/amp-* + for chan_info in channels_by_type[0]: + name = chan_info["native_channel_name"] + chan_info["sampling_rate"] = sr + chan_info["units"] = "uV" + chan_info["gain"] = 0.195 + chan_info["offset"] = -32768 * 0.195 + ordered_channels.append(chan_info) + data_dtype[0] += [(name, "uint16", BLOCK_SIZE)] + + # 1: RHD2000 auxiliary input channel stored in auxiliary.dat/aux-* + for chan_info in channels_by_type[1]: + name = chan_info["native_channel_name"] + chan_info["sampling_rate"] = sr / 4.0 + chan_info["units"] = "V" + chan_info["gain"] = 0.0000374 + chan_info["offset"] = 0.0 + ordered_channels.append(chan_info) + data_dtype[1] += [(name, "uint16", BLOCK_SIZE // 4)] + + # 2: RHD2000 supply voltage channel stored in supply.dat/vdd-* + for chan_info in channels_by_type[2]: + name = chan_info["native_channel_name"] + chan_info["sampling_rate"] = sr / BLOCK_SIZE + chan_info["units"] = "V" + chan_info["gain"] = 0.0000748 + chan_info["offset"] = 0.0 + ordered_channels.append(chan_info) + data_dtype[2] += [(name, "uint16",)] + + # temperature is not an official channel in the header + for i in range(global_info["num_temp_sensor_channels"]): + name = "temperature_{}".format(i) + chan_info = {"native_channel_name": name, "signal_type": 20} + chan_info["sampling_rate"] = sr / BLOCK_SIZE + chan_info["units"] = "Celsius" + chan_info["gain"] = 0.001 + chan_info["offset"] = 0.0 + ordered_channels.append(chan_info) + data_dtype += [(name, "int16",)] + + # 3: USB board ADC input channel stored in analogin.dat/board-ANALOG-* + for chan_info in channels_by_type[3]: + name = chan_info["native_channel_name"] + chan_info["sampling_rate"] = sr + chan_info["units"] = "V" + if global_info["eval_board_mode"] == 0: + chan_info["gain"] = 0.000050354 + chan_info["offset"] = 0.0 + elif global_info["eval_board_mode"] == 1: + chan_info["gain"] = 0.00015259 + chan_info["offset"] = -32768 * 0.00015259 + elif global_info["eval_board_mode"] == 13: + chan_info["gain"] = 0.0003125 + chan_info["offset"] = -32768 * 0.0003125 + ordered_channels.append(chan_info) + data_dtype[3] += [(name, "uint16", BLOCK_SIZE)] + + # 4: USB board digital input channel stored in digitalin.dat/board-DIGITAL-IN-* + # 5: USB board digital output channel stored in digitalout.dat/board-DIGITAL-OUT-* + for sig_type in [4, 5]: + # User can obtain digital channels from analog_signal_chunk and then process + # them themself + if len(channels_by_type[sig_type]) > 0: + name = {4: "DIGITAL-IN", 5: "DIGITAL-OUT"}[sig_type] + chan_info = channels_by_type[sig_type][0] + chan_info["native_channel_name"] = name # overwite to allow memmap to work + chan_info["sampling_rate"] = sr + chan_info["units"] = "TTL" # arbitrary units so I did TTL for the logic + chan_info["gain"] = 1.0 + chan_info["offset"] = 0.0 + ordered_channels.append(chan_info) + data_dtype[sig_type] += [(name, "uint16", BLOCK_SIZE)] + + if bool(global_info["notch_filter_mode"]) and version >= V("3.0"): + global_info["notch_filter_applied"] = True + else: + global_info["notch_filter_applied"] = False + + # because data_dtype has all possible dtypes we need to delete the dtypes which aren't contained in this dataset + # so if the data_dtype value is still an empty list we delete to ignore. + dtype_cleanup = [] + for key, value in data_dtype.items(): + if len(value)==0: + dtype_cleanup.append(key) + for key in dtype_cleanup: + _ = data_dtype.pop(key) + + return global_info, ordered_channels, data_dtype, BLOCK_SIZE diff --git a/neo/test/iotest/test_intanbinaryio.py b/neo/test/iotest/test_intanbinaryio.py new file mode 100644 index 000000000..fd2d13023 --- /dev/null +++ b/neo/test/iotest/test_intanbinaryio.py @@ -0,0 +1,19 @@ +import unittest + +from neo.io import IntanBinaryIO +from neo.test.iotest.common_io_test import BaseTestIO + + +class TestIntanBinaryIO(BaseTestIO, unittest.TestCase, ): + ioclass = IntanBinaryIO + entities_to_download = [ + 'intan' + ] + entities_to_test = [ + 'intan/intan_fpc_test_231117_052630/', + 'intan/intan_fps_test_231117_052500/', + ] + + +if __name__ == "__main__": + unittest.main() diff --git a/neo/test/rawiotest/test_intanbinaryrawio.py b/neo/test/rawiotest/test_intanbinaryrawio.py new file mode 100644 index 000000000..4bf4682e4 --- /dev/null +++ b/neo/test/rawiotest/test_intanbinaryrawio.py @@ -0,0 +1,20 @@ +import unittest + +from neo.rawio.intanbinaryrawio import IntanBinaryRawIO + +from neo.test.rawiotest.common_rawio_test import BaseTestRawIO + + +class TestIntanBinaryRawIO(BaseTestRawIO, unittest.TestCase, ): + rawioclass = IntanBinaryRawIO + entities_to_download = [ + 'intan' + ] + entities_to_test = [ + 'intan/intan_fpc_test_231117_052630/', + 'intan/intan_fps_test_231117_052500/', + ] + + +if __name__ == "__main__": + unittest.main()