Skip to content

Commit

Permalink
Add missing ReturnData fields in writeReturnData
Browse files Browse the repository at this point in the history 
Closes AMICI-dev#2591.
  • Loading branch information
@dweindl
dweindl committed Nov 26, 2024
1 parent c75909d commit 1430970
Show file tree
Hide file tree
Showing 3 changed files with 147 additions and 2 deletions.
12 changes: 12 additions & 0 deletions include/amici/hdf5.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ class ReturnData;
class ExpData;
class Model;
class Solver;
struct LogItem;

namespace hdf5 {

Expand Down Expand Up @@ -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<amici::LogItem> const& logItems,
std::string const& hdf5Location
);

/**
* @brief Create the given group and possibly parents.
* @param file HDF5 file to write to
Expand Down
4 changes: 2 additions & 2 deletions include/amici/rdata.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<LogItem> messages;
Expand All @@ -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`) */
Expand Down
133 changes: 133 additions & 0 deletions src/hdf5.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include <amici/hdf5.h>

#include <amici/edata.h>
#include <amici/logging.h>
#include <amici/model.h>
#include <amici/rdata.h>
#include <amici/solver.h>
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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> int_buffer(1);

int_buffer[0] = gsl::narrow<int>(rdata.newton_maxsteps);
H5LTset_attribute_int(
file.getId(), hdf5Location.c_str(), "newton_maxsteps",
int_buffer.data(), 1
);

int_buffer[0] = static_cast<int>(rdata.o2mode);
H5LTset_attribute_int(
file.getId(), hdf5Location.c_str(), "o2mode", int_buffer.data(), 1
);

int_buffer[0] = static_cast<int>(rdata.sensi);
H5LTset_attribute_int(
file.getId(), hdf5Location.c_str(), "sensi", int_buffer.data(), 1
);

int_buffer[0] = static_cast<int>(rdata.sensi_meth);
H5LTset_attribute_int(
file.getId(), hdf5Location.c_str(), "sensi_meth", int_buffer.data(), 1
);

int_buffer[0] = static_cast<int>(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<int>(rdata.pscale[i]);
createAndWriteInt1DDataset(file, hdf5Location + "/pscale", int_buffer);
}
writeLogItemsToHDF5(file, rdata.messages, hdf5Location + "/messages");

writeReturnDataDiagnosis(rdata, file, hdf5Location + "/diagnosis");
}

Expand Down Expand Up @@ -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<amici::LogItem> 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<LogItem> to std::vector<LogItemCStr>
std::vector<LogItemCStr> buffer(logItems.size());
for (size_t i = 0; i < logItems.size(); ++i) {
buffer[i].severity = static_cast<int>(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(
Expand Down Expand Up @@ -753,6 +883,7 @@ void createAndWriteDouble2DDataset(
const H5::H5File& file, std::string const& datasetName,
gsl::span<double const> 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(
Expand All @@ -765,6 +896,7 @@ void createAndWriteInt2DDataset(
H5::H5File const& file, std::string const& datasetName,
gsl::span<int const> 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(
Expand All @@ -777,6 +909,7 @@ void createAndWriteDouble3DDataset(
H5::H5File const& file, std::string const& datasetName,
gsl::span<double const> 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(
Expand Down

0 comments on commit 1430970

Please sign in to comment.