From 3088b4a60baa0ee6121b2f0366cb708173ac6194 Mon Sep 17 00:00:00 2001 From: Kyler Brown Date: Thu, 17 Aug 2017 13:57:23 -0500 Subject: [PATCH 01/37] Fix arf converter (#16) * fixing arf converter * Warn user of FileExistsError --- bark/io/arf2bark.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/bark/io/arf2bark.py b/bark/io/arf2bark.py index c7203c1..39ea85d 100644 --- a/bark/io/arf2bark.py +++ b/bark/io/arf2bark.py @@ -10,7 +10,8 @@ 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', @@ -20,7 +21,7 @@ def _parse_args(raw_args): help='timezone for data, tz database format (default is "America/Chicago")', default=None) 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 +31,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): 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: @@ -112,7 +110,7 @@ def transfer_dset(ds_name, ds, e_path, verbose=False): 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) if __name__ == '__main__': _main() From ebe8a42d51a69a2822320ea712ca6215ebc85153 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Thu, 17 Aug 2017 19:50:33 -0500 Subject: [PATCH 02/37] Fix bark-for-each (#17) * Allow entries to be relative paths * Glob command before running it --- bark/tools/barkforeach.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/bark/tools/barkforeach.py b/bark/tools/barkforeach.py index c4a0e92..ec3fc1f 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,29 @@ 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() for ename in entry_list: + os.chdir(base_dir) os.chdir(ename) if verbose: print('Working on ' + ename) - subprocess.run(cmd.split()) + expanded_cmd = glob_command(cmd) + subprocess.run(expanded_cmd) + +def glob_command(cmd): + expanded_cmd = [] + for token in cmd.split(): + if not(token[0] == '"' and token[-1] == '"'): + g = glob.glob(token) + if g: + expanded_cmd.extend(g) + else: + expanded_cmd.append(token) + else: + expanded_cmd.append(token) + return expanded_cmd def _main(): parsed_args = _parse_args(sys.argv[1:]) From f0201fd9b2e625d3935236a993de17d35f69943f Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Fri, 18 Aug 2017 21:45:22 -0500 Subject: [PATCH 03/37] Minor refactor (#21) * UNITS no longer needed * Make columns constructors always return dicts * Move pandas import to top * Make read_events mirror read_sampled * Make default meta a constant * Comment backwards-compatibility code * Regularize path treatment * Add bark-convert-spyking to to README * Rearrange installation instructions * Fix doc typos --- bark/bark.py | 66 +++++++++++++++++++++++--------------------- docs/CONTRIBUTING.md | 9 +++--- docs/README.md | 9 +++--- 3 files changed, 45 insertions(+), 39 deletions(-) diff --git a/bark/bark.py b/bark/bark.py index b135869..18336e7 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), @@ -168,18 +167,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 +194,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 +204,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 +232,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 +241,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 +278,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 +296,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 +319,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 +351,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 +362,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 +426,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: diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md index f78b863..5bd4db2 100644 --- a/docs/CONTRIBUTING.md +++ b/docs/CONTRIBUTING.md @@ -5,7 +5,8 @@ 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 ``` @@ -43,10 +44,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..d8afc6f 100644 --- a/docs/README.md +++ b/docs/README.md @@ -56,14 +56,14 @@ This repository contains: The python interface requires Python 3.5+. Installation with [Conda](http://conda.pydata.org/miniconda.html) is recommended. - git clone https://github.com/kylerbrown/bark - cd bark - git clone https://github.com/kylerbrown/resin cd resin pip install . cd .. - + + git clone https://github.com/margoliashlab/bark + cd bark + pip install -r requirements.txt pip install . @@ -113,6 +113,7 @@ 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 - `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 From 13ef4e8178b190d6593d5fd3909484ccdcea4df5 Mon Sep 17 00:00:00 2001 From: Kyler Brown Date: Thu, 24 Aug 2017 00:58:27 -0500 Subject: [PATCH 04/37] datspike bug fixes (#23) --- bark/tools/datspike.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) 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) From 751524ab554d8f8423a7553f1cdb59432d6eb5c3 Mon Sep 17 00:00:00 2001 From: Kyler Brown Date: Tue, 5 Sep 2017 13:34:34 -0500 Subject: [PATCH 05/37] Fix #24 (#25) --- bark/__init__.py | 1 + 1 file changed, 1 insertion(+) 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 From 3c44dbf083969c47635f03f9800a88faa655f60d Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 6 Sep 2017 15:20:06 -0500 Subject: [PATCH 06/37] Fix empty tle (#26) * Only create top-level entry if it won't be empty * Fix typo * Remove unnecessary function --- bark/io/arf2bark.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/bark/io/arf2bark.py b/bark/io/arf2bark.py index 39ea85d..6f2df4d 100644 --- a/bark/io/arf2bark.py +++ b/bark/io/arf2bark.py @@ -60,10 +60,13 @@ def arf2bark(arf_file, root_path, 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' + @@ -82,6 +85,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) @@ -105,8 +112,7 @@ 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:]) From add8eb6ab42a35dcdedc57ddb12f686c254e9172 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Thu, 7 Sep 2017 16:24:56 -0500 Subject: [PATCH 07/37] Privatize LazyDict (#28) --- bark/bark.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/bark/bark.py b/bark/bark.py index 18336e7..69dc6da 100644 --- a/bark/bark.py +++ b/bark/bark.py @@ -43,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 @@ -67,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] @@ -445,7 +445,7 @@ def read_entry(name, meta=DEFAULT_META): # 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) From 700bf793686689ab84369a1461e8382fca72ec70 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Sun, 10 Sep 2017 16:20:39 -0500 Subject: [PATCH 08/37] bark-for-each new features (#29) * Add progress to verbose output * Improve quoted string handling * Expand user & home directory in command --- bark/tools/barkforeach.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/bark/tools/barkforeach.py b/bark/tools/barkforeach.py index ec3fc1f..644e04d 100644 --- a/bark/tools/barkforeach.py +++ b/bark/tools/barkforeach.py @@ -18,25 +18,33 @@ def _parse_args(raw_args): 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) + 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 = [] - for token in cmd.split(): - if not(token[0] == '"' and token[-1] == '"'): - g = glob.glob(token) - if g: - expanded_cmd.extend(g) - else: - expanded_cmd.append(token) - else: - expanded_cmd.append(token) + 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(): From 9d3a9723ce72b0ae14cd0483ee2a73c9061bb4f6 Mon Sep 17 00:00:00 2001 From: Annali95 Date: Fri, 15 Sep 2017 14:54:00 -0500 Subject: [PATCH 09/37] Update README.md --- docs/README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/README.md b/docs/README.md index d8afc6f..00c15e1 100644 --- a/docs/README.md +++ b/docs/README.md @@ -106,6 +106,7 @@ There are many external tools for processing CSV files, including [pandas](http: 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 From 77f9d48285f9ba1048b124219c205fdea7aae542 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Mon, 2 Oct 2017 12:45:57 -0500 Subject: [PATCH 10/37] Add closers (#32) * Add close methods for Root and Entry * Write closing tests * Rename variables --- bark/bark.py | 17 +++++++++++++++++ tests/test_bark.py | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+) diff --git a/bark/bark.py b/bark/bark.py index 69dc6da..4462085 100644 --- a/bark/bark.py +++ b/bark/bark.py @@ -77,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(): @@ -99,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(): diff --git a/tests/test_bark.py b/tests/test_bark.py index 88b795e..9fdad97 100644 --- a/tests/test_bark.py +++ b/tests/test_bark.py @@ -121,3 +121,44 @@ 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')) From 9d687dcf801ab9418c5c9e384a16aef4967ecd83 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Tue, 10 Oct 2017 15:33:00 -0500 Subject: [PATCH 11/37] Split once (#34) * Add one_cut option to dat-split * Write tests for datchunk * Remove unnecessary import --- bark/tools/barkutils.py | 25 +++++++++++++++++-------- tests/test_bark.py | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+), 8 deletions(-) diff --git a/bark/tools/barkutils.py b/bark/tools/barkutils.py index 5a1ac98..ece0bfd 100644 --- a/bark/tools/barkutils.py +++ b/bark/tools/barkutils.py @@ -298,18 +298,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/tests/test_bark.py b/tests/test_bark.py index 9fdad97..2b842e8 100644 --- a/tests/test_bark.py +++ b/tests/test_bark.py @@ -162,3 +162,42 @@ def test_closing(tmpdir): 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() From 25ddb3277170e4933724daf13d0de3ffa7c4e80e Mon Sep 17 00:00:00 2001 From: Kyler Brown Date: Fri, 20 Oct 2017 14:01:47 -0500 Subject: [PATCH 12/37] better short description, fixing spec link (#35) --- docs/README.md | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/README.md b/docs/README.md index 00c15e1..5c82d10 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 From 37c0534afa013cdb6e9dfbee36370eb96166e306 Mon Sep 17 00:00:00 2001 From: Kyler Brown Date: Wed, 13 Dec 2017 13:05:06 -0600 Subject: [PATCH 13/37] Update rhd2bark.py (#37) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 👍 --- bark/io/rhd/rhd2bark.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bark/io/rhd/rhd2bark.py b/bark/io/rhd/rhd2bark.py index e27daec..ed5a12b 100644 --- a/bark/io/rhd/rhd2bark.py +++ b/bark/io/rhd/rhd2bark.py @@ -18,7 +18,7 @@ def bark_rhd_to_entry(): 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', From c144c19ebe321666cfb34668bc401daecdf5cf18 Mon Sep 17 00:00:00 2001 From: DipanshuSehjal Date: Tue, 19 Dec 2017 14:36:10 -0600 Subject: [PATCH 14/37] Convert wav to dat format (#36) * Convert wav to dat format * Update and rename dat-from-wav.py to datfromwav.py * Create datfromwav.py * Update datfromwav.py * Delete datfromwav.py * test for datfromwav * pytest for datfromwav and update to datfromwav.py * update changes to test_dat_from_wav --- bark/io/datfromwav.py | 40 ++++++++++++++++++++++++++++++ tests/test_datfromwav.py | 53 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 bark/io/datfromwav.py create mode 100644 tests/test_datfromwav.py 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/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' From ffcb8b686a980033e2e10edadd3ac92388eb5b40 Mon Sep 17 00:00:00 2001 From: Annali95 Date: Tue, 27 Mar 2018 17:24:20 -0400 Subject: [PATCH 15/37] Psg-view added (#20) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Psg-view added * fix some problems * Error fixed * modularize the bark label view Doesn’t change the functionality of bark. Only modularize the class including: Osc_Plot Spec_Plot Minimap_Plot * Setup added * Remove the model_label_view * Replace code with bark.write_metadata * fix some bugs * add some functions * add pyqt5 to requirement.txt * remove print * Enable arrows to control the keys * Input dialog added * update * add next page * Add color to psg_view * Update README.md * Update README.md --- bark/tools/psg_view.py | 1105 ++++++++++++++++++++++++++++++++++++++++ docs/README.md | 2 + requirements.txt | 1 + setup.py | 2 + 4 files changed, 1110 insertions(+) create mode 100644 bark/tools/psg_view.py 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/README.md b/docs/README.md index 5c82d10..50279c7 100644 --- a/docs/README.md +++ b/docs/README.md @@ -59,6 +59,8 @@ This repository contains: The python interface requires Python 3.5+. Installation with [Conda](http://conda.pydata.org/miniconda.html) is recommended. +If any error happens, please make sure your matplotlib version is 2.0.2 and the pyqt version is 5.6.0. + git clone https://github.com/kylerbrown/resin cd resin pip install . diff --git a/requirements.txt b/requirements.txt index 3968082..a7d4ab8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,3 +5,4 @@ PyYAML scipy arrow==0.10.0 dataset==0.8.0 +pyqt5 diff --git a/setup.py b/setup.py index a997797..c4fca58 100644 --- a/setup.py +++ b/setup.py @@ -48,5 +48,7 @@ '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' + ] }) From b7fcbd770c481bcb74a39422e2704b94a9fbc705 Mon Sep 17 00:00:00 2001 From: Kyler Brown Date: Wed, 2 May 2018 13:11:11 -0500 Subject: [PATCH 16/37] Forced QT5 backend for OS X (#40) * Forced QT5 backend for OS X * Remove debugging code * Add comment to OS X-specific backend --- bark/tools/labelview.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/bark/tools/labelview.py b/bark/tools/labelview.py index 4685ee0..38fb3a0 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,14 @@ 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. Qt5 does capture them correctly. + import matplotlib + matplotlib.use('Qt5Agg') + +import matplotlib.pyplot as plt + help_string = ''' Pressing any number or letter (uppercase or lowercase) will mark a segment. From a5e6d0e5797a62836b0cd0b1021c014662508e4f Mon Sep 17 00:00:00 2001 From: Annali95 Date: Wed, 2 May 2018 13:40:22 -0500 Subject: [PATCH 17/37] B plot (#38) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Psg-view added * fix some problems * Error fixed * modularize the bark label view Doesn’t change the functionality of bark. Only modularize the class including: Osc_Plot Spec_Plot Minimap_Plot * Setup added * Remove the model_label_view * Replace code with bark.write_metadata * fix some bugs * add some functions * b plot added * add set up bark-Bplot * delete psg_view * Fix function-to-script addressing --- bark/tools/B_PLot.py | 432 +++++++++++++++++++++++++++++++++++++++++++ setup.py | 4 +- 2 files changed, 434 insertions(+), 2 deletions(-) create mode 100644 bark/tools/B_PLot.py diff --git a/bark/tools/B_PLot.py b/bark/tools/B_PLot.py new file mode 100644 index 0000000..221377d --- /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.stop - self.window_size + self.stop = self.start + 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/setup.py b/setup.py index c4fca58..a06e9af 100644 --- a/setup.py +++ b/setup.py @@ -48,7 +48,7 @@ '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-psg-view=bark.tools.psg_view:_main', + 'bark-bplot=bark.tools.B_PLot:_run', ] }) From f602e6af09dccc22a96a9753f7776fb40bf43a73 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 2 May 2018 14:10:15 -0500 Subject: [PATCH 18/37] Fix out-of-bounds bug when viewing end of file (#42) --- bark/tools/B_PLot.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bark/tools/B_PLot.py b/bark/tools/B_PLot.py index 221377d..5988c46 100644 --- a/bark/tools/B_PLot.py +++ b/bark/tools/B_PLot.py @@ -345,8 +345,8 @@ 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.stop - self.window_size - self.stop = self.start + self.start = self.max_time - self.window_size + self.stop = self.max_time self.update_plot_data() def dec_i(self): From b53108fff5221177d0523d7ba28ef5da277ea009 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 2 May 2018 15:56:18 -0500 Subject: [PATCH 19/37] Refactor RHD I/O (#41) * Convert EOLs * Set up magic number removal * Alias long names * Combine numpy fromfile calls into one * Update attribution --- bark/io/rhd/read_one_data_block.py | 103 +++++++++++++++++------------ 1 file changed, 59 insertions(+), 44 deletions(-) diff --git a/bark/io/rhd/read_one_data_block.py b/bark/io/rhd/read_one_data_block.py index cd5f887..19141c7 100644 --- a/bark/io/rhd/read_one_data_block.py +++ b/bark/io/rhd/read_one_data_block.py @@ -1,44 +1,59 @@ -#! /bin/env python -# -# Michael Gibson 23 April 2015 - -import sys, struct -import numpy as np - -def read_one_data_block(data, header, indices, fid): - """Reads one 60-sample data block from fid into data, at the location indicated by indices.""" - - # In version 1.2, we moved from saving timestamps as unsigned - # integers to signed integers to accommodate negative (adjusted) - # timestamps for pretrigger data[' - if (header['version']['major'] == 1 and header['version']['minor'] >= 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))) - +#! /bin/env python +# +# Michael Gibson 23 April 2015 +# Graham Fetterman April 2018 + +import sys, struct +import numpy as np + +AMP_SAMPLES = 60 +AUX_SAMPLES = 15 +SUPPLY_SAMPLES = 1 +TEMP_SAMPLES = 1 +ADC_SAMPLES = 60 + +def read_one_data_block(data, header, indices, fid): + """Reads one 60-sample data block from fid into data, at the location indicated by indices.""" + + # In version 1.2, we moved from saving timestamps as unsigned + # integers to signed integers to accommodate negative (adjusted) + # timestamps for pretrigger data[' + if (header['version']['major'] == 1 and header['version']['minor'] >= 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))) + + num_amp_chan = header['num_amplifier_channels'] + num_aux_chan = header['num_aux_input_channels'] + num_supply_chan = header['num_supply_voltage_channels'] + num_temp_chan = header['num_temp_sensor_channels'] + num_adc_chan = header['num_board_adc_channels'] + num_samples = (num_amp_chan * AMP_SAMPLES + + num_aux_chan * AUX_SAMPLES + + num_supply_chan * SUPPLY_SAMPLES + + num_temp_chan * TEMP_SAMPLES + + num_adc_chan * ADC_SAMPLES) + + tmp = np.fromfile(fid, dtype='uint16', count=num_samples) + start = 0 + if num_amp_chan: + data['amplifier_data'][range(num_amp_chan), indices['amplifier']:(indices['amplifier'] + AMP_SAMPLES)] = tmp[start:(num_amp_chan * AMP_SAMPLES)].reshape(num_amp_chan, AMP_SAMPLES) + start += num_amp_chan * AMP_SAMPLES + if num_aux_chan: + data['aux_input_data'][range(num_aux_chan), indices['aux_input']:(indices['aux_input'] + AUX_SAMPLES)] = tmp[start:(start + num_aux_chan * AUX_SAMPLES)].reshape(num_aux_chan, AUX_SAMPLES) + start += num_aux_chan * AUX_SAMPLES + if num_supply_chan: + data['supply_voltage_data'][range(num_supply_chan), indices['supply_voltage']:(indices['supply_voltage'] + SUPPLY_SAMPLES)] = tmp[start:(start + num_supply_chan * SUPPLY_SAMPLES)].reshape(num_supply_chan, SUPPLY_SAMPLES) + start += num_supply_chan * SUPPLY_SAMPLES + if num_temp_chan: + data['temp_sensor_data'][range(num_temp_chan), indices['supply_voltage']:(indices['supply_voltage'] + TEMP_SAMPLES)] = tmp[start:(start + num_temp_chan * TEMP_SAMPLES)].reshape(num_temp_chan, TEMP_SAMPLES) + start += num_temp_chan * TEMP_SAMPLES + if num_adc_chan: + data['board_adc_data'][range(num_adc_chan), indices['board_adc']:(indices['board_adc'] + ADC_SAMPLES)] = tmp[start:(start + num_adc_chan * ADC_SAMPLES)].reshape(num_adc_chan, ADC_SAMPLES) + + 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))) + From b84c5298488029dda9c1b1823022fecd25d639ba Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Thu, 10 May 2018 18:46:32 -0500 Subject: [PATCH 20/37] Fix Spyking Circus bark conversion (#45) * Replace spykingcircus matlab reader with phy reader * Improve argparse help text * Update setup.py --- bark/io/spykingcircus.py | 149 +++++++++++++-------------------------- setup.py | 2 +- 2 files changed, 50 insertions(+), 101 deletions(-) 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/setup.py b/setup.py index a06e9af..daf2849 100644 --- a/setup.py +++ b/setup.py @@ -27,7 +27,7 @@ '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-db=bark.io.db:_run', 'dat-decimate=bark.tools.barkutils:rb_decimate', 'dat-resample=bark.tools.barkutils:rb_resample', From f6a964dab40a668628b10b29d2b59e784a293f88 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Thu, 10 May 2018 18:46:54 -0500 Subject: [PATCH 21/37] dat-ref speedup (#46) * Copy dataset with shutil * Get # channels and samples from data.shape * Reduce number of calls to np.mean() * Reduce number of calls to np.median() * Fix variable name bugs * Make output file argument required --- bark/tools/datref.py | 59 ++++++++++++++++++++++++++++---------------- 1 file changed, 38 insertions(+), 21 deletions(-) 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) From 896d396411e94871bb6fa18fc3f9796e01d3a75a Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Thu, 17 May 2018 15:37:34 -0500 Subject: [PATCH 22/37] Further RHD I/O speedup (#43) * Add more constants * Change function signature and docstring * Fold time, digital IO into read; read >=1 dbs * Rename file to reflect function * Update import * Grab dbs in chunks * Update init import * Add legacy option * Remove debugging --- bark/io/rhd/__init__.py | 2 +- bark/io/rhd/legacy_load_intan_rhd_format.py | 262 ++++++++++ bark/io/rhd/legacy_read_one_data_block.py | 44 ++ bark/io/rhd/load_intan_rhd_format.py | 530 ++++++++++---------- bark/io/rhd/read_data_blocks.py | 73 +++ bark/io/rhd/read_one_data_block.py | 59 --- bark/io/rhd/rhd2bark.py | 15 +- 7 files changed, 661 insertions(+), 324 deletions(-) create mode 100644 bark/io/rhd/legacy_load_intan_rhd_format.py create mode 100644 bark/io/rhd/legacy_read_one_data_block.py create mode 100644 bark/io/rhd/read_data_blocks.py delete mode 100644 bark/io/rhd/read_one_data_block.py diff --git a/bark/io/rhd/__init__.py b/bark/io/rhd/__init__.py index cc2b333..806f8d7 100644 --- a/bark/io/rhd/__init__.py +++ b/bark/io/rhd/__init__.py @@ -1,5 +1,5 @@ from .read_header import read_header from .get_bytes_per_data_block import get_bytes_per_data_block -from .read_one_data_block import read_one_data_block +from .read_data_blocks import read_data_blocks from .notch_filter import notch_filter from .data_to_result import data_to_result 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..656c8d3 --- /dev/null +++ b/bark/io/rhd/legacy_load_intan_rhd_format.py @@ -0,0 +1,262 @@ +#! /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.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 + +# 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 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..cd5f887 --- /dev/null +++ b/bark/io/rhd/legacy_read_one_data_block.py @@ -0,0 +1,44 @@ +#! /bin/env python +# +# Michael Gibson 23 April 2015 + +import sys, struct +import numpy as np + +def read_one_data_block(data, header, indices, fid): + """Reads one 60-sample data block from fid into data, at the location indicated by indices.""" + + # In version 1.2, we moved from saving timestamps as unsigned + # integers to signed integers to accommodate negative (adjusted) + # timestamps for pretrigger data[' + if (header['version']['major'] == 1 and header['version']['minor'] >= 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/load_intan_rhd_format.py b/bark/io/rhd/load_intan_rhd_format.py index 039bacf..b26c519 100644 --- a/bark/io/rhd/load_intan_rhd_format.py +++ b/bark/io/rhd/load_intan_rhd_format.py @@ -1,262 +1,268 @@ -#! /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 + +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_data_blocks import read_data_blocks +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 + chunk_size = 1000 # data blocks to read at once; diminishing returns above 1000 + chunks, remainder = divmod(num_data_blocks, chunk_size) + for i in range(chunks): + read_data_blocks(data, header, indices, fid, datablocks_per_chunk=chunk_size) + + # Increment indices + indices['amplifier'] += 60 * chunk_size + indices['aux_input'] += 15 * chunk_size + indices['supply_voltage'] += 1 * chunk_size + indices['board_adc'] += 60 * chunk_size + indices['board_dig_in'] += 60 * chunk_size + indices['board_dig_out'] += 60 * chunk_size + + fraction_done = 100 * (1.0 * i / num_data_blocks) + if fraction_done >= percent_done: + percent_done = percent_done + print_increment + if remainder: + read_data_blocks(data, header, indices, fid, datablocks_per_chunk=remainder) + # 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 diff --git a/bark/io/rhd/read_data_blocks.py b/bark/io/rhd/read_data_blocks.py new file mode 100644 index 0000000..befdda0 --- /dev/null +++ b/bark/io/rhd/read_data_blocks.py @@ -0,0 +1,73 @@ +#! /bin/env python +# +# Michael Gibson 23 April 2015 +# Graham Fetterman April 2018 + +import sys, struct +import numpy as np + +NAMES = ['time', 'amp', 'aux', 'supply', 'temp', 'adc', 'digin', 'digout'] + +TIME_SAMPLES = 60 +AMP_SAMPLES = 60 +AUX_SAMPLES = 15 +SUPPLY_SAMPLES = 1 +TEMP_SAMPLES = 1 +ADC_SAMPLES = 60 +DIGIN_SAMPLES = 60 +DIGOUT_SAMPLES = 60 +ALL_SAMPLES = [TIME_SAMPLES, AMP_SAMPLES, AUX_SAMPLES, SUPPLY_SAMPLES, TEMP_SAMPLES, ADC_SAMPLES, DIGIN_SAMPLES, DIGOUT_SAMPLES] + +TIME_DTYPE_NEW = 'i4' +TIME_DTYPE_OLD = 'u4' +AMP_DTYPE = 'u2' +AUX_DTYPE = 'u2' +SUPPLY_DTYPE = 'u2' +TEMP_DTYPE = 'u2' +ADC_DTYPE = 'u2' +DIGIN_DTYPE = 'u2' +DIGOUT_DTYPE = 'u2' +ALL_DTYPES = [TIME_DTYPE_NEW, AMP_DTYPE, AUX_DTYPE, SUPPLY_DTYPE, TEMP_DTYPE, ADC_DTYPE, DIGIN_DTYPE, DIGOUT_DTYPE] + +def read_data_blocks(data, header, indices, fid, datablocks_per_chunk=1): + """Reads a number of 60-sample data blocks from fid into data, at the location indicated by indices.""" + + # 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'] == 1 and header['version']['minor'] >= 2) or (header['version']['major'] > 1): + ALL_DTYPES[0] = TIME_DTYPE_NEW + else: + ALL_DTYPES[0] = TIME_DTYPE_OLD + + time_chan = 1 + amp_chan = header['num_amplifier_channels'] + aux_chan = header['num_aux_input_channels'] + supply_chan = header['num_supply_voltage_channels'] + temp_chan = header['num_temp_sensor_channels'] + adc_chan = header['num_board_adc_channels'] + digin_chan = header['num_board_dig_in_channels'] + digout_chan = header['num_board_dig_out_channels'] + all_chans = [time_chan, amp_chan, aux_chan, supply_chan, temp_chan, adc_chan, digin_chan, digout_chan] + + # create a structured dtype for one datablock + db_dtype = [(name, (dt, (chans * samples))) for name,dt,chans,samples in zip(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): + data['t_amplifier'][(indices['amplifier'] + idx * TIME_SAMPLES):(indices['amplifier'] + (idx + 1) * TIME_SAMPLES)] = db['time'] + if amp_chan: + data['amplifier_data'][range(amp_chan), (indices['amplifier'] + idx * AMP_SAMPLES):(indices['amplifier'] + (idx + 1) * AMP_SAMPLES)] = db['amp'].reshape(amp_chan, AMP_SAMPLES) + if aux_chan: + data['aux_input_data'][range(aux_chan), (indices['aux_input'] + idx * AUX_SAMPLES):(indices['aux_input'] + (idx + 1) * AUX_SAMPLES)] = db['aux'].reshape(aux_chan, AUX_SAMPLES) + if supply_chan: + data['supply_voltage_data'][range(supply_chan), (indices['supply_voltage'] + idx * SUPPLY_SAMPLES):(indices['supply_voltage'] + (idx + 1) * SUPPLY_SAMPLES)] = db['supply'].reshape(supply_chan, SUPPLY_SAMPLES) + if temp_chan: + data['temp_sensor_data'][range(temp_chan), (indices['supply_voltage'] + idx * TEMP_SAMPLES):(indices['supply_voltage'] + (idx + 1) * TEMP_SAMPLES)] = db['temp'].reshape(temp_chan, TEMP_SAMPLES) + if adc_chan: + data['board_adc_data'][range(adc_chan), (indices['board_adc'] + idx * ADC_SAMPLES):(indices['board_adc'] + (idx + 1) * ADC_SAMPLES)] = db['adc'].reshape(adc_chan, ADC_SAMPLES) + if digin_chan: + data['board_dig_in_raw'][(indices['board_dig_in'] + idx * DIGIN_SAMPLES):(indices['board_dig_in'] + (idx + 1) * DIGIN_SAMPLES)] = db['digin'].reshape(digin_chan, DIGIN_SAMPLES) + if digout_chan: + data['board_dig_out_raw'][(indices['board_dig_out'] + idx * DIGOUT_SAMPLES):(indices['board_dig_out'] + (idx + 1) * DIGOUT_SAMPLES)] = db['digout'].reshape(digout_chan, DIGOUT_SAMPLES) diff --git a/bark/io/rhd/read_one_data_block.py b/bark/io/rhd/read_one_data_block.py deleted file mode 100644 index 19141c7..0000000 --- a/bark/io/rhd/read_one_data_block.py +++ /dev/null @@ -1,59 +0,0 @@ -#! /bin/env python -# -# Michael Gibson 23 April 2015 -# Graham Fetterman April 2018 - -import sys, struct -import numpy as np - -AMP_SAMPLES = 60 -AUX_SAMPLES = 15 -SUPPLY_SAMPLES = 1 -TEMP_SAMPLES = 1 -ADC_SAMPLES = 60 - -def read_one_data_block(data, header, indices, fid): - """Reads one 60-sample data block from fid into data, at the location indicated by indices.""" - - # In version 1.2, we moved from saving timestamps as unsigned - # integers to signed integers to accommodate negative (adjusted) - # timestamps for pretrigger data[' - if (header['version']['major'] == 1 and header['version']['minor'] >= 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))) - - num_amp_chan = header['num_amplifier_channels'] - num_aux_chan = header['num_aux_input_channels'] - num_supply_chan = header['num_supply_voltage_channels'] - num_temp_chan = header['num_temp_sensor_channels'] - num_adc_chan = header['num_board_adc_channels'] - num_samples = (num_amp_chan * AMP_SAMPLES + - num_aux_chan * AUX_SAMPLES + - num_supply_chan * SUPPLY_SAMPLES + - num_temp_chan * TEMP_SAMPLES + - num_adc_chan * ADC_SAMPLES) - - tmp = np.fromfile(fid, dtype='uint16', count=num_samples) - start = 0 - if num_amp_chan: - data['amplifier_data'][range(num_amp_chan), indices['amplifier']:(indices['amplifier'] + AMP_SAMPLES)] = tmp[start:(num_amp_chan * AMP_SAMPLES)].reshape(num_amp_chan, AMP_SAMPLES) - start += num_amp_chan * AMP_SAMPLES - if num_aux_chan: - data['aux_input_data'][range(num_aux_chan), indices['aux_input']:(indices['aux_input'] + AUX_SAMPLES)] = tmp[start:(start + num_aux_chan * AUX_SAMPLES)].reshape(num_aux_chan, AUX_SAMPLES) - start += num_aux_chan * AUX_SAMPLES - if num_supply_chan: - data['supply_voltage_data'][range(num_supply_chan), indices['supply_voltage']:(indices['supply_voltage'] + SUPPLY_SAMPLES)] = tmp[start:(start + num_supply_chan * SUPPLY_SAMPLES)].reshape(num_supply_chan, SUPPLY_SAMPLES) - start += num_supply_chan * SUPPLY_SAMPLES - if num_temp_chan: - data['temp_sensor_data'][range(num_temp_chan), indices['supply_voltage']:(indices['supply_voltage'] + TEMP_SAMPLES)] = tmp[start:(start + num_temp_chan * TEMP_SAMPLES)].reshape(num_temp_chan, TEMP_SAMPLES) - start += num_temp_chan * TEMP_SAMPLES - if num_adc_chan: - data['board_adc_data'][range(num_adc_chan), indices['board_adc']:(indices['board_adc'] + ADC_SAMPLES)] = tmp[start:(start + num_adc_chan * ADC_SAMPLES)].reshape(num_adc_chan, ADC_SAMPLES) - - 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 ed5a12b..e16d09a 100644 --- a/bark/io/rhd/rhd2bark.py +++ b/bark/io/rhd/rhd2bark.py @@ -3,7 +3,8 @@ import arrow 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 +import bark.io.rhd.legacy_load_intan_rhd_format from bark import create_entry, write_metadata @@ -48,11 +49,16 @@ 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") args = p.parse_args() attrs = dict(args.keyvalues) if args.keyvalues else {} check_exists(args.rhdfiles) rhds_to_entry(args.rhdfiles, args.out, args.timestamp, args.parents, - args.maxgaps, args.timestamp, **attrs) + args.maxgaps, args.timestamp, legacy=args.legacy, **attrs) def rhd_filename_to_timestamp(fname, timezone): @@ -139,10 +145,15 @@ def rhds_to_entry(rhd_paths, parents=False, max_gaps=10, timezone='America/Chicago', + legacy=False, **attrs): """ Converts a temporally contiguous list of .rhd files to a bark entry. """ + if legacy: + read_data = bark.io.rhd.legacy_load_intan_rhd_format.read_data + else: + read_data = bark.io.rhd.load_intan_rhd_format.read_data if not timestamp: timestamp = rhd_filename_to_timestamp(rhd_paths[0], timezone) else: From f349a4590dfb961d7b4395d552c36662e0716924 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 25 Jul 2018 14:30:39 -0500 Subject: [PATCH 23/37] Stream RHD I/O (#51) * Convert large RHD files in chunks * Warn user of removal of notch filter if applied * Re-enable legacy functionality * Remove now-unnecessary indices * Refactor data memory preallocation * Remove notch filter application * Reformat refactor function * Move preallocation to read_data_blocks * Refactor rhd2bark * Fix bug * Remove unnecessary imports & tidy * Add and update docstrings * Reformat version check * Account for copies of data in memory restriction * Refactor uint->int conversion to use less memory * Avoid duplicating data in memory * Remove unnecessary parameters --- bark/io/rhd/load_intan_rhd_format.py | 283 +++++++++++---------------- bark/io/rhd/read_data_blocks.py | 74 +++++-- bark/io/rhd/rhd2bark.py | 175 ++++++++++------- 3 files changed, 284 insertions(+), 248 deletions(-) diff --git a/bark/io/rhd/load_intan_rhd_format.py b/bark/io/rhd/load_intan_rhd_format.py index b26c519..19eec15 100644 --- a/bark/io/rhd/load_intan_rhd_format.py +++ b/bark/io/rhd/load_intan_rhd_format.py @@ -2,16 +2,16 @@ # # Michael Gibson 17 July 2015 # Kyler Brown December 2016 +# Graham Fetterman July 2018 from __future__ import absolute_import, division, unicode_literals, print_function -import sys, struct, math, os, time +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 -from bark.io.rhd.notch_filter import notch_filter +from bark.io.rhd.read_data_blocks import read_data_blocks, preallocate_memory, AMP_SAMPLES from bark.io.rhd.data_to_result import data_to_result # constants @@ -23,20 +23,33 @@ ADC_BIT_VOLTS_0 = 50.353e-6 TEMP_BIT_CELCIUS = 0.01 - -def read_data(filename, no_floats=False): +def read_data(filename, no_floats=False, max_memory=0): """Reads Intan Technologies RHD2000 data file generated by evaluation board GUI. - Data are returned in a dictionary, for future extensibility. + 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 + """ - tic = time.time() 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'])) + print('Found {} amplifier channel{}.'.format(header[ 'num_amplifier_channels'], plural(header['num_amplifier_channels']))) print('Found {} auxiliary input channel{}.'.format(header[ @@ -72,14 +85,7 @@ def read_data(filename, no_floats=False): 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'] + record_time = AMP_SAMPLES * num_data_blocks / header['sample_rate'] if data_present: print( @@ -91,170 +97,113 @@ def read_data(filename, no_floats=False): 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 - chunk_size = 1000 # data blocks to read at once; diminishing returns above 1000 + # 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 i in range(chunks): - read_data_blocks(data, header, indices, fid, datablocks_per_chunk=chunk_size) - - # Increment indices - indices['amplifier'] += 60 * chunk_size - indices['aux_input'] += 15 * chunk_size - indices['supply_voltage'] += 1 * chunk_size - indices['board_adc'] += 60 * chunk_size - indices['board_dig_in'] += 60 * chunk_size - indices['board_dig_out'] += 60 * chunk_size - - fraction_done = 100 * (1.0 * i / num_data_blocks) - if fraction_done >= percent_done: - percent_done = percent_done + print_increment + read_data_blocks(data, header, fid, datablocks_per_chunk=chunk_size) + yield check_data_and_reformat(header, data, no_floats) if remainder: - read_data_blocks(data, header, indices, fid, datablocks_per_chunk=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.') - -# Close data file. + else: + yield data_to_result(header, {}, data_present) + # 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) +def check_data_and_reformat(header, data, no_floats): + """Performs some cleanup on data and builds a dictionary to return. -# 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 + 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 - if header['eval_board_mode'] == 1: - extras['ADC_input_bit_volts'] = ADC_BIT_VOLTS_1 + 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']): + 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 terabyte+ level. + extras['amplifier_bit_microvolts'] = 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'] = 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 + extras['ADC_input_bit_volts'] = ADC_BIT_VOLTS_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: - data = [] - -# Move variables to result struct. - result = data_to_result(header, data, data_present) + # 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'] + + # Move variables to result struct. + result = data_to_result(header, data, data_present=True) result.update(extras) return result @@ -262,7 +211,3 @@ def read_data(filename, no_floats=False): def plural(n): return '' if n == 1 else 's' - -if __name__ == '__main__': - a = read_data(sys.argv[1]) - #print a diff --git a/bark/io/rhd/read_data_blocks.py b/bark/io/rhd/read_data_blocks.py index befdda0..e862dbf 100644 --- a/bark/io/rhd/read_data_blocks.py +++ b/bark/io/rhd/read_data_blocks.py @@ -3,7 +3,6 @@ # Michael Gibson 23 April 2015 # Graham Fetterman April 2018 -import sys, struct import numpy as np NAMES = ['time', 'amp', 'aux', 'supply', 'temp', 'adc', 'digin', 'digout'] @@ -29,13 +28,63 @@ DIGOUT_DTYPE = 'u2' ALL_DTYPES = [TIME_DTYPE_NEW, AMP_DTYPE, AUX_DTYPE, SUPPLY_DTYPE, TEMP_DTYPE, ADC_DTYPE, DIGIN_DTYPE, DIGOUT_DTYPE] -def read_data_blocks(data, header, indices, fid, datablocks_per_chunk=1): - """Reads a number of 60-sample data blocks from fid into data, at the location indicated by indices.""" +def preallocate_memory(header, num_datablocks): + """Preallocates space for a chunk of data. + Args: + header (dict): metadata + num_datablocks (int): size of chunk in datablocks + + Returns: + dict of numpy arrays: array sizes and types depend on Intan specs + """ + data = {} + if ((header['version']['major'], header['version']['minor']) >= (1, 2)): + time_dt = np.int + else: + time_dt = np.uint + data['t_amplifier'] = np.zeros(AMP_SAMPLES * num_datablocks, dtype=time_dt) + data['amplifier_data'] = np.zeros([header['num_amplifier_channels'], + AMP_SAMPLES * num_datablocks], + dtype=np.uint16) + data['aux_input_data'] = np.zeros([header['num_aux_input_channels'], + AUX_SAMPLES * num_datablocks], + dtype=np.uint16) + data['supply_voltage_data'] = np.zeros([header['num_supply_voltage_channels'], + SUPPLY_SAMPLES * num_datablocks], + dtype=np.uint16) + data['temp_sensor_data'] = np.zeros([header['num_temp_sensor_channels'], + TEMP_SAMPLES * num_datablocks], + dtype=np.uint16) + data['board_adc_data'] = np.zeros([header['num_board_adc_channels'], + ADC_SAMPLES * num_datablocks], + dtype=np.uint16) + data['board_dig_in_data'] = np.zeros([header['num_board_dig_in_channels'], + DIGIN_SAMPLES * num_datablocks], + dtype=np.uint) + data['board_dig_in_raw'] = np.zeros(DIGIN_SAMPLES * num_datablocks, + dtype=np.uint) + data['board_dig_out_data'] = np.zeros([header['num_board_dig_out_channels'], + DIGOUT_SAMPLES * num_datablocks], + dtype=np.uint) + data['board_dig_out_raw'] = np.zeros(DIGOUT_SAMPLES * num_datablocks, + dtype=np.uint) + return data + +def read_data_blocks(data, header, fid, datablocks_per_chunk=1): + """Reads a number of 60-sample 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 + """ # 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'] == 1 and header['version']['minor'] >= 2) or (header['version']['major'] > 1): + if ((header['version']['major'], header['version']['minor']) >= (1, 2)): ALL_DTYPES[0] = TIME_DTYPE_NEW else: ALL_DTYPES[0] = TIME_DTYPE_OLD @@ -56,18 +105,19 @@ def read_data_blocks(data, header, indices, fid, datablocks_per_chunk=1): chunk = np.fromfile(fid, dtype=np.dtype(db_dtype), count=datablocks_per_chunk) for idx,db in enumerate(chunk): - data['t_amplifier'][(indices['amplifier'] + idx * TIME_SAMPLES):(indices['amplifier'] + (idx + 1) * TIME_SAMPLES)] = db['time'] + data['t_amplifier'][(idx * TIME_SAMPLES):((idx + 1) * TIME_SAMPLES)] = db['time'] if amp_chan: - data['amplifier_data'][range(amp_chan), (indices['amplifier'] + idx * AMP_SAMPLES):(indices['amplifier'] + (idx + 1) * AMP_SAMPLES)] = db['amp'].reshape(amp_chan, AMP_SAMPLES) + data['amplifier_data'][range(amp_chan), (idx * AMP_SAMPLES):((idx + 1) * AMP_SAMPLES)] = db['amp'].reshape(amp_chan, AMP_SAMPLES) if aux_chan: - data['aux_input_data'][range(aux_chan), (indices['aux_input'] + idx * AUX_SAMPLES):(indices['aux_input'] + (idx + 1) * AUX_SAMPLES)] = db['aux'].reshape(aux_chan, AUX_SAMPLES) + data['aux_input_data'][range(aux_chan), (idx * AUX_SAMPLES):((idx + 1) * AUX_SAMPLES)] = db['aux'].reshape(aux_chan, AUX_SAMPLES) if supply_chan: - data['supply_voltage_data'][range(supply_chan), (indices['supply_voltage'] + idx * SUPPLY_SAMPLES):(indices['supply_voltage'] + (idx + 1) * SUPPLY_SAMPLES)] = db['supply'].reshape(supply_chan, SUPPLY_SAMPLES) + data['supply_voltage_data'][range(supply_chan), (idx * SUPPLY_SAMPLES):((idx + 1) * SUPPLY_SAMPLES)] = db['supply'].reshape(supply_chan, SUPPLY_SAMPLES) if temp_chan: - data['temp_sensor_data'][range(temp_chan), (indices['supply_voltage'] + idx * TEMP_SAMPLES):(indices['supply_voltage'] + (idx + 1) * TEMP_SAMPLES)] = db['temp'].reshape(temp_chan, TEMP_SAMPLES) + data['temp_sensor_data'][range(temp_chan), (idx * TEMP_SAMPLES):((idx + 1) * TEMP_SAMPLES)] = db['temp'].reshape(temp_chan, TEMP_SAMPLES) if adc_chan: - data['board_adc_data'][range(adc_chan), (indices['board_adc'] + idx * ADC_SAMPLES):(indices['board_adc'] + (idx + 1) * ADC_SAMPLES)] = db['adc'].reshape(adc_chan, ADC_SAMPLES) + data['board_adc_data'][range(adc_chan), (idx * ADC_SAMPLES):((idx + 1) * ADC_SAMPLES)] = db['adc'].reshape(adc_chan, ADC_SAMPLES) if digin_chan: - data['board_dig_in_raw'][(indices['board_dig_in'] + idx * DIGIN_SAMPLES):(indices['board_dig_in'] + (idx + 1) * DIGIN_SAMPLES)] = db['digin'].reshape(digin_chan, DIGIN_SAMPLES) + data['board_dig_in_raw'][(idx * DIGIN_SAMPLES):((idx + 1) * DIGIN_SAMPLES)] = db['digin'].reshape(digin_chan, DIGIN_SAMPLES) if digout_chan: - data['board_dig_out_raw'][(indices['board_dig_out'] + idx * DIGOUT_SAMPLES):(indices['board_dig_out'] + (idx + 1) * DIGOUT_SAMPLES)] = db['digout'].reshape(digout_chan, DIGOUT_SAMPLES) + data['board_dig_out_raw'][(idx * DIGOUT_SAMPLES):((idx + 1) * DIGOUT_SAMPLES)] = db['digout'].reshape(digout_chan, DIGOUT_SAMPLES) + diff --git a/bark/io/rhd/rhd2bark.py b/bark/io/rhd/rhd2bark.py index e16d09a..069d73a 100644 --- a/bark/io/rhd/rhd2bark.py +++ b/bark/io/rhd/rhd2bark.py @@ -1,12 +1,14 @@ import sys import os.path import arrow +import itertools from dateutil import tz import numpy as np -import bark.io.rhd.load_intan_rhd_format -import bark.io.rhd.legacy_load_intan_rhd_format +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 @@ -54,12 +56,31 @@ def bark_rhd_to_entry(): "--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, legacy=args.legacy, **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( @@ -123,13 +144,14 @@ 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): @@ -146,68 +168,87 @@ def rhds_to_entry(rhd_paths, 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. """ - if legacy: - read_data = bark.io.rhd.legacy_load_intan_rhd_format.read_data - else: - read_data = bark.io.rhd.load_intan_rhd_format.read_data if not timestamp: 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']) 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 cur,old in zip((adc_channels, amp_channels), + (adc_chan_names(first), amp_chan_names(first))): + if cur != old: + msg = '{} has channels {}\n{} has channels {}' + raise ValueError(msg.format(rhd_file, cur, 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)) + From a0aea55505bd7bb94303ca726abb865ab678007b Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 25 Jul 2018 14:33:08 -0500 Subject: [PATCH 24/37] Mangle colliding attribute names in arf2bark (#49) * Add check for create_entry positional arguments * Add command-line option for mangle prefix --- bark/io/arf2bark.py | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/bark/io/arf2bark.py b/bark/io/arf2bark.py index 6f2df4d..e5e0de6 100644 --- a/bark/io/arf2bark.py +++ b/bark/io/arf2bark.py @@ -8,6 +8,8 @@ import numpy import collections as coll +ENTRY_PREFIX = 'entry' + def _parse_args(raw_args): desc = 'Unspool an HDF5 ARF file into a Bark tree.' epi = 'Fails if bark_root already exists.' @@ -20,6 +22,10 @@ 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('bark_root', help='location of new bark root') return parser.parse_args(raw_args) @@ -31,7 +37,7 @@ 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_path, timezone, verbose): +def arf2bark(arf_file, root_path, timezone, verbose, mangle_prefix=ENTRY_PREFIX): with arf.open_file(arf_file, 'r') as af: os.mkdir(root_path) root = bark.Root(root_path) @@ -43,6 +49,21 @@ def arf2bark(arf_file, root_path, 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) @@ -116,7 +137,7 @@ def transfer_dset(ds_name, ds, e_path, verbose=False): def _main(): args = _parse_args(sys.argv[1:]) - arf2bark(args.arf_file, args.bark_root, args.timezone, args.verbose) + arf2bark(args.arf_file, args.bark_root, args.timezone, args.verbose, args.mangle_prefix) if __name__ == '__main__': _main() From 07b466b55db25f6c135a54ab63e8577b24d88a51 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Sat, 25 Aug 2018 16:38:55 -0500 Subject: [PATCH 25/37] Add git instructions (#52) --- docs/CONTRIBUTING.md | 2 + docs/git_workflow.md | 102 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 104 insertions(+) create mode 100644 docs/git_workflow.md diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md index 5bd4db2..c0ee8a9 100644 --- a/docs/CONTRIBUTING.md +++ b/docs/CONTRIBUTING.md @@ -10,6 +10,8 @@ $ 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.) ``` 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. From 3615aa2b153cfb2a61884b1db8c4c3b648fc371f Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Sat, 25 Aug 2018 16:39:44 -0500 Subject: [PATCH 26/37] Write time bins even if amplifiers are absent (#55) --- bark/io/rhd/data_to_result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bark/io/rhd/data_to_result.py b/bark/io/rhd/data_to_result.py index 5906dbf..7f0456d 100644 --- a/bark/io/rhd/data_to_result.py +++ b/bark/io/rhd/data_to_result.py @@ -9,11 +9,11 @@ def data_to_result(header, data, data_present): result['notes'] = header['notes'] result['frequency_parameters'] = header['frequency_parameters'] + 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['t_amplifier'] = data['t_amplifier'] result['spike_triggers'] = header['spike_triggers'] if header['num_aux_input_channels'] > 0: From 74a5df9908c6589e658801a437bd95bcde695e09 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Tue, 11 Dec 2018 11:22:42 -0600 Subject: [PATCH 27/37] Bugfix for issue 56 (#57) * Preserve specified order when using --col-attr * Ignore channels without attribute being selected on --- bark/tools/barkutils.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/bark/tools/barkutils.py b/bark/tools/barkutils.py index ece0bfd..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) From c1ad7e05dc6fb211243f8ff7affb43e79e992d4b Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 13 Feb 2019 18:51:12 -0600 Subject: [PATCH 28/37] Add .mda I/O (#58) * Add .mda I/O * Update readme with new scripts * Add attribution note * Remove script-style __name__ guard --- bark/io/mda.py | 217 +++++++++++++++++++++++++++++++++++++++++++++++++ docs/README.md | 2 + setup.py | 2 + 3 files changed, 221 insertions(+) create mode 100644 bark/io/mda.py diff --git a/bark/io/mda.py b/bark/io/mda.py new file mode 100644 index 0000000..4b702ec --- /dev/null +++ b/bark/io/mda.py @@ -0,0 +1,217 @@ +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(' Date: Wed, 13 Mar 2019 16:23:15 -0500 Subject: [PATCH 29/37] .mda transpose bugfix (#60) * Correct .mda serialization bug * Avoid transposing if possible * Remove confusing information in docstring --- bark/io/mda.py | 50 +++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 41 insertions(+), 9 deletions(-) diff --git a/bark/io/mda.py b/bark/io/mda.py index 4b702ec..3305cd1 100644 --- a/bark/io/mda.py +++ b/bark/io/mda.py @@ -59,21 +59,39 @@ def read_mda_header(filename): raise ValueError('invalid datatype code: {}'.format(dt_code)) return MdaHeader(datatype, dimensions, uses_64bit_dims) -def read_mda_data(filename): +def read_mda_data(filename, chan_axis='columns'): """Read an .mda file's data. Args: filename (str): path to the .mda file to read + chan_axis ('rows' or 'columns'): which axis of the array returned should + correspond to channels. To get Bark-formatted data out, 'columns' + is appropriate. Returns: numpy.ndarray: Numpy array containing the data + + Raises: + ValueError: if `chan_axis` is neither 'rows' nor 'columns' """ hdr = read_mda_header(filename) + # a proper .mda file is written with channels as rows and serialization in + # 'F' or "Fortran" format. To get channels as columns (i.e., the transpose + # of the .mda standard) you read in 'C' format and reverse the dimensions + if chan_axis == 'rows': + order = 'F' + array_shape = ndr.dims + elif chan_axis == 'columns': + order = 'C' + array_shape = tuple(reversed(hdr.dims)) + else: + raise ValueError('chan_axis must be either "rows" or "columns"') with open(filename, 'rb') as mda_file: mda_file.seek(hdr.header_size) - return np.fromfile(mda_file, dtype=hdr.dt).reshape(hdr.dims, order='F') + return np.fromfile(mda_file, dtype=hdr.dt).reshape(array_shape, + order=order) -def write_mda_file(filename, data, dtype=None): +def write_mda_file(filename, data, dtype=None, chan_axis='columns'): """Write an array to a .mda file (including a header). The only departure from the .mda spec is that the dimension sizes are @@ -87,13 +105,27 @@ def write_mda_file(filename, data, dtype=None): data (ndarray): data to write dtype (str, or None): numpy-compatible datatype string; if `None`, the data is written in its current datatype + chan_axis ('rows' or 'columns'): which axis of `data` corresponds to + channels. For Bark datasets, 'columns' is the standard. Returns: None Raises: - ValueError: if `dtype` is not supported + ValueError: if `dtype` is not supported, or if `chan_axis` is neither + 'rows' nor 'columns' """ + # 'order' and 'array_shape' are chosen so that a proper .mda file is written + # this corresponds to channels as rows and serialization in 'F' or "Fortran" + # format + if chan_axis == 'rows': + order = 'F' + array_shape = data.shape + elif chan_axis == 'columns': + order = 'C' + array_shape = tuple(reversed(data.shape)) + else: + raise ValueError('chan_axis must be either "rows" or "columns"') if dtype is None: dtype = str(data.dtype) bytes_per_entry = np.dtype(dtype).itemsize @@ -106,12 +138,12 @@ def write_mda_file(filename, data, dtype=None): # to save some needless complexity, always use 64-bit ints for dim sizes ndim_code = -1 * data.ndim mda_file.write(struct.pack(' Date: Wed, 10 Apr 2019 16:16:34 -0500 Subject: [PATCH 30/37] Add col/row to spike_times_dataframe_from_array() (#62) --- bark/io/mda.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/bark/io/mda.py b/bark/io/mda.py index 3305cd1..aa41e65 100644 --- a/bark/io/mda.py +++ b/bark/io/mda.py @@ -80,7 +80,7 @@ def read_mda_data(filename, chan_axis='columns'): # of the .mda standard) you read in 'C' format and reverse the dimensions if chan_axis == 'rows': order = 'F' - array_shape = ndr.dims + array_shape = hdr.dims elif chan_axis == 'columns': order = 'C' array_shape = tuple(reversed(hdr.dims)) @@ -145,7 +145,11 @@ def write_mda_file(filename, data, dtype=None, chan_axis='columns'): data = data.transpose() data.tofile(mda_file) -def spike_times_dataframe_from_mda(mda_hdr, mda_data, sampling_rate, keep=None): +def spike_times_dataframe_from_array(mda_hdr, + mda_data, + sampling_rate, + chan_axis='columns', + keep=None): """Converts data from an .mda file to a Pandas DataFrame usable by Bark. Args: @@ -154,6 +158,8 @@ def spike_times_dataframe_from_mda(mda_hdr, mda_data, sampling_rate, keep=None): MountainSort output file 'firings.mda' sampling_rate (number): the sampling rate corresponding to the cluster times in the .mda data + chan_axis('rows' or 'columns'): which axis of `mda_data` corresponds to + channels. For Bark-formatted datasets, 'columns' is standard. keep (iterable or None): the clusters to keep (all others are dropped); if None, all clusters are kept @@ -161,10 +167,17 @@ def spike_times_dataframe_from_mda(mda_hdr, mda_data, sampling_rate, keep=None): pandas.DataFrame: with columns 'amplitude', 'name', and 'start' """ df = pd.DataFrame(columns=['center_channel', 'amplitude', 'name', 'start']) - df['center_channel'] = mda_data[0] # amplitude is not currently provided by MountainSort, so it'll be NaNs - df['name'] = mda_data[2].astype('int16') - df['start'] = mda_data[1] / sampling_rate + if chan_axis == 'rows': + df['center_channel'] = mda_data[0] + df['name'] = mda_data[2].astype('int16') + df['start'] = mda_data[1] / sampling_rate + elif chan_axis == 'columns': + df['center_channel'] = mda_data[:, 0] + df['name'] = mda_data[:, 2].astype('int16') + df['start'] = mda_data[:, 1] / sampling_rate + else: + raise ValueError('chan_axis must be either "rows" or "columns"') if keep is None: keep = df['name'].unique() return df[df['name'].isin(keep)] @@ -222,7 +235,7 @@ def bark_event_ds_from_mda(mda_file, out_file, sampling_rate, keep=None): """ hdr = read_mda_header(mda_file) data = read_mda_data(mda_file) - df = spike_times_dataframe_from_mda(hdr, data, sampling_rate, keep) + df = spike_times_dataframe_from_array(hdr, data, sampling_rate, keep=keep) metadata = bark_metadata_from_df(hdr, df, sampling_rate) return bark.write_events(out_file, df, **metadata) From c00226bda8d147f691c87354de5754381156aa8d Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 28 Jul 2021 17:39:26 -0500 Subject: [PATCH 31/37] Replace StopIteration on exhaustion with return (#65) --- bark/stream.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From d6dfba585b0cd4fe7a9c195dfee98a4474f6e621 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 28 Jul 2021 17:55:18 -0500 Subject: [PATCH 32/37] Remove divide-by-zero from tests (#67) Close #66 --- tests/test_stream.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 From 3c10ab1170dc3fc3b07c8acb36538ab70c143b69 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 28 Jul 2021 21:25:25 -0500 Subject: [PATCH 33/37] Replace deprecated matplotlib method (#69) `bark-label-view` now runs successfully under matplotlib up to 3.4.2. Fixes #68 --- bark/tools/labelview.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/bark/tools/labelview.py b/bark/tools/labelview.py index 38fb3a0..2d35165 100644 --- a/bark/tools/labelview.py +++ b/bark/tools/labelview.py @@ -208,13 +208,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), @@ -231,7 +231,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, From a6239b6f07503d987825bd0bb3b02c31017fc62d Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Wed, 28 Jul 2021 22:00:42 -0500 Subject: [PATCH 34/37] Update installation instructions (#70) This incorporates information gained and fixes produced during #69 , #67 , #65 , #61 , #53 . --- docs/README.md | 33 +++++++++++++++++++++++---------- requirements.txt | 7 +++---- 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/docs/README.md b/docs/README.md index 1bf00dc..7aacc78 100644 --- a/docs/README.md +++ b/docs/README.md @@ -59,13 +59,6 @@ This repository contains: The python interface requires Python 3.5+. Installation with [Conda](http://conda.pydata.org/miniconda.html) is recommended. -If any error happens, please make sure your matplotlib version is 2.0.2 and the pyqt version is 5.6.0. - - git clone https://github.com/kylerbrown/resin - cd resin - pip install . - cd .. - git clone https://github.com/margoliashlab/bark cd bark @@ -75,9 +68,29 @@ If any error happens, please make sure your matplotlib version is 2.0.2 and the # 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. + +This simple installation supports the main bark library, most of the conversion scripts, +and most of the 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 (if you don't intend to use them, +you don't need to worry about this): + +* `bark-label-view` (for hand-labeling audio data), requires: + * Matplotlib (>=2.0) + * the spectral analysis library [`resin`](https://github.com/margoliashlab/resin) +* `bark-psg-view` (for hand-scoring PSG data), requires: + * Matplotlib (2.0.2) + * PyQt (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 diff --git a/requirements.txt b/requirements.txt index a7d4ab8..e963261 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,8 +1,7 @@ +numpy<1.17 +scipy +pandas<1.3 pytest -pandas -numpy PyYAML -scipy arrow==0.10.0 dataset==0.8.0 -pyqt5 From 9892b8f629d1a3e97eb5bc16efb0d4257ecec8c8 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Fri, 30 Jul 2021 18:54:17 -0500 Subject: [PATCH 35/37] Enable bark-label-view split/add on OS X (#71) Fix #44 . `Ctrl+click` is difficult to get matplotlib backends to recognize on Macs. Some don't register any key presses during mouse clicks, and others don't register `Ctrl+click`. Recognition of `Cmd+click` is inconsistent as well. For this reason, **for Mac only**, the control sequences for splitting an interval and adding a new interval in `bark-label-view` have been changed to `Shift+click`. They remain `Ctrl+click` for linux and windows. In the process of investigating this, I discovered that the built-in `TkAgg` backend sometimes doesn't recognize `Shift+click` the first time it is produced; it does recognize subsequent attempts. This can be a pain point in some workflows, so I added a check for the PyQt5 backend `Qt5Agg`, which doesn't have this problem, falling back on `TkAgg` if PyQt5 isn't installed. It's a large library to have as a dependency for just this feature, so I have labeled it "optional" in the installation section. As I note there, `bark-label-view` is still usable under the `TkAgg` backend, just occasionally slightly less smooth. --- bark/tools/labelview.py | 38 +++++++++++++++++++++++++++++--------- docs/README.md | 14 +++++++------- 2 files changed, 36 insertions(+), 16 deletions(-) diff --git a/bark/tools/labelview.py b/bark/tools/labelview.py index 2d35165..66e7c87 100644 --- a/bark/tools/labelview.py +++ b/bark/tools/labelview.py @@ -12,12 +12,32 @@ 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. Qt5 does capture them correctly. + # 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 - matplotlib.use('Qt5Agg') + 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' -import matplotlib.pyplot as plt help_string = ''' @@ -39,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. @@ -364,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/docs/README.md b/docs/README.md index 7aacc78..488a5b5 100644 --- a/docs/README.md +++ b/docs/README.md @@ -69,19 +69,21 @@ The python interface requires Python 3.5+. Installation with [Conda](http://cond # optional tests pytest -v -This simple installation supports the main bark library, most of the conversion scripts, -and most of the command-line data manipulation tools. Exceptions are noted below. +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 (if you don't intend to use them, -you don't need to worry about 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) - * PyQt (5.6.0) + * 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` @@ -121,8 +123,6 @@ 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. From 26dbbf7b4f622ec747bfdc30a0d914ef19b06593 Mon Sep 17 00:00:00 2001 From: Graham Fetterman Date: Mon, 2 Aug 2021 16:16:03 -0500 Subject: [PATCH 36/37] Support Intan GUI v3 (#72) The transition from version 1 to version 3 of the Intan GUI brought two principal changes to `.rhd` files: 1. The version 2+ GUI writes an additional header qstring field (`digital_reference_channel`) to reflect new functionality (the GUI allows software referencing before writing data to disk). 2. Files now contain 128 samples per "data block", rather than the 60 of version 1. Another small, but annoying, change is that ports are now 1-indexed, where before they were 0-indexed. These changes have been incorporated into the `.rhd`-to-bark conversion code in `bark-convert-rhd` (along with some code cleanup). They also make it seem prudent to begin noting the Intan GUI version that produced a given set of data; that information is now included in the entry metadata produced by `bark-convert-rhd`. These changes apply equally to the "legacy" code that more closely follows that provided by Intan (found on their [downloads page](https://intantech.com/downloads.html?tabSelect=Software) and to the streaming version I wrote a couple of years ago to reduce intolerable conversion times on large datasets. However, the legacy/streaming split has not been relevant so far (and they have been verified to produce identical outputs across a wide range of files). The time may come soon to eliminate the legacy version entirely. --- bark/io/rhd/constants.py | 35 ++ bark/io/rhd/data_to_result.py | 111 ++--- bark/io/rhd/get_bytes_per_data_block.py | 94 ++-- bark/io/rhd/legacy_load_intan_rhd_format.py | 505 ++++++++++---------- bark/io/rhd/legacy_read_one_data_block.py | 149 ++++-- bark/io/rhd/load_intan_rhd_format.py | 178 +++---- bark/io/rhd/qstring.py | 74 ++- bark/io/rhd/read_data_blocks.py | 204 +++++--- bark/io/rhd/read_header.py | 331 +++++++------ bark/io/rhd/rhd2bark.py | 29 +- 10 files changed, 956 insertions(+), 754 deletions(-) create mode 100644 bark/io/rhd/constants.py diff --git a/bark/io/rhd/constants.py b/bark/io/rhd/constants.py new file mode 100644 index 0000000..e7f4cfd --- /dev/null +++ b/bark/io/rhd/constants.py @@ -0,0 +1,35 @@ +import numpy as np + +LAST_TESTED_MAJOR_VERSION = 3 + +# Datatypes for different channel types + +TIMESTAMP_DTYPE_LE_V1_1 = np.dtype('uint32') +TIMESTAMP_DTYPE_GE_V1_2 = np.dtype('int32') + +AMPLIFIER_DTYPE = np.dtype('uint16') + +AUXILIARY_DTYPE = np.dtype('uint16') + +SUPPLY_DTYPE = np.dtype('uint16') + +TEMP_DTYPE = np.dtype('uint16') + +ADC_DTYPE = np.dtype('uint16') + +DIG_IN_DTYPE = np.dtype('uint16') +DIG_OUT_DTYPE = np.dtype('uint16') + +# Voltage scaling for different channel types + +AMPLIFIER_BIT_MICROVOLTS = 0.195 + +AUX_BIT_VOLTS = 37.4e-6 + +SUPPLY_BIT_VOLTS = 74.8e-6 + +ADC_BIT_VOLTS_MODE_0 = 50.354e-6 +ADC_BIT_VOLTS_MODE_1 = 152.59e-6 +ADC_BIT_VOLTS_MODE_13 = 312.5e-6 + +TEMP_BIT_CELCIUS = 0.01 diff --git a/bark/io/rhd/data_to_result.py b/bark/io/rhd/data_to_result.py index 7f0456d..d4dd536 100644 --- a/bark/io/rhd/data_to_result.py +++ b/bark/io/rhd/data_to_result.py @@ -1,54 +1,57 @@ -#! /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['notes'] = header['notes'] - result['frequency_parameters'] = header['frequency_parameters'] - - 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 +#! /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 index 656c8d3..c640b49 100644 --- a/bark/io/rhd/legacy_load_intan_rhd_format.py +++ b/bark/io/rhd/legacy_load_intan_rhd_format.py @@ -1,262 +1,243 @@ -#! /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.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 - -# 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 +# 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 index cd5f887..cf09a84 100644 --- a/bark/io/rhd/legacy_read_one_data_block.py +++ b/bark/io/rhd/legacy_read_one_data_block.py @@ -1,44 +1,105 @@ -#! /bin/env python -# -# Michael Gibson 23 April 2015 - -import sys, struct -import numpy as np - -def read_one_data_block(data, header, indices, fid): - """Reads one 60-sample data block from fid into data, at the location indicated by indices.""" - - # In version 1.2, we moved from saving timestamps as unsigned - # integers to signed integers to accommodate negative (adjusted) - # timestamps for pretrigger data[' - if (header['version']['major'] == 1 and header['version']['minor'] >= 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))) - +#! /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 19eec15..16a85f7 100644 --- a/bark/io/rhd/load_intan_rhd_format.py +++ b/bark/io/rhd/load_intan_rhd_format.py @@ -2,26 +2,19 @@ # # Michael Gibson 17 July 2015 # Kyler Brown December 2016 -# Graham Fetterman July 2018 - -from __future__ import absolute_import, division, unicode_literals, print_function +# 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, AMP_SAMPLES +from bark.io.rhd.read_data_blocks import read_data_blocks, preallocate_memory from bark.io.rhd.data_to_result import data_to_result -# constants -AMPLIFIER_BIT_MICROVOLTS = 0.195 +from . import constants as const + 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, max_memory=0): """Reads Intan Technologies RHD2000 data file generated by evaluation board GUI. @@ -47,27 +40,21 @@ def read_data(filename, no_floats=False, max_memory=0): 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.' + msg = ('Warning: a notch filter ({}Hz) was applied in the GUI, ' + + 'but has not been applied here.') print(msg.format(header['notch_filter_frequency'])) - 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']))) + 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. @@ -80,21 +67,29 @@ def read_data(filename, no_floats=False, max_memory=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') + 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) - record_time = AMP_SAMPLES * num_data_blocks / header['sample_rate'] + 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. Amplifiers were sampled at {:0.2f} kS/s.'.format( - record_time, header['sample_rate'] / 1000)) + 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( - 'Header file contains no data. Amplifiers were sampled at {:0.2f} kS/s.'.format( - header['sample_rate'] / 1000)) + 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 @@ -102,10 +97,11 @@ def read_data(filename, no_floats=False, max_memory=0): # 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) + 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 i in range(chunks): + 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: @@ -133,74 +129,87 @@ def check_data_and_reformat(header, data, no_floats): 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']): - 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) + 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']): - 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) + 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'] = AMPLIFIER_BIT_MICROVOLTS + 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'] = AUX_BIT_VOLTS - extras['supply_bit_volts'] = SUPPLY_BIT_VOLTS - extras['temp_bit_celcius'] = TEMP_BIT_CELCIUS + 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'] = ADC_BIT_VOLTS_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'] = ADC_BIT_VOLTS_0 + 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) + 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. - 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 + 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: - data['board_adc_data'] = np.multiply( - ADC_BIT_VOLTS_1, (data['board_adc_data'].astype(np.int32) - - UINT16_BIT_OFFSET)) # units = volts + 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( - 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 + 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.not_equal(data['t_amplifier'][1:] - data[ - 't_amplifier'][:-1], 1)) + num_gaps = np.sum(np.diff(data['t_amplifier']) != 1) if num_gaps != 0: - print( - 'Warning: {0} gaps in timestamp data found. Time scale will not be uniform!'.format( - num_gaps)) + 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'][range(0, len(data[ - 't_amplifier']), 4)] - data['t_supply_voltage'] = data['t_amplifier'][range(0, len(data[ - 't_amplifier']), 60)] + 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'] - data['t_temp_sensor'] = data['t_supply_voltage'] # Move variables to result struct. result = data_to_result(header, data, data_present=True) @@ -210,4 +219,3 @@ def check_data_and_reformat(header, data, no_floats): 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_dt = np.int + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 else: - time_dt = np.uint - data['t_amplifier'] = np.zeros(AMP_SAMPLES * num_datablocks, dtype=time_dt) + 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'], - AMP_SAMPLES * num_datablocks], - dtype=np.uint16) + num_samples * num_datablocks], + dtype=const.AMPLIFIER_DTYPE) data['aux_input_data'] = np.zeros([header['num_aux_input_channels'], - AUX_SAMPLES * num_datablocks], - dtype=np.uint16) + (num_samples // 4) * num_datablocks], + dtype=const.AUXILIARY_DTYPE) data['supply_voltage_data'] = np.zeros([header['num_supply_voltage_channels'], - SUPPLY_SAMPLES * num_datablocks], - dtype=np.uint16) + 1 * num_datablocks], + dtype=const.SUPPLY_DTYPE) data['temp_sensor_data'] = np.zeros([header['num_temp_sensor_channels'], - TEMP_SAMPLES * num_datablocks], - dtype=np.uint16) + 1 * num_datablocks], + dtype=const.TEMP_DTYPE) data['board_adc_data'] = np.zeros([header['num_board_adc_channels'], - ADC_SAMPLES * num_datablocks], - dtype=np.uint16) + num_samples * num_datablocks], + dtype=const.ADC_DTYPE) data['board_dig_in_data'] = np.zeros([header['num_board_dig_in_channels'], - DIGIN_SAMPLES * num_datablocks], - dtype=np.uint) - data['board_dig_in_raw'] = np.zeros(DIGIN_SAMPLES * num_datablocks, - dtype=np.uint) + 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'], - DIGOUT_SAMPLES * num_datablocks], - dtype=np.uint) - data['board_dig_out_raw'] = np.zeros(DIGOUT_SAMPLES * num_datablocks, - dtype=np.uint) + 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 60-sample data blocks from fid into data. + """Reads a number of data blocks from fid into data. Args: data (dict of numpy arrays): having the same format as the return value @@ -81,43 +63,117 @@ def read_data_blocks(data, header, fid, datablocks_per_chunk=1): 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)): - ALL_DTYPES[0] = TIME_DTYPE_NEW + time_dtype = const.TIMESTAMP_DTYPE_GE_V1_2 else: - ALL_DTYPES[0] = TIME_DTYPE_OLD + time_dtype = const.TIMESTAMP_DTYPE_LE_V1_1 - time_chan = 1 - amp_chan = header['num_amplifier_channels'] - aux_chan = header['num_aux_input_channels'] - supply_chan = header['num_supply_voltage_channels'] - temp_chan = header['num_temp_sensor_channels'] - adc_chan = header['num_board_adc_channels'] - digin_chan = header['num_board_dig_in_channels'] - digout_chan = header['num_board_dig_out_channels'] - all_chans = [time_chan, amp_chan, aux_chan, supply_chan, temp_chan, adc_chan, digin_chan, digout_chan] + 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(NAMES, ALL_DTYPES, all_chans, ALL_SAMPLES)] + 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): - data['t_amplifier'][(idx * TIME_SAMPLES):((idx + 1) * TIME_SAMPLES)] = db['time'] - if amp_chan: - data['amplifier_data'][range(amp_chan), (idx * AMP_SAMPLES):((idx + 1) * AMP_SAMPLES)] = db['amp'].reshape(amp_chan, AMP_SAMPLES) - if aux_chan: - data['aux_input_data'][range(aux_chan), (idx * AUX_SAMPLES):((idx + 1) * AUX_SAMPLES)] = db['aux'].reshape(aux_chan, AUX_SAMPLES) - if supply_chan: - data['supply_voltage_data'][range(supply_chan), (idx * SUPPLY_SAMPLES):((idx + 1) * SUPPLY_SAMPLES)] = db['supply'].reshape(supply_chan, SUPPLY_SAMPLES) - if temp_chan: - data['temp_sensor_data'][range(temp_chan), (idx * TEMP_SAMPLES):((idx + 1) * TEMP_SAMPLES)] = db['temp'].reshape(temp_chan, TEMP_SAMPLES) - if adc_chan: - data['board_adc_data'][range(adc_chan), (idx * ADC_SAMPLES):((idx + 1) * ADC_SAMPLES)] = db['adc'].reshape(adc_chan, ADC_SAMPLES) - if digin_chan: - data['board_dig_in_raw'][(idx * DIGIN_SAMPLES):((idx + 1) * DIGIN_SAMPLES)] = db['digin'].reshape(digin_chan, DIGIN_SAMPLES) - if digout_chan: - data['board_dig_out_raw'][(idx * DIGOUT_SAMPLES):((idx + 1) * DIGOUT_SAMPLES)] = db['digout'].reshape(digout_chan, DIGOUT_SAMPLES) + 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(' Date: Tue, 17 Aug 2021 20:18:28 -0500 Subject: [PATCH 37/37] Update python version requirement --- docs/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/README.md b/docs/README.md index 488a5b5..bca899e 100644 --- a/docs/README.md +++ b/docs/README.md @@ -57,7 +57,7 @@ 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/margoliashlab/bark cd bark