diff --git a/bark/__init__.py b/bark/__init__.py index 093cb7b..1e55310 100644 --- a/bark/__init__.py +++ b/bark/__init__.py @@ -1,3 +1,4 @@ from __future__ import absolute_import from bark.bark import * from bark.bark import __version__ +from bark import stream diff --git a/bark/bark.py b/bark/bark.py index b135869..4462085 100644 --- a/bark/bark.py +++ b/bark/bark.py @@ -13,9 +13,11 @@ import arrow import yaml import numpy as np +import pandas as pd import functools as ft BUFFER_SIZE = 10000 +DEFAULT_META = '.meta.yaml' spec_version = "0.2" __version__ = "0.2" @@ -28,9 +30,6 @@ bark: %s """ % (__version__) -_Units = namedtuple('_Units', ['TIME_UNITS']) -UNITS = _Units(TIME_UNITS=('s', 'samples')) - _pairs = ((None, None), ('UNDEFINED', 0), ('ACOUSTIC', 1), ('EXTRAC_HP', 2), ('EXTRAC_LF', 3), ('EXTRAC_EEG', 4), ('INTRAC_CC', 5), ('INTRAC_VC', 6), ('EVENT', 1000), ('SPIKET', 1001), @@ -44,7 +43,7 @@ -class LazyDict(dict): +class _LazyDict(dict): """Allows lazy loading of data and memoizing the result. If value is a function, evaluates and replaces with the return @@ -68,7 +67,7 @@ def __init__(self, path): # entries are lazily loaded by creating a dictionary # with the entry name and a function, that when called # loads the data. See the custom LazyDict data structure - self.entries = LazyDict({os.path.split(x)[-1]: ft.partial(read_entry, name=x) for x in subdirs}) + self.entries = _LazyDict({os.path.split(x)[-1]: ft.partial(read_entry, name=x) for x in subdirs}) def __getitem__(self, item): return self.entries[item] @@ -78,6 +77,15 @@ def __len__(self): def __contains__(self, item): return self.entries.__contains__(item) + + def close(self): + for e_name in self.entries: + entry = self.entries.get(e_name) + if not callable(entry): + entry.close() + p = entry.path + self.entries[e_name] = ft.partial(read_entry, name=p) + del entry class Entry(): @@ -100,6 +108,14 @@ def __contains__(self, item): def __lt__(self, other): return self.timestamp < other.timestamp + + def close(self): + for ds_name in self.datasets: + dataset = self.datasets.get(ds_name) + if not callable(dataset): + p = dataset.path + self.datasets[ds_name] = ft.partial(read_dataset, fname=p) + del dataset class Data(): @@ -168,18 +184,20 @@ def event_columns(dataframe, columns=None): register with `dataframe`'s columns Returns: - dict, or None: Columns dictionary for `dataframe` + dict: Columns dictionary for `dataframe` """ if columns is None: - return template_columns(dataframe.columns) - for fieldkey in columns: - if fieldkey not in dataframe.columns: - del columns[fieldkey] - for col in dataframe.columns: - if col not in columns: - columns[col] = {'units': None} - if 'units' not in columns[col]: - columns[col]['units'] = None + columns = template_columns(dataframe.columns) + else: + for fieldkey in columns: + if fieldkey not in dataframe.columns: + del columns[fieldkey] + for col in dataframe.columns: + if col not in columns: + columns[col] = {'units': None} + if 'units' not in columns[col]: + columns[col]['units'] = None + return columns def sampled_columns(data, columns=None): @@ -193,7 +211,7 @@ def sampled_columns(data, columns=None): register with shape of `data` Returns: - dict, or None: Columns dictionary for `data` + dict: Columns dictionary for `data` Raises: ValueError: if the keys in `columns` don't match up with `data` @@ -203,16 +221,18 @@ def sampled_columns(data, columns=None): else: n_channels = data.shape[1] if columns is None: - return template_columns(range(n_channels)) - if len(columns) != n_channels: - raise ValueError( - 'the columns attribute does not match the number of columns') - for i in range(n_channels): - if i not in columns: + columns = template_columns(range(n_channels)) + else: + if len(columns) != n_channels: raise ValueError( - 'the columns attribute is missing column {}'.format(i)) - if 'units' not in columns[i]: - columns[i]['units'] = None + 'the columns attribute does not match the number of columns') + for i in range(n_channels): + if i not in columns: + raise ValueError( + 'the columns attribute is missing column {}'.format(i)) + if 'units' not in columns[i]: + columns[i]['units'] = None + return columns def write_sampled(datfile, data, sampling_rate, **params): @@ -229,6 +249,7 @@ def write_sampled(datfile, data, sampling_rate, **params): Returns: SampledData: sampled dataset containing `data` """ + path = os.path.abspath(datfile) if 'columns' not in params: params['columns'] = sampled_columns(data) params["dtype"] = data.dtype.str @@ -237,7 +258,7 @@ def write_sampled(datfile, data, sampling_rate, **params): mdata[:] = data[:] write_metadata(datfile, sampling_rate=sampling_rate, **params) params['sampling_rate'] = sampling_rate - return SampledData(mdata, datfile, params) + return SampledData(mdata, path, params) def read_sampled(datfile, mode="r"): @@ -274,7 +295,6 @@ def write_events(eventsfile, data, **params): Returns: EventData: event dataset containing `data` """ - import pandas as pd if 'columns' not in params: params['columns'] = event_columns(data) if data.empty and not list(data.columns): @@ -293,10 +313,10 @@ def read_events(eventsfile): Returns: EventData: event dataset containing `eventsfile`'s data """ - import pandas as pd - data = pd.read_csv(eventsfile).fillna('') + path = os.path.abspath(eventsfile) params = read_metadata(eventsfile) - return EventData(data, eventsfile, params) + data = pd.read_csv(eventsfile).fillna('') + return EventData(data, path, params) def read_dataset(fname): @@ -316,7 +336,7 @@ def read_dataset(fname): return dset -def read_metadata(path, meta='.meta.yaml'): +def read_metadata(path, meta=DEFAULT_META): """Loads metadata for a dataset. Args: @@ -348,7 +368,7 @@ def read_metadata(path, meta='.meta.yaml'): raise FileNotFoundError(m.format(path, metafile)) -def write_metadata(path, meta='.meta.yaml', **params): +def write_metadata(path, meta=DEFAULT_META, **params): """Writes metadata for a dataset. Args: @@ -359,6 +379,7 @@ def write_metadata(path, meta='.meta.yaml', **params): **params: all other keyword arguments are treated as dataset attributes, and added to the meta file """ + # the two following checks are for backwards-compatibility if 'n_channels' in params: del params['n_channels'] if 'n_samples' in params: @@ -422,7 +443,7 @@ def create_entry(name, timestamp, parents=False, **attributes): return read_entry(name) -def read_entry(name, meta=".meta.yaml"): +def read_entry(name, meta=DEFAULT_META): """Reads a Bark Entry from a directory. Args: @@ -441,7 +462,7 @@ def read_entry(name, meta=".meta.yaml"): # datasets are lazily loaded by creating a dictionary # with the dataset name and a function, that when called # loads the data. See the custom LazyDict data structure - datasets = LazyDict({name: ft.partial(read_dataset, fname=full_name) + datasets = _LazyDict({name: ft.partial(read_dataset, fname=full_name) for name, full_name in zip(dset_names, dset_full_names)}) return Entry(datasets, path, attrs) diff --git a/bark/io/arf2bark.py b/bark/io/arf2bark.py index c7203c1..e5e0de6 100644 --- a/bark/io/arf2bark.py +++ b/bark/io/arf2bark.py @@ -8,9 +8,12 @@ import numpy import collections as coll +ENTRY_PREFIX = 'entry' + def _parse_args(raw_args): desc = 'Unspool an HDF5 ARF file into a Bark tree.' - parser = argparse.ArgumentParser(description=desc) + epi = 'Fails if bark_root already exists.' + parser = argparse.ArgumentParser(description=desc, epilog=epi) parser.add_argument('-v', '--verbose', help='increase output verbosity', @@ -19,8 +22,12 @@ def _parse_args(raw_args): '--timezone', help='timezone for data, tz database format (default is "America/Chicago")', default=None) + parser.add_argument('-m', + '--mangle-prefix', + help='prefix for names which collide with bark function arguments (default: {})'.format(ENTRY_PREFIX), + default=ENTRY_PREFIX) parser.add_argument('arf_file', help='ARF file to convert') - parser.add_argument('root_parent', help='directory in which to place the Bark Root') + parser.add_argument('bark_root', help='location of new bark root') return parser.parse_args(raw_args) def copy_attrs(attrs): @@ -30,11 +37,8 @@ def copy_attrs(attrs): na.update({k: v.decode() for k,v in na.items() if isinstance(v, bytes)}) return na -def arf2bark(arf_file, root_parent, timezone, verbose): +def arf2bark(arf_file, root_path, timezone, verbose, mangle_prefix=ENTRY_PREFIX): with arf.open_file(arf_file, 'r') as af: - # root - root_dirname = os.path.splitext(arf_file)[0] - root_path = os.path.join(os.path.abspath(root_parent), root_dirname) os.mkdir(root_path) root = bark.Root(root_path) if verbose: @@ -45,6 +49,21 @@ def arf2bark(arf_file, root_parent, timezone, verbose): if isinstance(entry, h5py.Group): # entries entry_path = os.path.join(root_path, ename) entry_attrs = copy_attrs(entry.attrs) + for pos_arg in ('name', 'parents'): + # along with 'timestamp' below, these are positional arguments to create_entry + # for now, I prefer hard-coding them over messing with runtime introspection + new_name = pos_arg + while new_name in entry_attrs: + new_name = '{}_{}'.format(mangle_prefix, new_name) + try: + entry_attrs[new_name] = entry_attrs.pop(pos_arg) + except KeyError: + pass + else: + if verbose: + print('Renamed attribute {} of entry {} to {}'.format(pos_arg, + ename, + new_name)) timestamp = entry_attrs.pop('timestamp') if timezone: timestamp = bark.convert_timestamp(timestamp, timezone) @@ -62,10 +81,13 @@ def arf2bark(arf_file, root_parent, timezone, verbose): else: transfer_dset(ds_name, dataset, entry_path, verbose) elif isinstance(entry, h5py.Dataset): # top-level datasets - if tle is None: - path = os.path.join(root_path, 'top_level') - tle = bark.create_entry(path, 0, parents=False).path - transfer_dset(ename, entry, tle, verbose) + if arf.is_time_series(entry) or arf.is_marked_pointproc(entry): + if tle is None: + path = os.path.join(root_path, 'top_level') + tle = bark.create_entry(path, 0, parents=False).path + transfer_dset(ename, entry, tle, verbose) + else: + unknown_ds_warning(ename) # and skip, w/o creating TLE if found_trigin: print('Warning: found datasets named "trig_in". Jill-created ' + '"trig_in" datasets segfault when read, so these datasets' + @@ -84,6 +106,10 @@ def build_columns(units, column_names=None): d.update({k: None for k,v in d.items() if (k == 'units' and v == '')}) return cols +def unknown_ds_warning(ds_name): + print('Warning: unknown dataset type - neither time series nor point' + + ' process. Skipping dataset ' + ds_name) + def transfer_dset(ds_name, ds, e_path, verbose=False): ds_attrs = copy_attrs(ds.attrs) units = ds_attrs.pop('units', None) @@ -107,12 +133,11 @@ def transfer_dset(ds_name, ds, e_path, verbose=False): if verbose: print('Created event dataset: ' + ds_path) else: - print('Warning: unknown dataset type - neither time series nor point' + - ' process. Skipping dataset ' + ds_name) + unknown_ds_warning(ds_name) def _main(): args = _parse_args(sys.argv[1:]) - arf2bark(args.arf_file, args.root_parent, args.timezone, args.verbose) + arf2bark(args.arf_file, args.bark_root, args.timezone, args.verbose, args.mangle_prefix) if __name__ == '__main__': _main() diff --git a/bark/io/datfromwav.py b/bark/io/datfromwav.py new file mode 100644 index 0000000..14dabfb --- /dev/null +++ b/bark/io/datfromwav.py @@ -0,0 +1,40 @@ +from scipy.io import wavfile +import bark, os + + +def dat_from_wav(wav, barkname, **attrs): + rate, data = wavfile.read(wav) + return bark.write_sampled(barkname, data, rate,**attrs) + + +def _main(): + ''' Function for getting commandline args.''' + + import argparse + + p = argparse.ArgumentParser(description=''' + converts wav file to bark format + ''') + p.add_argument('wav', help='path to wav file') + p.add_argument('out', help="path to bark file") + p.add_argument("-a", + "--attributes", + action='append', + type=lambda kv: kv.split("="), + dest='keyvalues', + help="extra metadata in the form of KEY=VALUE") + + args = p.parse_args() + + if args.keyvalues: + dat_from_wav(args.wav, + args.out, + **dict(args.keyvalues)) + else: + dat_from_wav(args.wav, args.out) + + + + +if __name__ == '__main__': + _main() diff --git a/bark/io/mda.py b/bark/io/mda.py new file mode 100644 index 0000000..aa41e65 --- /dev/null +++ b/bark/io/mda.py @@ -0,0 +1,262 @@ +import argparse +import bark +import numpy as np +import pandas as pd +import struct +import sys + +# the primary I/O functions below are adapted from those in mountainlab_pytools +# (https://github.com/magland/mountainlab_pytools) + +_DATATYPE_FROM_DT_CODE = {-2: 'uint8', + -3: 'float32', + -4: 'int16', + -5: 'int32', + -6: 'uint16', + -7: 'float64', + -8: 'uint32'} +_DT_CODE_FROM_DATATYPE = {v: k for k, v in _DATATYPE_FROM_DT_CODE.items()} + +class MdaHeader: + def __init__(self, dtype, dimensions, uses_64bit_dims=False): + self.uses64bitdims = uses_64bit_dims + self.dt_code = _DT_CODE_FROM_DATATYPE[dtype] + self.dt = dtype + self.num_bytes_per_entry = np.dtype(dtype).itemsize + self.num_dims = len(dimensions) + self.dimprod = np.prod(dimensions) + self.dims = dimensions + self.header_size = 3 * 4 + self.num_dims * (8 if uses_64bit_dims else 4) + +def read_mda_header(filename): + """Read an .mda file's header information. + + Args: + filename (str): path to the .mda file to read + + Returns: + MdaHeader: object containing header information + + Raises: + ValueError: if the header indicates dimension count < 1, or if the + datatype code is not recognized + """ + with open(filename, 'rb') as mda_file: + # read first 3 header items to determine header size + dt_code, _, num_dims = struct.unpack(' 0: - result['amplifier_channels'] = header['amplifier_channels'] - if data_present: - result['amplifier_data'] = data['amplifier_data'] - result['t_amplifier'] = data['t_amplifier'] - result['spike_triggers'] = header['spike_triggers'] - - if header['num_aux_input_channels'] > 0: - result['aux_input_channels'] = header['aux_input_channels'] - if data_present: - result['aux_input_data'] = data['aux_input_data'] - result['t_aux_input'] = data['t_aux_input'] - - if header['num_supply_voltage_channels'] > 0: - result['supply_voltage_channels'] = header['supply_voltage_channels'] - if data_present: - result['supply_voltage_data'] = data['supply_voltage_data'] - result['t_supply_voltage'] = data['t_supply_voltage'] - - if header['num_board_adc_channels'] > 0: - result['board_adc_channels'] = header['board_adc_channels'] - if data_present: - result['board_adc_data'] = data['board_adc_data'] - result['t_board_adc'] = data['t_board_adc'] - - if header['num_board_dig_in_channels'] > 0: - result['board_dig_in_channels'] = header['board_dig_in_channels'] - if data_present: - result['board_dig_in_data'] = data['board_dig_in_data'] - result['t_dig'] = data['t_dig'] - - if header['num_board_dig_out_channels'] > 0: - result['board_dig_out_channels'] = header['board_dig_out_channels'] - if data_present: - result['board_dig_out_data'] = data['board_dig_out_data'] - result['t_dig'] = data['t_dig'] - - if header['num_temp_sensor_channels'] > 0: - if data_present: - result['temp_sensor_data'] = data['temp_sensor_data'] - result['t_temp_sensor'] = data['t_temp_sensor'] - - return result +#! /bin/env python +# +# Michael Gibson 27 April 2015 + +def data_to_result(header, data, data_present): + """Moves the header and data (if present) into a common object.""" + + result = {} + result['version'] = header['version'] + result['notes'] = header['notes'] + result['frequency_parameters'] = header['frequency_parameters'] + if 'reference_channel' in header: + result['digital_reference_channel'] = header['reference_channel'] + + result['t_amplifier'] = data['t_amplifier'] + if header['num_amplifier_channels'] > 0: + result['amplifier_channels'] = header['amplifier_channels'] + if data_present: + result['amplifier_data'] = data['amplifier_data'] + result['spike_triggers'] = header['spike_triggers'] + + if header['num_aux_input_channels'] > 0: + result['aux_input_channels'] = header['aux_input_channels'] + if data_present: + result['aux_input_data'] = data['aux_input_data'] + result['t_aux_input'] = data['t_aux_input'] + + if header['num_supply_voltage_channels'] > 0: + result['supply_voltage_channels'] = header['supply_voltage_channels'] + if data_present: + result['supply_voltage_data'] = data['supply_voltage_data'] + result['t_supply_voltage'] = data['t_supply_voltage'] + + if header['num_board_adc_channels'] > 0: + result['board_adc_channels'] = header['board_adc_channels'] + if data_present: + result['board_adc_data'] = data['board_adc_data'] + result['t_board_adc'] = data['t_board_adc'] + + if header['num_board_dig_in_channels'] > 0: + result['board_dig_in_channels'] = header['board_dig_in_channels'] + if data_present: + result['board_dig_in_data'] = data['board_dig_in_data'] + result['t_dig'] = data['t_dig'] + + if header['num_board_dig_out_channels'] > 0: + result['board_dig_out_channels'] = header['board_dig_out_channels'] + if data_present: + result['board_dig_out_data'] = data['board_dig_out_data'] + result['t_dig'] = data['t_dig'] + + if header['num_temp_sensor_channels'] > 0: + if data_present: + result['temp_sensor_data'] = data['temp_sensor_data'] + result['t_temp_sensor'] = data['t_temp_sensor'] + + return result diff --git a/bark/io/rhd/get_bytes_per_data_block.py b/bark/io/rhd/get_bytes_per_data_block.py index 332e5b4..655065e 100644 --- a/bark/io/rhd/get_bytes_per_data_block.py +++ b/bark/io/rhd/get_bytes_per_data_block.py @@ -1,35 +1,59 @@ -#! /bin/env python -# -# Michael Gibson 23 April 2015 - - -def get_bytes_per_data_block(header): - """Calculates the number of bytes in each 60-sample datablock.""" - - # Each data block contains 60 amplifier samples. - bytes_per_block = 60 * 4 # timestamp data - bytes_per_block = bytes_per_block + 60 * 2 * header['num_amplifier_channels'] - - # Auxiliary inputs are sampled 4x slower than amplifiers - bytes_per_block = bytes_per_block + 15 * 2 * header['num_aux_input_channels'] - - # Supply voltage is sampled 60x slower than amplifiers - bytes_per_block = bytes_per_block + 1 * 2 * header['num_supply_voltage_channels'] - - # Board analog inputs are sampled at same rate as amplifiers - bytes_per_block = bytes_per_block + 60 * 2 * header['num_board_adc_channels'] - - # Board digital inputs are sampled at same rate as amplifiers - if header['num_board_dig_in_channels'] > 0: - bytes_per_block = bytes_per_block + 60 * 2 - - # Board digital outputs are sampled at same rate as amplifiers - if header['num_board_dig_out_channels'] > 0: - bytes_per_block = bytes_per_block + 60 * 2 - - # Temp sensor is sampled 60x slower than amplifiers - if header['num_temp_sensor_channels'] > 0: - bytes_per_block = bytes_per_block + 1 * 2 * header['num_temp_sensor_channels'] - - return bytes_per_block - \ No newline at end of file +#! /bin/env python +# +# Michael Gibson 23 April 2015 +# 2018 changes by Adrian Foy merged in 2021 by Graham Fetterman + +from . import constants as const + +def get_bytes_per_data_block(header): + """Calculates the number of bytes in each datablock.""" + + num_samples = header['num_samples_per_data_block'] + bytes_per_block = 0 + + # Timebase (with version-specific data type) + if (header['version']['major'], header['version']['minor']) >= (1, 2): + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 + else: + time_dtype = const.TIMESTAMP_DTYPE_LE_V1_1 + bytes_per_block += num_samples * time_dtype.itemsize + + # Each data block contains a version-specific number of amplifier samples. + bytes_per_block += (num_samples * + const.AMPLIFIER_DTYPE.itemsize * + header['num_amplifier_channels']) + + # Auxiliary inputs are sampled 4x slower than amplifiers. + bytes_per_block += ((num_samples / 4) * + const.AUXILIARY_DTYPE.itemsize * + header['num_aux_input_channels']) + + # Supply voltage is sampled only once per data block (i.e., 60x or 128x + # slower than amplifiers). + bytes_per_block += (1 * + const.SUPPLY_DTYPE.itemsize + * header['num_supply_voltage_channels']) + + # Board analog inputs are sampled at same rate as amplifiers. + bytes_per_block += (num_samples * + const.ADC_DTYPE.itemsize * + header['num_board_adc_channels']) + + # Board digital inputs are sampled at same rate as amplifiers, and packed + # together. + if header['num_board_dig_in_channels'] > 0: + bytes_per_block += num_samples * const.DIG_IN_DTYPE.itemsize + + # Board digital outputs are sampled at same rate as amplifiers, and packed + # together. + if header['num_board_dig_out_channels'] > 0: + bytes_per_block += num_samples * const.DIG_OUT_DTYPE.itemsize + + # Temperature is sampled only once per data block (i.e., 60x or 128x slower + # than amplifiers). + if header['num_temp_sensor_channels'] > 0: + bytes_per_block += (1 * + const.TEMP_DTYPE.itemsize * + header['num_temp_sensor_channels']) + + return bytes_per_block diff --git a/bark/io/rhd/legacy_load_intan_rhd_format.py b/bark/io/rhd/legacy_load_intan_rhd_format.py new file mode 100644 index 0000000..c640b49 --- /dev/null +++ b/bark/io/rhd/legacy_load_intan_rhd_format.py @@ -0,0 +1,243 @@ +#! /bin/env python +# +# Michael Gibson 17 July 2015 +# Kyler Brown December 2016 +# 2021 changes by Adrian Foy merged in 2021 by Graham Fetterman + +import os +import numpy as np + +from bark.io.rhd.read_header import read_header +from bark.io.rhd.get_bytes_per_data_block import get_bytes_per_data_block +from bark.io.rhd.legacy_read_one_data_block import read_one_data_block +from bark.io.rhd.notch_filter import notch_filter +from bark.io.rhd.data_to_result import data_to_result + +from . import constants as const + +UINT16_BIT_OFFSET = int(2**15) + +def read_data(filename, + no_floats=False, + digital_io_data_dtype=np.uint): + """Reads Intan RHD2000 data file generated by evaluation board GUI. + + Data are returned in a dictionary, for future extensibility. + + Args: + filename (str): file to read data from + no_floats (bool): whether to expand sampled data outputs (including + amplifiers and ADCs) from int16 to float64 (defaults + to False, i.e., do the expansion). + digital_io_data_dtype (numpy dtype): what dtype the digital I/O data is + stored as (default: np.uint) + Returns: + dict: all data and metadata in the file + """ + + fid = open(filename, 'rb') + filesize = os.path.getsize(filename) + + header = read_header(fid) + + if header['notch_filter_frequency'] > 0: + msg = ('Warning: a notch filter ({}Hz) was applied in the GUI, ' + + 'but has not been applied here.') + print(msg.format(header['notch_filter_frequency'])) + + channel_reporting = [('amplifier', 'num_amplifier_channels'), + ('auxiliary input', 'num_aux_input_channels'), + ('supply voltage', 'num_supply_voltage_channels'), + ('board ADC', 'num_board_adc_channels'), + ('board digital input', 'num_board_dig_in_channels'), + ('board digital output', 'num_board_dig_out_channels'), + ('temperature sensor', 'num_temp_sensor_channels')] + for channel_name, channel_count_id in channel_reporting: + print('Found {} {} channel{}.'.format(header[channel_count_id], + channel_name, + plural(header[channel_count_id]))) + print('') + + # Determine how many samples the data file contains. + bytes_per_block = get_bytes_per_data_block(header) + + # How many data blocks remain in this file? + data_present = False + bytes_remaining = filesize - fid.tell() + if bytes_remaining > 0: + data_present = True + + if bytes_remaining % bytes_per_block != 0: + msg = ('Something is wrong with file size: ' + + 'should have a whole number of data blocks') + raise ValueError(msg) + + num_data_blocks = int(bytes_remaining / bytes_per_block) + + num_samples = header['num_samples_per_data_block'] + num_amplifier_samples = num_samples * num_data_blocks + num_aux_input_samples = (num_samples // 4) * num_data_blocks + num_supply_voltage_samples = 1 * num_data_blocks + num_board_adc_samples = num_samples * num_data_blocks + num_board_dig_in_samples = num_samples * num_data_blocks + num_board_dig_out_samples = num_samples * num_data_blocks + + record_time = num_amplifier_samples / header['sample_rate'] + + if data_present: + print('File contains {:0.3f} seconds of data.'.format(record_time), + 'Amplifiers were sampled at', + '{:0.2f} kHz.'.format(header['sample_rate'] / 1000)) + else: + print('File contains no data. Amplifiers were sampled at', + '{:0.2f} kHz.'.format(header['sample_rate'] / 1000)) + + if data_present: + # Pre-allocate memory for data. + data = {} + if (header['version']['major'], header['version']['minor']) >= (1, 2): + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 + else: + time_dtype = const.TIMESTAMP_DTYPE_LE_V1_1 + data['t_amplifier'] = np.zeros(num_amplifier_samples, dtype=time_dtype) + + data['amplifier_data'] = np.zeros([header['num_amplifier_channels'], + num_amplifier_samples], + dtype=const.AMPLIFIER_DTYPE) + data['aux_input_data'] = np.zeros([header['num_aux_input_channels'], + num_aux_input_samples], + dtype=const.AUXILIARY_DTYPE) + data['supply_voltage_data'] = np.zeros([header['num_supply_voltage_channels'], + num_supply_voltage_samples], + dtype=const.SUPPLY_DTYPE) + data['temp_sensor_data'] = np.zeros([header['num_temp_sensor_channels'], + num_supply_voltage_samples], + dtype=const.TEMP_DTYPE) + data['board_adc_data'] = np.zeros([header['num_board_adc_channels'], + num_board_adc_samples], + dtype=const.ADC_DTYPE) + data['board_dig_in_data'] = np.zeros([header['num_board_dig_in_channels'], + num_board_dig_in_samples], + dtype=digital_io_data_dtype) + data['board_dig_in_raw'] = np.zeros(num_board_dig_in_samples, + dtype=const.DIG_IN_DTYPE) + data['board_dig_out_data'] = np.zeros([header['num_board_dig_out_channels'], + num_board_dig_out_samples], + dtype=digital_io_data_dtype) + data['board_dig_out_raw'] = np.zeros(num_board_dig_out_samples, + dtype=const.DIG_OUT_DTYPE) + + # Initialize indices used in looping + indices = {} + indices['amplifier'] = 0 + indices['aux_input'] = 0 + indices['supply_voltage'] = 0 + indices['board_adc'] = 0 + indices['board_dig_in'] = 0 + indices['board_dig_out'] = 0 + + for i in range(num_data_blocks): + read_one_data_block(data, header, indices, fid) + + # Increment indices + indices['amplifier'] += num_samples + indices['aux_input'] += num_samples // 4 + indices['supply_voltage'] += 1 + indices['board_adc'] += num_samples + indices['board_dig_in'] += num_samples + indices['board_dig_out'] += num_samples + + # Make sure we have read exactly the right amount of data. + bytes_remaining = filesize - fid.tell() + if bytes_remaining != 0: + raise Exception('Error: End of file not reached.') + + # Close data file. + fid.close() + + extras = {} # extra parameters + if data_present: + + # Extract digital input channels to separate variables. + for i in range(header['num_board_dig_in_channels']): + mask = 1 << header['board_dig_in_channels'][i]['native_order'] + masked_bits = np.bitwise_and(data['board_dig_in_raw'], mask) + data['board_dig_in_data'][i,:] = masked_bits.astype(np.bool) + + # Extract digital output channels to separate variables. + for i in range(header['num_board_dig_out_channels']): + mask = 1 << header['board_dig_out_channels'][i]['native_order'] + masked_bits = np.bitwise_and(data['board_dig_out_raw'], mask) + data['board_dig_out_data'][i,:] = masked_bits.astype(np.bool) + + if no_floats: + # record the bit voltage scaling level but do not apply to the data + # converting to floats increases size 4x, which makes a big + # difference at the terabyte+ level. + extras['amplifier_bit_microvolts'] = const.AMPLIFIER_BIT_MICROVOLTS + data['amplifier_data'] = (data['amplifier_data'].astype(np.int32) - + UINT16_BIT_OFFSET).astype(np.int16) + extras['aux_bit_volts'] = const.AUX_BIT_VOLTS + extras['supply_bit_volts'] = const.SUPPLY_BIT_VOLTS + extras['temp_bit_celcius'] = const.TEMP_BIT_CELCIUS + + if header['eval_board_mode'] == 1: + extras['ADC_input_bit_volts'] = const.ADC_BIT_VOLTS_MODE_1 + elif header['eval_board_mode'] == 13: + extra['ADC_input_bit_volts'] = const.ADC_BIT_VOLTS_MODE_13 + else: + extras['ADC_input_bit_volts'] = const.ADC_BIT_VOLTS_MODE_0 + data['board_adc_data'] = (data['board_adc_data'].astype(np.int32) - + UINT16_BIT_OFFSET).astype(np.int16) + else: + # Scale voltage levels appropriately. + offset_amp_data = (data['amplifier_data'].astype(np.int32) - + UINT16_BIT_OFFSET) + data['amplifier_data'] = np.multiply(const.AMPLIFIER_BIT_MICROVOLTS, + offset_amp_data) # units of microvolts + data['aux_input_data'] = np.multiply(const.AUX_BIT_VOLTS, + data['aux_input_data']) # units of volts + data['supply_voltage_data'] = np.multiply(const.SUPPLY_BIT_VOLTS, + data['supply_voltage_data']) # units of volts + if header['eval_board_mode'] == 1: + offset_adc_data = (data['board_adc_data'].astype(np.int32) - + UINT16_BIT_OFFSET) + data['board_adc_data'] = np.multiply(const.ADC_BIT_VOLTS_MODE_1, + offset_adc_data) # units of volts + elif header['eval_board_mode'] == 13: + offset_adc_data = (data['board_adc_data'].astype(np.int32) - + UINT16_BIT_OFFSET) + data['board_adc_data'] = np.multiply(const.ADC_BIT_VOLTS_MODE_13, + offset_adc_data) # units of volts + else: + data['board_adc_data'] = np.multiply(const.ADC_BIT_VOLTS_MODE_0, + data['board_adc_data']) # units of volts + data['temp_sensor_data'] = np.multiply(const.TEMP_BIT_CELCIUS, + data['temp_sensor_data']) # units of degrees C + + # Check for gaps in timestamps. + num_gaps = np.sum(np.diff(data['t_amplifier']) != 1) + if num_gaps != 0: + print('Warning: {} gaps in timestamp data found.'.format(num_gaps), + 'Time scale will not be uniform!') + + # Scale time steps (units of seconds). + data['t_amplifier'] = data['t_amplifier'] / header['sample_rate'] + data['t_aux_input'] = data['t_amplifier'][0:len(data['t_amplifier']):4] + per_block = data['t_amplifier'][0:len(data['t_amplifier']):num_samples] + data['t_supply_voltage'] = per_block + data['t_temp_sensor'] = per_block + data['t_board_adc'] = data['t_amplifier'] + data['t_dig'] = data['t_amplifier'] + + else: + data = [] + + # Move variables to result struct. + result = data_to_result(header, data, data_present) + result.update(extras) + return result + + +def plural(n): + return '' if n == 1 else 's' diff --git a/bark/io/rhd/legacy_read_one_data_block.py b/bark/io/rhd/legacy_read_one_data_block.py new file mode 100644 index 0000000..cf09a84 --- /dev/null +++ b/bark/io/rhd/legacy_read_one_data_block.py @@ -0,0 +1,105 @@ +#! /bin/env python +# +# Michael Gibson 23 April 2015 +# 2018 changes by Adrian Foy merged in 2021 by Graham Fetterman + +import struct +import numpy as np + +from . import constants as const + +def read_one_data_block(data, header, indices, fid): + """Reads one data block from fid into data, at location per indices.""" + + num_samples = header['num_samples_per_data_block'] + + # Timebase + + # Prior to version 1.2, timestamps are unsigned integers. + # From version 1.2 onwards, timestamps are signed integers, to accommodate + # negative (adjusted) timestamps for pretrigger data. + if (header['version']['major'], header['version']['minor']) >= (1, 2): + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 + else: + time_dtype = const.TIMESTAMP_DTYPE_LE_V1_1 + start = indices['amplifier'] + stop = indices['amplifier'] + num_samples + byte_layout = '<' + time_dtype.char * num_samples + byte_count = time_dtype.itemsize * num_samples + values = np.array(struct.unpack(byte_layout, fid.read(byte_count))) + data['t_amplifier'][start:stop] = values + + # Amplifier channels + + if header['num_amplifier_channels'] > 0: + num_channels = header['num_amplifier_channels'] + start = indices['amplifier'] + stop = indices['amplifier'] + num_samples + num_values = num_samples * num_channels + values = np.fromfile(fid, dtype=const.AMPLIFIER_DTYPE, count=num_values) + values = values.reshape(num_channels, num_samples) + data['amplifier_data'][:,start:stop] = values + + # Auxiliary input channels + + if header['num_aux_input_channels'] > 0: + num_channels = header['num_aux_input_channels'] + start = indices['aux_input'] + stop = indices['aux_input'] + num_samples // 4 + num_values = (num_samples // 4) * num_channels + values = np.fromfile(fid, dtype=const.AUXILIARY_DTYPE, count=num_values) + values = values.reshape(num_channels, num_samples // 4) + data['aux_input_data'][:,start:stop] = values + + # Supply voltage channels + + if header['num_supply_voltage_channels'] > 0: + num_channels = header['num_supply_voltage_channels'] + start = indices['supply_voltage'] + stop = indices['supply_voltage'] + 1 + num_values = 1 * num_channels + values = np.fromfile(fid, dtype=const.SUPPLY_DTYPE, count=num_values) + values = values.reshape(num_channels, 1) + data['supply_voltage_data'][:,start:stop] = values + + # Temperature sensor channels + + if header['num_temp_sensor_channels'] > 0: + num_channels = header['num_temp_sensor_channels'] + start = indices['supply_voltage'] + stop = indices['supply_voltage'] + 1 + num_values = 1 * num_channels + values = np.fromfile(fid, dtype=const.TEMP_DTYPE, count=num_values) + values = values.reshape(num_channels, 1) + data['temp_sensor_data'][:,start:stop] = values + + # Board ADC channels + + if header['num_board_adc_channels'] > 0: + num_channels = header['num_board_adc_channels'] + start = indices['board_adc'] + stop = indices['board_adc'] + num_samples + num_values = num_samples * num_channels + values = np.fromfile(fid, dtype=const.ADC_DTYPE, count=num_values) + values = values.reshape(num_channels, num_samples) + data['board_adc_data'][:,start:stop] = values + + # Board digital input channels (packed together) + + if header['num_board_dig_in_channels'] > 0: + start = indices['board_dig_in'] + stop = indices['board_dig_in'] + num_samples + byte_layout = '<' + const.DIG_IN_DTYPE.char * num_samples + byte_count = const.DIG_IN_DTYPE.itemsize * num_samples + values = np.array(struct.unpack(byte_layout, fid.read(byte_count))) + data['board_dig_in_raw'][start:stop] = values + + # Board digital output channels (packed together) + + if header['num_board_dig_out_channels'] > 0: + start = indices['board_dig_out'] + stop = indices['board_dig_out'] + num_samples + byte_layout = '<' + const.DIG_OUT_DTYPE.char * num_samples + byte_count = const.DIG_OUT_DTYPE.itemsize * num_samples + values = np.array(struct.unpack(byte_layout, fid.read(byte_count))) + data['board_dig_out_raw'][start:stop] = values diff --git a/bark/io/rhd/load_intan_rhd_format.py b/bark/io/rhd/load_intan_rhd_format.py index 039bacf..16a85f7 100644 --- a/bark/io/rhd/load_intan_rhd_format.py +++ b/bark/io/rhd/load_intan_rhd_format.py @@ -1,262 +1,221 @@ -#! /bin/env python -# -# Michael Gibson 17 July 2015 -# Kyler Brown December 2016 - -from __future__ import absolute_import, division, unicode_literals, print_function - -import sys, struct, math, os, time -import numpy as np - -from bark.io.rhd.read_header import read_header -from bark.io.rhd.get_bytes_per_data_block import get_bytes_per_data_block -from bark.io.rhd.read_one_data_block import read_one_data_block -from bark.io.rhd.notch_filter import notch_filter -from bark.io.rhd.data_to_result import data_to_result - -# constants -AMPLIFIER_BIT_MICROVOLTS = 0.195 -UINT16_BIT_OFFSET = int(2**15) -AUX_BIT_VOLTS = 37.4e-6 -SUPPLY_BIT_VOLTS = 74.8e-6 -ADC_BIT_VOLTS_1 = 152.59e-6 -ADC_BIT_VOLTS_0 = 50.353e-6 -TEMP_BIT_CELCIUS = 0.01 - - -def read_data(filename, no_floats=False): - """Reads Intan Technologies RHD2000 data file generated by evaluation board GUI. - Data are returned in a dictionary, for future extensibility. - """ - - tic = time.time() - fid = open(filename, 'rb') - filesize = os.path.getsize(filename) - - header = read_header(fid) - - print('Found {} amplifier channel{}.'.format(header[ - 'num_amplifier_channels'], plural(header['num_amplifier_channels']))) - print('Found {} auxiliary input channel{}.'.format(header[ - 'num_aux_input_channels'], plural(header['num_aux_input_channels']))) - print('Found {} supply voltage channel{}.'.format(header[ - 'num_supply_voltage_channels'], plural(header[ - 'num_supply_voltage_channels']))) - print('Found {} board ADC channel{}.'.format(header[ - 'num_board_adc_channels'], plural(header['num_board_adc_channels']))) - print('Found {} board digital input channel{}.'.format(header[ - 'num_board_dig_in_channels'], plural(header[ - 'num_board_dig_in_channels']))) - print('Found {} board digital output channel{}.'.format(header[ - 'num_board_dig_out_channels'], plural(header[ - 'num_board_dig_out_channels']))) - print('Found {} temperature sensors channel{}.'.format(header[ - 'num_temp_sensor_channels'], plural(header[ - 'num_temp_sensor_channels']))) - print('') - - # Determine how many samples the data file contains. - bytes_per_block = get_bytes_per_data_block(header) - - # How many data blocks remain in this file? - data_present = False - bytes_remaining = filesize - fid.tell() - if bytes_remaining > 0: - data_present = True - - if bytes_remaining % bytes_per_block != 0: - raise Exception( - 'Something is wrong with file size : should have a whole number of data blocks') - - num_data_blocks = int(bytes_remaining / bytes_per_block) - - num_amplifier_samples = 60 * num_data_blocks - num_aux_input_samples = 15 * num_data_blocks - num_supply_voltage_samples = 1 * num_data_blocks - num_board_adc_samples = 60 * num_data_blocks - num_board_dig_in_samples = 60 * num_data_blocks - num_board_dig_out_samples = 60 * num_data_blocks - - record_time = num_amplifier_samples / header['sample_rate'] - - if data_present: - print( - 'File contains {:0.3f} seconds of data. Amplifiers were sampled at {:0.2f} kS/s.'.format( - record_time, header['sample_rate'] / 1000)) - else: - print( - 'Header file contains no data. Amplifiers were sampled at {:0.2f} kS/s.'.format( - header['sample_rate'] / 1000)) - - if data_present: - # Pre-allocate memory for data. - data = {} - if (header['version']['major'] == 1 and - header['version']['minor'] >= 2) or ( - header['version']['major'] > 1): - data['t_amplifier'] = np.zeros(num_amplifier_samples, dtype=np.int) - else: - data['t_amplifier'] = np.zeros(num_amplifier_samples, - dtype=np.uint) - - data['amplifier_data'] = np.zeros( - [header['num_amplifier_channels'], num_amplifier_samples], - dtype=np.uint16) - data['aux_input_data'] = np.zeros( - [header['num_aux_input_channels'], num_aux_input_samples], - dtype=np.uint16) - data['supply_voltage_data'] = np.zeros( - [header['num_supply_voltage_channels'], - num_supply_voltage_samples], - dtype=np.uint16) - data['temp_sensor_data'] = np.zeros( - [header['num_temp_sensor_channels'], num_supply_voltage_samples], - dtype=np.uint16) - data['board_adc_data'] = np.zeros( - [header['num_board_adc_channels'], num_board_adc_samples], - dtype=np.uint16) - data['board_dig_in_data'] = np.zeros( - [header['num_board_dig_in_channels'], num_board_dig_in_samples], - dtype=np.uint) - data['board_dig_in_raw'] = np.zeros(num_board_dig_in_samples, - dtype=np.uint) - data['board_dig_out_data'] = np.zeros( - [header['num_board_dig_out_channels'], num_board_dig_out_samples], - dtype=np.uint) - data['board_dig_out_raw'] = np.zeros(num_board_dig_out_samples, - dtype=np.uint) - - # Initialize indices used in looping - indices = {} - indices['amplifier'] = 0 - indices['aux_input'] = 0 - indices['supply_voltage'] = 0 - indices['board_adc'] = 0 - indices['board_dig_in'] = 0 - indices['board_dig_out'] = 0 - - print_increment = 10 - percent_done = print_increment - for i in range(num_data_blocks): - read_one_data_block(data, header, indices, fid) - - # Increment indices - indices['amplifier'] += 60 - indices['aux_input'] += 15 - indices['supply_voltage'] += 1 - indices['board_adc'] += 60 - indices['board_dig_in'] += 60 - indices['board_dig_out'] += 60 - - fraction_done = 100 * (1.0 * i / num_data_blocks) - if fraction_done >= percent_done: - percent_done = percent_done + print_increment - # Make sure we have read exactly the right amount of data. - bytes_remaining = filesize - fid.tell() - if bytes_remaining != 0: - raise Exception('Error: End of file not reached.') - -# Close data file. - fid.close() - - extras = {} # dictionary for extra parameters - if (data_present): - - # Extract digital input channels to separate variables. - for i in range(header['num_board_dig_in_channels']): - data['board_dig_in_data'][i, :] = np.not_equal( - np.bitwise_and(data['board_dig_in_raw'], ( - 1 << header['board_dig_in_channels'][i]['native_order'])), - 0) - -# Extract digital output channels to separate variables. - for i in range(header['num_board_dig_out_channels']): - data['board_dig_out_data'][i, :] = np.not_equal( - np.bitwise_and(data['board_dig_out_raw'], ( - 1 << header['board_dig_out_channels'][i]['native_order'])), - 0) - if no_floats: - # record the bit voltage scaling level but do not apply to the data - # converting to floats increases size 4x, which makes a big difference at the terrabyte+ level. - extras['amplifier_bit_microvolts'] = AMPLIFIER_BIT_MICROVOLTS - data['amplifier_data'] = (data['amplifier_data'].astype(np.int32) - - UINT16_BIT_OFFSET).astype(np.int16) - extras['aux_bit_volts'] = AUX_BIT_VOLTS - extras['supply_bit_volts'] = SUPPLY_BIT_VOLTS - extras['temp_bit_celcius'] = TEMP_BIT_CELCIUS - - if header['eval_board_mode'] == 1: - extras['ADC_input_bit_volts'] = ADC_BIT_VOLTS_1 - - else: - extras['ADC_input_bit_volts'] = ADC_BIT_VOLTS_0 - data['board_adc_data'] = (data['board_adc_data'].astype(np.int32) - - UINT16_BIT_OFFSET).astype(np.int16) - else: - # Scale voltage levels appropriately. - data['amplifier_data'] = np.multiply(AMPLIFIER_BIT_MICROVOLTS, ( - data['amplifier_data'].astype(np.int32) - UINT16_BIT_OFFSET) - ) # units = microvolts - data['aux_input_data'] = np.multiply( - AUX_BIT_VOLTS, data['aux_input_data']) # units = volts - data['supply_voltage_data'] = np.multiply( - SUPPLY_BIT_VOLTS, data['supply_voltage_data']) # units = volts - if header['eval_board_mode'] == 1: - data['board_adc_data'] = np.multiply( - ADC_BIT_VOLTS_1, (data['board_adc_data'].astype(np.int32) - - UINT16_BIT_OFFSET)) # units = volts - else: - data['board_adc_data'] = np.multiply( - ADC_BIT_VOLTS_0, data['board_adc_data']) # units = volts - data['temp_sensor_data'] = np.multiply( - TEMP_BIT_CELCIUS, data['temp_sensor_data']) # units = deg C - -# Check for gaps in timestamps. - num_gaps = np.sum(np.not_equal(data['t_amplifier'][1:] - data[ - 't_amplifier'][:-1], 1)) - if num_gaps != 0: - print( - 'Warning: {0} gaps in timestamp data found. Time scale will not be uniform!'.format( - num_gaps)) - -# Scale time steps (units = seconds). - data['t_amplifier'] = data['t_amplifier'] / header['sample_rate'] - data['t_aux_input'] = data['t_amplifier'][range(0, len(data[ - 't_amplifier']), 4)] - data['t_supply_voltage'] = data['t_amplifier'][range(0, len(data[ - 't_amplifier']), 60)] - data['t_board_adc'] = data['t_amplifier'] - data['t_dig'] = data['t_amplifier'] - data['t_temp_sensor'] = data['t_supply_voltage'] - - # If the software notch filter was selected during the recording, apply the - # same notch filter to amplifier data here. - if header['notch_filter_frequency'] > 0: - print('Applying notch filter...') - - print_increment = 10 - percent_done = print_increment - for i in range(header['num_amplifier_channels']): - data['amplifier_data'][i, :] = notch_filter( - data['amplifier_data'][i, :], header['sample_rate'], - header['notch_filter_frequency'], 10) - - fraction_done = 100 * (i / header['num_amplifier_channels']) - if fraction_done >= percent_done: - percent_done += print_increment - else: - data = [] - -# Move variables to result struct. - result = data_to_result(header, data, data_present) - result.update(extras) - return result - - -def plural(n): - return '' if n == 1 else 's' - - -if __name__ == '__main__': - a = read_data(sys.argv[1]) - #print a +#! /bin/env python +# +# Michael Gibson 17 July 2015 +# Kyler Brown December 2016 +# Graham Fetterman July 2018, 2021 + +import os +import numpy as np + +from bark.io.rhd.read_header import read_header +from bark.io.rhd.get_bytes_per_data_block import get_bytes_per_data_block +from bark.io.rhd.read_data_blocks import read_data_blocks, preallocate_memory +from bark.io.rhd.data_to_result import data_to_result + +from . import constants as const + +UINT16_BIT_OFFSET = int(2**15) + +def read_data(filename, no_floats=False, max_memory=0): + """Reads Intan Technologies RHD2000 data file generated by evaluation board GUI. + + Data are yielded in a dictionary, for future extensibility. + + A file's contents are split into chunks governed by max_memory. + + Args: + filename (str): .rhd file name + no_floats (bool): whether to avoid converting 16-bit integer values + to 32-bit floats (not converting saves disk space) + max_memory (int): size of chunks to split file's data into (in bytes) + + Yields: + dict: containing data fields and some metadata + + """ + + fid = open(filename, 'rb') + filesize = os.path.getsize(filename) + + header = read_header(fid) + + if header['notch_filter_frequency'] > 0: + msg = ('Warning: a notch filter ({}Hz) was applied in the GUI, ' + + 'but has not been applied here.') + print(msg.format(header['notch_filter_frequency'])) + + channel_reporting = [('amplifier', 'num_amplifier_channels'), + ('auxiliary input', 'num_aux_input_channels'), + ('supply voltage', 'num_supply_voltage_channels'), + ('board ADC', 'num_board_adc_channels'), + ('board digital input', 'num_board_dig_in_channels'), + ('board digital output', 'num_board_dig_out_channels'), + ('temperature sensor', 'num_temp_sensor_channels')] + for channel_name, channel_count_id in channel_reporting: + print('Found {} {} channel{}.'.format(header[channel_count_id], + channel_name, + plural(header[channel_count_id]))) + print('') + + # Determine how many samples the data file contains. + bytes_per_block = get_bytes_per_data_block(header) + + # How many data blocks remain in this file? + data_present = False + bytes_remaining = filesize - fid.tell() + if bytes_remaining > 0: + data_present = True + + if bytes_remaining % bytes_per_block != 0: + msg = ('Something is wrong with file size: ' + + 'should have a whole number of data blocks') + raise ValueError(msg) + + num_data_blocks = int(bytes_remaining / bytes_per_block) + + num_samples = header['num_samples_per_data_block'] + num_amplifier_samples = num_samples * num_data_blocks + num_aux_input_samples = (num_samples // 4) * num_data_blocks + num_supply_voltage_samples = 1 * num_data_blocks + num_board_adc_samples = num_samples * num_data_blocks + num_board_dig_in_samples = num_samples * num_data_blocks + num_board_dig_out_samples = num_samples * num_data_blocks + + record_time = num_amplifier_samples / header['sample_rate'] + + if data_present: + print('File contains {:0.3f} seconds of data.'.format(record_time), + 'Amplifiers were sampled at', + '{:0.2f} kHz.'.format(header['sample_rate'] / 1000)) + else: + print('File contains no data. Amplifiers were sampled at', + '{:0.2f} kHz.'.format(header['sample_rate'] / 1000)) + + if data_present: + # chunk_size governs how many datablocks are read in and then written + # to file at once + # minimum is 1 datablock, maximum is every datablock in the file + # two copies of the chunk are in memory at once for some operations + # (reformatting, changing dtypes), so max_memory is divided by 2 + chunk_size = min(max(int(0.5 * max_memory / bytes_per_block), 1), + num_data_blocks) + chunks, remainder = divmod(num_data_blocks, chunk_size) + data = preallocate_memory(header, chunk_size) + for _ in range(chunks): + read_data_blocks(data, header, fid, datablocks_per_chunk=chunk_size) + yield check_data_and_reformat(header, data, no_floats) + if remainder: + data = preallocate_memory(header, remainder) + read_data_blocks(data, header, fid, datablocks_per_chunk=remainder) + yield check_data_and_reformat(header, data, no_floats) + # Make sure we have read exactly the right amount of data. + bytes_remaining = filesize - fid.tell() + if bytes_remaining != 0: + raise Exception('Error: End of file not reached.') + else: + yield data_to_result(header, {}, data_present) + # Close data file. + fid.close() + +def check_data_and_reformat(header, data, no_floats): + """Performs some cleanup on data and builds a dictionary to return. + + Args: + header (dict): metadata for the data + data (dict): different data fields are contained in numpy arrays + no_floats (bool): whether to expand 16-bit ints into 32-bit floats + + Returns: + dict: combining data and some metadata + """ + extras = {} # dictionary for extra parameters + + # Extract digital input channels to separate variables. + for i in range(header['num_board_dig_in_channels']): + mask = 1 << header['board_dig_in_channels'][i]['native_order'] + masked_bits = np.bitwise_and(data['board_dig_in_raw'], mask) + data['board_dig_in_data'][i,:] = masked_bits.astype(np.bool) + + # Extract digital output channels to separate variables. + for i in range(header['num_board_dig_out_channels']): + mask = 1 << header['board_dig_out_channels'][i]['native_order'] + masked_bits = np.bitwise_and(data['board_dig_out_raw'], mask) + data['board_dig_out_data'][i,:] = masked_bits.astype(np.bool) + if no_floats: + # record the bit voltage scaling level but do not apply to the data + # converting to floats increases size 4x, which makes a big + # difference at the terabyte+ level. + extras['amplifier_bit_microvolts'] = const.AMPLIFIER_BIT_MICROVOLTS + # numpy doesn't do over/underflow checks, so this actually works as intended + np.subtract(data['amplifier_data'], + UINT16_BIT_OFFSET, + data['amplifier_data'], + casting='unsafe') + data['amplifier_data'] = data['amplifier_data'].astype(np.int16, + copy=False) + extras['aux_bit_volts'] = const.AUX_BIT_VOLTS + extras['supply_bit_volts'] = const.SUPPLY_BIT_VOLTS + extras['temp_bit_celcius'] = const.TEMP_BIT_CELCIUS + + if header['eval_board_mode'] == 1: + extras['ADC_input_bit_volts'] = const.ADC_BIT_VOLTS_MODE_1 + elif header['eval_board_mode'] == 13: + extra['ADC_input_bit_volts'] = const.ADC_BIT_VOLTS_MODE_13 + else: + extras['ADC_input_bit_volts'] = const.ADC_BIT_VOLTS_MODE_0 + # underflow is intentional here as well + np.subtract(data['board_adc_data'], + UINT16_BIT_OFFSET, + data['board_adc_data'], + casting='unsafe') + data['board_adc_data'] = data['board_adc_data'].astype(np.int16, + copy=False) + else: + # Scale voltage levels appropriately. + offset_amp_data = (data['amplifier_data'].astype(np.int32) - + UINT16_BIT_OFFSET) + data['amplifier_data'] = np.multiply(const.AMPLIFIER_BIT_MICROVOLTS, + offset_amp_data) # units of microvolts + data['aux_input_data'] = np.multiply(const.AUX_BIT_VOLTS, + data['aux_input_data']) # units of volts + data['supply_voltage_data'] = np.multiply(const.SUPPLY_BIT_VOLTS, + data['supply_voltage_data']) # units of volts + if header['eval_board_mode'] == 1: + offset_adc_data = (data['board_adc_data'].astype(np.int32) - + UINT16_BIT_OFFSET) + data['board_adc_data'] = np.multiply(const.ADC_BIT_VOLTS_MODE_1, + offset_adc_data) # units of volts + elif header['eval_board_mode'] == 13: + offset_adc_data = (data['board_adc_data'].astype(np.int32) - + UINT16_BIT_OFFSET) + data['board_adc_data'] = np.multiply(const.ADC_BIT_VOLTS_MODE_13, + offset_adc_data) # units of volts + else: + data['board_adc_data'] = np.multiply(const.ADC_BIT_VOLTS_MODE_0, + data['board_adc_data']) # units of volts + data['temp_sensor_data'] = np.multiply(const.TEMP_BIT_CELCIUS, + data['temp_sensor_data']) # units of degrees C + + # Check for gaps in timestamps. + num_gaps = np.sum(np.diff(data['t_amplifier']) != 1) + if num_gaps != 0: + print('Warning: {} gaps in timestamp data found.'.format(num_gaps), + 'Time scale will not be uniform!') + + # Scale time steps (units = seconds). + data['t_amplifier'] = data['t_amplifier'] / header['sample_rate'] + data['t_aux_input'] = data['t_amplifier'][0:len(data['t_amplifier']):4] + per_block = data['t_amplifier'][0:len(data['t_amplifier']):header['num_samples_per_data_block']] + data['t_supply_voltage'] = per_block + data['t_temp_sensor'] = per_block + data['t_board_adc'] = data['t_amplifier'] + data['t_dig'] = data['t_amplifier'] + + # Move variables to result struct. + result = data_to_result(header, data, data_present=True) + result.update(extras) + return result + + +def plural(n): + return '' if n == 1 else 's' diff --git a/bark/io/rhd/qstring.py b/bark/io/rhd/qstring.py index e7e21b2..3a546c1 100644 --- a/bark/io/rhd/qstring.py +++ b/bark/io/rhd/qstring.py @@ -1,41 +1,33 @@ -#! /bin/env python -# -# Michael Gibson 23 April 2015 - - -import sys, struct, os - -def read_qstring(fid): - """Read Qt style QString. - - The first 32-bit unsigned number indicates the length of the string (in bytes). - If this number equals 0xFFFFFFFF, the string is null. - - Strings are stored as unicode. - """ - - length, = struct.unpack(' (os.fstat(fid.fileno()).st_size - fid.tell() + 1) : - print(length) - raise Exception('Length too long.') - - # convert length from bytes to 16-bit Unicode words - length = int(length / 2) - - data = [] - for i in range(0, length): - c, = struct.unpack('= (3,0): - a = ''.join([chr(c) for c in data]) - else: - a = ''.join([unichr(c) for c in data]) - - return a - -if __name__ == '__main__': - a=read_qstring(open(sys.argv[1], 'rb')) - print(a) +#! /bin/env python +# +# Michael Gibson 23 April 2015 + + +import struct, os + +def read_qstring(fid): + """Read Qt style QString. + + The first 32-bit unsigned number indicates the length of the string (in bytes). + If this number equals 0xFFFFFFFF, the string is null. + + Strings are stored as unicode. + """ + + length, = struct.unpack(' (os.fstat(fid.fileno()).st_size - fid.tell() + 1) : + print(length) + raise Exception('Length too long.') + + # convert length from bytes to 16-bit Unicode words + length = int(length / 2) + + data = [] + for i in range(0, length): + c, = struct.unpack('= (1, 2)): + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 + else: + time_dtype = const.TIMESTAMP_DTYPE_LE_V1_1 + num_samples = header['num_samples_per_data_block'] + data['t_amplifier'] = np.zeros(num_samples * num_datablocks, dtype=time_dtype) + data['amplifier_data'] = np.zeros([header['num_amplifier_channels'], + num_samples * num_datablocks], + dtype=const.AMPLIFIER_DTYPE) + data['aux_input_data'] = np.zeros([header['num_aux_input_channels'], + (num_samples // 4) * num_datablocks], + dtype=const.AUXILIARY_DTYPE) + data['supply_voltage_data'] = np.zeros([header['num_supply_voltage_channels'], + 1 * num_datablocks], + dtype=const.SUPPLY_DTYPE) + data['temp_sensor_data'] = np.zeros([header['num_temp_sensor_channels'], + 1 * num_datablocks], + dtype=const.TEMP_DTYPE) + data['board_adc_data'] = np.zeros([header['num_board_adc_channels'], + num_samples * num_datablocks], + dtype=const.ADC_DTYPE) + data['board_dig_in_data'] = np.zeros([header['num_board_dig_in_channels'], + num_samples * num_datablocks], + dtype=digital_io_data_dtype) + data['board_dig_in_raw'] = np.zeros(num_samples * num_datablocks, + dtype=const.DIG_IN_DTYPE) + data['board_dig_out_data'] = np.zeros([header['num_board_dig_out_channels'], + num_samples * num_datablocks], + dtype=digital_io_data_dtype) + data['board_dig_out_raw'] = np.zeros(num_samples * num_datablocks, + dtype=const.DIG_OUT_DTYPE) + return data + +def read_data_blocks(data, header, fid, datablocks_per_chunk=1): + """Reads a number of data blocks from fid into data. + + Args: + data (dict of numpy arrays): having the same format as the return value + of preallocate_memory() above + header (dict): metadata + fid (file object): file to read from + datablocks_per_chunk (int): how many datablocks to read in + """ + all_names = ['time', + 'amp', + 'aux', + 'supply', + 'temp', + 'adc', + 'digin', + 'digout'] + + # In version 1.2, Intan moved from saving timestamps as unsigned + # integers to signed integers to accommodate negative (adjusted) + # timestamps for pretrigger data. + if ((header['version']['major'], header['version']['minor']) >= (1, 2)): + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 + else: + time_dtype = const.TIMESTAMP_DTYPE_LE_V1_1 + + all_dtypes = [time_dtype, + const.AMPLIFIER_DTYPE, + const.AUXILIARY_DTYPE, + const.SUPPLY_DTYPE, + const.TEMP_DTYPE, + const.ADC_DTYPE, + const.DIG_IN_DTYPE, + const.DIG_OUT_DTYPE] + + num_time_chans = 1 + all_chans = [num_time_chans, + header['num_amplifier_channels'], + header['num_aux_input_channels'], + header['num_supply_voltage_channels'], + header['num_temp_sensor_channels'], + header['num_board_adc_channels'], + header['num_board_dig_in_channels'], + header['num_board_dig_out_channels']] + + num_samples = header['num_samples_per_data_block'] + all_samples = [num_samples, + num_samples, + num_samples // 4, + 1, + 1, + num_samples, + num_samples, + num_samples] + + # create a structured dtype for one datablock + db_dtype = [(name, (dt, chans * samples)) + for name, dt, chans, samples + in zip(all_names, all_dtypes, all_chans, all_samples)] + + chunk = np.fromfile(fid, dtype=np.dtype(db_dtype), count=datablocks_per_chunk) + + for idx, db in enumerate(chunk): + + # Timebase + start = idx * num_samples + stop = (idx + 1) * num_samples + data['t_amplifier'][start:stop] = db['time'] + + # Amplifier channels + if header['num_amplifier_channels']: + start = idx * num_samples + stop = (idx + 1) * num_samples + values = db['amp'].reshape(header['num_amplifier_channels'], + num_samples) + data['amplifier_data'][:,start:stop] = values + + # Auxiliary channels + if header['num_aux_input_channels']: + start = idx * (num_samples // 4) + stop = (idx + 1) * (num_samples // 4) + values = db['aux'].reshape(header['num_aux_input_channels'], + num_samples // 4) + data['aux_input_data'][:,start:stop] = values + + # Supply voltage channels + if header['num_supply_voltage_channels'] > 0: + start = idx * 1 + stop = (idx + 1) * 1 + values = db['supply'].reshape(header['num_supply_voltage_channels'], + 1) + data['supply_voltage_data'][:,start:stop] = values + + # Temperature sensor channels + if header['num_temp_sensor_channels'] > 0: + start = idx * 1 + stop = (idx + 1) * 1 + values = db['temp'].reshape(header['num_temp_sensor_channels'], 1) + data['temp_sensor_data'][:,start:stop] = values + + # Board ADC channels + if header['num_board_adc_channels'] > 0: + start = idx * num_samples + stop = (idx + 1) * num_samples + values = db['adc'].reshape(header['num_board_adc_channels'], + num_samples) + data['board_adc_data'][:,start:stop] = values + + # Board digital input channels (packed together) + if header['num_board_dig_in_channels'] > 0: + start = idx * num_samples + stop = (idx + 1) * num_samples + values = db['digin'].reshape(header['num_board_dig_in_channels'], + num_samples) + data['board_dig_in_raw'][start:stop] = values + + # Board digital output channels (packed together) + if header['num_board_dig_out_channels'] > 0: + start = idx * num_samples + stop = (idx + 1) * num_samples + values = db['digout'].reshape(header['num_board_dig_out_channels'], + num_samples) + data['board_dig_out_raw'][start:stop] = values diff --git a/bark/io/rhd/read_header.py b/bark/io/rhd/read_header.py index de1c7a4..45edd51 100644 --- a/bark/io/rhd/read_header.py +++ b/bark/io/rhd/read_header.py @@ -1,145 +1,186 @@ -#! /bin/env python -# -# Michael Gibson 23 April 2015 - -import sys, struct -from .qstring import read_qstring - - -def read_header(fid): - """Reads the Intan File Format header from the given file.""" - - # Check 'magic number' at beginning of file to make sure this is an Intan - # Technologies RHD2000 data file. - magic_number, = struct.unpack('= 1) or ( - version['major'] > 1): - header['num_temp_sensor_channels'], = struct.unpack('= 3)) or (version['major'] > 1): - header['eval_board_mode'], = struct.unpack(' 0) and (signal_group_enabled > 0): - for signal_channel in range(0, signal_group_num_channels): - new_channel = {'port_name': signal_group_name, - 'port_prefix': signal_group_prefix, - 'port_number': signal_group} - new_channel['native_channel_name'] = read_qstring(fid) - new_channel['custom_channel_name'] = read_qstring(fid) - (new_channel['native_order'], new_channel['custom_order'], - signal_type, channel_enabled, new_channel['chip_channel'], - new_channel['board_stream']) = struct.unpack(' LAST_TESTED_MAJOR_VERSION: + print('Warning: this converter has only been tested up to major', + 'version {}.\n'.format(LAST_TESTED_MAJOR_VERSION)) + + freq = {} + + # Read information of sampling rate and amplifier frequency settings. + header['sample_rate'], = struct.unpack('= (1, 1): + (header['num_temp_sensor_channels'],) = struct.unpack('= (1, 3): + (header['eval_board_mode'],) = struct.unpack(' 1: + header['reference_channel'] = read_qstring(fid) + + # Place frequency-related information in data structure. + # (Note: much of this structure is set above) + freq['amplifier_sample_rate'] = header['sample_rate'] + freq['aux_input_sample_rate'] = header['sample_rate'] / 4 + freq['supply_voltage_sample_rate'] = (header['sample_rate'] / + header['num_samples_per_data_block']) + freq['board_adc_sample_rate'] = header['sample_rate'] + freq['board_dig_in_sample_rate'] = header['sample_rate'] + + header['frequency_parameters'] = freq + + # Create lists for each type of data channel. + header['spike_triggers'] = [] + header['amplifier_channels'] = [] + header['aux_input_channels'] = [] + header['supply_voltage_channels'] = [] + header['board_adc_channels'] = [] + header['board_dig_in_channels'] = [] + header['board_dig_out_channels'] = [] + + # Read signal summary from data file header. + + (number_of_signal_groups,) = struct.unpack(' 0 and signal_group_enabled > 0: + for signal_channel in range(0, signal_group_num_channels): + new_channel = {'port_name': signal_group_name, + 'port_prefix': signal_group_prefix, + 'port_number': signal_group} + new_channel['native_channel_name'] = read_qstring(fid) + new_channel['custom_channel_name'] = read_qstring(fid) + (new_channel['native_order'], + new_channel['custom_order'], + signal_type, + channel_enabled, + new_channel['chip_channel'], + new_channel['board_stream']) = struct.unpack('= 2) or (header['version']['major'] > 1): - data['t_amplifier'][indices['amplifier']:(indices['amplifier']+60)] = np.array(struct.unpack('<' + 'i' *60, fid.read(240))) - else: - data['t_amplifier'][indices['amplifier']:(indices['amplifier']+60)] = np.array(struct.unpack('<' + 'I' *60, fid.read(240))) - - if header['num_amplifier_channels'] > 0: - tmp = np.fromfile(fid, dtype='uint16', count=60 * header['num_amplifier_channels']) - data['amplifier_data'][range(header['num_amplifier_channels']), indices['amplifier']:(indices['amplifier']+60)] = tmp.reshape(header['num_amplifier_channels'], 60) - - if header['num_aux_input_channels'] > 0: - tmp = np.fromfile(fid, dtype='uint16', count=15 * header['num_aux_input_channels']) - data['aux_input_data'][range(header['num_aux_input_channels']), indices['aux_input']:(indices['aux_input']+15)] = tmp.reshape(header['num_aux_input_channels'], 15) - - if header['num_supply_voltage_channels'] > 0: - tmp = np.fromfile(fid, dtype='uint16', count=1 * header['num_supply_voltage_channels']) - data['supply_voltage_data'][range(header['num_supply_voltage_channels']), indices['supply_voltage']:(indices['supply_voltage']+1)] = tmp.reshape(header['num_supply_voltage_channels'], 1) - - if header['num_temp_sensor_channels'] > 0: - tmp = np.fromfile(fid, dtype='uint16', count=1 * header['num_temp_sensor_channels']) - data['temp_sensor_data'][range(header['num_temp_sensor_channels']), indices['supply_voltage']:(indices['supply_voltage']+1)] = tmp.reshape(header['num_temp_sensor_channels'], 1) - - if header['num_board_adc_channels'] > 0: - tmp = np.fromfile(fid, dtype='uint16', count=60 * header['num_board_adc_channels']) - data['board_adc_data'][range(header['num_board_adc_channels']), indices['board_adc']:(indices['board_adc']+60)] = tmp.reshape(header['num_board_adc_channels'], 60) - - if header['num_board_dig_in_channels'] > 0: - data['board_dig_in_raw'][indices['board_dig_in']:(indices['board_dig_in']+60)] = np.array(struct.unpack('<' + 'H' *60, fid.read(120))) - - if header['num_board_dig_out_channels'] > 0: - data['board_dig_out_raw'][indices['board_dig_out']:(indices['board_dig_out']+60)] = np.array(struct.unpack('<' + 'H' *60, fid.read(120))) - diff --git a/bark/io/rhd/rhd2bark.py b/bark/io/rhd/rhd2bark.py index e27daec..d10e1ee 100644 --- a/bark/io/rhd/rhd2bark.py +++ b/bark/io/rhd/rhd2bark.py @@ -1,24 +1,27 @@ import sys import os.path import arrow +import itertools from dateutil import tz import numpy as np -from bark.io.rhd.load_intan_rhd_format import read_data +import bark.io.rhd.load_intan_rhd_format as lirf +import bark.io.rhd.legacy_load_intan_rhd_format as legacy_lirf from bark import create_entry, write_metadata +DEFAULT_MAX_MEM = '3GB' def bark_rhd_to_entry(): import argparse default_max_gaps = 10 p = argparse.ArgumentParser( - description="""Create a Bark entry from RHD files + description="""Create a Bark entry from RHD files. RHD files should be contiguous in time. An error is raised if the RHD files do not all have the same channels recorded. """) p.add_argument("rhdfiles", help="RHD file(s) to convert", nargs="+") - p.add_argument("-o", "--out", help="name of bark entry") + p.add_argument("-o", "--out", help="name of bark entry", required=True) p.add_argument("-a", "--attributes", action='append', @@ -28,7 +31,7 @@ def bark_rhd_to_entry(): p.add_argument( "-t", "--timestamp", - help="""format: YYYY-MM-DD or YYYY-MM-DDTHH:MM:SS.S, if left unspecified + help="""format: YYYY-MM-DD or YYYY-MM-DDTHH:MM:SS.S; if left unspecified the timestamp will be inferred from the filename of the first RHD file.""") p.add_argument('--timezone', @@ -38,7 +41,7 @@ def bark_rhd_to_entry(): p.add_argument( "-p", "--parents", - help="No error if already exists, new meta-data written, \ + help="No error if entry already exists, new meta-data written, \ and datasets will be overwritten.", action="store_true") p.add_argument( @@ -48,12 +51,35 @@ def bark_rhd_to_entry(): .format(default_max_gaps), type=int, default=default_max_gaps) + p.add_argument( + "-l", + "--legacy", + help="Use original Intan-derived code (slower)", + action="store_true") + p.add_argument( + "-m", + "--max-memory", + help='Rough maximum memory usage (default: "{}")'.format(DEFAULT_MAX_MEM), + default=DEFAULT_MAX_MEM) args = p.parse_args() attrs = dict(args.keyvalues) if args.keyvalues else {} check_exists(args.rhdfiles) + max_mem = max_mem_from_string(args.max_memory) rhds_to_entry(args.rhdfiles, args.out, args.timestamp, args.parents, - args.maxgaps, args.timestamp, **attrs) - + args.maxgaps, args.timestamp, legacy=args.legacy, max_mem=max_mem, **attrs) + +def max_mem_from_string(s): + mem_dict = {'kb': 10, 'mb': 20, 'gb': 30} + if len(s) < 2: + raise ValueError('{} is not a valid amount of memory'.format(s)) + s = ''.join(s.split()).lower() + if s[-2:] in mem_dict: + pwr = mem_dict[s[-2:]] + elif s[-1] == 'b': + pwr = 0 + else: + raise ValueError('{} is not a valid amount of memory'.format(s)) + return int(float(s[:-2])) * (2 ** pwr) def rhd_filename_to_timestamp(fname, timezone): return arrow.get(fname, 'YYMMDD_HHmmss').replace( @@ -96,6 +122,8 @@ def amplifier_metadata(result, dsetname): sampling_rate=result['frequency_parameters'][ 'amplifier_sample_rate'], ) attrs.update(result['frequency_parameters']) + if 'digital_reference_channel' in result: + attrs['digital_reference_channel'] = result['digital_reference_channel'] columns = {i: chan_attrs for i, chan_attrs in enumerate(result['amplifier_channels'])} for k in columns: @@ -117,20 +145,19 @@ def not_implemented_warnings(result): print("TEMP SENSOR DATA CONVERSION NOT YET IMPLEMENTED") -def check_timestamp_gaps(data, max_gaps): - num_gaps = np.sum(~np.isclose( - np.diff(data['t_amplifier']), 1. / data['frequency_parameters'][ - 'amplifier_sample_rate'])) - if num_gaps > max_gaps: - raise Exception("{} data gaps exceeds maximum limit {}".format( - num_gaps, max_gaps)) +def count_timestamp_gaps(data, last_chunk_last_timestamp): + epsilon = 1. / data['frequency_parameters']['amplifier_sample_rate'] + num_gaps = 0 + if not np.isclose(data['t_amplifier'][0] - last_chunk_last_timestamp, + epsilon): + num_gaps += 1 + num_gaps += np.sum(~np.isclose(np.diff(data['t_amplifier']), epsilon)) + return num_gaps def check_exists(rhd_paths): - for filepath in rhd_paths: - if not os.path.exists(filepath): - print("file {} does not exist".format(filepath)) - sys.exit(0) + if not all(os.path.exists(filepath) for filepath in rhd_paths): + raise FileNotFoundError(filepath) def rhds_to_entry(rhd_paths, @@ -139,6 +166,8 @@ def rhds_to_entry(rhd_paths, parents=False, max_gaps=10, timezone='America/Chicago', + legacy=False, + max_mem=DEFAULT_MAX_MEM, **attrs): """ Converts a temporally contiguous list of .rhd files to a bark entry. @@ -147,56 +176,80 @@ def rhds_to_entry(rhd_paths, timestamp = rhd_filename_to_timestamp(rhd_paths[0], timezone) else: timestamp = input_string_to_timestamp(timestamp, timezone) - # extract data and metadata from first file - print(rhd_paths[0]) - result = read_data(rhd_paths[0], no_floats=True) - not_implemented_warnings(result) - check_timestamp_gaps(result, max_gaps) - # make entry - entry_attrs = result['notes'] - attrs.update(entry_attrs) + # process first file and create entry and datasets as needed + first, rest = data_feed(rhd_paths[0], legacy, max_mem) + attrs.update(first['notes']) + attrs['intan_gui_version'] = '{}.{}'.format(first['version']['major'], + first['version']['minor']) create_entry(entry_name, timestamp, parents, **attrs) - # make datasets - board_channels = adc_chan_names(result) - if board_channels: - dsetname = os.path.join(entry_name, 'board_adc.dat') - board_adc_metadata(result, dsetname) - with open(dsetname, 'wb') as fp: - fp.write(result['board_adc_data'].T.tobytes()) - - amplifier_channels = amp_chan_names(result) - if amplifier_channels: - dsetname = os.path.join(entry_name, 'amplifier.dat') - amplifier_metadata(result, dsetname) - with open(dsetname, 'wb') as fp: - fp.write(result['amplifier_data'].T.tobytes()) - - # now that the metadata has been written (and data from the first file) - # write data for the remainder of the files - for rhdfile in rhd_paths[1:]: - print(rhdfile) - result = read_data(rhdfile, no_floats=True) - not_implemented_warnings(result) - check_timestamp_gaps(result, max_gaps) - cur_board_channels = adc_chan_names(result) - cur_amplifier_channels = amp_chan_names(result) - - # check that the same channels are being recorded - if board_channels != cur_board_channels: - raise ValueError("""{} has channels {} - {} has channels {} """.format( - rhdfile, cur_board_channels, rhd_paths[0], board_channels)) - if amplifier_channels != cur_amplifier_channels: - raise ValueError("""{} has channels {} - {} has channels {}""" - .format(rhdfile, cur_amplifier_channels, - rhd_paths[0], amplifier_channels)) - # write data - if cur_board_channels: - dsetname = os.path.join(entry_name, 'board_adc.dat') - with open(dsetname, 'ab') as fp: - fp.write(result['board_adc_data'].T.tobytes()) - if cur_amplifier_channels: - dsetname = os.path.join(entry_name, 'amplifier.dat') - with open(dsetname, 'ab') as fp: - fp.write(result['amplifier_data'].T.tobytes()) + adc_channels = adc_chan_names(first) + amp_channels = amp_chan_names(first) + adc_dset_name = None + amp_dset_name = None + if adc_channels: + adc_dset_name = os.path.join(entry_name, 'board_adc.dat') + board_adc_metadata(first, adc_dset_name) + open(adc_dset_name, 'wb').close() + if amp_channels: + amp_dset_name = os.path.join(entry_name, 'amplifier.dat') + amplifier_metadata(first, amp_dset_name) + open(amp_dset_name, 'wb').close() + write_data_feed(first, rest, adc_dset_name, amp_dset_name, max_gaps) + # process the rest of the files + for rhd_file in rhd_paths[1:]: + first, rest = data_feed(rhd_file, legacy, max_mem) + # check that channels are all the same as first file + for curr, old in zip((adc_channels, amp_channels), + (adc_chan_names(first), amp_chan_names(first))): + if curr != old: + msg = '{} has channels {}\n{} has channels {}' + raise ValueError(msg.format(rhd_file, curr, rhd_paths[0], old)) + write_data_feed(first, rest, adc_dset_name, amp_dset_name, max_gaps) + +def data_feed(rhd_file, legacy, max_memory): + """Set up a stream to feed data from rhd_file in chunks. + + Args: + rhd_file (str): filename + legacy (bool): whether to use legacy Intan-provided code + max_memory (number): memory size of chunks, in bytes + + Returns: + tuple(dict, iterable): first chunk in stream, plus rest of stream + """ + print(rhd_file) + if legacy: + # the legacy code reads the entire file's contents into memory at once, + # so there's nothing left after the first chunk + return (legacy_lirf.read_data(rhd_file, no_floats=True), []) + else: + feed = lirf.read_data(rhd_file, no_floats=True, max_memory=max_memory) + return (next(feed), feed) + +def write_data_feed(first, rest, adc_fn, amp_fn, max_gaps): + """Write a data feed to disk. + + Args: + first (dict): first chunk in the data feed + rest (iterable): rest of the data feed + adc_fn (str or None): filename of the ADC bark dataset + amp_fn (str or None): filename of the amplifier bark dataset + max_gaps (int): maximum number of "non-small" gaps in the timestamps + that will be tolerated + """ + not_implemented_warnings(first) + last_timestamp = first['t_amplifier'][0] + timestamp_gaps = 0 + for data_chunk in itertools.chain([first], rest): + timestamp_gaps += count_timestamp_gaps(data_chunk, last_timestamp) + last_timestamp = data_chunk['t_amplifier'][-1] + if adc_fn: + with open(adc_fn, 'ab') as fp: + data_chunk['board_adc_data'].T.tofile(fp) + if amp_fn: + with open(amp_fn, 'ab') as fp: + data_chunk['amplifier_data'].T.tofile(fp) + if timestamp_gaps > max_gaps: + msg = '{} timestamp gaps in data exceed maximum limit {}' + raise Exception(msg.format(timestamp_gaps, max_gaps)) + diff --git a/bark/io/spykingcircus.py b/bark/io/spykingcircus.py index 819928c..ec11881 100644 --- a/bark/io/spykingcircus.py +++ b/bark/io/spykingcircus.py @@ -1,109 +1,58 @@ -import argparse +# original by KJB + import bark -import collections -import h5py import numpy as np -import os -import pandas -import sys - -SC_GRADES_DICT = {0.0: 'O', - 1.0: 'F', - 2.0: 'E', - 3.0: 'D', - 4.0: 'C', - 5.0: 'B', - 6.0: 'A'} -SC_TEMPLATE_PREFIX = 'temp_' - -SpikeEvent = collections.namedtuple('SpikeEvent', ['name', 'time', 'amplitude']) - -def get_sc_path(entry_fn, dataset, sc_suffix, sc_filename): - sc_dir = os.path.splitext(dataset)[0] - if sc_suffix is None: - sc_suffix = '' - else: - sc_suffix = '-' + sc_suffix - fn = sc_dir + '.spyc.' + sc_filename + sc_suffix + '.hdf5' - return os.path.join(entry_fn, sc_dir, fn) - -def unique_temp_name(tn): - return tn[len(SC_TEMPLATE_PREFIX):] - -def long_temp_name(stn): - return SC_TEMPLATE_PREFIX + stn - -def extract_sc(entry_fn, dataset, sc_suffix, out_fn): - sr = bark.read_metadata(os.path.join(entry_fn, dataset))['sampling_rate'] - # determine file names - results_path = get_sc_path(entry_fn, dataset, sc_suffix, 'result') - templates_path = get_sc_path(entry_fn, dataset, sc_suffix, 'templates') - # extract times and amplitudes - with h5py.File(results_path, 'r') as rf: - cluster_times = {unique_temp_name(name): np.array(indices).astype(float) / sr - for name,indices in rf['spiketimes'].items()} - cluster_amplitudes = {unique_temp_name(name): np.array(amplitudes) - for name,amplitudes in rf['amplitudes'].items()} - cluster_names = sorted(cluster_times.keys(), key=int) - event_list = [] - for n in cluster_names: - event_list.extend([SpikeEvent(n, time[0], amp[0]) - for time,amp in zip(cluster_times[n], cluster_amplitudes[n])]) - event_list.sort(key=lambda se: se.time) - # extract grades and center pad - with h5py.File(templates_path, 'r') as tf: - cluster_grades = [SC_GRADES_DICT[tag[0]] for tag in tf['tagged']] - cluster_grades = {n: cluster_grades[idx] for idx,n in enumerate(cluster_names)} - NUM_TEMPLATES = int(tf['temp_shape'][2][0] / 2) - NUM_CHANNELS = int(tf['temp_shape'][0][0]) - NUM_SAMPLES = int(tf['temp_shape'][1][0]) - CHAN_BY_SAMPLE = NUM_CHANNELS * NUM_SAMPLES - full_templates = {} - for t in range(NUM_TEMPLATES): - y_vals = tf['temp_y'][0] == t - x_vals = tf['temp_x'][:,y_vals][0].astype(int) - reconst = np.zeros(CHAN_BY_SAMPLE) - for loc in x_vals: - reconst[loc] = tf['temp_data'][:,loc][0] - reshaped = reconst.reshape((NUM_CHANNELS, -1)) - full_templates[t] = np.copy(reshaped) - center_channel = {} - for t in full_templates: - # note that this assumes negative-going spikes - min_across_channels = list(np.amin(full_templates[t], axis=1)) - total_min = min(min_across_channels) - center_channel[str(t)] = min_across_channels.index(total_min) - # write times and amplitudes to event dataset +import pandas as pd +import os.path + + +def create_data(guifolder, sampling_rate): + times = np.load(os.path.join(guifolder, 'spike_times.npy')) / sampling_rate + amplitude = np.load(os.path.join(guifolder, 'amplitudes.npy')) + name = np.load(os.path.join(guifolder, 'spike_clusters.npy')) + positions = np.load(os.path.join(guifolder, 'channel_positions.npy') + ) # not used? might be needed for missing channels + data = pd.DataFrame({'name': name, 'start': times, 'amplitude': amplitude}) + return data + + +def create_metadata(guifolder): + cluster_groups = pd.read_csv( + os.path.join(guifolder, 'cluster_groups.csv'), + sep='\t') + templates = np.load(guifolder + '/templates.npy') + channel = np.argmax(np.argmax(np.abs(templates), 1), 1) + template_dict = {int(name): {'score': cluster_groups.group[i], + 'center_channel': int(channel[i])} + for i, name in enumerate(cluster_groups.cluster_id)} attrs = {'columns': {'start': {'units': 's'}, 'name': {'units': None}, 'amplitude': {'units': None}}, 'datatype': 1001, - 'sampling_rate': sr, - 'templates': {name: {'score': cluster_grades[name], - 'sc_name': long_temp_name(name), - 'center_channel': center_channel[name]} - for name in cluster_names}} - return bark.write_events(os.path.join(entry_fn, out_fn), - pandas.DataFrame({'start': [event.time for event in event_list], - 'name': [event.name for event in event_list], - 'amplitude': [event.amplitude for event in event_list]}), - **attrs) + 'templates': template_dict, + 'filetype': 'csv', + 'creator': 'phy', } + return attrs + + +def main(): + import argparse + p = argparse.ArgumentParser(description=''' + Convert Spyking Circus PHY GUI output to Bark event dataset. + ''') + p.add_argument('phydir', help='directory containing PHY GUI output files') + p.add_argument('out', help='name of output event dataset') + p.add_argument('-r', + '--rate', + required=True, + type=float, + help='sampling rate of original data') + args = p.parse_args() + data = create_data(args.phydir, args.rate) + attrs = create_metadata(args.phydir) + bark.write_events(args.out, data, **attrs) -def _parse_args(raw_args): - desc = 'Extract Spyking Circus spike-sorting info into a Bark-readable form.' - epi = 'Assumes standard Spyking Circus naming conventions for output files.\n' - epi += 'Currently only supports extraction of spike times, amplitudes, and quality tags.' - parser = argparse.ArgumentParser(description=desc, epilog=epi) - dflt = 'sc_template_times.csv' - parser.add_argument('-o', '--out', help='output filename (default: {})'.format(dflt), default=dflt) - parser.add_argument('-s', '--suffix', help='Spyking Circus suffix (if applicable)') - parser.add_argument('entry', help='Bark entry path') - parser.add_argument('dataset', help='sampled dataset name to generate Spyking Circus filenames') - return parser.parse_args(raw_args) - -def _main(): - parsed_args = _parse_args(sys.argv[1:]) - extract_sc(parsed_args.entry, parsed_args.dataset, parsed_args.suffix, parsed_args.out) if __name__ == '__main__': - _main() + main() + diff --git a/bark/stream.py b/bark/stream.py index e01ca3b..59b0426 100644 --- a/bark/stream.py +++ b/bark/stream.py @@ -13,9 +13,9 @@ def array_iterator(data, chunksize): try: result = data[index:index + chunksize] except IndexError: - raise StopIteration + return if result.shape[0] == 0: - raise StopIteration + return yield result index += chunksize diff --git a/bark/tools/B_PLot.py b/bark/tools/B_PLot.py new file mode 100644 index 0000000..5988c46 --- /dev/null +++ b/bark/tools/B_PLot.py @@ -0,0 +1,432 @@ +import os +import sys +import string +import yaml +import numpy as np +from scipy.signal import spectrogram +import matplotlib.pyplot as plt +import bark +from bark.io.eventops import (OpStack, write_stack, read_stack, Update, Merge, + Split, Delete, New) +import warnings +warnings.filterwarnings('ignore') # suppress matplotlib warnings +from bark.tools.spectral import BarkSpectra + + +help_string = ''' + +Shortcuts +--------- +any letter or number annotate segment +up arrow zoom out +down arrow zoom in +up left move left +down right move right +click on map to move through the graphic +drag the graphic to move the graphic + +''' + +zoom_size = 50000 + +# kill all the shorcuts +def kill_shortcuts(plt): + plt.rcParams['keymap.all_axes'] = '' + plt.rcParams['keymap.back'] = '' + plt.rcParams['keymap.forward'] = '' + plt.rcParams['keymap.fullscreen'] = '' + plt.rcParams['keymap.grid'] = '' + plt.rcParams['keymap.home'] = '' + plt.rcParams['keymap.pan'] = '' + #plt.rcParams['keymap.quit'] = '' + plt.rcParams['keymap.save'] = '' + plt.rcParams['keymap.xscale'] = '' + plt.rcParams['keymap.yscale'] = '' + plt.rcParams['keymap.zoom'] = '' + + +class Plot: + def __init__(self, + ax, + sampled): + + self.ax = ax + self.data = sampled.data.ravel() + self.sr = sampled.sampling_rate + + def update_x_axis(self,start,stop): + self.ax.set_xlim(start,stop) + + def clear_plot(self): + self.ax.cla() + + ''' + class Osc_Plot + + parameter: + + ax : axis object to plot spectrogram on + sampled : the data read from .dat file including sound data and sample rate + opstack : the opstack stack contains label information + + ''' +class Osc_Plot(Plot): + def __init__(self, + ax, + sampled): + Plot.__init__(self,ax,sampled) + self.N_points = 35000 + self.ax.set_axis_bgcolor('k') + self.ax.tick_params(axis='x', + which='both', + bottom='off', + top='off', + labelbottom='off') + self.osc_line, = self.ax.plot( + np.arange(self.N_points), + np.zeros(self.N_points), + color='gray') + self.ax.figure.tight_layout() + + def update_oscillogram(self,buf_start,buf_stop): + self.selected_boundary = None + self.buffer_start_samp = buf_start + self.buffer_stop_samp = buf_stop + self.buf_start = self.buffer_start_samp / self.sr + self.buf_stop = self.buffer_stop_samp / self.sr + self.update_oscillo() + + + def update_oscillo(self): + + x = self.data[self.buffer_start_samp:self.buffer_stop_samp] + t = np.arange(len(x)) / self.sr + self.buf_start + if len(x) > 10000: + t_interp = np.linspace(self.buf_start, self.buf_stop, 10000) + x_interp = np.interp(t_interp, t, x) + else: + t_interp = t + x_interp = x + self.osc_line.set_data(t_interp, x_interp) + self.ax.set_xlim(self.buf_start, self.buf_stop) + self.ax.set_ylim(min(x), max(x)) + + ''' + class Spec_Plot + + parameter: + + ax : axis object to plot spectrogram on + sampled : the data read from .dat file including sound data and sample rate + + ''' + +class Spec_Plot(Plot): + def __init__(self, + ax, + sampled): + Plot.__init__(self,ax,sampled) + self.ax.set_axis_bgcolor('k') + + def update_spectrogram(self,start,stop): + self.ax.clear() + self.plot_spectrogram(self.data, + self.sr, + start, + stop, + ax=self.ax) + self.ax.set_xlim(start, stop) + + def plot_spectrogram(self,data, + sr, + start, + stop, + ms_nfft=15, + ax=None, + lowfreq=300, + highfreq=8000, + n_tapers=2, + NW=1.5, + derivative=True, + window=('kaiser', 8), + **kwargs): + ''' + data : a vector of samples, first sample starts at time = 0 + sr : sampling rate + start : start time to slice data, units: seconds + stop : stop time to slice data, units: seconds + ms_nfft : width of fourier transform in milliseconds + ax : axis object to plot spectrogram on. + lowfreq : lowest frequency to plot + highfreq : highest frequency to plot + n_tapers : Number of tapers to use in a custom multi-taper Fourier + transform estimate + NW : multi-taper bandwidth parameter for custom multi-taper Fourier + transform estimate increasing this value reduces side-band + ripple, decreasing sharpens peaks + derivative: if True, plots the spectral derivative, SAP style + + ''' + nfft = int(ms_nfft / 1000. * sr) + start_samp = int(start * sr) - nfft // 2 + if start_samp < 0: + start_samp = 0 + stop_samp = int(stop * sr) - nfft // 2 + x = data[start_samp:stop_samp] + + # determine overlap based on screen size. + # We don't need more points than pixels + pixels = 1000 + samples_per_pixel = int((stop - start) * sr / pixels) + noverlap = max(nfft - samples_per_pixel, 0) + + from matplotlib import colors + + spa = BarkSpectra(sr, + NFFT=nfft, + noverlap=noverlap, + data_window=int(0.01 * sr), + n_tapers=n_tapers, + NW=NW, + freq_range=(lowfreq, highfreq)) + spa.signal(x) + pxx, f, t, thresh = spa.spectrogram(ax=ax, derivative=derivative) + + # calculate the parameter for the plot + freq_mask = (f > lowfreq) & (f < highfreq) + fsub = f[freq_mask] + Sxxsub = pxx[freq_mask, :] + t += start + + # plot the spectrogram + if derivative: + image = ax.pcolorfast(t, + fsub, + Sxxsub, + cmap='inferno', + norm=colors.SymLogNorm(linthresh=thresh)) + else: + image = ax.pcolorfast(t, + fsub, + Sxxsub, + cmap='inferno', + norm=colors.LogNorm(vmin=thresh)) + + plt.sca(ax) + plt.ylim(lowfreq, highfreq) + return image + + ''' + class Minimap_Plot + + parameter: + + ax : axis object to plot spectrogram on + opstack : the opstack stack contains label information + + ''' + +class Minimap_Plot: + + def __init__(self, + ax, + max_time): + self.max_time = max_time + self.ax = ax + self.ax.set_axis_bgcolor('k') + self.ax.tick_params(axis='y', + which='both', + left='off', + right='off', + labelleft='off') + self.ax.set_xlim(0, self.max_time) + self.current = self.ax.axvline(color='r') + def update_minimap(self, x_Data): + self.current.set_xdata((x_Data, x_Data)) + + + + + +class SegmentReviewer: + def __init__(self, + osc_ax, + spec_ax, + map_ax, + sampled): + self.canvas = osc_ax.get_figure().canvas + + self.osc = Osc_Plot(ax = osc_ax, sampled = sampled) + + self.spec = Spec_Plot(ax = spec_ax, sampled = sampled) + + self.data = sampled.data.ravel() + self.sr = sampled.sampling_rate + self.max_time =int(round(len(self.data)/self.sr)) + self.N_points = 300000 + self.window_size = self.N_points/self.sr + self.start = 0 + self.stop = self.window_size + self.press_flag = 0 + zoom_size = self.sr*5 + self.map = Minimap_Plot(ax = map_ax, max_time = self.max_time) + + self.update_plot_data() + + + def update_plot_data(self): + 'updates plot data on all three axes' + sr = self.sr + start_samp = int(self.start * sr) + stop_samp = int(self.stop * sr) + syl_samps = stop_samp - start_samp + self.buffer_start_samp = start_samp - (self.N_points - syl_samps) // 2 + if self.buffer_start_samp < 0: + self.buffer_start_samp = 0 + self.buffer_stop_samp = self.buffer_start_samp + self.N_points + + if self.buffer_stop_samp >= self.data.shape[0]: + self.buffer_stop_samp = self.data.shape[0] - 1 + + self.buf_start = self.buffer_start_samp / sr + self.buf_stop = self.buffer_stop_samp / sr + + self.spec.update_spectrogram(self.buf_start,self.buf_stop) + self.osc.update_oscillogram(self.buffer_start_samp,self.buffer_stop_samp) + self.map.update_minimap((self.start+self.stop)/2) + + self.canvas.draw() + + def connect(self): + 'creates all the event connections' + self.cid_key_press = self.canvas.mpl_connect('key_press_event', + self.on_key_press) + self.cid_mouse_press = self.canvas.mpl_connect('button_press_event', + self.on_mouse_press) + self.cid_mouse_motion = self.canvas.mpl_connect('motion_notify_event', + self.on_mouse_motion) + self.cid_mouse_release = self.canvas.mpl_connect( + 'button_release_event', self.on_mouse_release) + + def on_mouse_press(self, event): + if event.inaxes == self.map.ax: + self.start = event.xdata - self.window_size/2 + self.stop = self.start + self.window_size + if self.start < 0: + self.start = 0 + self.stop = self.window_size + if self.stop > self.max_time: + self.start = self.max_time - self.window_size + self.stop = self.max_time + self.update_plot_data() + if event.inaxes in (self.osc.ax, self.spec.ax): + self.press_x = event.xdata + self.press_flag = 1 + self.canvas.draw() + + def on_mouse_motion(self, event): + self.canvas.draw() + + def on_mouse_release(self, event): + if event.inaxes in (self.osc.ax, self.spec.ax) and self.press_flag == 1: + self.start -= (event.xdata - self.press_x) + self.stop = self.start + self.window_size + if self.stop > self.max_time: + self.start = self.stop - self.window_size + self.stop = self.start + if self.start < 0: + self.start = 0 + self.stop = self.window_size + self.update_plot_data() + self.press_flag = 0 + self.canvas.draw() + + def inc_i(self): + self.start += self.window_size/2 + self.stop = self.start + self.window_size + if self.stop > self.max_time: + self.start = self.max_time - self.window_size + self.stop = self.max_time + self.update_plot_data() + + def dec_i(self): + self.start -= self.window_size/2 + self.stop = self.start - self.window_size + if self.start < 0: + self.start = 0 + self.stop = self.window_size + self.update_plot_data() + + def on_key_press(self, event): + # print('you pressed ', event.key) + if event.key in ('pagedown', ' ', 'right'): + self.inc_i() + elif event.key in ('pageup', 'backspace', 'left'): + self.dec_i() + elif event.key in ('ctrl+i', 'down'): + if self.N_points > zoom_size: + self.N_points -= zoom_size + self.stop = self.start + self.window_size + self.window_size = self.N_points/self.sr + self.update_plot_data() + elif event.key in ('ctrl+o', 'up'): + if self.N_points - zoom_size < self.max_time * self.sr : + self.N_points += zoom_size + self.window_size = self.N_points/self.sr + self.stop = self.start + self.window_size + self.update_plot_data() + + + def __enter__(self): + return self + + def __exit__(self, *args): + return 0 + + +def build_shortcut_map(mapfile=None): + allkeys = string.digits + string.ascii_letters + shortcut_map = {x: x for x in allkeys} + # load keys from file + if mapfile: + custom = {str(key): value + for key, value in yaml.load(open(mapfile, 'r')).items()} + print('custom keymaps:', custom) + shortcut_map.update(custom) + return shortcut_map + + + + +def main(datfile): + + kill_shortcuts(plt) + sampled = bark.read_sampled(datfile) + # assert len(sampled.attrs['columns']) == 1 + plt.figure() + # Oscillogram and Spectrogram get + osc_ax = plt.subplot2grid((7, 1), (0, 0), rowspan=3) + spec_ax = plt.subplot2grid((7, 1), (3, 0), rowspan=3, sharex=osc_ax) + map_ax = plt.subplot2grid((7, 1), (6, 0), rowspan=1) + + # Segement review is a context manager to ensure a save prompt + # on exit. see SegmentReviewer.__exit__ + with SegmentReviewer(osc_ax, spec_ax, map_ax, sampled) as reviewer: + reviewer.connect() + plt.show(block=True) + + +def _run(): + import argparse + + p = argparse.ArgumentParser(description=''' + Review and annotate segments + ''') + p.add_argument('dat', help='name of a sampled dataset') + + args = p.parse_args() + main(args.dat) + + +if __name__ == '__main__': + _run() diff --git a/bark/tools/barkforeach.py b/bark/tools/barkforeach.py index c4a0e92..644e04d 100644 --- a/bark/tools/barkforeach.py +++ b/bark/tools/barkforeach.py @@ -1,4 +1,5 @@ import argparse +import glob import subprocess import os import sys @@ -14,12 +15,37 @@ def _parse_args(raw_args): parser.add_argument('entries', nargs='+', help='entries') return parser.parse_args(raw_args) -def bark_for_each(cmd, entry_list, verbose): +def bark_for_each(cmd, entry_list, verbose, base_dir=None): + if base_dir is None: + base_dir = os.getcwd() + ct = 1 + total = len(entry_list) for ename in entry_list: + os.chdir(base_dir) os.chdir(ename) if verbose: - print('Working on ' + ename) - subprocess.run(cmd.split()) + print('Working on {} ({} of {})'.format(ename, ct, total)) + expanded_cmd = glob_command(cmd) + subprocess.run(expanded_cmd) + ct += 1 + +def glob_command(cmd): + expanded_cmd = [] + quote_split = cmd.split('"') + if len(quote_split) % 2 == 0: + raise ValueError('cannot parse command: un-paired quotation marks') + for idx,chunk in enumerate(quote_split): + if idx % 2 == 0: # if the chunk was not enclosed in quotes + for token in chunk.split(): + token = os.path.expanduser(token) + g = glob.glob(token) + if g: + expanded_cmd.extend(g) + else: + expanded_cmd.append(token) + else: # if the chunk was enclosed in quotes + expanded_cmd.append('"' + chunk + '"') + return expanded_cmd def _main(): parsed_args = _parse_args(sys.argv[1:]) diff --git a/bark/tools/barkutils.py b/bark/tools/barkutils.py index 5a1ac98..8857c31 100644 --- a/bark/tools/barkutils.py +++ b/bark/tools/barkutils.py @@ -167,9 +167,10 @@ def rb_select(): stream = bark.read_sampled(fname).toStream() if col_attr: columns = stream.attrs['columns'] - channels = [i - for i in range(len(columns)) - if columns[i][col_attr] in channels] + rev_attr = {col[col_attr]: idx + for idx, col in columns.items() + if col_attr in col} # so you can tag only some channels + channels = [rev_attr[c] for c in channels] else: channels = [int(c) for c in channels] stream[channels].write(outfname) @@ -298,18 +299,27 @@ def _datchunk(): p.add_argument("--seconds", help="specify seconds instead of samples", action='store_true') + p.add_argument("--onecut", + help="only perform the first cut", + action="store_true") args = p.parse_args() - datchunk(args.dat, args.stride, args.seconds) + datchunk(args.dat, args.stride, args.seconds, args.onecut) -def datchunk(dat, stride, use_seconds): +def datchunk(dat, stride, use_seconds, one_cut): + def write_chunk(chunk, attrs, i): + filename = "{}-chunk-{}.dat".format(basename, i) + attrs['offset'] = stride * i + bark.write_sampled(filename, chunk, **attrs) attrs = bark.read_metadata(dat) - sr = attrs['sampling_rate'] if use_seconds: - stride = stride * sr + stride = stride * attrs['sampling_rate'] stride = int(stride) basename = os.path.splitext(dat)[0] - for i, chunk in enumerate(stream.read(dat, chunksize=stride)): - filename = "{}-chunk-{}.dat".format(basename, i) - attrs['offset'] = stride * i - bark.write_sampled(filename, chunk, **attrs) + if one_cut: + sds = bark.read_sampled(dat) + write_chunk(sds.data[:stride,:], attrs, 0) + write_chunk(sds.data[stride:,:], attrs, 1) + else: + for i, chunk in enumerate(stream.read(dat, chunksize=stride)): + write_chunk(chunk, attrs, i) diff --git a/bark/tools/datref.py b/bark/tools/datref.py index b3c6a94..15c5bc7 100644 --- a/bark/tools/datref.py +++ b/bark/tools/datref.py @@ -1,26 +1,25 @@ import numpy as np import bark +import shutil BUF = bark.BUFFER_SIZE #COEF_EST_MAX_SIZE = BUF * 500 def datref(datfile, outfile): - dataset = bark.read_sampled(datfile) - data, params = dataset.data, dataset.attrs - outparams = params.copy() - bark.write_sampled(outfile, data, outparams) + shutil.copyfile(datfile, outfile) + shutil.copyfile(datfile + '.meta.yaml', outfile + '.meta.yaml') outdset = bark.read_sampled(outfile, 'r+') out = outdset.data # determine reference coefficient - n_channels = len(params["columns"]) - coefs = np.zeros((n_channels, len(range(0, len(out), BUF)))) + n_samples, n_channels = out.shape + coefs = np.zeros((n_channels, len(range(0, n_samples, BUF)))) power = np.zeros_like(coefs) - for ith, i in enumerate(range(0, len(out), BUF)): + for ith, i in enumerate(range(0, n_samples, BUF)): + total_mean = np.mean(out[i:i + BUF, :], axis=1) for c in range(n_channels): - refs = np.delete(data[i:i + BUF, :], c, axis=1) # remove col c - ref = np.mean(refs, axis=1) - x = data[i:i + BUF, c] + x = out[i:i + BUF, c] + # this way we avoid re-calculating the entire mean for each channel + ref = (total_mean * n_channels - x) / (n_channels - 1) coefs[c, ith] = np.dot(x, ref) / np.dot(ref, ref) - best_C = np.zeros(n_channels) for c in range(n_channels): c_coefs = coefs[c, :] @@ -29,16 +28,34 @@ def datref(datfile, outfile): best_C[c] = np.nanmean(c_coefs[mask]) print("best reference coefficients: {}".format(best_C)) for i, c in enumerate(best_C): - outparams['columns'][i]['reference_coefficient'] = float(c) - for i in range(0, len(out), BUF): + outdset.attrs['columns'][i]['reference_coefficient'] = float(c) + # we want to avoid re-calculating the median from scratch for each channel + # unfortunately, the "new median after removing an element" calculation + # is less succinct than for the mean + if n_channels % 2 == 0: + median_idx = [int(n_channels / 2) - 1, int(n_channels / 2)] + idx_smaller = [median_idx[0] + 1] # new median if elt removed < median + idx_equal = [median_idx[0]] # new median if elt removed == median + idx_greater = [median_idx[0]] # new median if elt removed > median + else: + median_idx = [int(n_channels / 2)] + idx_smaller = [median_idx[0], median_idx[0] + 1] + idx_equal = [median_idx[0] - 1, median_idx[0] + 1] + idx_greater = [median_idx[0] - 1, median_idx[0]] + for i in range(0, n_samples, BUF): + sorted_buffer = np.sort(out[i:i + BUF, :], axis=1) + total_medians = np.mean(sorted_buffer[:, median_idx], axis=1) + new_med_smaller = np.mean(sorted_buffer[:, idx_smaller], axis=1) + new_med_equal = np.mean(sorted_buffer[:, idx_equal], axis=1) + new_med_greater = np.mean(sorted_buffer[:, idx_greater], axis=1) for c in range(n_channels): - refs = np.delete(data[i:i + BUF, :], c, axis=1) # remove col c - ref = np.mean(refs, axis=1) - x = data[i:i + BUF, c] - out[i:i + BUF, c] = data[i:i + BUF, c] - best_C[c] * np.median( - refs, - axis=1) - bark.write_metadata(outfile, **outparams) + less = np.less(out[i:i + BUF, c], total_medians) + equal = np.equal(out[i:i + BUF, c], total_medians) + greater = np.greater(out[i:i + BUF, c], total_medians) + out[i:i + BUF, c][less] = out[i:i + BUF, c][less] - best_C[c] * new_med_smaller[less] + out[i:i + BUF, c][equal] = out[i:i + BUF, c][equal] - best_C[c] * new_med_equal[equal] + out[i:i + BUF, c][greater] = out[i:i + BUF, c][greater] - best_C[c] * new_med_greater[greater] + bark.write_metadata(outfile, **outdset.attrs) def main(): @@ -47,7 +64,7 @@ def main(): References each channel from the median of all the others """) p.add_argument("dat", help="dat file") - p.add_argument("-o", "--out", help="name of output dat file") + p.add_argument("-o", "--out", help="name of output dat file", required=True) opt = p.parse_args() datref(opt.dat, opt.out) diff --git a/bark/tools/datspike.py b/bark/tools/datspike.py index 44863a0..9a919b1 100644 --- a/bark/tools/datspike.py +++ b/bark/tools/datspike.py @@ -1,5 +1,6 @@ import numpy as np import bark +from bark import stream from scipy.signal import argrelextrema default_order = 5 @@ -12,9 +13,9 @@ def thres_extrema(x, y, thresh): def compute_std(dat): - s = bark.stream.read(dat) + s = stream.read(dat) std = np.zeros(len(s.attrs['columns'])) - for i, x in enumerate(bark.stream.read(dat)): + for i, x in enumerate(stream.read(dat)): std += np.std(x, 0) return std / (i + 1) @@ -51,7 +52,7 @@ def main(dat, csv, thresh, is_std, order=default_order, min_dist=0): n_channels = bark.read_sampled(dat).data.shape[1] threshs = np.ones(n_channels) * thresh print('thresholds:', threshs) - s = bark.stream.read(dat) + s = stream.read(dat) pad_len = order with open(csv, 'w') as fp: fp.write('channel,start\n') @@ -60,7 +61,7 @@ def main(dat, csv, thresh, is_std, order=default_order, min_dist=0): bark.write_metadata(csv, datatype=1000, columns={'channel': {'units': None}, - 'start': {'units', 's'}}, + 'start': {'units': 's'}}, thresholds=threshs, order=order, source=dat) diff --git a/bark/tools/labelview.py b/bark/tools/labelview.py index 4685ee0..66e7c87 100644 --- a/bark/tools/labelview.py +++ b/bark/tools/labelview.py @@ -4,7 +4,6 @@ import yaml import numpy as np from scipy.signal import spectrogram -import matplotlib.pyplot as plt import bark from bark.io.eventops import (OpStack, write_stack, read_stack, Update, Merge, Split, Delete, New) @@ -12,6 +11,34 @@ warnings.filterwarnings('ignore') # suppress matplotlib warnings from bark.tools.spectral import BarkSpectra +if sys.platform == 'darwin': + # Keystrokes aren't correctly captured by many matplotlib backends on + # Mac OS X, including the native Cocoa backend. + # Both Qt5 and Tk capture them (mostly) correctly. Qt5 is a slightly + # better experience, but Tk is fine - and Tk is available out-of-the-box. + import matplotlib + try: + matplotlib.use('Qt5Agg') + import matplotlib.pyplot as plt + except ImportError: # PyQt not installed + matplotlib.use('TkAgg') + import matplotlib.pyplot as plt +else: # default backends on linux and windows are fine + import matplotlib.pyplot as plt + +# Use of control, windows, or command is tricky for cross-platform +# compatibility, as they're the keys most likely to be treated differently +# by different OSes and GUI frameworks. +# This manifests here relating to detecting keydowns during a mouse click. +# Control and command recognition is dodgy, but shift is fine. +# This is a bit of a kludge, but it preserves existing behavior on linux +# and windows. +if sys.platform == 'darwin': + click_meta_char = 'shift' +else: + click_meta_char = 'control' + + help_string = ''' Pressing any number or letter (uppercase or lowercase) will mark a segment. @@ -32,8 +59,8 @@ ctrl+w close click on segment boundary move boundary -ctrl+click inside a segment split segment -ctrl+click outside a segment new segment (TODO) +ctrl+click inside a segment split segment (shift+click on Mac) +ctrl+click outside a segment new segment (shift+click on Mac) click on segment boundaries to adjust them. @@ -201,13 +228,13 @@ def __init__(self, self.update_plot_data() def initialize_plots(self): - self.osc_ax.set_axis_bgcolor('k') + self.osc_ax.set_facecolor('k') self.osc_ax.tick_params(axis='x', which='both', bottom='off', top='off', labelbottom='off') - self.spec_ax.set_axis_bgcolor('k') + self.spec_ax.set_facecolor('k') self.osc_line, = self.osc_ax.plot( np.arange(self.N_points), np.zeros(self.N_points), @@ -224,7 +251,7 @@ def initialize_plots(self): def initialize_minimap(self): times, values = labels_to_scatter_coords(self.opstack.events) - self.map_ax.set_axis_bgcolor('k') + self.map_ax.set_facecolor('k') self.map_ax.scatter(times, values, c=values, @@ -357,18 +384,18 @@ def on_mouse_press(self, event): self.i = i self.update_plot_data() # sylable splitting - elif (event.key == 'control' and event.xdata > start_pos and + elif (event.key == click_meta_char and event.xdata > start_pos and event.xdata < stop_pos): self.opstack.push(Split(self.i, float(event.xdata))) self.update_plot_data() # new syllable before - elif event.key == 'control' and event.xdata < start_pos: + elif event.key == click_meta_char and event.xdata < start_pos: self.opstack.push(New(self.i, name='', start=float(event.xdata), stop=float(event.xdata) + .020)) self.update_plot_data() - elif event.key == 'control' and event.xdata > stop_pos: + elif event.key == click_meta_char and event.xdata > stop_pos: self.opstack.push(New(self.i + 1, name='', start=float(event.xdata), diff --git a/bark/tools/psg_view.py b/bark/tools/psg_view.py new file mode 100644 index 0000000..3745432 --- /dev/null +++ b/bark/tools/psg_view.py @@ -0,0 +1,1105 @@ +from PyQt5.QtWidgets import QMenu, QVBoxLayout, QSizePolicy, QMessageBox, QWidget, QPushButton +from PyQt5.QtGui import QIcon +from PyQt5 import QtGui + +from PyQt5 import QtCore, QtWidgets +from PyQt5.QtCore import Qt +from PyQt5.QtWidgets import QWidget, QInputDialog, QLineEdit, QFileDialog +from PyQt5.QtWidgets import QApplication, QMainWindow, QSizePolicy, QWidget, QComboBox, QLabel, QRadioButton, QCheckBox, QGridLayout, QLineEdit, QScrollArea + + +from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas +from matplotlib.figure import Figure +import matplotlib.pyplot as plt + +import os +import sys +import string +import yaml +import numpy as np + +import bark +from bark.io.eventops import (OpStack, write_stack, read_stack, Update, Merge, + Split, Delete, New) +import warnings +warnings.filterwarnings('ignore') # suppress matplotlib warnings + + +help_string = ''' +Pressing any number or letter (uppercase or lowercase) will mark a segment. + +Shortcuts +--------- +any letter or number annotate segment +ctrl+s saves the annotation data +ctrl+h prints this message +down arrow zoom out +up arrow zoom in +right next segment +left previous segment + +ctrl+z undo last operation +ctrl+y redo +ctrl+w close + +click on the right of the current label next segment +click on the left of the current label previous segment + +The top panel is a map of all label locations. +Click on a label to travel to that location. +Open file->shortcut to set label's shortcut +On close, an operation file and the final event file will be written. +Do not kill from terminal unless you want to prevent a save. + +To create custom label, input the label name in the right toolbar. + + +''' + +color = 'yellow' +color_label = 'white' +color_map = plt.get_cmap('hsv') +fontsize = 12 +default_gap = 3 +psg_scale = 2 #change the scale of the psg_viewer +number_shortcut = {'1':'1', '2':'2','3':'3', '4':'4', '5':'5', '6':'6', '7':'7', '8':'8', '9':'9', '0':'0'} + + +# kill all the shorcuts +def kill_shortcuts(plt): + plt.rcParams['keymap.all_axes'] = '' + plt.rcParams['keymap.back'] = '' + plt.rcParams['keymap.forward'] = '' + plt.rcParams['keymap.fullscreen'] = '' + plt.rcParams['keymap.grid'] = '' + plt.rcParams['keymap.home'] = '' + plt.rcParams['keymap.pan'] = '' + plt.rcParams['keymap.save'] = '' + plt.rcParams['keymap.xscale'] = '' + plt.rcParams['keymap.yscale'] = '' + plt.rcParams['keymap.zoom'] = '' + + +def labels_to_scatter_coords(labels): + times = [x['start'] for x in labels] + values = [] + for record in labels: + name = record['name'] + if not isinstance(name, str) or name == '': + v = 0 + elif name.isdigit(): + v = int(name) + elif name[0].isalpha(): + # alphabet in range 11-36 + v = 133 - ord(name[0].lower()) + else: + v = 37 + values.append(v) + return times, values + +def label_to_scatter(name): + if not isinstance(name, str) or name == '': + v = 0 + elif name.isdigit(): + v = int(name) + elif name[0].isalpha(): + # alphabet in range 11-36 + v = 133 - ord(name[0].lower()) + else: + v = 37 + return v + + +def nearest_label(labels, xdata): + return np.argmin(np.abs(xdata - np.array([x['start'] for x in labels]))) + + +def write_metadata(path, meta='.meta.yaml'): + import codecs + params = {'columns': {'name': {'units': 'null'}, + 'start': {'units': 's'}, + 'stop': {'units': 's'}}, + 'datatype': 2002} + bark.write_metadata(path, meta, **params) + + +def build_shortcut_map(mapfile=None): + allkeys = string.digits + string.ascii_letters + shortcut_map = {x: x for x in allkeys} + # load keys from file + if mapfile: + custom = {str(key): value + for key, value in yaml.load(open(mapfile, 'r')).items()} + print('custom keymaps:', custom) + shortcut_map.update(custom) + return shortcut_map + + +def to_seconds(dset): + 'TODO Converts bark EventData object to units of seconds.' + if 'offset' in dset.attrs and dset.attrs['offset'] != 0: + raise Exception('offsets are not yet supported in event file') + if dset.attrs['columns']['start']['units'] == 's': + pass + elif 'units' in dset.attrs and dset.attrs['units'] == 's': + pass + else: + raise Exception('only units of s are supported in event file') + return dset + + +def load_opstack(opsfile, labelfile, labeldata, use_ops): + load_ops = os.path.exists(opsfile) and use_ops + if load_ops: + opstack = read_stack(opsfile) + print('Reading operations from {}.'.format(opsfile)) + if len(opstack.original_events) != len(labeldata): + print("The number of segments in autosave file is incorrect.") + sys.exit(0) + for stack_event, true_event in zip(opstack.original_events, labeldata): + if (stack_event['name'] != true_event['name'] or + not np.allclose(stack_event['start'], true_event['start']) + or + not np.allclose(stack_event['stop'], true_event['stop'])): + print("Warning! Autosave:\n {}\n Original:\n{}" + .format(stack_event, true_event)) + else: + opstack = OpStack(labeldata) + return opstack + + + +def createlabel(name,start,end,interval): + """. + create a new empty label file with customize gap + + name: the name of label file + start: the start time + endL the end time + interval: the gap between label + + Returns: None + """ + import pandas as pd + data = [] + while start+interval < end: + x=start + y=start+interval + dict = {"start":x,"stop":y,"name":""} + data.append(dict) + start += interval + df = pd.DataFrame(data) + df.to_csv(name,index=False) + +def getfiles(): + file = FileDialog() + files = file.openFileNamesDialog() + if not files: + sys.exit(app.exec_()) + sampled = [bark.read_sampled(file) for file in files] + readonlylabelfile = file.openFileNameDialog() + if not readonlylabelfile: + import pandas as pd + origin_labels = pd.DataFrame() + else: + origin_labels = bark.read_events(readonlylabelfile).data + return files, sampled, sampled + + + +def readfiles(outfile=None, shortcutfile=None, use_ops=True): + """Read all files from the fileDialog and create files if those files are missing. + + If no .dat files, exit. + Auto find label file named with '[dat_name]_split.csv' + If not exist, create a new one with customize label and a .meta file + create opstack and outfiles + + Returns: origin_labels,trace_num, gap, sampled, opstack, shortcuts, outfile, labels.attrs, opsfile + """ + gap = 0 + file = FileDialog() + files = file.openFileNamesDialog() + if not files: + sys.exit(app.exec_()) + files.reverse() + sampled = [bark.read_sampled(file) for file in files] + readonlylabelfile = file.openFileNameDialog() + if not readonlylabelfile: + import pandas as pd + origin_labels = pd.DataFrame() + else: + origin_labels = bark.read_events(readonlylabelfile).data + trace_num = len(files) + dat = files[0] + labelfile = os.path.splitext(dat)[0] + '_split.csv' + exist = os.path.exists(labelfile) + kill_shortcuts(plt) + opsfile = labelfile + '.ops.json' + metadata = labelfile + '.meta.yaml' + if not os.path.exists(labelfile): + write_metadata(labelfile) + + if not os.path.exists(labelfile): + showDia = Input() + gap = int(showDia.showDialog()) + start = 0 + end = int(round(len(sampled[0].data)/sampled[0].attrs["sampling_rate"])) + trace_num = len(sampled) + createlabel(labelfile,start,end,gap) + + + labels = bark.read_events(labelfile) + labeldata = to_seconds(labels).data.to_dict('records') + if len(labeldata) == 0: + print('{} contains no intervals.'.format(labelfile)) + return + opstack = load_opstack(opsfile, labelfile, labeldata, use_ops) + if not gap: + if len(opstack.events) == 0: + print('opstack is empty. Please delete {}.'.format(opstack)) + return + gap = opstack.events[0]['stop'] - opstack.events[0]['start'] + + shortcuts = build_shortcut_map(shortcutfile) + #create a new outfile + if not outfile: + outfile = os.path.splitext(labelfile)[0] + '_edit.csv' + channelname = [] + import re + for name in files: + searchObj = re.search( r'(.*)/(.*).dat', name, re.M|re.I) + channelname.append(searchObj.group(2)) + + + return origin_labels,trace_num, channelname, gap, sampled, opstack, shortcuts, outfile, labels.attrs, opsfile + + +class Plot: + def __init__(self,ax, x_visible=True, y_visible=True, gap=default_gap): + self.ax = ax + self.ax.set_axis_bgcolor('k') + self.boundary_start = self.ax.axvline(color=color) + self.boundary_stop = self.ax.axvline(color=color) + self.gap = gap + self.label = self.ax.text(0,0,'',fontsize=fontsize,color=color) + self.label.set_visible(False) + self.ax.get_xaxis().set_visible(x_visible) + self.ax.get_yaxis().set_visible(y_visible) + + def update_x_axis(self,start,stop): + self.ax.set_xlim(start,stop) + + def clear_plot(self): + self.ax.cla() + + def update_one_label(self,start,stop,name): + # self.boundary_start.setcolor + self.boundary_start.set_xdata((start, start)) + self.boundary_stop.set_xdata((stop, stop)) + ymin, ymax = self.ax.get_ylim() + y = (ymin+ymax)/4 + self.label.set_text(name) + self.label.set_color(color_map(label_to_scatter(name)*7)) + self.label.set_x((start+stop)/2) + self.label.set_y(y) + self.label.set_visible(True) + + +class Psg_Plot(Plot): + + def __init__(self,ax,trace_num,data,sr,channelname, N_points, x_visible=True, y_visible=True, yrange=0.5, gap=default_gap): + + + Plot.__init__(self,ax,x_visible,y_visible,gap) + self.init = False + self.data = data + self.sr = sr + self.yrange = yrange + self.N_points = N_points + self.psg_line = [] + self.trace_num = trace_num + self.display = [] + self.channelname = channelname + self.scale = [] + + for i in range(self.trace_num): + line, = self.ax.plot( + np.arange(self.N_points), + np.zeros(self.N_points), + color='gray') + self.psg_line.append(line) + self.display.append(True) + self.scale.append(1) + + + def update_y_axis(self,yrange): + self.yrange = yrange + + def update_boundary(self,start,stop): + self.boundary_start.set_xdata((start, start)) + self.boundary_stop.set_xdata((stop, stop)) + + def get_yrange(self,buffer_start_samp,buffer_stop_samp): + x = self.data[0][buffer_start_samp:buffer_stop_samp] + return max(x) - min(x) + + def update_psgillograms(self,buffer_start_samp,buffer_stop_samp): + + offset = self.yrange/2 + buf_start = buffer_start_samp / self.sr + buf_stop = buffer_stop_samp / self.sr + + for i in range(self.trace_num): + if self.display[i]: + self.psg_line[i].set_visible(True) + else: + self.psg_line[i].set_visible(False) + + x = self.data[i][buffer_start_samp:buffer_stop_samp] + t = np.arange(len(x)) / self.sr + buf_start + if len(x) > 10000: + t_interp = np.linspace(buf_start, buf_stop, 10000) + x_interp = np.interp(t_interp, t, x) + else: + t_interp = t + x_interp = x + x_interp = list(map(lambda num: int(num*self.scale[i])+offset, x_interp)) + + self.psg_line[i].set_data(t_interp, x_interp) + offset += self.yrange + + self.update_x_axis(buf_start, buf_stop) + self.ax.yaxis.set_ticks(np.arange(0, offset, self.yrange/4)) + self.ax.set_ylim(0,offset) + self.set_y() + self.set_x() + + self.init = True + def set_x(self): + labels = [item.get_text() for item in self.ax.get_xticklabels()] + xtickslocs = self.ax.get_xticks() + + for i in range(len(labels)-1): + hour = int(xtickslocs[i])/3600 + temp = int(xtickslocs[i])% 3600 + minute = temp/60 + sec = temp%60 + labels[i] = "%.d : %.d : %.d" %(hour,minute,sec) + + self.ax.set_xticklabels(labels) + def set_y(self): + labels = [item.get_text() for item in self.ax.get_yticklabels()] + + for i in range(len(labels)-1): + if i%4 == 1: + a = -self.yrange/(4*self.scale[int(i/4)]) + labels[i] = "%.2f" % a + elif i%4 == 2: + labels[i] = self.channelname[int(i/4)] + elif i%4 == 3: + a = self.yrange/(4*self.scale[int(i/4)]) + labels[i] = "%.2f" % a + + else : + labels[i] = "" + + self.ax.set_yticklabels(labels) +''' +class Label_Plot + +parameter: + +ax : axis object to plot spectrogram on +x_visible : if x axis visible +y_visible: if y axis visible +gap: the gap between two label + +''' + +class Label_Plot(Plot): + def __init__(self,ax, x_visible=True, y_visible=True,gap=default_gap): + Plot.__init__(self,ax,x_visible,y_visible,gap) + self.labels = [self.ax.text(0,0,'',fontsize=fontsize,color=color_label) for _ in range(20)] + self.boundaries_start = [self.ax.axvline(color=color_label) for _ in range(20)] + self.boundaries_stop = [self.ax.axvline(color=color_label) for _ in range(20)] + def update_y_axis(self,yrange): + self.yrange = yrange + + def clear_plot(self): + for a,b,c in zip(self.labels,self.boundaries_start,self.boundaries_stop): + a.set_visible(False) + b.set_visible(False) + c.set_visible(False) + + + def update_labels(self,current,data,origin_data=False): + if origin_data == False: + # update the current creating labels + self.update_opstack_labels(current,data) + else: + # update the origin loaded labels + self.update_origin_labels(current,data) + + def update_opstack_labels(self,current,opstack,y_pos = 4): + xmin, xmax = self.ax.get_xlim() + + 'labels for current syl and two on either side' + for i in range(0,len(opstack.events)): + if(opstack.events[i]["start"]>=xmin): + break + start_i = i + while i >= 0 and i < len(opstack.events) and opstack.events[i]["stop"] < xmax: + label_i = i - start_i + text = self.labels[label_i] + start_line = self.boundaries_start[label_i] + start_line.set_visible(True) + stop_line = self.boundaries_stop[label_i] + stop_line.set_visible(True) + + start = opstack.events[i]['start'] + stop = opstack.events[i]['stop'] + x = (start+stop) / 2 + ymin, ymax = self.ax.get_ylim() + y = (ymin+ymax)/y_pos + name = opstack.events[i]['name'] + + if isinstance(name, str): + text.set_x(x) + text.set_visible(True) + text.set_text(name) + text.set_y(y) + text.set_color(color_map(label_to_scatter(name)*7)) + + else: + text.set_visible(False) + + start_line.set_xdata((start, start)) + stop_line.set_xdata((stop, stop)) + + if i == current: + # The start line overlap with the stop line of previous label. + # Fix this by putting start line and stop line on the top layer. + start_line_1 = self.boundaries_stop[label_i-1] + start_line_1.set_color(color) + stop_line.set_color(color) + else: + self.boundaries_stop[label_i].set_color(color_label) + self.boundaries_stop[label_i].set_color(color_label) + + i += 1 + + def update_origin_labels(self,current,origin_data): + self.ax.cla() + xmin, xmax = self.ax.get_xlim() + for i in range(0,len(origin_data)): + if(origin_data["start"][i]>=xmin): + break + while origin_data["stop"][i] < xmax: + name = origin_data["name"][i] + self.ax.axvline(x=origin_data["start"][i],color='pink') + self.ax.axvline(x=origin_data["stop"][i],color='pink') + ymin, ymax = self.ax.get_ylim() + pos_y = (ymin+ymax)/3 + pos_x = origin_data["start"][i] + (origin_data["stop"][i]-origin_data["start"][i])/2 + self.ax.text(pos_x, pos_y, name, fontsize=12, color = color_map(label_to_scatter(name)*7)) + i += 1 + + + +class PlotCanvas(FigureCanvas): + + def __init__(self, + origin_data, + trace_num, + gap, + sampled, + channelname, + opstack, + keymap, + outfile, + out_attrs, + opsfile=None, + parent=None, + width=5, + height=10, + dpi=100): + + fig = Figure(figsize =(width,height),dpi=dpi) + self.maxpoint = 20 + pos_1 = [0.1, 0.9, 0.8, 0.06] + pos_2 = [0.1, 0.82, 0.8, 0.06] + pos_psg = [0.1, 0.2, 0.8, 0.6] + pos_map = [0.1, 0.02, 0.8, 0.1] + + self.axes_1 = fig.add_axes(pos_1) + self.axes_2 = fig.add_axes(pos_2) + self.axes_4 = fig.add_axes(pos_psg) + self.axes = fig.add_axes(pos_map) + + self.origin_data = origin_data + self.trace_num = trace_num + self.data = [] + self.sr = 0 + self.gap = gap + self.yrange = 4000 + for dataset in sampled: + self.data.append(dataset.data.ravel()) + self.sr = dataset.sampling_rate + self.N_points = int(round(self.sr*self.gap*self.maxpoint)) + self.label_attrs = out_attrs + self.opstack = opstack + self.opsfile = opsfile + self.outfile = outfile + self.keymap = keymap + self.y_init = False + + + if opstack.ops: + self.label_index = opstack.ops[-1].index + else: + self.label_index = 0 + + self.psg_ax = Psg_Plot(ax= self.axes_4, + gap=gap, + trace_num=trace_num, + N_points = self.N_points, + data=self.data, + sr = self.sr, + channelname = channelname + ) + if not self.origin_data.empty: + self.label_ax = Label_Plot(ax=self.axes_1,x_visible=False, y_visible=False,gap=gap) + else : + self.label_ax = Label_Plot(ax=self.axes_1,x_visible=False, y_visible=False,gap=gap) + self.label_ax_2 = Label_Plot(ax=self.axes_2,x_visible=False, y_visible=False,gap=gap) + self.map_ax = self.axes + self.initialize_minimap() + + FigureCanvas.__init__(self, fig) + self.setParent(parent) + + FigureCanvas.setSizePolicy(self, + QSizePolicy.Expanding, + QSizePolicy.Expanding) + FigureCanvas.updateGeometry(self) + self.update_plot_data() + + + def update_plot_data(self): + + + if not self.opstack.events: + print('no segments') + plt.close("all") + return + + i = self.label_index + sr = self.sr + start = self.opstack.events[i]['start'] + start_samp = int(start * sr) + stop = self.opstack.events[i]['stop'] + stop_samp = int(stop * sr) + syl_samps = stop_samp - start_samp + buffer_start_samp = start_samp - (self.N_points - syl_samps) // 2 + + if buffer_start_samp < 0: + buffer_start_samp = 0 + + buffer_stop_samp = buffer_start_samp + self.N_points + + if buffer_stop_samp >= self.data[0].shape[0]: + buffer_stop_samp = self.data[0].shape[0] - 1 + + buf_start = buffer_start_samp / sr + buf_stop = buffer_stop_samp / sr + + name = self.opstack.events[i]['name'] + #fix me ax line covered by the grey line + + self.psg_ax.yrange = self.yrange + self.psg_ax.update_psgillograms(buffer_start_samp,buffer_stop_samp) + self.psg_ax.update_boundary(start,stop) + + if self.N_points > int(round(self.sr*self.gap*self.maxpoint)): + if not self.origin_data.empty: + self.label_ax.clear_plot() + self.label_ax_2.clear_plot() + self.label_ax_2.update_x_axis(buf_start,buf_stop) + self.label_ax_2.update_one_label(start,stop,name) + + else: + + if not self.origin_data.empty: + self.label_ax.update_x_axis(buf_start,buf_stop) + self.label_ax.update_labels(current = i,data = self.origin_data,origin_data=True) + + self.label_ax_2.update_x_axis(buf_start,buf_stop) + self.label_ax_2.update_labels(current = i,data = self.opstack,origin_data=False) + + self.update_minimap() + + if self.opstack.ops: + last_command = str(self.opstack.ops[-1]) + else: + last_command = 'none' + if i == 0: + self.map_ax.set_title('ctrl+h for help, prints to terminal') + else: + self.map_ax.set_title('Epoch {}/ {}'.format(i + 1, len( + self.opstack.events))) + + self.draw() + + + + def initialize_minimap(self): + times, values = labels_to_scatter_coords(self.opstack.events) + self.map_ax.set_axis_bgcolor('k') + self.map_ax.scatter(times, + values, + c=values, + vmin=0, + vmax=37, + cmap=plt.get_cmap('hsv'), + edgecolors='none') + self.map_ax.vlines(self.opstack.events[self.label_index]['start'], + -1, + 38, + zorder=0.5, + color='w', + linewidth=1) + self.map_ax.tick_params(axis='y', + which='both', + left='off', + right='off', + labelleft='off') + self.map_ax.set_ylim(-1, 38) + self.map_ax.get_xaxis().set_visible(False) + + + def update_minimap(self): + # If perfomance lags, may need to adjust plot elements instead of + # clearing everything and starting over. + self.map_ax.clear() + self.initialize_minimap() + + + + def connect(self): + 'creates all the event connections' + + self.cid_key_press = self.mpl_connect('key_press_event', + self.on_key_press) + self.cid_mouse_press = self.mpl_connect('button_press_event', + self.on_mouse_press) + + + def on_mouse_press(self, event): + start_pos = self.psg_ax.boundary_start.get_xdata()[0] + stop_pos = self.psg_ax.boundary_stop.get_xdata()[0] + + # jump to syllable from map click + if event.inaxes == self.map_ax: + i = nearest_label(self.opstack.events, float(event.xdata)) + self.label_index = i + self.update_plot_data() + + elif event.inaxes == self.axes_4: + if event.xdata < start_pos: + self.dec_i() + elif event.xdata > stop_pos: + self.inc_i() + + + + def inc_i(self): + 'Go to next syllable.' + if self.label_index < len(self.opstack.events) - 1: + self.label_index += 1 + self.update_plot_data() + def inc_page(self): + 'Go to next syllable.' + if self.label_index +int(self.N_points/int(self.sr * self.gap)) < len(self.opstack.events): + self.label_index += int(self.N_points/int(self.sr * self.gap)) + self.update_plot_data() + + + def dec_i(self): + 'Go to previous syllable.' + if self.label_index > 0: + self.label_index -= 1 + self.update_plot_data() + def dec_page(self): + 'Go to next syllable.' + if self.label_index - int(self.N_points/int(self.sr * self.gap)) > 0: + self.label_index -= int(self.N_points/int(self.sr * self.gap)) + self.update_plot_data() + + def on_key_press(self, event): + + if event.key() == Qt.Key_Right: + self.inc_i() + elif event.key() == Qt.Key_Left: + self.dec_i() + if event.key() == Qt.Key_Down: + self.zoom_in_x() + elif event.key() == Qt.Key_Up: + self.zoom_out_x() + elif event.key() <= Qt.Key_Z and event.key() >= Qt.Key_A: + if self.N_points > int(round(self.sr*self.gap*self.maxpoint)): + return + newlabel = chr(event.key()) + self.opstack.push(Update(self.label_index, 'name', newlabel)) + self.inc_i() + elif event.key() <= Qt.Key_9 and event.key() >= Qt.Key_0: + if self.N_points > int(round(self.sr*self.gap*self.maxpoint)): + return + newlabel = number_shortcut[chr(event.key())] + self.opstack.push(Update(self.label_index, 'name', newlabel)) + self.inc_i() + + def addlabel(self,str): + if self.N_points > int(round(self.sr*self.gap*self.maxpoint)): + return + self.opstack.push(Update(self.label_index, 'name', str)) + self.inc_i() + + def deletelabel(self): + if self.N_points > int(round(self.sr*self.gap*self.maxpoint)): + return + self.opstack.push(Update(self.label_index, 'name', "")) + self.dec_i() + + + def zoom_in_x(self): + if self.N_points < int(self.sr * self.gap)*120: + self.N_points += self.N_points + self.update_plot_data() + + def zoom_out_x(self): + if self.N_points >= int(self.sr * self.gap)*2: + self.N_points -= int(self.N_points/2); + self.update_plot_data() + + def zoom_in_y(self): + self.yrange *= 2 + self.update_plot_data() + + def zoom_out_y(self): + self.yrange /= 2 + self.update_plot_data() + + def delete(self): + self.opstack.push(Delete(self.label_index)) + if self.label_index >= len(self.opstack.events): + self.label_index = len(self.opstack.events) - 1 + self.update_plot_data() + + def redo(self): + if self.opstack.undo_ops: + self.opstack.redo() + self.label_index = self.opstack.ops[-1].index + self.update_plot_data() + + def undo(self): + if self.opstack.ops: + self.opstack.undo() + self.label_index = self.opstack.undo_ops[-1].index + self.update_plot_data() + + def save(self): + 'Writes out labels to file.' + from pandas import DataFrame + label_data = DataFrame(self.opstack.events) + bark.write_events(self.outfile, label_data, **self.label_attrs) + print(self.outfile, 'written') + if self.opsfile: + write_stack(self.opsfile, self.opstack) + print(self.opsfile, 'written') + + def __enter__(self): + return self + + def __exit__(self, *args): + self.save() + + + +class FileDialog(QWidget): + + def __init__(self): + super().__init__() + self.title = 'PyQt5 file dialogs - pythonspot.com' + self.left = 10 + self.top = 10 + self.width = 640 + self.height = 480 + self.initUI() + + def initUI(self): + self.setWindowTitle(self.title) + self.setGeometry(self.left, self.top, self.width, self.height) + self.show() + + def openFileNamesDialog(self): + options = QFileDialog.Options() + options |= QFileDialog.DontUseNativeDialog + files, _ = QFileDialog.getOpenFileNames(self,"Choose all the .dat files (One channel per file)", "","dat Files (*.dat);;Python Files (*.py)", options=options) + if files: + return files + else: + self.close() + + def openFileNameDialog(self): + options = QFileDialog.Options() + options |= QFileDialog.DontUseNativeDialog + fileName, _ = QFileDialog.getOpenFileName(self,"Choose a label file", "","csv Files (*.csv);;Python Files (*.py)", options=options) + if fileName: + return fileName + + +class Input(QWidget): + + def __init__(self): + super().__init__() + + def showDialog(self): + gap, ok = QInputDialog.getText(self, 'Input Dialog', + 'Enter the length of each epoch for labeling (in seconds):') + if ok: + return gap + +class Shortcut_map(QtWidgets.QDialog): + def __init__(self, parent=None): + super(Shortcut_map, self).__init__(parent) + self.d = QtWidgets.QDialog() + b1 = QPushButton("ok",self.d) + b1.move(40,540) + start_x = 20 + start_y = 20 + self.testf = [] + for i in range(0,10): + lbl = QLabel(str(i),self.d) + self.testf.append( QLineEdit(number_shortcut[str(i)],self.d)) + lbl.move(start_x, start_y) + start_y += 20 + self.testf[i].move(start_x, start_y) + start_y += 30 + b1.clicked.connect(self.changesetting) + self.d.setWindowTitle("Dialog") + self.d.setWindowModality(Qt.ApplicationModal) + self.d.exec_() + def changesetting(self): + for i in range(0,10): + number_shortcut[str(i)] = self.testf[i].text() + self.d.close() + + +class ScrollView (QtWidgets.QScrollArea): + def __init__(self, parent = None): + super(ScrollView, self).__init__(parent) + + def keyPressEvent (self, e): + super(ScrollView, self).keyPressEvent(e) + e.ignore() + +class App(QMainWindow): + + def __init__(self): + super().__init__() + + self.left = 0 + self.top = 20 + self.title = 'PsgView' + app = QtWidgets.QApplication(sys.argv) + rect = app.primaryScreen().availableGeometry() + self.width = rect.width() + self.height = rect.height() + self.widget = QWidget() + self.setCentralWidget(self.widget) + self.widget.setLayout(QVBoxLayout()) + self.widget.layout().setContentsMargins(0,0,0,0) + self.widget.layout().setSpacing(0) + self.initUI() + self.initToolbar() + self.initToolbox() + self.show() + + def initToolbar(self): + + self.file_menu = QtWidgets.QMenu('&File', self) + # Quit + self.file_menu.addAction('&Quit', self.fileQuit, + QtCore.Qt.CTRL + QtCore.Qt.Key_Q) + self.menuBar().addMenu(self.file_menu) + self.file_menu.addAction('&Save', self.reviewer.save, QtCore.Qt.CTRL + QtCore.Qt.Key_S) + self.file_menu.addAction('&Help', self.help, QtCore.Qt.CTRL + QtCore.Qt.Key_H) + self.file_menu.addAction('&Redo', self.reviewer.redo, QtCore.Qt.CTRL + QtCore.Qt.Key_Y) + self.file_menu.addAction('&Undo', self.reviewer.undo, QtCore.Qt.CTRL + QtCore.Qt.Key_Z) + self.file_menu.addAction('&shortcut', self.key_map) + + + # control + self.control_menu = QtWidgets.QMenu('&Control', self) + self.menuBar().addSeparator() + self.menuBar().addMenu(self.control_menu) + + self.control_menu.addAction('&zoom in x', self.reviewer.zoom_in_x, QtCore.Qt.CTRL + QtCore.Qt.Key_I) + self.control_menu.addAction('&zoom out x',self.reviewer.zoom_out_x, QtCore.Qt.CTRL + QtCore.Qt.Key_O) + self.control_menu.addAction('&zoom in y', self.reviewer.zoom_in_y, QtCore.Qt.CTRL + QtCore.Qt.Key_W) + self.control_menu.addAction('&zoom out y',self.reviewer.zoom_out_y, QtCore.Qt.CTRL + QtCore.Qt.Key_E) + self.control_menu.addAction('&next',self.reviewer.inc_i, QtCore.Qt.CTRL + QtCore.Qt.Key_J) + self.control_menu.addAction('&previous',self.reviewer.dec_i,QtCore.Qt.CTRL + QtCore.Qt.Key_F) + self.control_menu.addAction('&nextpage',self.reviewer.inc_page,QtCore.Qt.Key_PageDown) + self.control_menu.addAction('&previouspage',self.reviewer.dec_page,QtCore.Qt.Key_PageUp) + self.control_menu.addAction('&delete',self.reviewer.deletelabel,QtCore.Qt.Key_Backspace) + + + def initUI(self): + + self.setWindowTitle(self.title) + self.setGeometry(self.left, self.top, self.width, self.height) + origin_labels,trace_num, channelname, gap, sampled, opstack, shortcuts, outfile, attrs, opsfile = readfiles() + height = trace_num*psg_scale + 3 + self.width = self.width/100 + self.trace_num = trace_num + self.reviewer = PlotCanvas(origin_labels,trace_num, gap, sampled, channelname, + opstack, shortcuts, outfile, attrs, opsfile, parent = self, width=self.width, height=height) + self.reviewer.connect() + self.reviewer.move(20,0) + self.scroll = ScrollView(self.widget) + self.scroll.setWidget(self.reviewer) + self.widget.layout().addWidget(self.scroll) + + self.checkboxes = [] + self.scale_button_in = [] + self.scale_button_out = [] + def initToolbox(self): + trace_num = self.trace_num + lbl = QLabel(self) + lbl.setText("Customize Label") + self.label = QLineEdit(self) + label_button = QPushButton("Add Label",self) + + lbl_2 = QLabel(self) + lbl_2.setText("Number of Labels to display") + self.label_2 = QLineEdit(self) + label_2_button = QPushButton("Ok",self) + label_2_button.setToolTip('Change the number of labels to display') + start_x = (self.width-1) * 100 -50 + lbl.move(start_x, 30) + self.label.move(start_x, 60) + label_button.move(start_x,90) + label_button.clicked.connect(self.add_label) + + lbl_2.move(start_x, 120) + self.label_2.move(start_x, 150) + label_2_button.move(start_x,180) + label_2_button.clicked.connect(self.change_label) + + for i in range(0,trace_num): + box = QCheckBox('Track'+ str(i+1), self) + box.move(start_x,310+25*i) + box.ind = trace_num-1-i + box.stateChanged.connect(self.state_changed) + box.setChecked(True) + self.checkboxes.append(box) + button_1 =QPushButton("+",self) + button_1.setFixedWidth(20) + button_1.setFixedHeight(20) + + button_2 =QPushButton("-",self) + button_2.setFixedWidth(20) + button_2.setFixedHeight(20) + + button_1.move(start_x+80,315+25*i) + button_2.move(start_x+110,315+25*i) + button_1.ind = trace_num-1-i + button_2.ind = trace_num-1-i + button_1.clicked.connect(self.scale_unit_in) + button_2.clicked.connect(self.scale_unit_out) + + self.scale_button_in.append(button_1) + self.scale_button_out.append(button_2) + def keyPressEvent(self,e): + self.reviewer.on_key_press(e) + + def add_label(self): + text = self.label.text() + self.reviewer.addlabel(text) + + def scale_unit_in(self): + target = self.sender() + index = target.ind + self.reviewer.psg_ax.scale[index] += 2 + self.reviewer.update_plot_data() + + def scale_unit_out(self): + target = self.sender() + index = target.ind + if self.reviewer.psg_ax.scale[index] - 2 > 0: + self.reviewer.psg_ax.scale[index] -= 2 + self.reviewer.update_plot_data() + + def change_label(self): + number = self.label_2.text() + if number.isdigit(): + self.reviewer.N_points = int(round(self.reviewer.sr*self.reviewer.gap*int(number))) + self.reviewer.update_plot_data() + + def change_label_max(self): + number = int(self.label_3.text()) + self.reviewer.maxpoint = number + self.reviewer.update_plot_data() + + + def change_ylim(self): + inputNumber = self.ylim.text() + if inputNumber.isdigit(): + self.reviewer.yrange = float(inputNumber) + self.reviewer.update_plot_data() + + def change_gap(self): + inputNumber = self.gap_input.text() + if inputNumber.isdigit(): + self.reviewer.gap = int(inputNumber) + self.reviewer.update_plot_data() + else: + info = "Please select a number, `{0}` isn't valid!" + + def fileQuit(self): + self.close() + + def key_map(self): + w = Shortcut_map() + w.show() + + def about(self): + QtWidgets.QMessageBox.about(self, "About", + """This is bark psg viewer""" + ) + def help(self): + QtWidgets.QMessageBox.about(self, "Help", help_string) + + def state_changed(self, state): + target = self.sender() + index = target.ind + if state == Qt.Checked: + self.reviewer.psg_ax.display[index] = True + self.reviewer.psg_ax.scale[index] += 2 + else: + self.reviewer.psg_ax.display[index] = False + + self.reviewer.update_plot_data() + + +def _main(): + app = QApplication(sys.argv) + ex = App() + sys.exit(app.exec_()) + +if __name__ == '__main__': + _main() + + diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md index f78b863..c0ee8a9 100644 --- a/docs/CONTRIBUTING.md +++ b/docs/CONTRIBUTING.md @@ -5,10 +5,13 @@ Clone the repo locally, create a branch for the feature or bugfix, and start work! ``` -$ git clone https://github.com/margoliashlab/bark dir/to/use +$ git clone https://github.com/margoliashlab/bark +$ cd bark $ git checkout -b awesome-feature ``` +For more help with `git`, see the [git workflow suggestions](git_workflow.md). + Before you submit a pull request, make sure the tests work. (If it makes sense to, write tests for your new feature, too.) ``` @@ -43,10 +46,10 @@ Documentation lives in two places in the repository, depending on what it is. Documentation on the usage of modules, functions, classes, and scripts lives in their docstrings. These should be formatted according to the [Google Python docstring guidelines](http://sphinxcontrib-napoleon.readthedocs.io/en/latest/example_google.html). They are collected via Sphinx's `autodoc` extension and parsed by the `napoleon` extension. -Documentation on the repository as a whole, including installation, example usage and workflows, and other information not directly related to the code itself, lives in the `docs` folder, where it is built by Sphinx. It may be written in either [reStructuredText](http://docutils.sourceforge.net/docs/user/rst/quickref.html) or [(CommonMark-compatible) Markdown](http://commonmark.org/help/) (which is quite similar, but not quite identical, to Github-flavored Markdown. +Documentation on the repository as a whole, including installation, example usage and workflows, and other information not directly related to the code itself, lives in the `docs` folder, where it is built by Sphinx. It may be written in either [reStructuredText](http://docutils.sourceforge.net/docs/user/rst/quickref.html) or [(CommonMark-compatible) Markdown](http://commonmark.org/help/) (which is similar, but not quite identical, to Github-flavored Markdown). -The Sphinx-built documentation can be viewed online at [Read the Docs](bark.readthedocs.io/). +The Sphinx-built documentation can be viewed online at [Read the Docs](http://bark.readthedocs.io/). ## Additional notes -This file is modeled on [Atom's contributing guidelines](https://github.com/atom/atom/blob/master/CONTRIBUTING.md). \ No newline at end of file +This file is modeled on [Atom's contributing guidelines](https://github.com/atom/atom/blob/master/CONTRIBUTING.md). diff --git a/docs/README.md b/docs/README.md index 15f5f81..bca899e 100644 --- a/docs/README.md +++ b/docs/README.md @@ -1,5 +1,8 @@ # Bark -Bark is a standard for electrophysiology data. +Bark is: +1. a standard for time-series data, and a python implementation for reading and writing bark formatted data. +2. A python module for signal processing on larger-than-memory data sets. +3. A set of command-line tools for building data processing pipelines. [![Build Status](https://travis-ci.org/kylerbrown/bark.svg?branch=master)](https://travis-ci.org/kylerbrown/bark) @@ -13,13 +16,13 @@ Version: 0.2 By emphasizing filesystem directories, plain text files and a common binary array format, Bark makes it easy to use both large external projects and simple command-line utilities. -Bark's [small specification](specification.md) and Python implementation are easy to use in custom tools. +Bark's [small specification](../specification.md) and Python implementation are easy to use in custom tools. These tools can be chained together using GNU Make to build data pipelines. ## Why use Bark? -Bark takes the architecture of ARF and replaces HDF5 with common data storage formats, the advantages of this approach are: +Inspired by ARF, Bark uses a hierarchy of common data storage formats. The advantages of this approach are: - Use standard Unix tools to explore your data (cd, ls, grep, find, mv) - Build robust data processing pipelines with shell scripting or @@ -54,25 +57,42 @@ This repository contains: ## Installation -The python interface requires Python 3.5+. Installation with [Conda](http://conda.pydata.org/miniconda.html) is recommended. +The python interface runs under Python 3.5 through 3.8. Installation with [Conda](http://conda.pydata.org/miniconda.html) is recommended. - git clone https://github.com/kylerbrown/bark + git clone https://github.com/margoliashlab/bark cd bark - git clone https://github.com/kylerbrown/resin - cd resin - pip install . - cd .. - pip install -r requirements.txt pip install . # optional tests pytest -v - - -You'll also probably want to install [Neuroscope](http://neurosuite.sourceforge.net/). [Sox](http://sox.sourceforge.net/sox.html) is also useful. + +These installation instructions cover the main bark library and almost all of the conversion +scripts and command-line data manipulation tools. Exceptions are noted below. + +The requirements file omits dependencies for a few optional graphical tools included in this +repository. Their additional requirements are as follows, and are not shared across them. +If you don't intend to use one, you can ignore its requirements. + +* `bark-label-view` (for hand-labeling audio data), requires: + * Matplotlib (>=2.0) + * the spectral analysis library [`resin`](https://github.com/margoliashlab/resin) + * (optional) PyQt5 (provides a slightly better experience, but `bark-label-view` is + perfectly usable without it) +* `bark-psg-view` (for hand-scoring PSG data), requires: + * Matplotlib (2.0.2) + * PyQt5 (5.6.0) +* `bark-scope` opens a sampled data file in [neuroscope](http://neurosuite.sourceforge.net/). + It obviously requires an installation of neuroscope. + Note for MacOS users: you need to link the installed neuroscope to where `bark-scope` + expects to find it: + `$ ln -s /Applications/neuroscope.app/Contents/MacOS/neuroscope /usr/local/bin/neuroscope` + +Finally, [Sox](http://sox.sourceforge.net/sox.html) is also extremely useful for working +with audio data. One conversion routine, `dat-to-audio`, is a wrapper around Sox, and thus +requires it to be installed. ## Shell Commands @@ -103,9 +123,8 @@ There are many external tools for processing CSV files, including [pandas](http: ### Visualizations - `bark-scope` -- opens a sampled data file in [neuroscope](http://neurosuite.sourceforge.net/). (Requires an installation of neuroscope) -Note for MacOS users: run this command in the terminal: -`$ ln -s /Applications/neuroscope.app/Contents/MacOS/neuroscope /usr/local/bin/neuroscope` - `bark-label-view` -- Annotate or review events in relation to a sampled dataset, such as birdsong syllable labels on a microphone recording. +- `bark-psg-view` -- Annotate or review on mutiply channels of .dat files. ### Conversion @@ -113,6 +132,8 @@ Note for MacOS users: run this command in the terminal: - `bark-convert-rhd` -- converts [Intan](http://intantech.com/) .rhd files to datasets in a Bark entry - `bark-convert-openephys` -- converts a folder of [Open-Ephys](http://www.open-ephys.org/) .kwd files to datasets in a Bark entry - `bark-convert-arf` -- converts an ARF file to entries in a Bark Root +- `bark-convert-spyking` -- converts [Spyking Circus](https://spyking-circus.readthedocs.io/en/latest/) spike-sorted event data to a Bark event dataset +- `bark-convert-mountainsort` -- converts [MountainSort](https://github.com/flatironinstitute/mountainlab-js) spike-sorted data to a Bark event dataset - `csv-from-waveclus` -- converts a [wave_clus](https://github.com/csn-le/wave_clus) spike time file to a CSV - `csv-from-textgrid` -- converts a [praat](http://www.fon.hum.uva.nl/praat/) TextGrid file to a CSV - `csv-from-lbl` -- converts an [aplot](https://github.com/melizalab/aplot) [lbl](https://github.com/kylerbrown/lbl) file to a CSV @@ -120,6 +141,7 @@ Note for MacOS users: run this command in the terminal: - `dat-to-wave-clus` -- convert a sampled dataset to a [wave_clus](https://github.com/csn-le/wave_clus) compatible Matlab file - `dat-to-audio` -- convert a sampled dataset to an audio file. Uses [SOX](http://sox.sourceforge.net/) under the hood, and so it can convert to any file type SOX supports. +- `dat-to-mda` -- convert a Bark sampled dataset to a [MountainSort](https://github.com/flatironinstitute/mountainlab-js)-compatible `.mda` file ### Control Flow diff --git a/docs/git_workflow.md b/docs/git_workflow.md new file mode 100644 index 0000000..f3eb9aa --- /dev/null +++ b/docs/git_workflow.md @@ -0,0 +1,102 @@ +# Margoliash lab git and github workflow and guidelines + +## General principles + +1. All development should take place in dedicated branches. +2. The **only** commits made to your `master` branch should be merges from other forks (presumably, mostly from `margoliashlab`). +3. Write and run tests. Some code (particularly I/O) can be hard or impossible to properly test; omitting tests for such code should be a careful and deliberate decision. +4. Follow up on pull requests - you aren't done once you hit "Create pull request". There may be multiple rounds of review and requested or debated changes. +5. All pull requests should be reviewed by at least one person. This really helps cut down on bugs. + +## Instructions for common operations + +All of these instructions are written using the command-line `git` tools. + +### when you begin work on a new feature or bugfix of an existing repository + +We'll pretend the repository is named `testrepo`. + +1. Clone the repo onto your local machine and `cd` into it. + ``` + $ git clone https://github.com/margoliashlab/testrepo.git + $ cd testrepo + ``` +2. Create a new branch for your feature or bugfix. We'll call our new branch "testbranch". + ``` + $ git branch testbranch + $ git checkout testbranch + ``` + Or you can combine both of those operations into a single command: + ``` + $ git checkout -b testbranch + ``` +3. Begin working. + +### update your local copy of the repo + +There may have been changes made to the margoliashlab `master` branch while you were working. To make sure your local copy is up-to-date: + +1. Make sure you're on your local `master` branch. + ``` + $ git checkout master + ``` +2. Merge new commits from the margoliashlab `master` branch into your local `master` branch. + ``` + $ git pull margoliashlab master + ``` + +If you've kept your repo and branch organization tidy, by following this guide, this procedure should be painless and automatic. + +### incorporate changes from your local `master` branch into your new branch + +Once you've merged changes from margoliashlab/master into your local master, you also need to incorporate changes into your new branch. This will involve a rebase, rather than a merge. + +1. Make sure you're on your new branch. + ``` + $ git checkout testbranch + ``` +2. Rebase your branch onto the new local master branch. + ``` + $ git rebase master + ``` + +If `git` determines that this operation can be done automatically, it will do so, and you'll be done. However, if changes you've made in your new branch conflict with changes made to your `master` branch, you'll need to resolve the conflicts manually, by opening the files `git` mentions and deciding which version to keep (the rebase will have inserted some information into the files to help you decide). Once you resolve a round of conflicts, you need to commit your changes, and then continue the rebase with +``` +$ git rebase --continue +``` +This will continue the rebase already in progress. It's possible that the rebase will encounter further conflicts. Resolve them, commit the changes, and continue the rebase as already described until you're done. + +### when you're done writing your new branch + +Let's say you've finished writing a new feature, or you've successfully fixed a bug. Now you want to merge your changes into the margoliashlab `master` branch, so other lab members can use it. + +1. Make sure the tests run. + ``` + $ pytest -v + ``` + The Travis CI service will run the tests when you submit a pull request, but making sure before you submit is a good idea. +2. If there have been changes to the margoliashlab `master` branch while you were working, you may want to incorporate them into your branch before you submit a pull request. However, if none of the files you modified were changed on the margoliashlab `master` branch, you probably won't need to do this - `git` is smart enough to figure that out, based on commit timestamps. If you do need to update your local copy of the repo, follow the instructions elsewhere in this guide to merge them into your `master` branch and rebase your new branch on top of them. +3. Update your repo on Github. (This command assumes your Github copy of this repo is configured to be a remote called "origin".) + ``` + $ git push origin testbranch + ``` + Note that you're pushing to `testbranch` on Github, **not** `master`! If `testbranch` doesn't already exist on Github, `git` will create it. +4. Submit a pull request via the Github website. Make sure you're choosing the correct base and head forks and branches. It's also a very good idea to examine the commits and files that Github thinks should be included in the pull request **before** you hit the "Create pull request" button. This will often clue you in to errors you may be about to make - comparing across the wrong fork or branch, or forgetting to update your `.gitignore` to exclude some file that you needed during development but doesn't belong in the repo itself. +5. Whatever you do, **do not approve your own pull request**. Get somebody else to do it, after they look over the changes you made. This is not something that only newbies have to do - getting a second pair of eyes to look over your work is useful for developers of every level of experience. If there's someone you think is particularly appropriate to review your pull request, you can request that they do so when you submit it. +6. Address any questions, concerns, or requested changes to your pull request that the reviewer brings up. Unless you do this, your work won't be incorporated into the margoliashlab repo, and other lab members won't be able to use it. + +### post-merge cleanup + +Once your pull request has been approved and merged (likely via a squash-merge) into the margoliashlab repo, there are a couple of little bits of tidying-up to do: + +1. If you're done working on your new feature or bugfix, you can delete the branch from your local repository. + ``` + $ git checkout master + $ git branch -D testbranch + ``` + Don't worry about losing work - it's already a part of the margoliashlab `master` branch. And definitely **do not** merge your branch back into your local `master`. + You can also delete the branch on your Github repo. You can either do this on the Github website, or on the command line: + ``` + $ git push origin --delete testbranch + ``` +2. Merge the changes to margoliashlab `master` into your local `master`, following the instructions elsewhere in this guide. This now brings your new code into your local `master` branch. Doing it this way, rather than merging your new branch into your local `master` directly, will save you a lot of headaches later. Trust me on this. diff --git a/requirements.txt b/requirements.txt index 3968082..e963261 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ +numpy<1.17 +scipy +pandas<1.3 pytest -pandas -numpy PyYAML -scipy arrow==0.10.0 dataset==0.8.0 diff --git a/setup.py b/setup.py index a997797..bf98c10 100644 --- a/setup.py +++ b/setup.py @@ -27,7 +27,8 @@ 'bark-convert-rhd=bark.io.rhd.rhd2bark:bark_rhd_to_entry', 'bark-convert-openephys=bark.io.openephys.kwik2dat:kwd_to_entry', 'bark-convert-arf=bark.io.arf2bark:_main', - 'bark-convert-spyking=bark.io.spykingcircus:_main', + 'bark-convert-spyking=bark.io.spykingcircus:main', + 'bark-convert-mountainsort=bark.io.mda:bark_from_mountainsort', 'bark-db=bark.io.db:_run', 'dat-decimate=bark.tools.barkutils:rb_decimate', 'dat-resample=bark.tools.barkutils:rb_resample', @@ -45,8 +46,11 @@ 'dat-split=bark.tools.barkutils:_datchunk', 'dat-to-audio=bark.tools.barkutils:rb_to_audio', 'dat-to-wave-clus=bark.tools.barkutils:rb_to_wave_clus', + 'dat-to-mda=bark.io.mda:mda_from_bark_sampled', 'bark-label-view=bark.tools.labelview:_run', 'bark-for-each=bark.tools.barkforeach:_main', 'bark-rasters=bark.tools.genrasters:_main', + 'bark-psg-view=bark.tools.psg_view:_main', + 'bark-bplot=bark.tools.B_PLot:_run', ] }) diff --git a/tests/test_bark.py b/tests/test_bark.py index 88b795e..2b842e8 100644 --- a/tests/test_bark.py +++ b/tests/test_bark.py @@ -121,3 +121,83 @@ def test_datatypes(): assert bark.DATATYPES.code_to_name[1] == 'ACOUSTIC' assert bark.DATATYPES.code_to_name[2002] == 'COMPONENTL' assert bark.DATATYPES.code_to_name[None] is None + +def test_closing(tmpdir): + # setup + ds_name = 'test_sampled.dat' + entry1_path = os.path.join(tmpdir.strpath, "entry1") + dtime = arrow.get("2020-01-02T03:04:05+06:00").datetime + entry1 = bark.create_entry(entry1_path, dtime, food="pizza") + entry2_path = os.path.join(tmpdir.strpath, "entry2") + dtime = arrow.get("2020-01-10T03:04:05+06:00").datetime + entry2 = bark.create_entry(entry2_path, dtime, food="burritos") + data = np.zeros((10,3), dtype='int16') + params = dict(sampling_rate=30000, units="mV", unit_scale=0.025, extra="barley") + dset_path = os.path.join(entry1_path, ds_name) + dset = bark.write_sampled(datfile=dset_path, data=data, **params) + del entry1, entry2, dset + r = bark.read_root(tmpdir.strpath) + # initial checking + assert len(r.entries) == 2 + for ename in r.entries: + assert callable(r.entries.get(ename)) + # load entry1 + entry1 = r.entries['entry1'] + assert isinstance(r.entries.get('entry1'), bark.Entry) + assert callable(r.entries.get('entry2')) + # load sampled dataset + assert callable(entry1.datasets.get(ds_name)) + ds1 = entry1.datasets[ds_name] + assert not callable(entry1.datasets.get(ds_name)) + assert isinstance(ds1, bark.SampledData) + # close entry + del ds1 + assert not callable(entry1.datasets.get(ds_name)) + assert isinstance(entry1.datasets.get(ds_name), bark.SampledData) + entry1.close() + assert callable(entry1.datasets.get(ds_name)) + # close root + del entry1 + assert not callable(r.entries.get('entry1')) + assert isinstance(r.entries.get('entry1'), bark.Entry) + r.close() + assert callable(r.entries.get('entry1')) + +def test_datchunk(tmpdir): + from bark.tools import barkutils + CHUNK = 350 + TOTAL_SIZE = 1000 + data = np.array([range(TOTAL_SIZE), + range(TOTAL_SIZE, 2 * TOTAL_SIZE)]).transpose() + params = dict(sampling_rate=30000, units="mV", unit_scale=0.025, + extra="barley") + dset = bark.write_sampled(os.path.join(tmpdir.strpath, "test.dat"), data=data, **params) + barkutils.datchunk(dset.path, CHUNK, use_seconds=False, one_cut=True) + first_fn = os.path.join(tmpdir.strpath, "test-chunk-0.dat") + second_fn = os.path.join(tmpdir.strpath, "test-chunk-1.dat") + assert os.path.exists(first_fn) + assert os.path.exists(second_fn) + first = bark.read_sampled(first_fn) + second = bark.read_sampled(second_fn) + assert (first.data == dset.data[:CHUNK,:]).all() + assert first.attrs.pop('offset') == 0 + assert first.attrs == dset.attrs + assert (second.data == dset.data[CHUNK:TOTAL_SIZE,:]).all() + assert second.attrs.pop('offset') == CHUNK + assert second.attrs == dset.attrs + del first, second + os.remove(first_fn) + os.remove(second_fn) + assert not os.path.exists(first_fn) + assert not os.path.exists(second_fn) + barkutils.datchunk(dset.path, CHUNK, use_seconds=False, one_cut=False) + third_fn = os.path.join(tmpdir.strpath, "test-chunk-2.dat") + assert os.path.exists(first_fn) + assert os.path.exists(second_fn) + assert os.path.exists(third_fn) + first = bark.read_sampled(first_fn) + second = bark.read_sampled(second_fn) + third = bark.read_sampled(third_fn) + assert (first.data == dset.data[:CHUNK,:]).all() + assert (second.data == dset.data[CHUNK:2*CHUNK,:]).all() + assert (third.data == dset.data[2*CHUNK:,:]).all() diff --git a/tests/test_datfromwav.py b/tests/test_datfromwav.py new file mode 100644 index 0000000..1348d9a --- /dev/null +++ b/tests/test_datfromwav.py @@ -0,0 +1,53 @@ +import bark, pytest, numpy as np, os.path +from scipy.io import wavfile +from bark.io import datfromwav as dfw + + + +'''Tests for testing successful: +1. creation of wav files +2. creation of .dat and metadata file(.dat.metadata.yaml) + + Assuring: +1. sampling rates of wav and dat file are same +2. Data (numpy array) and dtype are same. +3. Attributes(if any) in .dat files are same as provided at creation of .dat file.''' + +def test_dat_from_wav_int(tmpdir): + + # dtype: int array + dir_path = tmpdir.strpath + fname_wav = os.path.join(dir_path, 'test.wav') + arr_int = np.random.randint(1000, size=(1,100)) + fname_dat = os.path.join(dir_path, 'test.dat') + _helper_test(fname_wav, fname_dat, dir_path, arr_int) + + +def test_dat_from_wav_float(tmpdir): + + # dtype: float array + dir_path = tmpdir.strpath + fname_wav = os.path.join(dir_path, 'test.wav') + arr_float = np.random.uniform(low=0.0, high=1000, size=(1,100)) + fname_dat = os.path.join(dir_path, 'test.dat') + _helper_test(fname_wav, fname_dat, dir_path, arr_float) + + +def _helper_test(fname_wav, fname_dat, dir_path, data): + rate = 48000 + #1. Create wav file + wavfile.write(fname_wav, rate, data) + assert os.path.exists(fname_wav), 'File Not Found: test.wav' + + #2. Generate .dat and .dat.meta.yaml file + attrs = {"name": "hello bark", "project": "bark"} + dat_file = dfw.dat_from_wav(fname_wav, fname_dat, **attrs) + assert os.path.exists(fname_dat), 'File Not Found: test.dat' + assert os.path.exists(os.path.join(dir_path, 'test.dat.meta.yaml')), 'File Not Found: test.dat.meta.yaml' + + #3. Compare data, dtype, rate in .dat file + assert np.array_equal(data, bark.read_sampled(fname_dat).data), 'Data in .wav and .dat files does not match' + assert data.dtype==bark.read_sampled(fname_dat).data.dtype, 'dtypes does not match' + assert rate == bark.read_sampled(fname_dat).sampling_rate, 'Sampling rates does not match' + assert 'hello bark' == bark.read_sampled(fname_dat).attrs['name'], 'name attribute does not match' + assert 'bark' == bark.read_sampled(fname_dat).attrs['project'], 'project attribute does not match' diff --git a/tests/test_stream.py b/tests/test_stream.py index 9ee78d7..c7b5bfd 100644 --- a/tests/test_stream.py +++ b/tests/test_stream.py @@ -8,10 +8,10 @@ eq = np.allclose -data1 = np.arange(30).reshape(10, 3) -data2 = np.arange(500).reshape(100, 5) -data3 = np.arange(1000).reshape(500, 2) -data4 = np.arange(11111).reshape(-1, 1) +data1 = np.arange(1, 31).reshape(10, 3) +data2 = np.arange(1, 501).reshape(100, 5) +data3 = np.arange(1, 1001).reshape(500, 2) +data4 = np.arange(1, 11112).reshape(-1, 1) def dummyf(x): return x