Harden SED-ML uniform time-course execution

VERSION.txt

@@ -1 +1 @@
-1.20261009.1
+1.20261009.2

src/api/libopencor/seduniformtimecourse.h

@@ -30,6 +30,7 @@ namespace libOpenCOR {
 class LIBOPENCOR_EXPORT SedUniformTimeCourse: public SedSimulation
 {
     friend class SedInstanceTask;
+    friend class SedSimulation;
 
 public:
     /**

src/sed/sedinstance.cpp

@@ -20,6 +20,7 @@ limitations under the License.
 
 #include "libopencor/seddocument.h"
 
+#include <atomic>
 #include <chrono>
 #include <exception>
 #include <memory>
@@ -200,13 +201,32 @@ bool SedInstance::Impl::startRun()
 
     mRunning.store(true, std::memory_order_release);
 
-    mRunFuture = std::async(std::launch::async, [this]() {
-        const auto result = run();
+    // Start our run in a separate thread.
+    // Note #1: we must be flagged as not running anymore once our run is done, even if run() throws an exception (it
+    //          reports a failure as an issue, but it might still throw, e.g., std::bad_alloc when restoring our
+    //          issues), hence we use a guard to do so. Otherwise, we would be stuck in RUNNING.
+    // Note #2: std::async() may throw an exception (e.g., std::system_error if no thread could be created), in which
+    //          case we must also be flagged as not running anymore. We have no way to trigger such an exception in our
+    //          tests, hence we ignore our try...catch statement during code coverage.
 
+#ifndef CODE_COVERAGE_ENABLED
+    try {
+#endif
+        mRunFuture = std::async(std::launch::async, [this]() {
+            auto resetRunning = [](std::atomic<bool> *pRunning) {
+                pRunning->store(false, std::memory_order_release);
+            };
+            const std::unique_ptr<std::atomic<bool>, decltype(resetRunning)> runningGuard {&mRunning, resetRunning};
+
+            return run();
+        });
+#ifndef CODE_COVERAGE_ENABLED
+    } catch (...) {
         mRunning.store(false, std::memory_order_release);
 
-        return result;
-    });
+        throw;
+    }
+#endif
 
     return true;
 }

src/sed/sedinstancetask.cpp

@@ -27,8 +27,10 @@ limitations under the License.
 #include <cmath>
 #include <cstdint>
 #include <condition_variable>
+#include <limits>
 #include <memory>
 #include <mutex>
+#include <stdexcept>
 #include <unordered_set>
 
 namespace libOpenCOR {
@@ -434,14 +436,17 @@ void SedInstanceTask::Impl::run(double pVoiStart, double pVoiEnd, double pVoiInt
     };
 
     // Compute the differential model.
+    // Note: when tracking our results, we take exactly as many steps as our results can hold (see run() below).
 
     auto *odeSolverPimpl {mOdeSolver->pimpl()};
+    const auto lastVoiCounter {pTrackResults ? mResults.resultsSize - 1 : 0};
     size_t voiCounter {0};
+    auto voiEndReached {pTrackResults ? (lastVoiCounter == 0) : fuzzyCompare(mVoi, pVoiEnd)};
 
     const auto computeRates = mRuntime->computeRates();
     const auto computeVariablesForDifferentialModel = mRuntime->computeVariablesForDifferentialModel();
 
-    while (!fuzzyCompare(mVoi, pVoiEnd)) {
+    while (!voiEndReached) {
         // Check whether a pause or stop has been requested.
 
         const auto runControl = mRunControl->load(std::memory_order_relaxed);
@@ -468,9 +473,19 @@ void SedInstanceTask::Impl::run(double pVoiStart, double pVoiEnd, double pVoiInt
             return;
         }
 
+        // Determine the point that we want to reach.
+        // Note: our last point is always exactly pVoiEnd. Indeed, pVoiStart + n * pVoiInterval may differ slightly from
+        //       pVoiEnd due to rounding errors, which would otherwise result in an extra step being taken and, when
+        //       tracking our results, its results being written past the end of our results arrays.
+
+        auto voiTarget {pVoiStart + static_cast<double>(++voiCounter) * pVoiInterval};
+
+        voiEndReached = pTrackResults ? (voiCounter == lastVoiCounter) : (voiTarget >= pVoiEnd);
+        voiTarget = voiEndReached ? pVoiEnd : std::min(voiTarget, pVoiEnd);
+
         // Update our model's state.
 
-        if (!odeSolverPimpl->solve(mVoi, std::min(pVoiStart + static_cast<double>(++voiCounter) * pVoiInterval, pVoiEnd))) {
+        if (!odeSolverPimpl->solve(mVoi, voiTarget)) {
             // Note: an NLA system that could not be solved is a likely reason for our ODE solver to have failed (see
             //       rhsFunction() in solvercvode.cpp).
 
@@ -526,25 +541,36 @@ double SedInstanceTask::Impl::run()
     auto startTime {std::chrono::high_resolution_clock::now()};
 
     // Reset our progress counters.
-    // Note: we retrieve our number of steps only once since it is used in several places below.
-
-    const auto *sedUniformTimeCoursePimpl {mDifferentialModel ? mSedUniformTimeCourse->pimpl() : nullptr};
-    const auto numberOfSteps {mDifferentialModel ? sedUniformTimeCoursePimpl->mNumberOfSteps : 1};
-    const auto totalSteps {(numberOfSteps > 0) ? static_cast<size_t>(numberOfSteps) : 0};
 
     mCompletedSteps.store(0, std::memory_order_relaxed);
-    mTotalSteps.store(totalSteps, std::memory_order_relaxed);
+    mTotalSteps.store(0, std::memory_order_relaxed);
 
-    // Make sure that our number of steps is valid.
-    // Note: our number of steps was validated when we were instantiated (see SedSimulation::Impl::isValid()), but it
-    //       may have been changed since then.
+    // Make sure that our simulation is valid.
+    // Note: our simulation was validated when we were instantiated (see SedSimulation::Impl::isValid()), but its times
+    //       and/or number of steps may have been changed since then.
 
-    if (numberOfSteps <= 0) {
-        addIssue(Issue::Type::ERROR, sedUniformTimeCoursePimpl->invalidNumberOfStepsError(mModel, numberOfSteps), "Simulation");
+    const auto *sedUniformTimeCoursePimpl {mDifferentialModel ? mSedUniformTimeCourse->pimpl() : nullptr};
 
-        return 0.0;
+    if (mDifferentialModel) {
+        const auto errors {sedUniformTimeCoursePimpl->validationErrors(mModel)};
+
+        if (!errors.empty()) {
+            for (const auto &error : errors) {
+                addIssue(Issue::Type::ERROR, error, "Simulation");
+            }
+
+            return 0.0;
+        }
     }
 
+    // Set our total number of steps.
+    // Note: we retrieve our number of steps only once since it is used in several places below.
+
+    const auto numberOfSteps {mDifferentialModel ? sedUniformTimeCoursePimpl->mNumberOfSteps : 1};
+    const auto totalSteps {static_cast<size_t>(numberOfSteps)};
+
+    mTotalSteps.store(totalSteps, std::memory_order_relaxed);
+
     // (Re)initialise our model.
     // Note: reinitialise our model because we initialised it when we created the instance task.
 
@@ -566,29 +592,48 @@ double SedInstanceTask::Impl::run()
             }
         }
 
-        // Initialise our results structure.
-        // Note: we release our previous results first (to limit our peak memory usage) and only then allocate our new
-        //       results, which we do into a local structure that we move into ours once all of its arrays have been
-        //       allocated. This means that should an allocation fail (e.g., std::bad_alloc), our results structure
-        //       would be left empty rather than with an inconsistent results size and arrays of mismatched sizes (which
-        //       would result in out-of-bounds accesses when retrieving our results).
+        // Initialise our results structure, if needed.
+        // Note #1: we only (re)allocate our results if their size has changed. This avoids needless allocations, but
+        //          it also means that our results remain where they are in memory, which matters since our Python
+        //          bindings return zero-copy NumPy arrays.
+        // Note #2: when (re)allocating our results, we release our previous results first (to limit our peak memory
+        //          usage) and only then allocate our new results, which we do into a local structure that we move into
+        //          ours once all of its arrays have been allocated. This means that should an allocation fail (e.g.,
+        //          std::bad_alloc), our results structure would be left empty rather than with an inconsistent results
+        //          size and arrays of mismatched sizes (which would result in out-of-bounds accesses when retrieving
+        //          our results).
+        // Note #3: on a 32-bit platform (e.g., WebAssembly), the size of an array may overflow (e.g., with many
+        //          constants and a large number of steps), in which case we would allocate an array that is too small
+        //          for our results, hence we check for it. We cannot trigger such an overflow on a 64-bit platform, so
+        //          we ignore our check during code coverage.
 
         const auto resultsSize {totalSteps + 1};
 
-        mResults = {};
+        if (mResults.resultsSize != resultsSize) {
+            mResults = {};
+
+            auto arraySize = [resultsSize](size_t pCount) {
+#ifndef CODE_COVERAGE_ENABLED
+                if ((pCount != 0) && (resultsSize > (std::numeric_limits<size_t>::max() / pCount))) {
+                    throw std::length_error("the results are too large to be allocated");
+                }
+#endif
 
-        SedInstanceTaskResults results;
+                return pCount * resultsSize;
+            };
+            SedInstanceTaskResults results;
 
-        results.resultsSize = resultsSize;
+            results.resultsSize = resultsSize;
 
-        results.voi.resize(resultsSize);
-        results.states.resize(mStateCount * resultsSize);
-        results.rates.resize(mStateCount * resultsSize);
-        results.constants.resize(mConstantCount * resultsSize);
-        results.computedConstants.resize(mComputedConstantCount * resultsSize);
-        results.algebraicVariables.resize(mAlgebraicVariableCount * resultsSize);
+            results.voi.resize(resultsSize);
+            results.states.resize(arraySize(mStateCount));
+            results.rates.resize(arraySize(mStateCount));
+            results.constants.resize(arraySize(mConstantCount));
+            results.computedConstants.resize(arraySize(mComputedConstantCount));
+            results.algebraicVariables.resize(arraySize(mAlgebraicVariableCount));
 
-        mResults = std::move(results);
+            mResults = std::move(results);
+        }
 
         // Run our simulation from the output start time to the output end time, tracking our results.
 

src/sed/sedsimulation.cpp

@@ -22,8 +22,6 @@ limitations under the License.
 #include "solvernla_p.h"
 #include "solverode_p.h"
 
-#include "utils.h"
-
 namespace libOpenCOR {
 
 static constexpr auto ID_PREFIX {"simulation"};
@@ -33,25 +31,6 @@ SedSimulation::Impl::Impl(const SedDocumentPtr &pDocument)
 {
 }
 
-std::string SedSimulation::Impl::invalidNumberOfStepsError(const SedModelPtr &pModel, int pNumberOfSteps) const
-{
-    const auto &modelId = pModel->pimpl()->mId;
-    const auto numberOfSteps {toString(pNumberOfSteps)};
-    std::string res;
-
-    res.reserve(mId.size() + modelId.size() + numberOfSteps.size() + 112); // NOLINT
-
-    res += "Simulation '";
-    res += mId;
-    res += "' is to be used with model '";
-    res += modelId;
-    res += "' which requires a strictly positive number of steps but ";
-    res += numberOfSteps;
-    res += " is provided.";
-
-    return res;
-}
-
 bool SedSimulation::Impl::isValid(const SedModelPtr &pModel, const SedUniformTimeCoursePtr &pUniformTimeCourse)
 {
     auto modelType {pModel->pimpl()->mFile->pimpl()->mCellmlFile->type()};
@@ -77,8 +56,7 @@ bool SedSimulation::Impl::isValid(const SedModelPtr &pModel, const SedUniformTim
 
     removeAllIssues();
 
-    // Make sure that we are a uniform time course simulation with a strictly positive number of steps if our model is
-    // an ODE/DAE model.
+    // Make sure that we are a valid uniform time course simulation if our model is an ODE/DAE model.
     //---GRY--- WE DON'T CURRENTLY SUPPORT OTHER TYPES OF SIMULATION FOR ODE/DAE MODELS.
 
     if ((modelType == libcellml::AnalyserModel::Type::ODE)
@@ -95,8 +73,10 @@ bool SedSimulation::Impl::isValid(const SedModelPtr &pModel, const SedUniformTim
             error += "' which (currently) requires a uniform time course simulation.";
 
             addError(error);
-        } else if (pUniformTimeCourse->numberOfSteps() <= 0) {
-            addError(invalidNumberOfStepsError(pModel, pUniformTimeCourse->numberOfSteps()));
+        } else {
+            for (const auto &error : pUniformTimeCourse->pimpl()->validationErrors(pModel)) {
+                addError(error);
+            }
         }
     }
 

src/sed/sedsimulation_p.h

@@ -30,8 +30,6 @@ class SedSimulation::Impl: public SedBase::Impl
 
     explicit Impl(const SedDocumentPtr &pDocument);
 
-    std::string invalidNumberOfStepsError(const SedModelPtr &pModel, int pNumberOfSteps) const;
-
     bool isValid(const SedModelPtr &pModel, const SedUniformTimeCoursePtr &pUniformTimeCourse);
 
     const SolverOdePtr &odeSolver() const;

src/sed/seduniformtimecourse.cpp

@@ -18,6 +18,9 @@ limitations under the License.
 
 #include "utils.h"
 
+#include "libopencor/sedmodel.h"
+
+#include <cmath>
 #include <memory>
 
 namespace libOpenCOR {
@@ -27,6 +30,61 @@ SedUniformTimeCourse::Impl::Impl(const SedDocumentPtr &pDocument)
 {
 }
 
+Strings SedUniformTimeCourse::Impl::validationErrors(const SedModelPtr &pModel) const
+{
+    // Make sure that our times are finite and such that initialTime <= outputStartTime < outputEndTime, and that our
+    // number of steps is strictly positive.
+    // Note: we are validated both when instantiating a document (see SedSimulation::Impl::isValid()) and when running
+    //       an instance task (see SedInstanceTask::Impl::run()) since our times and number of steps may have been
+    //       changed in between.
+
+    const auto &modelId = pModel->id();
+    Strings res;
+
+    if (!std::isfinite(mInitialTime) || !std::isfinite(mOutputStartTime) || !std::isfinite(mOutputEndTime)
+        || (mInitialTime > mOutputStartTime) || (mOutputStartTime >= mOutputEndTime)) {
+        const auto initialTime {toString(mInitialTime)};
+        const auto outputStartTime {toString(mOutputStartTime)};
+        const auto outputEndTime {toString(mOutputEndTime)};
+        std::string error;
+
+        error.reserve(mId.size() + modelId.size() + initialTime.size() + outputStartTime.size() + outputEndTime.size() + 152); // NOLINT
+
+        error += "Simulation '";
+        error += mId;
+        error += "' is to be used with model '";
+        error += modelId;
+        error += "' which requires finite times such that initialTime <= outputStartTime < outputEndTime but ";
+        error += initialTime;
+        error += ", ";
+        error += outputStartTime;
+        error += ", and ";
+        error += outputEndTime;
+        error += " are provided.";
+
+        res.push_back(error);
+    }
+
+    if (mNumberOfSteps <= 0) {
+        const auto numberOfSteps {toString(mNumberOfSteps)};
+        std::string error;
+
+        error.reserve(mId.size() + modelId.size() + numberOfSteps.size() + 112); // NOLINT
+
+        error += "Simulation '";
+        error += mId;
+        error += "' is to be used with model '";
+        error += modelId;
+        error += "' which requires a strictly positive number of steps but ";
+        error += numberOfSteps;
+        error += " is provided.";
+
+        res.push_back(error);
+    }
+
+    return res;
+}
+
 double SedUniformTimeCourse::Impl::initialTime() const
 {
     return mInitialTime;

src/sed/seduniformtimecourse_p.h

@@ -32,6 +32,8 @@ class SedUniformTimeCourse::Impl: public SedSimulation::Impl
 
     explicit Impl(const SedDocumentPtr &pDocument);
 
+    Strings validationErrors(const SedModelPtr &pModel) const;
+
     double initialTime() const;
     void setInitialTime(double pInitialTime);