Timeseries

include/amici/amici.h

@@ -38,6 +38,20 @@ std::vector<std::unique_ptr<ReturnData>> runAmiciSimulations(
     Model const& model, bool failfast, int num_threads
 );
 
+/**
+ * @brief Same as runAmiciSimulation, but for a timeseries of ExpData instances.
+ *
+ * @param solver Solver instance
+ * @param edatas experimental data objects
+ * @param model model specification object
+ * @param failfast flag to allow early termination
+ * @return vector of pointers to return data objects
+ */
+std::vector<std::unique_ptr<ReturnData>> runAmiciSimulationsTimeseries(
+    Solver const& solver, std::vector<ExpData*> const& edatas,
+    Model const& model, bool failfast
+);
+
 /**
  * @brief Get the string representation of the given simulation status code
  * (see ReturnData::status).

include/amici/backwardproblem.h

@@ -70,6 +70,12 @@ class BackwardProblem {
      */
     AmiVector const& getAdjointState() const { return xB_; }
 
+    /**
+     * @brief Accessor for dxB
+     * @return xB
+     */
+    AmiVector const& getAdjointStateDerivative() const { return dxB_; }
+
     /**
      * @brief Accessor for xQB
      * @return xQB

include/amici/edata.h

@@ -398,10 +398,15 @@ class ExpData : public SimulationParameters {
     void clear_observations();
 
     /**
-     * @brief Arbitrary (not necessarily unique) identifier.
+     * @brief Unique identifier.
      */
     std::string id;
 
+    /**
+     * @brief Unique Identifier of the ExpData to use for initialization.
+     */
+    std::string initialization_id;
+
   protected:
     /**
      * @brief resizes observedData, observedDataStdDev, observedEvents and
@@ -440,6 +445,12 @@ class ExpData : public SimulationParameters {
         std::vector<realtype> const& input, char const* fieldname
     ) const;
 
+    /**
+     * @brief check whether this instance implements only equilibration
+     * without simulation
+     */
+    bool isEquilibration() const;
+
     /** @brief number of observables */
     int nytrue_{0};
 
@@ -480,6 +491,7 @@ inline bool operator==(ExpData const& lhs, ExpData const& rhs) {
     return *dynamic_cast<SimulationParameters const*>(&lhs)
                == *dynamic_cast<SimulationParameters const*>(&rhs)
            && lhs.id == rhs.id && lhs.nytrue_ == rhs.nytrue_
+           && lhs.initialization_id == rhs.initialization_id
            && lhs.nztrue_ == rhs.nztrue_ && lhs.nmaxevent_ == rhs.nmaxevent_
            && is_equal(lhs.observed_data_, rhs.observed_data_)
            && is_equal(lhs.observed_data_std_dev_, rhs.observed_data_std_dev_)

include/amici/rdata.h

@@ -95,10 +95,15 @@ class ReturnData : public ModelDimensions {
         Model& model, Solver const& solver, ExpData const* edata
     );
     /**
-     * @brief Arbitrary (not necessarily unique) identifier.
+     * @brief Unique identifier.
      */
     std::string id;
 
+    /**
+     * @brief Unique Identifier of the ExpData to use for initialization.
+     */
+    std::string initialization_id;
+
     /**
      * timepoints (shape `nt`)
      */

python/sdist/amici/swig_wrappers.py

@@ -73,10 +73,35 @@ def runAmiciSimulation(
     return numpy.ReturnDataView(rdata)
 
 
+class ExpDataTimeseries:
+    """
+    Prototype for timeseries data storage.
+    """
+
+    _exp_datas: list[amici.ExpData]
+
+    def __len__(self):
+        return len(self._exp_datas)
+
+    def __init__(self, exp_datas: list[amici.ExpData]):
+        self._exp_datas = exp_datas
+
+    def __getitem__(self, key):
+        return next(edata for edata in self._exp_datas if edata.id == key)
+
+    def __setitem__(self, key, value):
+        self._exp_datas = [
+            edata if edata.id != key else value for edata in self._exp_datas
+        ]
+
+    def get_ptr_vector(self):
+        return amici_swig.ExpDataPtrVector(self._exp_datas)
+
+
 def runAmiciSimulations(
     model: AmiciModel,
     solver: AmiciSolver,
-    edata_list: AmiciExpDataVector,
+    edata_list: AmiciExpDataVector | ExpDataTimeseries,
     failfast: bool = True,
     num_threads: int = 1,
 ) -> list["numpy.ReturnDataView"]:
@@ -105,14 +130,23 @@ def runAmiciSimulations(
             stacklevel=1,
         )
 
-    edata_ptr_vector = amici_swig.ExpDataPtrVector(edata_list)
-    rdata_ptr_list = amici_swig.runAmiciSimulations(
-        _get_ptr(solver),
-        edata_ptr_vector,
-        _get_ptr(model),
-        failfast,
-        num_threads,
-    )
+    if isinstance(edata_list, ExpDataTimeseries):
+        edata_ptr_vector = edata_list.get_ptr_vector()
+        rdata_ptr_list = amici_swig.runAmiciSimulationsTimeseries(
+            _get_ptr(solver),
+            edata_ptr_vector,
+            _get_ptr(model),
+            failfast,
+        )
+    else:
+        edata_ptr_vector = amici_swig.ExpDataPtrVector(edata_list)
+        rdata_ptr_list = amici_swig.runAmiciSimulations(
+            _get_ptr(solver),
+            edata_ptr_vector,
+            _get_ptr(model),
+            failfast,
+            num_threads,
+        )
     for rdata in rdata_ptr_list:
         _log_simulation(rdata)
         if solver.getReturnDataReportingMode() == amici.RDataReporting.full: