From 60e8d1b4e067ffb381d0031fd76aaaa676fe971b Mon Sep 17 00:00:00 2001 From: Matthias Koenig Date: Mon, 2 Nov 2026 20:16:19 +0100 Subject: [PATCH] Evaluate the reaction rates once per output row of the selections Selecting a reaction rate evaluated all reaction rates of the model for that one value, so an output row with the reaction rates selected was quadratic in the number of reactions. getSelectedValues evaluates them once per row and reads the selected rates from them. --- source/rrRoadRunner.cpp | 31 ++++++++++++++++--- source/rrRoadRunner.h | 16 +++++++++- test/cxx_api_tests/RoadRunnerAPITests.h | 25 +++++++++++++++ .../RoadRunnerAPITestsWithLLJit.cpp | 4 +++ .../RoadRunnerAPITestsWithMCJit.cpp | 4 +++ 5 files changed, 75 insertions(+), 5 deletions(-) diff --git a/source/rrRoadRunner.cpp b/source/rrRoadRunner.cpp index e26ec034ef..33104c041d 100644 --- a/source/rrRoadRunner.cpp +++ b/source/rrRoadRunner.cpp @@ -229,6 +229,11 @@ namespace rr { std::vector mSelectionList; + /** + * the reaction rates of an output row, see getSelectedReactionRates + */ + std::vector selectedReactionRates; + /** * ModelGenerator obtained from the factory */ @@ -1410,16 +1415,31 @@ namespace rr { return dResult; } - double RoadRunner::getNthSelectedOutput(size_t index, double currentTime) { + double RoadRunner::getNthSelectedOutput(size_t index, double currentTime, + const double *reactionRates) { const SelectionRecord &record = impl->mSelectionList[index]; if (record.selectionType == SelectionRecord::TIME) { return currentTime; + } else if (reactionRates && record.selectionType == SelectionRecord::REACTION_RATE) { + return reactionRates[record.index]; } else { return getValue(record); } } + const double *RoadRunner::getSelectedReactionRates() { + for (const SelectionRecord &record: impl->mSelectionList) { + if (record.selectionType == SelectionRecord::REACTION_RATE) { + std::vector &rates = impl->selectedReactionRates; + rates.resize(impl->model->getNumReactions()); + impl->model->getReactionRates(rates.size(), nullptr, rates.data()); + return rates.data(); + } + } + return nullptr; + } + int RoadRunner::getTimeRowIndex() { for (u_int j = 0; j < impl->mSelectionList.size(); j++) { const SelectionRecord &record = impl->mSelectionList[j]; @@ -1432,8 +1452,9 @@ namespace rr { } void RoadRunner::getSelectedValues(ls::DoubleMatrix &results, int nRow, double currentTime) { + const double *reactionRates = getSelectedReactionRates(); for (u_int j = 0; j < impl->mSelectionList.size(); j++) { - double out = getNthSelectedOutput(j, currentTime); + double out = getNthSelectedOutput(j, currentTime, reactionRates); results(nRow, j) = out; } } @@ -1443,9 +1464,10 @@ namespace rr { assert(results.size() == impl->mSelectionList.size() && "given std::vector and selection list different size"); + const double *reactionRates = getSelectedReactionRates(); size_t size = results.size(); for (size_t i = 0; i < size; ++i) { - results[i] = getNthSelectedOutput(i, currentTime); + results[i] = getNthSelectedOutput(i, currentTime, reactionRates); } } @@ -5222,8 +5244,9 @@ namespace rr { std::vector result; result.resize(impl->mSelectionList.size()); + const double *reactionRates = getSelectedReactionRates(); for (int i = 0; i < impl->mSelectionList.size(); i++) { - result[i] = getNthSelectedOutput(i, impl->model->getTime()); + result[i] = getNthSelectedOutput(i, impl->model->getTime(), reactionRates); } return result; } diff --git a/source/rrRoadRunner.h b/source/rrRoadRunner.h index 3fb94a6ab4..a847aa00e9 100644 --- a/source/rrRoadRunner.h +++ b/source/rrRoadRunner.h @@ -1906,7 +1906,21 @@ namespace rr { bool populateResult(); - double getNthSelectedOutput(size_t index, double currentTime); + /** + * the value of the n'th selection, a reaction rate from reactionRates + * if given, see getSelectedReactionRates + */ + double getNthSelectedOutput(size_t index, double currentTime, + const double *reactionRates = nullptr); + + /** + * Evaluates the reaction rates once for the selections of an output + * row and returns them, or nullptr if no reaction rate is selected. + * The rate of one reaction is evaluated with all reaction rates of the + * model, so reading every selected rate on its own made an output row + * quadratic in the number of reactions. + */ + const double *getSelectedReactionRates(); bool isParameterUsed(const std::string &sid); diff --git a/test/cxx_api_tests/RoadRunnerAPITests.h b/test/cxx_api_tests/RoadRunnerAPITests.h index d79c8801b8..a6b91a88a9 100644 --- a/test/cxx_api_tests/RoadRunnerAPITests.h +++ b/test/cxx_api_tests/RoadRunnerAPITests.h @@ -588,6 +588,31 @@ class RoadRunnerAPITests : public RoadRunnerTest { delete testModel; } + void getSelectedReactionRates() { + // the reaction rates of every output row are the ones of its state, + // also selected twice and between other selections + TestModel *testModel = TestModelFactory("SimpleFlux"); + RoadRunner rr(testModel->str()); + rr.setSelections({"time", "_J1", "[S1]", "_J0", "[S2]", "_J1"}); + const ls::DoubleMatrix *result = rr.simulate(0, 10, 11); + ASSERT_EQ(result->numRows(), 11); + for (int i = 0; i < result->numRows(); ++i) { + const double s1 = (*result)(i, 2); + const double s2 = (*result)(i, 4); + EXPECT_DOUBLE_EQ((*result)(i, 1), 0.01 * s2); + EXPECT_DOUBLE_EQ((*result)(i, 3), 0.1 * s1); + EXPECT_DOUBLE_EQ((*result)(i, 5), 0.01 * s2); + } + // the last row is the state of the model at the end + std::vector x = rr.getSelectedValues(); + ASSERT_EQ(x.size(), 6); + EXPECT_EQ(x[1], rr.getValue("_J1")); + EXPECT_EQ(x[3], rr.getValue("_J0")); + EXPECT_EQ(x[5], rr.getValue("_J1")); + EXPECT_EQ(x[3], (*result)(10, 3)); + delete testModel; + } + void getIndependentFloatingSpeciesIds() { TestModel *testModel = TestModelFactory("SimpleFlux"); RoadRunner rr(testModel->str()); diff --git a/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp b/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp index 343f2b2659..f06e109093 100644 --- a/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp +++ b/test/cxx_api_tests/RoadRunnerAPITestsWithLLJit.cpp @@ -194,6 +194,10 @@ TEST_F(RoadRunnerAPITestsWithLLJit, getSelectedValues){ getSelectedValues(); } +TEST_F(RoadRunnerAPITestsWithLLJit, getSelectedReactionRates){ + getSelectedReactionRates(); +} + TEST_F(RoadRunnerAPITestsWithLLJit, getIndependentFloatingSpeciesIds){ getIndependentFloatingSpeciesIds(); } diff --git a/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp b/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp index 464fc96e77..546bd09199 100644 --- a/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp +++ b/test/cxx_api_tests/RoadRunnerAPITestsWithMCJit.cpp @@ -194,6 +194,10 @@ TEST_F(RoadRunnerAPITestsWithMCJit, getSelectedValues){ getSelectedValues(); } +TEST_F(RoadRunnerAPITestsWithMCJit, getSelectedReactionRates){ + getSelectedReactionRates(); +} + TEST_F(RoadRunnerAPITestsWithMCJit, getIndependentFloatingSpeciesIds){ getIndependentFloatingSpeciesIds(); }