sx_rdata_

dweindl · Sep 25, 2024 · 7077a81 · 7077a81
1 parent cc17c5e
commit 7077a81
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 21 deletions.
diff --git a/include/amici/model.h b/include/amici/model.h
@@ -1432,14 +1432,15 @@ class Model : public AbstractModel, public ModelDimensions {
     /**
      * @brief Expand conservation law for state sensitivities.
      * @param sx_rdata Output buffer for state variables sensitivities with
-     * conservation laws expanded (stored in `amici::ReturnData`).
+     * conservation laws expanded
+     * (stored in `amici::ReturnData` shape `nplist` x `nx`, row-major).
      * @param sx_solver State variables sensitivities with conservation laws
      * applied (solver returns this)
      * @param x_solver State variables with conservation laws
      * applied (solver returns this)
      */
     void fsx_rdata(
-        AmiVectorArray& sx_rdata, AmiVectorArray const& sx_solver,
+        gsl::span<realtype> sx_rdata, AmiVectorArray const& sx_solver,
         AmiVector const& x_solver
     );
 

diff --git a/include/amici/rdata.h b/include/amici/rdata.h
@@ -463,10 +463,6 @@ class ReturnData : public ModelDimensions {
     /** offset for sigma_residuals */
     realtype sigma_offset;
 
-    /** full sensitivity state vector array, including states eliminated from
-     * conservation laws */
-    AmiVectorArray sx_rdata_;
-
     /** array of number of found roots for a certain event type
      * (shape `ne`) */
     std::vector<int> nroots_;

diff --git a/src/model.cpp b/src/model.cpp
@@ -1998,17 +1998,17 @@ void Model::fx_rdata(gsl::span<realtype> x_rdata, AmiVector const& x) {
 }
 
 void Model::fsx_rdata(
-    AmiVectorArray& sx_rdata, AmiVectorArray const& sx,
+    gsl::span<realtype> sx_rdata, AmiVectorArray const& sx,
     AmiVector const& x_solver
 ) {
     realtype* stcl = nullptr;
     for (int ip = 0; ip < nplist(); ip++) {
         if (ncl() > 0)
             stcl = &state_.stotal_cl.at(plist(ip) * ncl());
         fsx_rdata(
-            sx_rdata.data(ip), sx.data(ip), stcl,
-            state_.unscaledParameters.data(), state_.fixedParameters.data(),
-            x_solver.data(), state_.total_cl.data(), plist(ip)
+            &sx_rdata[ip], sx.data(ip), stcl, state_.unscaledParameters.data(),
+            state_.fixedParameters.data(), x_solver.data(),
+            state_.total_cl.data(), plist(ip)
         );
     }
 }

diff --git a/src/rdata.cpp b/src/rdata.cpp
@@ -46,7 +46,6 @@ ReturnData::ReturnData(
     , rdata_reporting(rdrm)
     , sigma_res(sigma_res)
     , sigma_offset(sigma_offset)
-    , sx_rdata_(nx, nplist)
     , nroots_(ne) {
     switch (rdata_reporting) {
     case RDataReporting::full:
@@ -198,9 +197,7 @@ void ReturnData::processPreEquilibration(
         model.fx_rdata(x_ss, simulation_state.x);
     }
     if (!sx_ss.empty() && sensi >= SensitivityOrder::first) {
-        model.fsx_rdata(sx_rdata_, simulation_state.sx, simulation_state.x);
-        for (int ip = 0; ip < nplist; ip++)
-            writeSlice(sx_rdata_[ip], slice(sx_ss, ip, nx));
+        model.fsx_rdata(sx_ss, simulation_state.sx, simulation_state.x);
     }
     /* Get cpu time for pre-equilibration in milliseconds */
     preeq_cpu_time = preeq.getCPUTime();
@@ -254,9 +251,7 @@ void ReturnData::processForwardProblem(
     }
 
     if (!sx0.empty()) {
-        model.fsx_rdata(sx_rdata_, initialState.sx, initialState.x);
-        for (int ip = 0; ip < nplist; ip++)
-            writeSlice(sx_rdata_[ip], slice(sx0, ip, nx));
+        model.fsx_rdata(sx0, initialState.sx, initialState.x);
     }
 
     // process timepoint data
@@ -331,10 +326,7 @@ void ReturnData::getDataSensisFSA(
     ExpData const* edata
 ) {
     if (!sx.empty()) {
-        model.fsx_rdata(sx_rdata_, simulation_state.sx, simulation_state.x);
-        for (int ip = 0; ip < nplist; ip++) {
-            writeSlice(sx_rdata_[ip], slice(sx, it * nplist + ip, nx));
-        }
+        model.fsx_rdata(sx, simulation_state.sx, simulation_state.x);
     }
 
     if (!sy.empty()) {