From 5608c5a6e97b5eaa803e146caa07fa8ac3644263 Mon Sep 17 00:00:00 2001 From: Oliver Ruebel Date: Wed, 18 Sep 2024 23:31:36 -0700 Subject: [PATCH] Sync add_read with main branch (#101) * Merged main into add_read * Fix docs build for SpikeEventSeries * Fix code formatting * Fix segfault due to duplicate declaration of NWBFile.io parameter --- CMakeLists.txt | 2 + src/io/BaseIO.hpp | 9 ++ src/io/hdf5/HDF5IO.cpp | 10 ++ src/io/hdf5/HDF5IO.hpp | 8 ++ src/nwb/NWBFile.cpp | 142 ++++++++++++++++++---- src/nwb/NWBFile.hpp | 30 ++++- src/nwb/RecordingContainers.cpp | 16 +++ src/nwb/RecordingContainers.hpp | 18 ++- src/nwb/ecephys/SpikeEventSeries.cpp | 57 +++++++++ src/nwb/ecephys/SpikeEventSeries.hpp | 84 +++++++++++++ tests/examples/testWorkflowExamples.cpp | 3 + tests/examples/test_ecephys_data_read.cpp | 4 +- tests/testEcephys.cpp | 141 +++++++++++++++++++-- tests/testNWBFile.cpp | 84 ++++++++++++- tests/testRecordingWorkflow.cpp | 12 +- tests/testUtils.hpp | 16 ++- 16 files changed, 590 insertions(+), 46 deletions(-) create mode 100644 src/nwb/ecephys/SpikeEventSeries.cpp create mode 100644 src/nwb/ecephys/SpikeEventSeries.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 56b94e4c..5e5fa408 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,6 +28,7 @@ add_library( src/nwb/base/TimeSeries.cpp src/nwb/device/Device.cpp src/nwb/ecephys/ElectricalSeries.cpp + src/nwb/ecephys/SpikeEventSeries.cpp src/nwb/file/ElectrodeGroup.cpp src/nwb/file/ElectrodeTable.cpp src/nwb/hdmf/base/Container.cpp @@ -77,6 +78,7 @@ target_include_directories( ) target_compile_features(aqnwb_aqnwb PUBLIC cxx_std_17) +target_compile_definitions(aqnwb_aqnwb PUBLIC BOOST_NO_CXX98_FUNCTION_BASE) # ---- Additional libraries needed ---- find_package(HDF5 REQUIRED COMPONENTS CXX) diff --git a/src/io/BaseIO.hpp b/src/io/BaseIO.hpp index e2a11210..de5445a0 100644 --- a/src/io/BaseIO.hpp +++ b/src/io/BaseIO.hpp @@ -343,6 +343,14 @@ class BaseIO virtual std::unique_ptr getDataSet( const std::string& path) = 0; + /** + * @brief Checks whether a Dataset, Group, or Link already exists at the + * location in the file. + * @param path The location of the object in the file. + * @return Whether the object exists. + */ + virtual bool objectExists(const std::string& path) = 0; + /** * @brief Convenience function for creating NWB related attributes. * @param path The location of the object in the file. @@ -376,6 +384,7 @@ class BaseIO * @return The status of the operation. */ Status createTimestampsAttributes(const std::string& path); + /** * @brief Returns true if the file is open. * @return True if the file is open, false otherwise. diff --git a/src/io/hdf5/HDF5IO.cpp b/src/io/hdf5/HDF5IO.cpp index d2d56947..3ee98604 100644 --- a/src/io/hdf5/HDF5IO.cpp +++ b/src/io/hdf5/HDF5IO.cpp @@ -736,6 +736,16 @@ bool HDF5IO::canModifyObjects() return statusOK && !inSWMRMode; } +bool HDF5IO::objectExists(const std::string& path) +{ + htri_t exists = H5Lexists(file->getId(), path.c_str(), H5P_DEFAULT); + if (exists > 0) { + return true; + } else { + return false; + } +} + std::unique_ptr HDF5IO::getDataSet( const std::string& path) { diff --git a/src/io/hdf5/HDF5IO.hpp b/src/io/hdf5/HDF5IO.hpp index c0e76c28..e3d3cffb 100644 --- a/src/io/hdf5/HDF5IO.hpp +++ b/src/io/hdf5/HDF5IO.hpp @@ -286,6 +286,14 @@ class HDF5IO : public BaseIO std::unique_ptr getDataSet( const std::string& path) override; + /** + * @brief Checks whether a Dataset, Group, or Link already exists at the + * location in the file. + * @param path The location of the object in the file. + * @return Whether the object exists. + */ + bool objectExists(const std::string& path) override; + /** * @brief Returns the HDF5 type of object at a given path. * @param path The location in the file of the object. diff --git a/src/nwb/NWBFile.cpp b/src/nwb/NWBFile.cpp index 974fdeb5..b6667ec9 100644 --- a/src/nwb/NWBFile.cpp +++ b/src/nwb/NWBFile.cpp @@ -13,15 +13,18 @@ #include "io/BaseIO.hpp" #include "nwb/device/Device.hpp" #include "nwb/ecephys/ElectricalSeries.hpp" +#include "nwb/ecephys/SpikeEventSeries.hpp" #include "nwb/file/ElectrodeGroup.hpp" -#include "nwb/file/ElectrodeTable.hpp" #include "spec/core.hpp" #include "spec/hdmf_common.hpp" #include "spec/hdmf_experimental.hpp" using namespace AQNWB::NWB; -constexpr SizeType CHUNK_XSIZE = 2048; +constexpr SizeType CHUNK_XSIZE = + 2048; // TODO - replace these with io settings input +constexpr SizeType SPIKE_CHUNK_XSIZE = + 8; // TODO - replace with io settings input std::vector NWBFile::emptyContainerIndexes = {}; @@ -58,10 +61,8 @@ Status NWBFile::createFileStructure(const std::string& identifierText, if (!io->canModifyObjects()) { return Status::Failure; } - io->createCommonNWBAttributes("/", "core", "NWBFile", ""); io->createAttribute(AQNWB::SPEC::CORE::version, "/", "nwb_version"); - io->createGroup("/acquisition"); io->createGroup("/analysis"); io->createGroup("/processing"); @@ -94,12 +95,12 @@ Status NWBFile::createFileStructure(const std::string& identifierText, io->createStringDataSet("/session_start_time", time); io->createStringDataSet("/timestamps_reference_time", time); io->createStringDataSet("/identifier", identifierText); - return Status::Success; } Status NWBFile::createElectricalSeries( std::vector recordingArrays, + std::vector recordingNames, const IO::BaseDataType& dataType, RecordingContainers* recordingContainers, std::vector& containerIndexes) @@ -108,26 +109,43 @@ Status NWBFile::createElectricalSeries( return Status::Failure; } - // store all recorded data in the acquisition group - std::string rootPath = "/acquisition/"; + if (recordingNames.size() != recordingArrays.size()) { + return Status::Failure; + } - // Setup electrode table - ElectrodeTable elecTable = ElectrodeTable(io); - elecTable.initialize(); + // Setup electrode table if it was not yet created + bool electrodeTableCreated = + io->objectExists(ElectrodeTable::electrodeTablePath); + if (!electrodeTableCreated) { + elecTable = std::make_unique(io); + elecTable->initialize(); + + // Add electrode information to table (does not write to datasets yet) + for (const auto& channelVector : recordingArrays) { + elecTable->addElectrodes(channelVector); + } + } + + // Create datasets + for (size_t i = 0; i < recordingArrays.size(); ++i) { + const auto& channelVector = recordingArrays[i]; + const std::string& recordingName = recordingNames[i]; - // Create continuous datasets - for (const auto& channelVector : recordingArrays) { // Setup electrodes and devices std::string groupName = channelVector[0].groupName; std::string devicePath = "/general/devices/" + groupName; std::string electrodePath = "/general/extracellular_ephys/" + groupName; - std::string electricalSeriesPath = rootPath + groupName; + std::string electricalSeriesPath = acquisitionPath + "/" + recordingName; - Device device = Device(devicePath, io); - device.initialize("description", "unknown"); + // Check if device exists for groupName, create device and electrode group + // if not + if (!io->objectExists(devicePath)) { + Device device = Device(devicePath, io); + device.initialize("description", "unknown"); - ElectrodeGroup elecGroup = ElectrodeGroup(electrodePath, io); - elecGroup.initialize("description", "unknown", device); + ElectrodeGroup elecGroup = ElectrodeGroup(electrodePath, io); + elecGroup.initialize("description", "unknown", device); + } // Setup electrical series datasets auto electricalSeries = @@ -141,14 +159,94 @@ Status NWBFile::createElectricalSeries( SizeArray {CHUNK_XSIZE, 0}); recordingContainers->addContainer(std::move(electricalSeries)); containerIndexes.push_back(recordingContainers->containers.size() - 1); + } + + // write electrode information to datasets + // (requires that the ElectrodeGroup has been written) + if (!electrodeTableCreated) { + elecTable->finalize(); + } + + return Status::Success; +} - // Add electrode information to electrode table (does not write to datasets - // yet) - elecTable.addElectrodes(channelVector); +Status NWBFile::createSpikeEventSeries( + std::vector recordingArrays, + std::vector recordingNames, + const IO::BaseDataType& dataType, + RecordingContainers* recordingContainers, + std::vector& containerIndexes) +{ + if (!io->canModifyObjects()) { + return Status::Failure; + } + + if (recordingNames.size() != recordingArrays.size()) { + return Status::Failure; + } + + // Setup electrode table if it was not yet created + bool electrodeTableCreated = + io->objectExists(ElectrodeTable::electrodeTablePath); + if (!electrodeTableCreated) { + elecTable = std::make_unique(io); + elecTable->initialize(); + + // Add electrode information to table (does not write to datasets yet) + for (const auto& channelVector : recordingArrays) { + elecTable->addElectrodes(channelVector); + } + } + + // Create datasets + for (size_t i = 0; i < recordingArrays.size(); ++i) { + const auto& channelVector = recordingArrays[i]; + const std::string& recordingName = recordingNames[i]; + + // Setup electrodes and devices + std::string groupName = channelVector[0].groupName; + std::string devicePath = "/general/devices/" + groupName; + std::string electrodePath = "/general/extracellular_ephys/" + groupName; + std::string spikeEventSeriesPath = acquisitionPath + "/" + recordingName; + + // Check if device exists for groupName, create device and electrode group + // if not + if (!io->objectExists(devicePath)) { + Device device = Device(devicePath, io); + device.initialize("description", "unknown"); + + ElectrodeGroup elecGroup = ElectrodeGroup(electrodePath, io); + elecGroup.initialize("description", "unknown", device); + } + + // Setup Spike Event Series datasets + SizeArray dsetSize; + SizeArray chunkSize; + if (channelVector.size() == 1) { + dsetSize = SizeArray {0, 0}; + chunkSize = SizeArray {SPIKE_CHUNK_XSIZE, 1}; + } else { + dsetSize = SizeArray {0, channelVector.size(), 0}; + chunkSize = SizeArray {SPIKE_CHUNK_XSIZE, 1, 1}; + } + + auto spikeEventSeries = + std::make_unique(spikeEventSeriesPath, io); + spikeEventSeries->initialize( + dataType, + channelVector, + "Stores spike waveforms from an extracellular ephys recording", + dsetSize, + chunkSize); + recordingContainers->addContainer(std::move(spikeEventSeries)); + containerIndexes.push_back(recordingContainers->containers.size() - 1); } // write electrode information to datasets - elecTable.finalize(); + // (requires that the ElectrodeGroup has been written) + if (!electrodeTableCreated) { + elecTable->finalize(); + } return Status::Success; } @@ -160,7 +258,7 @@ void NWBFile::cacheSpecifications( const std::array, N>& specVariables) { - io->createGroup("/specifications/" + specPath + "/"); + io->createGroup("/specifications/" + specPath); io->createGroup("/specifications/" + specPath + "/" + versionNumber); for (const auto& [name, content] : specVariables) { diff --git a/src/nwb/NWBFile.hpp b/src/nwb/NWBFile.hpp index 0ab8c4af..02646777 100644 --- a/src/nwb/NWBFile.hpp +++ b/src/nwb/NWBFile.hpp @@ -11,6 +11,7 @@ #include "io/BaseIO.hpp" #include "nwb/RecordingContainers.hpp" #include "nwb/base/TimeSeries.hpp" +#include "nwb/file/ElectrodeTable.hpp" /*! * \namespace AQNWB::NWB @@ -73,14 +74,38 @@ class NWBFile : public Container * @param recordingArrays vector of ChannelVector indicating the electrodes to * record from. A separate ElectricalSeries will be * created for each ChannelVector. + * @param recordingNames vector indicating the names of the ElectricalSeries + * within the acquisition group + * @param dataType The data type of the elements in the data block. * @param recordingContainers The container to store the created TimeSeries. * @param containerIndexes The indexes of the containers added to * recordingContainers - * @param dataType The data type of the elements in the data block. * @return Status The status of the object creation operation. */ Status createElectricalSeries( std::vector recordingArrays, + std::vector recordingNames, + const IO::BaseDataType& dataType = IO::BaseDataType::I16, + RecordingContainers* recordingContainers = nullptr, + std::vector& containerIndexes = emptyContainerIndexes); + + /** + * @brief Create SpikeEventSeries objects to record data into. + * Created objects are stored in recordingContainers. + * @param recordingArrays vector of ChannelVector indicating the electrodes to + * record from. A separate ElectricalSeries will be + * created for each ChannelVector. + * @param recordingNames vector indicating the names of the SpikeEventSeries + * within the acquisition group + * @param dataType The data type of the elements in the data block. + * @param recordingContainers The container to store the created TimeSeries. + * @param containerIndexes The indexes of the containers added to + * recordingContainers + * @return Status The status of the object creation operation. + */ + Status createSpikeEventSeries( + std::vector recordingArrays, + std::vector recordingNames, const IO::BaseDataType& dataType = IO::BaseDataType::I16, RecordingContainers* recordingContainers = nullptr, std::vector& containerIndexes = emptyContainerIndexes); @@ -132,7 +157,10 @@ class NWBFile : public Container const std::array, N>& specVariables); + std::unique_ptr elecTable; + const std::string identifierText; static std::vector emptyContainerIndexes; + inline const static std::string acquisitionPath = "/acquisition"; }; } // namespace AQNWB::NWB \ No newline at end of file diff --git a/src/nwb/RecordingContainers.cpp b/src/nwb/RecordingContainers.cpp index d7464bca..4658d551 100644 --- a/src/nwb/RecordingContainers.cpp +++ b/src/nwb/RecordingContainers.cpp @@ -2,6 +2,7 @@ #include "nwb/RecordingContainers.hpp" #include "nwb/ecephys/ElectricalSeries.hpp" +#include "nwb/ecephys/SpikeEventSeries.hpp" #include "nwb/hdmf/base/Container.hpp" using namespace AQNWB::NWB; @@ -63,3 +64,18 @@ Status RecordingContainers::writeElectricalSeriesData( es->writeChannel(channel.localIndex, numSamples, data, timestamps); } + +Status RecordingContainers::writeSpikeEventData(const SizeType& containerInd, + const SizeType& numSamples, + const SizeType& numChannels, + const void* data, + const void* timestamps) +{ + SpikeEventSeries* ses = + dynamic_cast(getContainer(containerInd)); + + if (ses == nullptr) + return Status::Failure; + + ses->writeSpike(numSamples, numChannels, data, timestamps); +} diff --git a/src/nwb/RecordingContainers.hpp b/src/nwb/RecordingContainers.hpp index aa003086..f38d84c7 100644 --- a/src/nwb/RecordingContainers.hpp +++ b/src/nwb/RecordingContainers.hpp @@ -71,7 +71,7 @@ class RecordingContainers const void* timestamps); /** - * @brief Write ElectricalSereis data to a recordingContainer dataset. + * @brief Write ElectricalSeries data to a recordingContainer dataset. * @param containerInd The index of the electrical series dataset within the * electrical series group. * @param channel The channel index to use for writing timestamps. @@ -89,6 +89,22 @@ class RecordingContainers const void* data, const void* timestamps); + /** + * @brief Write SpikeEventSeries data to a recordingContainer dataset. + * @param containerInd The index of the SpikeEventSeries dataset within the + * SpikeEventSeries containers. + * @param numSamples Number of samples in the time for the single event. + * @param numChannels Number of channels in the time for the single event. + * @param data A pointer to the data block. + * @param timestamps A pointer to the timestamps block + * @return The status of the write operation. + */ + Status writeSpikeEventData(const SizeType& containerInd, + const SizeType& numSamples, + const SizeType& numChannels, + const void* data, + const void* timestamps); + std::vector> containers; std::string name; }; diff --git a/src/nwb/ecephys/SpikeEventSeries.cpp b/src/nwb/ecephys/SpikeEventSeries.cpp new file mode 100644 index 00000000..9fc646b5 --- /dev/null +++ b/src/nwb/ecephys/SpikeEventSeries.cpp @@ -0,0 +1,57 @@ +#include "nwb/ecephys/SpikeEventSeries.hpp" + +using namespace AQNWB::NWB; + +// SpikeEventSeries + +/** Constructor */ +SpikeEventSeries::SpikeEventSeries(const std::string& path, + std::shared_ptr io) + : ElectricalSeries(path, io) +{ +} + +/** Destructor */ +SpikeEventSeries::~SpikeEventSeries() {} + +void SpikeEventSeries::initialize(const IO::BaseDataType& dataType, + const Types::ChannelVector& channelVector, + const std::string& description, + const SizeArray& dsetSize, + const SizeArray& chunkSize, + const float& conversion, + const float& resolution, + const float& offset) +{ + ElectricalSeries::initialize(dataType, + channelVector, + description, + dsetSize, + chunkSize, + conversion, + resolution, + offset); + + this->eventsRecorded = 0; +} + +Status SpikeEventSeries::writeSpike(const SizeType& numSamples, + const SizeType& numChannels, + const void* data, + const void* timestamps) +{ + // get offsets and datashape + std::vector dataShape; + std::vector positionOffset; + if (numChannels == 1) { + dataShape = {1, numSamples}; + positionOffset = {this->eventsRecorded, 0}; + } else { + dataShape = {1, numChannels, numSamples}; + positionOffset = {this->eventsRecorded, 0, 0}; + } + this->eventsRecorded += 1; + + // write channel data + return writeData(dataShape, positionOffset, data, timestamps); +} \ No newline at end of file diff --git a/src/nwb/ecephys/SpikeEventSeries.hpp b/src/nwb/ecephys/SpikeEventSeries.hpp new file mode 100644 index 00000000..804db0ac --- /dev/null +++ b/src/nwb/ecephys/SpikeEventSeries.hpp @@ -0,0 +1,84 @@ +#pragma once + +#include + +#include "Channel.hpp" +#include "io/BaseIO.hpp" +#include "nwb/ecephys/ElectricalSeries.hpp" + +namespace AQNWB::NWB +{ +/** + * @brief Stores snapshots/snippets of recorded spike events (i.e., threshold + * crossings). + */ +class SpikeEventSeries : public ElectricalSeries +{ +public: + /** + * @brief Constructor. + * @param path The location of the SpikeEventSeries in the file. + * @param io A shared pointer to the IO object. + */ + SpikeEventSeries(const std::string& path, std::shared_ptr io); + + /** + * @brief Destructor + */ + ~SpikeEventSeries(); + + /** + * @brief Initializes the Electrical Series + * + * @param dataType The data type to use for storing the recorded voltage + * @param channelVector The electrodes to use for recording + * @param description The description of the SpikeEventSeries, should describe + * how events were detected. + * @param dsetSize Initial size of the main dataset. This must be a vector + * with two elements. The first element specifies the length + * in time and the second element must be equal to the + * length of channelVector + * @param chunkSize Chunk size to use. The number of elements must be two to + * specify the size of a chunk in the time and electrode + * dimension + * @param conversion Scalar to multiply each element in data to convert it to + * the specified ‘unit’ + * @param resolution Smallest meaningful difference between values in data, + * stored in the specified by unit + * @param offset Scalar to add to the data after scaling by ‘conversion’ to + * finalize its coercion to the specified ‘unit' + */ + void initialize(const IO::BaseDataType& dataType, + const Types::ChannelVector& channelVector, + const std::string& description, + const SizeArray& dsetSize, + const SizeArray& chunkSize, + const float& conversion = 1.0f, + const float& resolution = -1.0f, + const float& offset = 0.0f); + + /** + * @brief Write a single spike series event + * + * @param numSamples The number of samples in the event + * @param numChannels The number of channels in the event + * @param data The data of the event + * @param timestamps The timestamps of the event + */ + Status writeSpike(const SizeType& numSamples, + const SizeType& numChannels, + const void* data, + const void* timestamps); + +private: + /** + * @brief The neurodataType of the SpikeEventSeries. + */ + std::string neurodataType = "SpikeEventSeries"; + + /** + * @brief The number of events already written. + */ + SizeType eventsRecorded; +}; +} // namespace AQNWB::NWB diff --git a/tests/examples/testWorkflowExamples.cpp b/tests/examples/testWorkflowExamples.cpp index c41092f7..5bfff8fb 100644 --- a/tests/examples/testWorkflowExamples.cpp +++ b/tests/examples/testWorkflowExamples.cpp @@ -27,6 +27,8 @@ TEST_CASE("workflowExamples") std::vector mockRecordingArrays = getMockChannelArrays(); + std::vector mockChannelNames = + getMockChannelArrayNames("esdata"); std::vector> mockData = getMockData2D(numSamples, numChannels); std::vector mockTimestamps = getMockTimestamps(numSamples); @@ -49,6 +51,7 @@ TEST_CASE("workflowExamples") // [example_workflow_datasets_snippet] std::vector containerIndexes; nwbfile->createElectricalSeries(mockRecordingArrays, + mockChannelNames, BaseDataType::I16, recordingContainers.get(), containerIndexes); diff --git a/tests/examples/test_ecephys_data_read.cpp b/tests/examples/test_ecephys_data_read.cpp index 7f9a3def..9a959881 100644 --- a/tests/examples/test_ecephys_data_read.cpp +++ b/tests/examples/test_ecephys_data_read.cpp @@ -28,6 +28,8 @@ TEST_CASE("ElectricalSeriesReadExample", "[ecephys]") SizeType numChannels = 2; std::vector mockArrays = getMockChannelArrays(); BaseDataType dataType = BaseDataType::F32; + std::vector mockChannelNames = + getMockChannelArrayNames("esdata"); std::vector> mockData = getMockData2D(numSamples, numChannels); std::vector mockTimestamps = getMockTimestamps(numSamples, 1); @@ -59,7 +61,7 @@ TEST_CASE("ElectricalSeriesReadExample", "[ecephys]") // create a new ElectricalSeries Status resultCreate = nwbfile.createElectricalSeries( - mockArrays, dataType, recordingContainers.get()); + mockArrays, mockChannelNames, dataType, recordingContainers.get()); REQUIRE(resultCreate == Status::Success); // get the new ElectricalSeries diff --git a/tests/testEcephys.cpp b/tests/testEcephys.cpp index ded15228..194bea18 100644 --- a/tests/testEcephys.cpp +++ b/tests/testEcephys.cpp @@ -9,6 +9,7 @@ #include "io/hdf5/HDF5IO.hpp" #include "nwb/device/Device.hpp" #include "nwb/ecephys/ElectricalSeries.hpp" +#include "nwb/ecephys/SpikeEventSeries.hpp" #include "nwb/file/ElectrodeGroup.hpp" #include "nwb/file/ElectrodeTable.hpp" #include "testUtils.hpp" @@ -32,16 +33,6 @@ TEST_CASE("ElectricalSeries", "[ecephys]") std::string devicePath = "/device"; std::string electrodePath = "/elecgroup/"; - SECTION("test initialization") - { - // TODO - } - - SECTION("test linking to electrode table region") - { - // TODO - } - SECTION("test writing channels") { // setup io object @@ -159,3 +150,133 @@ TEST_CASE("ElectricalSeries", "[ecephys]") REQUIRE_THAT(dataOut[1], Catch::Matchers::Approx(mockData[1]).margin(1)); } } + +TEST_CASE("SpikeEventSeries", "[ecephys]") +{ + // setup recording info + SizeType numSamples = 32; + SizeType numEvents = 10; + std::string dataPath = "/sesdata"; + BaseDataType dataType = BaseDataType::F32; + std::vector mockTimestamps = getMockTimestamps(numEvents, 1); + std::string devicePath = "/device"; + std::string electrodePath = "/elecgroup/"; + + SECTION("test writing events - events x channels x samples") + { + // setup mock data + SizeType numChannels = 4; + std::vector mockArrays = + getMockChannelArrays(numChannels); + std::vector> mockData = + getMockData2D(numSamples * numChannels, numEvents); + + // setup io object + std::string path = getTestFilePath("SpikeEventSeries3D.h5"); + std::shared_ptr io = createIO("HDF5", path); + io->open(); + io->createGroup("/general"); + io->createGroup("/general/extracellular_ephys"); + + // setup electrode table, device, and electrode group + NWB::ElectrodeTable elecTable = NWB::ElectrodeTable(io); + elecTable.initialize(); + + // setup electrical series + NWB::SpikeEventSeries ses = NWB::SpikeEventSeries(dataPath, io); + ses.initialize(dataType, + mockArrays[0], + "no description", + SizeArray {0, numChannels, numSamples}, + SizeArray {8, 1, 1}); + + // write channel data + for (SizeType e = 0; e < numEvents; ++e) { + double timestamp = mockTimestamps[e]; + ses.writeSpike(numSamples, numChannels, mockData[e].data(), ×tamp); + } + io->close(); + + // Read data back from file + std::unique_ptr file = + std::make_unique(path, H5F_ACC_RDONLY); + std::unique_ptr dataset = + std::make_unique(file->openDataSet(dataPath + "/data")); + std::vector> dataOut( + numEvents, std::vector(numSamples * numChannels)); + float* buffer = new float[numEvents * numSamples * numChannels]; + + H5::DataSpace fSpace = dataset->getSpace(); + hsize_t dims[3]; + fSpace.getSimpleExtentDims(dims, NULL); + hsize_t memdims = dims[0] * dims[1] * dims[2]; + dataset->read(buffer, H5::PredType::NATIVE_FLOAT, fSpace, fSpace); + + for (SizeType i = 0; i < numEvents; ++i) { + for (SizeType j = 0; j < (numSamples * numChannels); ++j) { + dataOut[i][j] = buffer[i * (numSamples * numChannels) + j]; + } + } + delete[] buffer; + REQUIRE_THAT(dataOut[0], Catch::Matchers::Approx(mockData[0]).margin(1)); + REQUIRE_THAT(dataOut[1], Catch::Matchers::Approx(mockData[1]).margin(1)); + } + + SECTION("test writing events - events x samples") + { + // setup mock data + std::vector mockArrays = getMockChannelArrays(1); + std::vector> mockData = + getMockData2D(numSamples, numEvents); + + // setup io object + std::string path = getTestFilePath("SpikeEventSeries2D.h5"); + std::shared_ptr io = createIO("HDF5", path); + io->open(); + io->createGroup("/general"); + io->createGroup("/general/extracellular_ephys"); + + // setup electrode table, device, and electrode group + NWB::ElectrodeTable elecTable = NWB::ElectrodeTable(io); + elecTable.initialize(); + + // setup electrical series + NWB::SpikeEventSeries ses = NWB::SpikeEventSeries(dataPath, io); + ses.initialize(dataType, + mockArrays[0], + "no description", + SizeArray {0, numSamples}, + SizeArray {8, 1}); + + // write channel data + for (SizeType e = 0; e < numEvents; ++e) { + double timestamp = mockTimestamps[e]; + ses.writeSpike(numSamples, 1, mockData[e].data(), ×tamp); + } + io->close(); + + // Read data back from file + std::unique_ptr file = + std::make_unique(path, H5F_ACC_RDONLY); + std::unique_ptr dataset = + std::make_unique(file->openDataSet(dataPath + "/data")); + std::vector> dataOut(numEvents, + std::vector(numSamples)); + float* buffer = new float[numEvents * numSamples]; + + H5::DataSpace fSpace = dataset->getSpace(); + hsize_t dims[3]; + fSpace.getSimpleExtentDims(dims, NULL); + hsize_t memdims = dims[0] * dims[1] * dims[2]; + dataset->read(buffer, H5::PredType::NATIVE_FLOAT, fSpace, fSpace); + + for (SizeType i = 0; i < numEvents; ++i) { + for (SizeType j = 0; j < (numSamples); ++j) { + dataOut[i][j] = buffer[i * (numSamples) + j]; + } + } + delete[] buffer; + REQUIRE_THAT(dataOut[0], Catch::Matchers::Approx(mockData[0]).margin(1)); + REQUIRE_THAT(dataOut[1], Catch::Matchers::Approx(mockData[1]).margin(1)); + } +} \ No newline at end of file diff --git a/tests/testNWBFile.cpp b/tests/testNWBFile.cpp index 68b2cf18..bcf9123c 100644 --- a/tests/testNWBFile.cpp +++ b/tests/testNWBFile.cpp @@ -6,6 +6,7 @@ #include "nwb/NWBFile.hpp" #include "nwb/RecordingContainers.hpp" #include "nwb/base/TimeSeries.hpp" +#include "nwb/ecephys/SpikeEventSeries.hpp" #include "testUtils.hpp" using namespace AQNWB; @@ -15,7 +16,8 @@ TEST_CASE("saveNWBFile", "[nwb]") std::string filename = getTestFilePath("testSaveNWBFile.nwb"); // initialize nwbfile object and create base structure - NWB::NWBFile nwbfile(std::make_unique(filename)); + auto io = std::make_shared(filename); + NWB::NWBFile nwbfile(io); nwbfile.initialize(generateUuid()); nwbfile.finalize(); } @@ -32,10 +34,15 @@ TEST_CASE("createElectricalSeries", "[nwb]") // create Electrical Series std::vector mockArrays = getMockChannelArrays(1, 2); + std::vector mockChannelNames = + getMockChannelArrayNames("esdata"); std::unique_ptr recordingContainers = std::make_unique(); - Status resultCreate = nwbfile.createElectricalSeries( - mockArrays, BaseDataType::F32, recordingContainers.get()); + Status resultCreate = + nwbfile.createElectricalSeries(mockArrays, + mockChannelNames, + BaseDataType::F32, + recordingContainers.get()); REQUIRE(resultCreate == Status::Success); // start recording @@ -60,6 +67,71 @@ TEST_CASE("createElectricalSeries", "[nwb]") nwbfile.finalize(); } +TEST_CASE("createMultipleEcephysDatasets", "[nwb]") +{ + std::string filename = getTestFilePath("createESandSES.nwb"); + + // initialize nwbfile object and create base structure + std::shared_ptr io = std::make_shared(filename); + NWB::NWBFile nwbfile(io); + nwbfile.initialize(generateUuid()); + + // create Electrical Series + std::vector mockArrays = getMockChannelArrays(1, 2); + std::vector mockChannelNames = + getMockChannelArrayNames("esdata"); + std::unique_ptr recordingContainers = + std::make_unique(); + Status resultCreateES = + nwbfile.createElectricalSeries(mockArrays, + mockChannelNames, + BaseDataType::F32, + recordingContainers.get()); + REQUIRE(resultCreateES == Status::Success); + + // create SpikeEventSeries + SizeType numSamples = 5; + std::vector mockSpikeChannelNames = + getMockChannelArrayNames("spikedata"); + Status resultCreateSES = + nwbfile.createSpikeEventSeries(mockArrays, + mockSpikeChannelNames, + BaseDataType::F32, + recordingContainers.get()); + + // start recording + Status resultStart = io->startRecording(); + REQUIRE(resultStart == Status::Success); + + // write electrical series data + std::vector mockData = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f}; + std::vector mockTimestamps = {0.1, 0.2, 0.3, 0.4, 0.5}; + std::vector positionOffset = {0, 0}; + std::vector dataShape = {mockData.size(), 0}; + + NWB::TimeSeries* ts0 = + static_cast(recordingContainers->getContainer(0)); + ts0->writeData( + dataShape, positionOffset, mockData.data(), mockTimestamps.data()); + NWB::TimeSeries* ts1 = + static_cast(recordingContainers->getContainer(1)); + ts1->writeData( + dataShape, positionOffset, mockData.data(), mockTimestamps.data()); + + // write spike event series data + SizeType numEvents = 10; + NWB::SpikeEventSeries* ses0 = + static_cast(recordingContainers->getContainer(2)); + NWB::SpikeEventSeries* ses1 = + static_cast(recordingContainers->getContainer(3)); + for (SizeType i = 0; i < numEvents; ++i) { + ses0->writeSpike(numSamples, 1, mockData.data(), &mockTimestamps[0]); + ses1->writeSpike(numSamples, 1, mockData.data(), &mockTimestamps[0]); + } + + nwbfile.finalize(); +} + TEST_CASE("setCanModifyObjectsMode", "[nwb]") { std::string filename = getTestFilePath("testCanModifyObjectsMode.nwb"); @@ -80,8 +152,10 @@ TEST_CASE("setCanModifyObjectsMode", "[nwb]") // test that dataset creation fails after starting the recording std::vector mockArrays = getMockChannelArrays(1, 2); - Status resultCreatePostStart = - nwbfile.createElectricalSeries(mockArrays, BaseDataType::F32); + std::vector mockChannelNames = + getMockChannelArrayNames("esdata"); + Status resultCreatePostStart = nwbfile.createElectricalSeries( + mockArrays, mockChannelNames, BaseDataType::F32); REQUIRE(resultCreatePostStart == Status::Failure); // stop recording diff --git a/tests/testRecordingWorkflow.cpp b/tests/testRecordingWorkflow.cpp index 9651dc88..21142928 100644 --- a/tests/testRecordingWorkflow.cpp +++ b/tests/testRecordingWorkflow.cpp @@ -29,6 +29,8 @@ TEST_CASE("writeContinuousData", "[recording]") std::vector mockRecordingArrays = getMockChannelArrays(); + std::vector mockChannelNames = + getMockChannelArrayNames("esdata"); std::vector> mockData = getMockData2D(numSamples, numChannels); std::vector mockTimestamps = getMockTimestamps(numSamples); @@ -46,8 +48,10 @@ TEST_CASE("writeContinuousData", "[recording]") nwbfile->initialize(generateUuid()); // 4. create datasets and add to recording containers - nwbfile->createElectricalSeries( - mockRecordingArrays, BaseDataType::F32, recordingContainers.get()); + nwbfile->createElectricalSeries(mockRecordingArrays, + mockChannelNames, + BaseDataType::F32, + recordingContainers.get()); // 5. start the recording io->startRecording(); @@ -93,7 +97,7 @@ TEST_CASE("writeContinuousData", "[recording]") nwbfile->finalize(); // check contents of data - std::string dataPath = "/acquisition/array0/data"; + std::string dataPath = "/acquisition/esdata0/data"; std::unique_ptr file = std::make_unique(path, H5F_ACC_RDONLY); std::unique_ptr dataset = @@ -119,7 +123,7 @@ TEST_CASE("writeContinuousData", "[recording]") REQUIRE_THAT(dataOut[1], Catch::Matchers::Approx(mockData[1]).margin(1)); // check contents of timestamps - std::string timestampsPath = "/acquisition/array0/timestamps"; + std::string timestampsPath = "/acquisition/esdata0/timestamps"; std::unique_ptr tsDataset = std::make_unique(file->openDataSet(timestampsPath)); double* tsBuffer = new double[numSamples]; diff --git a/tests/testUtils.hpp b/tests/testUtils.hpp index 06778827..66c9c9ef 100644 --- a/tests/testUtils.hpp +++ b/tests/testUtils.hpp @@ -34,14 +34,16 @@ inline std::string getTestFilePath(std::string filename) } inline std::vector getMockChannelArrays( - SizeType numChannels = 2, SizeType numArrays = 2) + SizeType numChannels = 2, + SizeType numArrays = 2, + std::string groupName = "array") { std::vector arrays(numArrays); for (SizeType i = 0; i < numArrays; i++) { std::vector chGroup; for (SizeType j = 0; j < numChannels; j++) { Channel ch("ch" + std::to_string(j), - "array" + std::to_string(i), + groupName + std::to_string(i), i, j, i * numArrays + j); @@ -52,6 +54,16 @@ inline std::vector getMockChannelArrays( return arrays; } +inline std::vector getMockChannelArrayNames( + std::string baseName = "esdata", SizeType numArrays = 2) +{ + std::vector arrayNames(numArrays); + for (SizeType i = 0; i < numArrays; i++) { + arrayNames[i] = baseName + std::to_string(i); + } + return arrayNames; +} + inline std::vector getMockData1D(SizeType numSamples = 1000) { std::vector mockData(numSamples);