diff --git a/source/llvm/LLVMModelGenerator.cpp b/source/llvm/LLVMModelGenerator.cpp index 9e809aa1aa..aaed4bbafa 100644 --- a/source/llvm/LLVMModelGenerator.cpp +++ b/source/llvm/LLVMModelGenerator.cpp @@ -312,25 +312,46 @@ namespace rrllvm { std::vector newSymbols = newModel->getRateRuleSymbols(); + // Species whose value is derived from a conserved-moiety total (see + // LLVMExecutableModel::setFloatingSpeciesAmounts) must be transferred + // after every other floating species: setting one of them shifts the + // conserved total by the difference from its currently derived value, + // which is only correct once the other species in that conservation + // law hold their transferred values rather than the new model's + // freshly-generated init values. + auto transferFloatingSpecies = [&](int i, const std::string& id) { + double value = 0; + oldModel->getFloatingSpeciesAmounts(1, &i, &value); + try { + newModel->setValue(id, value); + } + catch (const exception& e) { + rrLog(Logger::LOG_WARNING) << "regenerateModel: failed to transfer current " + "value for floating species '" << id << "' (value " << value + << "): " << e.what(); + } + }; + + std::vector deferredConservedMoietySpecies; + for (int i = 0; i < oldModel->getNumFloatingSpecies(); i++) { std::string id = oldModel->getFloatingSpeciesId(i); int index = newModel->getFloatingSpeciesIndex(id); if (index >= 0) { // new model has this species - double value = 0; - oldModel->getFloatingSpeciesAmounts(1, &i, &value); - try { - newModel->setValue(id, value); - } - catch (const exception& e) { - rrLog(Logger::LOG_WARNING) << "regenerateModel: failed to transfer current " - "value for floating species '" << id << "' (value " << value - << "): " << e.what(); + if (newModel->symbols->isConservedMoietySpecies(id)) { + deferredConservedMoietySpecies.push_back(i); + continue; } + transferFloatingSpecies(i, id); } } + for (int i : deferredConservedMoietySpecies) { + transferFloatingSpecies(i, oldModel->getFloatingSpeciesId(i)); + } + for (int i = 0; i < oldModel->getNumBoundarySpecies(); i++) { std::string id = oldModel->getBoundarySpeciesId(i); diff --git a/test/model_analysis/model_analysis.cpp b/test/model_analysis/model_analysis.cpp index 52a32d1520..19e16b6d26 100644 --- a/test/model_analysis/model_analysis.cpp +++ b/test/model_analysis/model_analysis.cpp @@ -22,6 +22,29 @@ class ModelAnalysisTests : public RoadRunnerTest { }; +TEST_F(ModelAnalysisTests, issue1351_steadyStateAfterSimulate) { + // https://github.com/sys-bio/roadrunner/issues/1351 + // This model (ATP <-> ADP <-> AMP) has one conserved moiety, so + // ATP + ADP + AMP must equal the same total (2.1 + 1.5 + 0.33 = 3.93) + // whether it's read right after construction, after simulate(), or + // after steadyState() runs following a simulate(). + const double conservedTotal = 2.1 + 1.5 + 0.33; + + RoadRunner ssOnly((modelAnalysisModelsDir / "steady_state_after_sim.xml").string()); + ssOnly.steadyState(); + EXPECT_NEAR(ssOnly.getValue("ATP") + ssOnly.getValue("ADP") + ssOnly.getValue("AMP"), + conservedTotal, 1e-6); + + RoadRunner rr((modelAnalysisModelsDir / "steady_state_after_sim.xml").string()); + rr.simulate(0, 100, 200); + EXPECT_NEAR(rr.getValue("ATP") + rr.getValue("ADP") + rr.getValue("AMP"), + conservedTotal, 1e-6); + + rr.steadyState(); + EXPECT_NEAR(rr.getValue("ATP") + rr.getValue("ADP") + rr.getValue("AMP"), + conservedTotal, 1e-6); +} + // https://github.com/sys-bio/roadrunner/issues/1339 // SimulateOptions is a long-lived object (RoadRunner::self.simulateOpt) that // the C API mutates field-by-field across repeated simulate() calls. The diff --git a/test/models/ModelAnalysis/steady_state_after_sim.xml b/test/models/ModelAnalysis/steady_state_after_sim.xml new file mode 100644 index 0000000000..dff1dac757 --- /dev/null +++ b/test/models/ModelAnalysis/steady_state_after_sim.xml @@ -0,0 +1,87 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + kf + AMP + ATP + + + + kr + + + ADP + 2 + + + + + + + + + + + + + + + + + + kc + ATP + + + + + + + + + + + + + + + + kg + ADP + + + + + + +