Add id field to ExpData and ReturnData (#1622)

A ReturnData will have the same id as the ExpData used for simulation. Allows for more informative output and makes it easier to link ReturnDatas to specific simulation conditions.

include/amici/edata.h

@@ -385,6 +385,11 @@ class ExpData : public SimulationParameters {
      */
     const realtype *getObservedEventsStdDevPtr(int ie) const;
 
+    /**
+     * @brief Arbitrary (not necessarily unique) identifier.
+     */
+    std::string id;
+
   protected:
     /**
      * @brief resizes observedData, observedDataStdDev, observedEvents and

include/amici/hdf5.h

@@ -240,6 +240,17 @@ void createAndWriteDouble3DDataset(H5::H5File const &file,
                                    gsl::span<const double> buffer, hsize_t m,
                                    hsize_t n, hsize_t o);
 
+/**
+ * @brief Read string attribute from HDF5 object.
+ * @param file HDF5 file
+ * @param optionsObject Object to read attribute from
+ * @param attributeName Name of attribute to read
+ * @return Attribute value
+ */
+std::string getStringAttribute(H5::H5File const& file,
+                               std::string const& optionsObject,
+                               std::string const& attributeName);
+
 /**
  * @brief Read scalar native double attribute from HDF5 object.
  * @param file HDF5 file

include/amici/rdata.h

@@ -92,6 +92,10 @@ class ReturnData: public ModelDimensions {
                                   SteadystateProblem const *posteq,
                                   Model &model, Solver const &solver,
                                   ExpData const *edata);
+    /**
+     * @brief Arbitrary (not necessarily unique) identifier.
+     */
+    std::string id;
 
     /**
      * timepoints (shape `nt`)

include/amici/serialization.h

@@ -170,6 +170,7 @@ void serialize(Archive &ar, amici::SimulationParameters &s, const unsigned int /
 template <class Archive>
 void serialize(Archive &ar, amici::ReturnData &r, const unsigned int /*version*/) {
     ar &dynamic_cast<amici::ModelDimensions&>(r);
+    ar &r.id;
     ar &r.nx;
     ar &r.nxtrue;
     ar &r.nplist;

python/amici/numpy.py

@@ -48,6 +48,9 @@ def __getitem__(self, item: str) -> Union[np.ndarray, float]:
         if item in self._cache:
             return self._cache[item]
 
+        if item == 'id':
+            return getattr(self._swigptr, item)
+
         if item not in self._field_names:
             self.__missing__(item)
 

python/amici/petab_objective.py

@@ -528,7 +528,9 @@ def create_edata_for_condition(
 
     # create an ExpData object
     edata = amici.ExpData(amici_model)
-
+    edata.id = condition[SIMULATION_CONDITION_ID]
+    if condition.get(PREEQUILIBRATION_CONDITION_ID):
+        edata.id += "+" + condition.get(PREEQUILIBRATION_CONDITION_ID)
     ##########################################################################
     # enable initial parameters reinitialization
     species_in_condition_table = [

python/tests/test_sbml_import.py

@@ -160,8 +160,12 @@ def test_steadystate_simulation(model_steadystate_module):
     solver.setSensitivityOrder(amici.SensitivityOrder.first)
     rdata = amici.runAmiciSimulation(model, solver)
     edata = [amici.ExpData(rdata, 1, 0)]
+    edata[0].id = "some condition ID"
     rdata = amici.runAmiciSimulations(model, solver, edata)
 
+    assert rdata[0].status == amici.AMICI_SUCCESS
+    assert rdata[0].id == edata[0].id
+
     # check roundtripping of DataFrame conversion
     df_edata = amici.getDataObservablesAsDataFrame(model, edata)
     edata_reconstructed = amici.getEdataFromDataFrame(model, df_edata)

src/amici.cpp

@@ -112,6 +112,9 @@ AmiciApplication::runAmiciSimulation(Solver& solver,
 
     std::unique_ptr<ReturnData> rdata = std::make_unique<ReturnData>(solver,
                                                                      model);
+    if(edata) {
+        rdata->id = edata->id;
+    }
 
     std::unique_ptr<SteadystateProblem> preeq {};
     std::unique_ptr<ForwardProblem> fwd {};

src/edata.cpp

@@ -98,7 +98,9 @@ ExpData::ExpData(ReturnData const& rdata, std::vector<realtype> sigma_y,
             observed_events_.at(iz + rdata.nztrue * ie) = rdata.z.at(iz + rdata.nz * ie) + e(gen);
             observed_data_std_dev_.at(iz + rdata.nztrue * ie) = sigma;
         }
-        }
+    }
+
+    id = rdata.id;
 }
 
 void ExpData::setTimepoints(const std::vector<realtype> &ts) {

src/hdf5.cpp

@@ -112,6 +112,10 @@ std::unique_ptr<ExpData> readSimulationExpData(std::string const& hdf5Filename,
 
     auto edata = std::unique_ptr<ExpData>(new ExpData(model));
 
+    if(attributeExists(file, hdf5Root, "id")) {
+        edata->id = getStringAttribute(file, hdf5Root, "id");
+    }
+
     if (model.ny * model.nt() > 0) {
         if(locationExists(file,  hdf5Root + "/Y")) {
             auto my = getDoubleDataset2D(file, hdf5Root + "/Y", m, n);
@@ -192,6 +196,10 @@ void writeSimulationExpData(const ExpData &edata, H5::H5File const& file,
     if(!locationExists(file, hdf5Location))
         createGroup(file, hdf5Location);
 
+    H5LTset_attribute_string(file.getId(), hdf5Location.c_str(), "id",
+                             edata.id.c_str());
+
+
     if (edata.nt())
         createAndWriteDouble1DDataset(file, hdf5Location + "/ts",
                                       edata.getTimepoints());
@@ -245,6 +253,9 @@ void writeReturnData(const ReturnData &rdata, H5::H5File const& file, const std:
     if (!rdata.ts.empty())
         createAndWriteDouble1DDataset(file, hdf5Location + "/t", rdata.ts);
 
+    H5LTset_attribute_string(file.getId(), hdf5Location.c_str(), "id",
+                             rdata.id.c_str());
+
     H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
                              "llh", &rdata.llh, 1);
     H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
@@ -461,6 +472,34 @@ void writeReturnData(ReturnData const& rdata,
     writeReturnData(rdata, file, hdf5Location);
 }
 
+std::string getStringAttribute(H5::H5File const& file,
+                               std::string const& optionsObject,
+                               std::string const& attributeName) {
+    hsize_t dims;
+    H5T_class_t type_class;
+    size_t type_size;
+    auto status = H5LTget_attribute_info(file.getId(), optionsObject.c_str(),
+                                         attributeName.c_str(), &dims,
+                                         &type_class,&type_size);
+    if(status < 0) {
+        throw AmiException("Could get info for attribute %s for object %s.",
+                           attributeName.c_str(), optionsObject.c_str());
+    }
+    std::vector<char> value(type_size);
+    status = H5LTget_attribute_string(file.getId(), optionsObject.c_str(),
+                                      attributeName.c_str(), value.data());
+
+#ifdef AMI_HDF5_H_DEBUG
+    printf("%s: %s\n", attributeName.c_str(), value.data());
+#endif
+
+    if(status < 0)
+        throw AmiException("Attribute %s not found for object %s.",
+                           attributeName.c_str(), optionsObject.c_str());
+
+    return std::string(value.data());
+}
+
 double getDoubleScalarAttribute(H5::H5File const& file,
                                 std::string const& optionsObject,
                                 std::string const& attributeName) {

tests/cpp/unittests/testSerialization.cpp

@@ -11,6 +11,7 @@
 void
 checkReturnDataEqual(amici::ReturnData const& r, amici::ReturnData const& s)
 {
+    ASSERT_EQ(r.id, s.id);
     ASSERT_EQ(r.np, s.np);
     ASSERT_EQ(r.nk, s.nk);
     ASSERT_EQ(r.nx, s.nx);
@@ -227,7 +228,7 @@ TEST(ReturnDataSerializationTest, ToString)
         std::vector<int>(nz, 0));
 
     amici::ReturnData r(solver, m);
-
+    r.id = "some_id";
     std::string serialized = amici::serializeToString(r);
 
     checkReturnDataEqual(