diff --git a/include/amici/hdf5.h b/include/amici/hdf5.h index d8acf8881e..cbb5d93f74 100644 --- a/include/amici/hdf5.h +++ b/include/amici/hdf5.h @@ -26,6 +26,7 @@ class ReturnData; class ExpData; class Model; class Solver; +struct LogItem; namespace hdf5 { @@ -137,6 +138,17 @@ void writeReturnDataDiagnosis( std::string const& hdf5Location ); +/** + * @brief Write log message to HDF5 file + * @param file HDF5 file to write to + * @param logItems Log items to write + * @param hdf5Location Full dataset path inside the HDF5 file (will be created) + */ +void writeLogItemsToHDF5( + H5::H5File const& file, std::vector const& logItems, + std::string const& hdf5Location +); + /** * @brief Create the given group and possibly parents. * @param file HDF5 file to write to diff --git a/include/amici/rdata.h b/include/amici/rdata.h index df25f39923..793be9435a 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -453,7 +453,7 @@ class ReturnData : public ModelDimensions { /** boolean indicating whether residuals for standard deviations have been * added */ - bool sigma_res; + bool sigma_res{false}; /** log messages */ std::vector messages; @@ -463,7 +463,7 @@ class ReturnData : public ModelDimensions { protected: /** offset for sigma_residuals */ - realtype sigma_offset; + realtype sigma_offset{0.0}; /** array of number of found roots for a certain event type * (shape `ne`) */ diff --git a/src/hdf5.cpp b/src/hdf5.cpp index 8ee20c0dbe..c8d4ec4b66 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -412,6 +413,11 @@ void writeReturnData( file, hdf5Location + "/y", rdata.y, rdata.nt, rdata.ny ); + if (!rdata.w.empty()) + createAndWriteDouble2DDataset( + file, hdf5Location + "/w", rdata.w, rdata.nt, rdata.nw + ); + if (!rdata.z.empty()) createAndWriteDouble2DDataset( file, hdf5Location + "/z", rdata.z, rdata.nmaxevent, rdata.nz @@ -480,6 +486,51 @@ void writeReturnData( rdata.nplist, rdata.nz ); + // TODO currently unused + /* + if (!rdata.s2rz.empty()) + createAndWriteDouble4DDataset( + file, hdf5Location + "/s2rz", rdata.s2rz, rdata.nmaxevent, + rdata.nztrue, rdata.nplist, rdata.nplist + ); + */ + + std::vector int_buffer(1); + + int_buffer[0] = gsl::narrow(rdata.newton_maxsteps); + H5LTset_attribute_int( + file.getId(), hdf5Location.c_str(), "newton_maxsteps", + int_buffer.data(), 1 + ); + + int_buffer[0] = static_cast(rdata.o2mode); + H5LTset_attribute_int( + file.getId(), hdf5Location.c_str(), "o2mode", int_buffer.data(), 1 + ); + + int_buffer[0] = static_cast(rdata.sensi); + H5LTset_attribute_int( + file.getId(), hdf5Location.c_str(), "sensi", int_buffer.data(), 1 + ); + + int_buffer[0] = static_cast(rdata.sensi_meth); + H5LTset_attribute_int( + file.getId(), hdf5Location.c_str(), "sensi_meth", int_buffer.data(), 1 + ); + + int_buffer[0] = static_cast(rdata.rdata_reporting); + H5LTset_attribute_int( + file.getId(), hdf5Location.c_str(), "rdrm", int_buffer.data(), 1 + ); + + if (!rdata.pscale.empty()) { + int_buffer.resize(rdata.pscale.size()); + for (int i = 0; (unsigned)i < rdata.pscale.size(); i++) + int_buffer[i] = static_cast(rdata.pscale[i]); + createAndWriteInt1DDataset(file, hdf5Location + "/pscale", int_buffer); + } + writeLogItemsToHDF5(file, rdata.messages, hdf5Location + "/messages"); + writeReturnDataDiagnosis(rdata, file, hdf5Location + "/diagnosis"); } @@ -634,6 +685,85 @@ void writeReturnDataDiagnosis( createAndWriteDouble2DDataset( file, hdf5Location + "/J", rdata.J, rdata.nx, rdata.nx ); + + if (!rdata.x_ss.empty()) + createAndWriteDouble1DDataset(file, hdf5Location + "/x_ss", rdata.x_ss); + + if (!rdata.sx_ss.empty()) + createAndWriteDouble2DDataset( + file, hdf5Location + "/sx_ss", rdata.sx_ss, rdata.nplist, + rdata.nx_rdata + ); +} + +// work-around for macos segfaults, use struct without std::string +struct LogItemCStr { + int severity; + const char* identifier; + const char* message; +}; + +void writeLogItemsToHDF5( + H5::H5File const& file, std::vector const& logItems, + std::string const& hdf5Location +) { + if (logItems.empty()) + return; + + try { + hsize_t dims[1] = {logItems.size()}; + H5::DataSpace dataspace(1, dims); + + // works on Ubuntu, but segfaults on macos: + /* + // Create a compound datatype for the LogItem struct. + H5::CompType logItemType(sizeof(amici::LogItem)); + logItemType.insertMember( + "severity", HOFFSET(amici::LogItem, severity), + H5::PredType::NATIVE_INT + ); + auto vlstr_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + logItemType.insertMember( + "identifier", HOFFSET(amici::LogItem, identifier), vlstr_type + ); + logItemType.insertMember( + "message", HOFFSET(amici::LogItem, message), vlstr_type + ); + H5::DataSet dataset + = file.createDataSet(hdf5Location, logItemType, dataspace); + + dataset.write(logItems.data(), logItemType); + */ + + // ... therefore, as a workaround, we use a struct without std::string + H5::CompType logItemType(sizeof(LogItemCStr)); + logItemType.insertMember( + "severity", HOFFSET(LogItemCStr, severity), + H5::PredType::NATIVE_INT + ); + auto vlstr_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + logItemType.insertMember( + "identifier", HOFFSET(LogItemCStr, identifier), vlstr_type + ); + logItemType.insertMember( + "message", HOFFSET(LogItemCStr, message), vlstr_type + ); + H5::DataSet dataset + = file.createDataSet(hdf5Location, logItemType, dataspace); + + // Convert std::vector to std::vector + std::vector buffer(logItems.size()); + for (size_t i = 0; i < logItems.size(); ++i) { + buffer[i].severity = static_cast(logItems[i].severity); + buffer[i].identifier = logItems[i].identifier.c_str(); + buffer[i].message = logItems[i].message.c_str(); + } + + // Write the data to the dataset. + dataset.write(buffer.data(), logItemType); + } catch (H5::Exception& e) { + throw AmiException(e.getCDetailMsg()); + } } void writeReturnData( @@ -753,6 +883,7 @@ void createAndWriteDouble2DDataset( const H5::H5File& file, std::string const& datasetName, gsl::span buffer, hsize_t m, hsize_t n ) { + Expects(buffer.size() == m * n); hsize_t const adims[]{m, n}; H5::DataSpace dataspace(2, adims); auto dataset = file.createDataSet( @@ -765,6 +896,7 @@ void createAndWriteInt2DDataset( H5::H5File const& file, std::string const& datasetName, gsl::span buffer, hsize_t m, hsize_t n ) { + Expects(buffer.size() == m * n); hsize_t const adims[]{m, n}; H5::DataSpace dataspace(2, adims); auto dataset = file.createDataSet( @@ -777,6 +909,7 @@ void createAndWriteDouble3DDataset( H5::H5File const& file, std::string const& datasetName, gsl::span buffer, hsize_t m, hsize_t n, hsize_t o ) { + Expects(buffer.size() == m * n * o); hsize_t const adims[]{m, n, o}; H5::DataSpace dataspace(3, adims); auto dataset = file.createDataSet(