diff --git a/.github/workflows/update_save_files.yml b/.github/workflows/update_save_files.yml index ffd47b1eb2..03f14f08d0 100644 --- a/.github/workflows/update_save_files.yml +++ b/.github/workflows/update_save_files.yml @@ -31,7 +31,7 @@ jobs: libroadrunner_deps_owner: [ "sys-bio" ] libroadrunner_deps_repo: [ "libroadrunner-deps" ] libroadrunner_deps_name: [ "libroadrunner-deps" ] - libroadrunner_deps_release_version: [ "v2.2.22" ] + libroadrunner_deps_release_version: [ "v2.3.3" ] llvm_owner: [ "sys-bio" ] llvm_repo: [ "llvm-13.x" ] llvm_name: [ "llvm-13.x" ] diff --git a/CMakeLists.txt b/CMakeLists.txt index abad131c23..0e4c61dbfe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,8 +8,8 @@ cmake_minimum_required(VERSION 3.16) # Version information and include modules set(ROADRUNNER_VERSION_MAJOR 2) -set(ROADRUNNER_VERSION_MINOR 9) -set(ROADRUNNER_VERSION_PATCH 3) +set(ROADRUNNER_VERSION_MINOR 10) +set(ROADRUNNER_VERSION_PATCH 0) set(ROADRUNNER_VERSION "${ROADRUNNER_VERSION_MAJOR}.${ROADRUNNER_VERSION_MINOR}.${ROADRUNNER_VERSION_PATCH}") diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index e716866e33..7eebc5b520 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -256,6 +256,7 @@ if (BUILD_LLVM) llvm/ModelInitialValueSymbolResolver.cpp llvm/ModelResources.cpp llvm/Random.cpp + llvm/ResetRateRuleStoichCodeGen.cpp llvm/SBMLInitialValueSymbolResolver.cpp llvm/SBMLModelObjectCache.cpp llvm/SBMLSupportFunctions.cpp @@ -304,6 +305,7 @@ if (BUILD_LLVM) llvm/ModelInitialValueSymbolResolver.h llvm/ModelResources.h llvm/Random.h + llvm/ResetRateRuleStoichCodeGen.h llvm/SBMLInitialValueSymbolResolver.h llvm/SBMLModelObjectCache.h llvm/SBMLSupportFunctions.h diff --git a/source/llvm/EvalInitialConditionsCodeGen.cpp b/source/llvm/EvalInitialConditionsCodeGen.cpp index e1d9fae4b5..07d5c5a311 100644 --- a/source/llvm/EvalInitialConditionsCodeGen.cpp +++ b/source/llvm/EvalInitialConditionsCodeGen.cpp @@ -214,7 +214,16 @@ void EvalInitialConditionsCodeGen::codeGenStoichiometry( // in 1100 tests) if (!nz.id.empty() && dataSymbols.hasRateRule(nz.id)) { - modelDataBuilder.createRateRuleValueStore(nz.id, stoichValue); + // stoichValue is the net, CSR-signed value (negative for a + // reactant); the rate rule slot holds the reference's own + // magnitude, so undo that sign flip before storing it there. + Value *rateRuleSeed = stoichValue; + if (nz.type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant) + { + Value *negOne = ConstantFP::get(builder.getContext(), APFloat(-1.0)); + rateRuleSeed = builder.CreateFMul(negOne, stoichValue, "unneg_" + nz.id); + } + modelDataBuilder.createRateRuleValueStore(nz.id, rateRuleSeed); } Value *row = ConstantInt::get(Type::getInt32Ty(context), nz.row, true); diff --git a/source/llvm/EvalVolatileStoichCodeGen.cpp b/source/llvm/EvalVolatileStoichCodeGen.cpp index 4dadfb94ab..8ca9d51f3b 100644 --- a/source/llvm/EvalVolatileStoichCodeGen.cpp +++ b/source/llvm/EvalVolatileStoichCodeGen.cpp @@ -101,17 +101,23 @@ namespace rrllvm { "generating update code for non-constant species " "reference reactant " << r->getId(); - const StoichiometryMath *sm = r->getStoichiometryMath(); - if (!sm){ - rrLog(Logger::LOG_WARNING) << "No stoichiometry found for " - "species \"" << r->getId() << "\"" - " in reaction \"" << reaction->getName() << "\"" << std::endl; + Value *value = 0; + + if (dataSymbols.hasAssignmentRule(r->getId()) + || dataSymbols.hasRateRule(r->getId())) { + value = resolver.loadSymbolValue(r->getId()); + } else if (r->isSetStoichiometryMath()) { + const StoichiometryMath *sm = r->getStoichiometryMath(); + value = astCodeGen.codeGenDouble(sm->getMath()); + } else { + rrLog(Logger::LOG_WARNING) << "species reference " + << r->getId() << " has been determined to be " + "non-constant, but it has no rules or MathML, so" + " no update code will be generated"; continue; } -// assert(sm && "SmoichiometryMath variable sm is nullptr"); - - Value *value = astCodeGen.codeGenDouble(sm->getMath()); + assert(value && "value for species reference stoichiometry is 0"); // reactants are consumed, so they get a negative stoichiometry Value *negOne = ConstantFP::get(builder.getContext(), APFloat(-1.0)); diff --git a/source/llvm/Jit.cpp b/source/llvm/Jit.cpp index 8f77cb9d57..73e79b1420 100644 --- a/source/llvm/Jit.cpp +++ b/source/llvm/Jit.cpp @@ -135,6 +135,8 @@ namespace rrllvm { "getPiecewiseTrigger")); modelResources->evalVolatileStoichPtr = reinterpret_cast(lookupFunctionAddress( "evalVolatileStoich")); + modelResources->resetRateRuleStoichPtr = reinterpret_cast(lookupFunctionAddress( + "resetRateRuleStoich")); modelResources->evalConversionFactorPtr = reinterpret_cast(lookupFunctionAddress( "evalConversionFactor")); if (options & LoadSBMLOptions::READ_ONLY) { diff --git a/source/llvm/LLVMExecutableModel.cpp b/source/llvm/LLVMExecutableModel.cpp index 3a7c0985d0..d993df2033 100644 --- a/source/llvm/LLVMExecutableModel.cpp +++ b/source/llvm/LLVMExecutableModel.cpp @@ -184,6 +184,7 @@ LLVMExecutableModel::LLVMExecutableModel() : eventAssignPtr(0), getPiecewiseTriggerPtr(0), evalVolatileStoichPtr(0), + resetRateRuleStoichPtr(0), evalConversionFactorPtr(0), setBoundarySpeciesAmountPtr(0), setFloatingSpeciesAmountPtr(0), @@ -232,6 +233,7 @@ LLVMExecutableModel::LLVMExecutableModel( eventAssignPtr(modelResources->eventAssignPtr), getPiecewiseTriggerPtr(modelResources->getPiecewiseTriggerPtr), evalVolatileStoichPtr(modelResources->evalVolatileStoichPtr), + resetRateRuleStoichPtr(modelResources->resetRateRuleStoichPtr), evalConversionFactorPtr(modelResources->evalConversionFactorPtr), setBoundarySpeciesAmountPtr(modelResources->setBoundarySpeciesAmountPtr), setFloatingSpeciesAmountPtr(modelResources->setFloatingSpeciesAmountPtr), @@ -300,6 +302,7 @@ LLVMExecutableModel::LLVMExecutableModel(std::istream& in, uint modelGeneratorOp eventAssignPtr = resources->eventAssignPtr; getPiecewiseTriggerPtr = resources->getPiecewiseTriggerPtr; evalVolatileStoichPtr = resources->evalVolatileStoichPtr; + resetRateRuleStoichPtr = resources->resetRateRuleStoichPtr; evalConversionFactorPtr = resources->evalConversionFactorPtr; setBoundarySpeciesAmountPtr = resources->setBoundarySpeciesAmountPtr; setFloatingSpeciesAmountPtr = resources->setFloatingSpeciesAmountPtr; @@ -832,6 +835,14 @@ void LLVMExecutableModel::reset(int opt) resetOneType(opt, SelectionRecord::_GLOBAL_PARAMETER, modelData->numIndGlobalParameters, getNumGlobalParameters(), &LLVMExecutableModel::getGlobalParameterInitValues, &LLVMExecutableModel::setGlobalParameterValues, &LLVMModelDataSymbols::getGlobalParameterId, buffer, inits, initvals); //std::cout << this; + // Reset rate-rule-controlled stoichiometries, mirroring how rate-rule- + // controlled compartments/species/parameters are restored above under + // the RATE flag. + if ((opt & SelectionRecord::RATE) && !symbols->isConservedMoietyAnalysis()) + { + resetRateRuleStoichPtr(modelData); + } + // Whether were we forced to reset cms: bool reset_cm = false; @@ -1242,15 +1253,10 @@ void LLVMExecutableModel::getIds(int types, std::list &ids) } if (checkExact(SelectionRecord::STOICHIOMETRY, types)) { - for (size_t s = 0; s < getNumIndFloatingSpecies(); ++s) { - string sid = getFloatingSpeciesId(s); - for (size_t r = 0; r < getNumReactions(); ++r) { - if (getStoichiometry(s, r) != 0) - { - string rid = getReactionId(r); - ids.push_back(symbols->getStoichiometryIdFor(sid, rid)); - } - } + for (const LLVMModelDataSymbols::SpeciesReferenceInfo& entry : symbols->getStoichiometryList()) { + string sid = getFloatingSpeciesId(entry.row); + string rid = getReactionId(entry.column); + ids.push_back(symbols->getStoichiometryIdFor(sid, rid)); } } @@ -2369,6 +2375,13 @@ int LLVMExecutableModel::setStoichiometry(int index, double value) std::list::const_iterator stoichiometry = stoichiometryIndx.begin(); for (int i = 0; i < index; i++) ++stoichiometry; + + if (!stoichiometry->id.empty() && symbols->hasAssignmentRule(stoichiometry->id)) + throw LLVMException("Cannot set stoichiometry '" + stoichiometry->id + "': it is governed by an assignment rule"); + + if (!stoichiometry->id.empty() && symbols->hasRateRule(stoichiometry->id)) + modelData->rateRuleValuesAlias[symbols->getRateRuleIndex(stoichiometry->id)] = value; + if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Product) return setStoichiometry(stoichiometry->row, stoichiometry->column, value); else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant) @@ -2417,6 +2430,13 @@ int LLVMExecutableModel::setInitStoichiometry(int index, double value) std::list::const_iterator stoichiometry = stoichiometryIndx.begin(); for (int i = 0; i < index; i++) ++stoichiometry; + + if (!stoichiometry->id.empty() && symbols->hasAssignmentRule(stoichiometry->id)) + throw LLVMException("Cannot set stoichiometry '" + stoichiometry->id + "': it is governed by an assignment rule"); + + if (!stoichiometry->id.empty() && symbols->hasRateRule(stoichiometry->id)) + modelData->rateRuleValuesAlias[symbols->getRateRuleIndex(stoichiometry->id)] = value; + if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Product) return setInitStoichiometry(stoichiometry->row, stoichiometry->column, value); else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant) @@ -2444,6 +2464,13 @@ double LLVMExecutableModel::getStoichiometry(int index) if (index < 0) throw LLVMException("The stoichiometry index is not valid"); + + // volatile (rule-governed) entries are only resynced into the CSR at + // RHS-eval / event-root time during simulation; resync here too so a + // plain getValue() reflects the model's current time/state exactly, + // rather than whatever the last internal integrator step left behind. + evalVolatileStoichPtr(modelData); + std::list stoichiometryIndx = symbols->getStoichiometryList(); std::list::const_iterator stoichiometry = stoichiometryIndx.begin(); for (int i = 0; i < index; i++) diff --git a/source/llvm/LLVMExecutableModel.h b/source/llvm/LLVMExecutableModel.h index 9add0ce22f..f11486f0a3 100644 --- a/source/llvm/LLVMExecutableModel.h +++ b/source/llvm/LLVMExecutableModel.h @@ -34,6 +34,7 @@ #include "EventTriggerCodeGen.h" #include "GetPiecewiseTriggerCodeGen.h" #include "EvalVolatileStoichCodeGen.h" +#include "ResetRateRuleStoichCodeGen.h" #include "EvalConversionFactorCodeGen.h" #include "SetValuesCodeGen.h" #include "SetInitialValuesCodeGen.h" @@ -678,6 +679,7 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel EventAssignCodeGen::FunctionPtr eventAssignPtr; GetPiecewiseTriggerCodeGen::FunctionPtr getPiecewiseTriggerPtr; EvalVolatileStoichCodeGen::FunctionPtr evalVolatileStoichPtr; + ResetRateRuleStoichCodeGen::FunctionPtr resetRateRuleStoichPtr; EvalConversionFactorCodeGen::FunctionPtr evalConversionFactorPtr; // set model values externally. diff --git a/source/llvm/LLVMModelDataSymbols.cpp b/source/llvm/LLVMModelDataSymbols.cpp index abc8bcfc9b..849dfdeddb 100644 --- a/source/llvm/LLVMModelDataSymbols.cpp +++ b/source/llvm/LLVMModelDataSymbols.cpp @@ -168,7 +168,7 @@ LLVMModelDataSymbols::LLVMModelDataSymbols(const libsbml::Model *model, unsigned else if (rule->getTypeCode() == SBML_RATE_RULE) { size_t rri = rateRules.size(); - rateRules[rule->getId()] = rri; + rateRules[rule->getVariable()] = rri; } else if (rule->getTypeCode() == SBML_ALGEBRAIC_RULE) { @@ -267,14 +267,10 @@ LLVMModelDataSymbols::SymbolIndexType LLVMModelDataSymbols::getSymbolIndex( result = i->second; return EVENT; } - else + else if ((i = stoichiometryMap.find(name)) != stoichiometryMap.end()) { - for (int stoichIndex = 0 ; stoichIndex < stoichIds.size(); ++stoichIndex) { - if (stoichIds.at(stoichIndex) == name) { - result = stoichIndex; - return STOICHIOMETRY; - } - } + result = i->second; + return STOICHIOMETRY; } result = -1; @@ -376,12 +372,11 @@ size_t LLVMModelDataSymbols::getReactionSize() const int LLVMModelDataSymbols::getStoichiometryIndex(const std::string& id) const { - for (int i = 0; i < stoichIds.size(); ++i) + StringUIntMap::const_iterator i = stoichiometryMap.find(id); + if (i != stoichiometryMap.end()) { - if (stoichIds.at(i) == id) - return i; + return i->second; } - return -1; } @@ -1402,6 +1397,11 @@ void LLVMModelDataSymbols::initReactions(const libsbml::Model* model) // index of the just added Reactant speciesMap[speciesIdx] = stoichTypes.size() - 1; + if (r->isSetId() && r->getId().length() > 0) + { + stoichiometryMap[r->getId()] = stoichIds.size() - 1; + } + } else { @@ -1462,6 +1462,11 @@ void LLVMModelDataSymbols::initReactions(const libsbml::Model* model) // index of the just added Reactant speciesMap[speciesIdx] = stoichTypes.size() - 1; + if (p->isSetId() && p->getId().length() > 0) + { + stoichiometryMap[p->getId()] = stoichIds.size() - 1; + } + } else { @@ -1478,21 +1483,21 @@ void LLVMModelDataSymbols::initReactions(const libsbml::Model* model) } } } - for (size_t ns = 0; ns < namedStoichiometryIds.size(); ns++) { - const Rule* rule = model->getRule(namedStoichiometryIds[ns]); - if (rule != NULL) { - std::string msg = "The named stoichiometry '" + namedStoichiometryIds[ns] + "' (also called a speciesReference) has "; - if (rule->isAssignment()) { - msg += "an assignment rule"; - } - else { - assert(rule->isRate()); - msg += "a rate rule"; - } - msg += ", which means that the stoichiometry of its reaction varies in time. Variable stoichiometries are not supported by roadrunner."; - throw_llvm_exception(msg); - } - } + //for (size_t ns = 0; ns < namedStoichiometryIds.size(); ns++) { + // const Rule* rule = model->getRule(namedStoichiometryIds[ns]); + // if (rule != NULL) { + // std::string msg = "The named stoichiometry '" + namedStoichiometryIds[ns] + "' (also called a speciesReference) has "; + // if (rule->isAssignment()) { + // msg += "an assignment rule"; + // } + // else { + // assert(rule->isRate()); + // msg += "a rate rule"; + // } + // msg += ", which means that the stoichiometry of its reaction varies in time. Variable stoichiometries are not supported by roadrunner."; + // throw_llvm_exception(msg); + // } + //} } @@ -1960,6 +1965,8 @@ void LLVMModelDataSymbols::saveState(std::ostream& out) const rr::saveBinary(out, stoichTypes); + rr::saveBinary(out, stoichiometryMap); + rr::saveBinary(out, assignmentRules); rr::saveBinary(out, rateRules); rr::saveBinary(out, globalParameterRateRules); @@ -2010,6 +2017,8 @@ void LLVMModelDataSymbols::loadState(std::istream& in) rr::loadBinary(in, stoichTypes); + rr::loadBinary(in, stoichiometryMap); + rr::loadBinary(in, assignmentRules); rr::loadBinary(in, rateRules); diff --git a/source/llvm/LLVMModelDataSymbols.h b/source/llvm/LLVMModelDataSymbols.h index d6c117e965..c36a64c065 100644 --- a/source/llvm/LLVMModelDataSymbols.h +++ b/source/llvm/LLVMModelDataSymbols.h @@ -710,6 +710,13 @@ class LLVMModelDataSymbols std::vector stoichTypes; + /** + * maps a named species reference's id to its index in the + * stoichIds / stoichRowIndx / stoichColIndx / stoichTypes arrays, + * so lookups don't require scanning stoichIds linearly. + */ + StringUIntMap stoichiometryMap; + /** * the set of rule, these contain the variable name of the rule so that * we can quickly see if a symbol has an associated rule. diff --git a/source/llvm/LLVMModelGenerator.cpp b/source/llvm/LLVMModelGenerator.cpp index 5a868107db..121be0028a 100644 --- a/source/llvm/LLVMModelGenerator.cpp +++ b/source/llvm/LLVMModelGenerator.cpp @@ -106,6 +106,7 @@ namespace rrllvm { dst->eventAssignPtr = src->eventAssignPtr; dst->getPiecewiseTriggerPtr = src->getPiecewiseTriggerPtr; dst->evalVolatileStoichPtr = src->evalVolatileStoichPtr; + dst->resetRateRuleStoichPtr = src->resetRateRuleStoichPtr; dst->evalConversionFactorPtr = src->evalConversionFactorPtr; } @@ -126,6 +127,7 @@ namespace rrllvm { EventAssignCodeGen(context).createFunction(); GetPiecewiseTriggerCodeGen(context).createFunction(); EvalVolatileStoichCodeGen(context).createFunction(); + ResetRateRuleStoichCodeGen(context).createFunction(); EvalConversionFactorCodeGen(context).createFunction(); Function* setBoundarySpeciesAmountIR = nullptr; @@ -401,8 +403,18 @@ namespace rrllvm { int new_s = newModel->getFloatingSpeciesIndex(sID); if (new_s < 0) continue; if (new_s >= newModel->getNumIndFloatingSpecies()) continue; - double stoich = oldModel->getStoichiometry(s, r); - newModel->setStoichiometry(new_s, new_r, stoich); + int oldIndex = oldModel->getStoichiometryIndex(sID, rID); + int newIndex = newModel->getStoichiometryIndex(sID, rID); + if (oldIndex < 0 || newIndex < 0) continue; + try { + double stoich = oldModel->getStoichiometry(oldIndex); + newModel->setStoichiometry(newIndex, stoich); + } + catch (const exception& e) { + rrLog(Logger::LOG_WARNING) << "regenerateModel: failed to transfer current " + "stoichiometry value for species '" << sID << "' in reaction '" << rID + << "': " << e.what(); + } } } diff --git a/source/llvm/ModelResources.h b/source/llvm/ModelResources.h index 98b6ed1060..0d05613614 100644 --- a/source/llvm/ModelResources.h +++ b/source/llvm/ModelResources.h @@ -63,6 +63,7 @@ namespace rrllvm EventAssignCodeGen::FunctionPtr eventAssignPtr; GetPiecewiseTriggerCodeGen::FunctionPtr getPiecewiseTriggerPtr; EvalVolatileStoichCodeGen::FunctionPtr evalVolatileStoichPtr; + ResetRateRuleStoichCodeGen::FunctionPtr resetRateRuleStoichPtr; EvalConversionFactorCodeGen::FunctionPtr evalConversionFactorPtr; SetBoundarySpeciesAmountCodeGen::FunctionPtr setBoundarySpeciesAmountPtr; SetFloatingSpeciesAmountCodeGen::FunctionPtr setFloatingSpeciesAmountPtr; diff --git a/source/llvm/ResetRateRuleStoichCodeGen.cpp b/source/llvm/ResetRateRuleStoichCodeGen.cpp new file mode 100644 index 0000000000..f54c84d7d0 --- /dev/null +++ b/source/llvm/ResetRateRuleStoichCodeGen.cpp @@ -0,0 +1,80 @@ +/* + * ResetRateRuleStoichCodeGen.cpp + */ +#pragma hdrstop + +#include "ResetRateRuleStoichCodeGen.h" +#include "LLVMException.h" +#include "ASTNodeCodeGen.h" +#include "SBMLInitialValueSymbolResolver.h" +#include "rrLogger.h" + +#include +#include + +namespace rrllvm { + using namespace rr; + using namespace llvm; + using namespace libsbml; + + + const char *ResetRateRuleStoichCodeGen::FunctionName = "resetRateRuleStoich"; + + ResetRateRuleStoichCodeGen::ResetRateRuleStoichCodeGen( + const ModelGeneratorContext &mgc) : + CodeGenBase(mgc) { + } + + ResetRateRuleStoichCodeGen::~ResetRateRuleStoichCodeGen() { + } + + Value *ResetRateRuleStoichCodeGen::codeGen() { + Value *modelData = 0; + + codeGenVoidModelDataHeader(FunctionName, modelData); + + // Reads declared/SBML values, same source evalInitialConditions uses, + // so this stays consistent with whatever the model was last + // (re)compiled with (e.g. after a setValue("init(...)") call, which + // mutates the SBML and regenerates the model). + SBMLInitialValueSymbolResolver initialValueResolver(modelData, modelGenContext); + ModelDataIRBuilder mdbuilder(modelData, dataSymbols, builder); + ASTNodeCodeGen astCodeGen(builder, initialValueResolver, modelGenContext, modelData); + + std::list stoichEntries = + dataSymbols.getStoichiometryList(); + + for (std::list::iterator i = + stoichEntries.begin(); i != stoichEntries.end(); i++) + { + const LLVMModelDataSymbols::SpeciesReferenceInfo& nz = *i; + + if (nz.id.empty() || !dataSymbols.hasRateRule(nz.id)) + { + continue; + } + + const ASTNode *node = modelSymbols.createStoichiometryNode(nz.row, nz.column); + Value *stoichValue = astCodeGen.codeGenDouble(node); + delete node; + + // stoichValue is the net, CSR-signed value (negative for a + // reactant); the rate rule slot holds the reference's own + // magnitude, so undo that sign flip before storing it there. + Value *rateRuleSeed = stoichValue; + if (nz.type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant) + { + Value *negOne = ConstantFP::get(builder.getContext(), APFloat(-1.0)); + rateRuleSeed = builder.CreateFMul(negOne, stoichValue, "unneg_" + nz.id); + } + + mdbuilder.createRateRuleValueStore(nz.id, rateRuleSeed); + mdbuilder.createStoichiometryStore(nz.row, nz.column, stoichValue, nz.id); + } + + builder.CreateRetVoid(); + + return verifyFunction(); + } + +} /* namespace rrllvm */ diff --git a/source/llvm/ResetRateRuleStoichCodeGen.h b/source/llvm/ResetRateRuleStoichCodeGen.h new file mode 100644 index 0000000000..d9088c6af7 --- /dev/null +++ b/source/llvm/ResetRateRuleStoichCodeGen.h @@ -0,0 +1,39 @@ +/* + * ResetRateRuleStoichCodeGen.h + */ + +#ifndef RESETRATERULESTOICHCODEGEN_H_ +#define RESETRATERULESTOICHCODEGEN_H_ + +#include "CodeGenBase.h" +#include "ModelGeneratorContext.h" +#include "ModelDataIRBuilder.h" +#include + +namespace rrllvm +{ + + typedef void (*ResetRateRuleStoichCodeGen_FunctionPtr)(LLVMModelData*); + + /** @class ResetRateRuleStoichCodeGen + * Restores every rate-rule-governed named stoichiometry (and the rate rule's + * own state) to the value declared in the SBML, without requiring the full + * evalInitialConditions pass. This is what lets a plain reset() (TIME | RATE + * | FLOATING) restore rate-rule-controlled stoichiometries the same way it + * already restores rate-rule-controlled compartments, species, and + * parameters. + */ + class ResetRateRuleStoichCodeGen : + public CodeGenBase + { + public: + ResetRateRuleStoichCodeGen(const ModelGeneratorContext& mgc); + virtual ~ResetRateRuleStoichCodeGen(); + + llvm::Value* codeGen(); + + static const char* FunctionName; + }; + +} /* namespace rrllvm */ +#endif /* RESETRATERULESTOICHCODEGEN_H_ */ diff --git a/source/rrUtils.cpp b/source/rrUtils.cpp index 1740b00c6f..b2a94d409b 100644 --- a/source/rrUtils.cpp +++ b/source/rrUtils.cpp @@ -468,8 +468,8 @@ bool hasUnimplementedTags(const std::string& descriptionFileName, const string& badtags.push_back("FastReaction"); badtags.push_back("VolumeConcentrationRate"); //badtags.push_back("RateOf"); - badtags.push_back("AssignedVariableStoichiometry"); - badtags.push_back("AssignedConstantStoichiometry"); + //badtags.push_back("AssignedVariableStoichiometry"); + //badtags.push_back("AssignedConstantStoichiometry"); if (integrator == "rk4" || integrator == "rk45") { diff --git a/test/models/SBMLFeatures/named_stoic.xml b/test/models/SBMLFeatures/named_stoic.xml new file mode 100644 index 0000000000..025666d7f1 --- /dev/null +++ b/test/models/SBMLFeatures/named_stoic.xml @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + 0.5 + + + + + + diff --git a/test/models/SBMLFeatures/named_stoic_multi_reactant.xml b/test/models/SBMLFeatures/named_stoic_multi_reactant.xml new file mode 100644 index 0000000000..bbbbeec9a8 --- /dev/null +++ b/test/models/SBMLFeatures/named_stoic_multi_reactant.xml @@ -0,0 +1,35 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + k + A + + + + + + + diff --git a/test/models/SBMLFeatures/stoich_ar.xml b/test/models/SBMLFeatures/stoich_ar.xml new file mode 100644 index 0000000000..0ee3a94046 --- /dev/null +++ b/test/models/SBMLFeatures/stoich_ar.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + time + + + + + + + + + + + 0.5 + + + + + + diff --git a/test/models/SBMLFeatures/stoich_ar_const.xml b/test/models/SBMLFeatures/stoich_ar_const.xml new file mode 100644 index 0000000000..e08d0aea1a --- /dev/null +++ b/test/models/SBMLFeatures/stoich_ar_const.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + 1 + + + + + + + + + + + 0.5 + + + + + + diff --git a/test/models/SBMLFeatures/stoich_event.xml b/test/models/SBMLFeatures/stoich_event.xml new file mode 100644 index 0000000000..4b8281f9a9 --- /dev/null +++ b/test/models/SBMLFeatures/stoich_event.xml @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + + 0.5 + + + + + + + + + + + time + 1 + + + + + + + 1 + + + + + + + diff --git a/test/models/SBMLFeatures/stoich_rr.xml b/test/models/SBMLFeatures/stoich_rr.xml new file mode 100644 index 0000000000..f3315b756e --- /dev/null +++ b/test/models/SBMLFeatures/stoich_rr.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + 1 + + + + + + + + + + + 0.5 + + + + + + diff --git a/test/sbml_features/CMakeLists.txt b/test/sbml_features/CMakeLists.txt index e596e8283c..28dd86801d 100644 --- a/test/sbml_features/CMakeLists.txt +++ b/test/sbml_features/CMakeLists.txt @@ -11,6 +11,7 @@ add_test_executable( rateOf.cpp unsupported_packages.cpp valid_sbml.cpp + variable_stoich.cpp ) # make test_targets list global to all tests diff --git a/test/sbml_features/named_stoich.cpp b/test/sbml_features/named_stoich.cpp index 9b31aa36bc..45ea77a7d7 100644 --- a/test/sbml_features/named_stoich.cpp +++ b/test/sbml_features/named_stoich.cpp @@ -100,12 +100,73 @@ TEST_F(SBMLFeatures, named_stoich_values) { } -TEST_F(SBMLFeatures, no_variable_named_stoich) { +TEST_F(SBMLFeatures, add_rule_to_named_stoich) { + // "n" is the reactant B's stoichiometry (declared 1); "m" (C, declared 2) + // and "q" (D, declared 3) are unaffected. Kinetic law rate = (n+m+q)*A. rr::RoadRunner rr((SBMLFeaturesDir / "named_stoic_in_kinetic_law.xml").string()); - EXPECT_THROW(rr.addAssignmentRule("n", "5", true), rrllvm::LLVMException); + EXPECT_NO_THROW(rr.addAssignmentRule("n", "5", true)); + EXPECT_EQ(rr.getValue("n"), 5.0); + // rate = (5 + 2 + 3) * A(1) = 10 + EXPECT_NEAR(rr.getValue("J0"), 10.0, 1e-9); + // "q" is the product D's stoichiometry (declared 3). rr::RoadRunner rr2((SBMLFeaturesDir / "named_stoic_in_kinetic_law.xml").string()); - EXPECT_THROW(rr2.addRateRule("q", "5", true), rrllvm::LLVMException); + EXPECT_NO_THROW(rr2.addRateRule("q", "5", true)); + // rate rules haven't integrated yet at t=0; q should still be its declared value. + EXPECT_NEAR(rr2.getValue("q"), 3.0, 1e-9); + rr2.oneStep(0, 1.0); + // dq/dt = 5, so after 1 time unit q should have grown by ~5. + EXPECT_NEAR(rr2.getValue("q"), 8.0, 1e-6); +} + + +TEST_F(SBMLFeatures, add_rate_rule_to_reactant_stoich) { + // "n" is reactant B's stoichiometry (declared 1). Mirrors the "q" + // (product) rate-rule case in add_rule_to_named_stoich, but exercises the + // reactant sign-handling path: createStoichiometryNode negates reactant + // values for the CSR cell, and that same negated value must not leak into + // the rate-rule integration slot. + rr::RoadRunner rr((SBMLFeaturesDir / "named_stoic_in_kinetic_law.xml").string()); + EXPECT_NO_THROW(rr.addRateRule("n", "5", true)); + // rate rules haven't integrated yet at t=0; n should still be its declared value. + EXPECT_NEAR(rr.getValue("n"), 1.0, 1e-9); + rr.oneStep(0, 1.0); + // dn/dt = 5, so after 1 time unit n should have grown by ~5, mirroring q's 3->8. + EXPECT_NEAR(rr.getValue("n"), 6.0, 1e-6); +} + + +TEST_F(SBMLFeatures, add_assignment_rule_to_reactant_stoich_tracks_time) { + // "n" is reactant B's stoichiometry. An assignment rule that depends on + // time must be resynced into the stoichiometry matrix as the simulation + // progresses, not just read once at t=0. + rr::RoadRunner rr((SBMLFeaturesDir / "named_stoic_in_kinetic_law.xml").string()); + EXPECT_NO_THROW(rr.addAssignmentRule("n", "1 + time", true)); + EXPECT_NEAR(rr.getValue("n"), 1.0, 1e-9); + rr.oneStep(0, 2.0); + EXPECT_NEAR(rr.getValue("n"), 3.0, 1e-9); +} + + +TEST_F(SBMLFeatures, set_value_on_multi_reactant_stoich_before_simulate) { + // A appears twice as a reactant in J0, via two independently named, + // rule-free species references (r1=2, r2=3). Setting one before + // simulating must not disturb the other, and the combined (net) cell used + // in the reaction must reflect the sum of both. + rr::RoadRunner rr((SBMLFeaturesDir / "named_stoic_multi_reactant.xml").string()); + EXPECT_EQ(rr.getValue("r1"), 2.0); + EXPECT_EQ(rr.getValue("r2"), 3.0); + + EXPECT_NO_THROW(rr.setValue("r1", 5)); + EXPECT_EQ(rr.getValue("r1"), 5.0); + EXPECT_EQ(rr.getValue("r2"), 3.0); + + // combined cell should be r1 + r2 = 8 + EXPECT_NEAR(rr.getValue("stoich(A, J0)"), 8.0, 1e-9); + + // dA/dt = -(r1+r2)*k*A < 0, A should be decreasing + rr.oneStep(0, 0.01); + EXPECT_LT(rr.getValue("A"), 10.0); } @@ -127,6 +188,28 @@ TEST_F(SBMLFeatures, named_stoich_selectors) { } +TEST_F(SBMLFeatures, named_stoich_set_and_reset) { + RoadRunner rr((SBMLFeaturesDir / "named_stoic.xml").string()); + rr.getSimulateOptions().setDuration(10); + + const ls::DoubleMatrix run1 = *rr.simulate(); + + rr.reset(SelectionRecord::ALL); + rr.setValue("N", 2); + const ls::DoubleMatrix run2 = *rr.simulate(); + + rr.reset(SelectionRecord::ALL); + const ls::DoubleMatrix run3 = *rr.simulate(); + + ASSERT_EQ(run1.numRows(), run2.numRows()); + ASSERT_EQ(run1.numRows(), run3.numRows()); + for (int i = 0; i < run1.numRows(); i++) { + EXPECT_EQ(run1(i, 1), run3(i, 1)); + EXPECT_NEAR(run2(i, 1), 2 * run1(i, 1), 1e-12); + } +} + + TEST_F(SBMLFeatures, get_named_stoich_value_from_model) { rr::RoadRunner rr((SBMLFeaturesDir / "named_stoic_in_kinetic_law.xml").string()); ExecutableModel* em = rr.getModel(); diff --git a/test/sbml_features/variable_stoich.cpp b/test/sbml_features/variable_stoich.cpp new file mode 100644 index 0000000000..30360f51d3 --- /dev/null +++ b/test/sbml_features/variable_stoich.cpp @@ -0,0 +1,187 @@ +#include "gtest/gtest.h" +#include "rrRoadRunner.h" +#include "rrException.h" +#include "rrUtils.h" +#include "rrTestSuiteModelSimulation.h" +#include "sbml/SBMLTypes.h" +#include "sbml/SBMLReader.h" +#include "../test_util.h" +#include +#include +#include "RoadRunnerTest.h" +#include "llvm/LLVMException.h" + +using namespace testing; +using namespace rr; +using namespace std; +using std::filesystem::path; + +class SBMLFeatures : public RoadRunnerTest { +public: + path SBMLFeaturesDir = rrTestModelsDir_ / "SBMLFeatures"; + SBMLFeatures() = default; +}; + +// stoich_rr.xml: N governed by a rate rule, dN/dt = 1, N(0) = 0. dX/dt = N*0.5. +// stoich_ar.xml: N governed by an assignment rule, N = time. dX/dt = N*0.5. +// stoich_ar_const.xml: N governed by an assignment rule, N = 1. dX/dt = N*0.5 = 0.5. +// stoich_event.xml: N unruled, starts at 0, event at time>1 sets N = 1. dX/dt = N*0.5. + +TEST_F(SBMLFeatures, variable_stoich_rr_simulate) { + RoadRunner rr((SBMLFeaturesDir / "stoich_rr.xml").string()); + const ls::DoubleMatrix result = *rr.simulate(0.0, 2.0, 11); + for (int i = 0; i < result.numRows(); i++) { + EXPECT_NEAR(result(i, 1), 0.01 * i * i, 1e-6); + } +} + +TEST_F(SBMLFeatures, variable_stoich_ar_simulate) { + RoadRunner rr((SBMLFeaturesDir / "stoich_ar.xml").string()); + const ls::DoubleMatrix result = *rr.simulate(0.0, 2.0, 11); + for (int i = 0; i < result.numRows(); i++) { + EXPECT_NEAR(result(i, 1), 0.01 * i * i, 1e-6); + } +} + +TEST_F(SBMLFeatures, variable_stoich_ar_const_simulate) { + RoadRunner rr((SBMLFeaturesDir / "stoich_ar_const.xml").string()); + const ls::DoubleMatrix result = *rr.simulate(0.0, 2.0, 11); + for (int i = 0; i < result.numRows(); i++) { + EXPECT_NEAR(result(i, 1), 0.1 * i, 1e-6); + } +} + +TEST_F(SBMLFeatures, variable_stoich_event_simulate) { + RoadRunner rr((SBMLFeaturesDir / "stoich_event.xml").string()); + const ls::DoubleMatrix result = *rr.simulate(0.0, 2.0, 11); + for (int i = 0; i < result.numRows(); i++) { + double expected = (i <= 5) ? 0.0 : 0.1 * (i - 5); + EXPECT_NEAR(result(i, 1), expected, 1e-4); + } +} + + +TEST_F(SBMLFeatures, variable_stoich_rr_set_initial) { + RoadRunner rr((SBMLFeaturesDir / "stoich_rr.xml").string()); + rr.setValue("N", 2); + const ls::DoubleMatrix result = *rr.simulate(0.0, 2.0, 11); + for (int i = 0; i < result.numRows(); i++) { + EXPECT_NEAR(result(i, 1), 0.2 * i + 0.01 * i * i, 1e-6); + } +} + +TEST_F(SBMLFeatures, variable_stoich_ar_set_throws) { + RoadRunner rr((SBMLFeaturesDir / "stoich_ar.xml").string()); + EXPECT_THROW(rr.setValue("N", 2), rrllvm::LLVMException); +} + +TEST_F(SBMLFeatures, variable_stoich_ar_const_set_throws) { + RoadRunner rr((SBMLFeaturesDir / "stoich_ar_const.xml").string()); + EXPECT_THROW(rr.setValue("N", 2), rrllvm::LLVMException); +} + +TEST_F(SBMLFeatures, variable_stoich_event_set_initial) { + RoadRunner rr((SBMLFeaturesDir / "stoich_event.xml").string()); + rr.setValue("N", 2); + const ls::DoubleMatrix result = *rr.simulate(0.0, 2.0, 11); + for (int i = 0; i < result.numRows(); i++) { + double expected = (i <= 5) ? 0.2 * i : 0.5 + 0.1 * i; + EXPECT_NEAR(result(i, 1), expected, 1e-4); + } +} + + +// A plain reset (TIME | RATE | FLOATING) mirrors global parameter semantics: +// it only restores a stoichiometry if a rate rule governs it. stoich_rr's N +// does have a rate rule, so plain reset() restores it. +TEST_F(SBMLFeatures, variable_stoich_rr_reset_restores_initial) { + RoadRunner rr((SBMLFeaturesDir / "stoich_rr.xml").string()); + const ls::DoubleMatrix run1 = *rr.simulate(0.0, 2.0, 11); + + rr.reset(); + rr.setValue("init(N)", 2); + const ls::DoubleMatrix run2 = *rr.simulate(0.0, 2.0, 11); + + rr.reset(); + const ls::DoubleMatrix run3 = *rr.simulate(0.0, 2.0, 11); + + for (int i = 0; i < run1.numRows(); i++) { + EXPECT_NEAR(run1(i, 1), 0.01 * i * i, 1e-6); + EXPECT_NEAR(run2(i, 1), 0.2 * i + 0.01 * i * i, 1e-6); + EXPECT_NEAR(run3(i, 1), 0.2 * i + 0.01 * i * i, 1e-6); + } +} + +// An explicit full reset (SelectionRecord::ALL, which includes SBML_INITIALIZE) +// restores a stoichiometry to whatever init(N) is CURRENTLY configured to -- +// not back to the original SBML-declared value. Since init(N) was changed to 2 +// before the second reset, run3 should match run2 (N0=2), not run1 (N0=0). +TEST_F(SBMLFeatures, variable_stoich_rr_full_reset_restores_initial) { + RoadRunner rr((SBMLFeaturesDir / "stoich_rr.xml").string()); + const ls::DoubleMatrix run1 = *rr.simulate(0.0, 2.0, 11); + + rr.reset(SelectionRecord::ALL); + rr.setValue("init(N)", 2); + const ls::DoubleMatrix run2 = *rr.simulate(0.0, 2.0, 11); + + rr.reset(SelectionRecord::ALL); + const ls::DoubleMatrix run3 = *rr.simulate(0.0, 2.0, 11); + + for (int i = 0; i < run1.numRows(); i++) { + EXPECT_NEAR(run3(i, 1), 0.2 * i + 0.01 * i * i, 1e-6); + } +} + +// stoich_event's N has no rate rule, so a plain reset() correctly leaves it +// alone (mirrors how a plain, non-rate-rule global parameter isn't restored +// by default reset() either). Only an explicit full reset restores it, to +// whatever init(N) is currently configured to -- see run2's formula below. +TEST_F(SBMLFeatures, variable_stoich_event_full_reset_restores_initial) { + RoadRunner rr((SBMLFeaturesDir / "stoich_event.xml").string()); + const ls::DoubleMatrix run1 = *rr.simulate(0.0, 2.0, 11); + + rr.reset(SelectionRecord::ALL); + rr.setValue("init(N)", 2); + const ls::DoubleMatrix run2 = *rr.simulate(0.0, 2.0, 11); + + rr.reset(SelectionRecord::ALL); + const ls::DoubleMatrix run3 = *rr.simulate(0.0, 2.0, 11); + + for (int i = 0; i < run1.numRows(); i++) { + double expected2 = (i <= 5) ? 0.2 * i : 0.5 + 0.1 * i; + EXPECT_NEAR(run3(i, 1), expected2, 1e-4); + } +} + + +TEST_F(SBMLFeatures, variable_stoich_rr_selection_type) { + RoadRunner rr((SBMLFeaturesDir / "stoich_rr.xml").string()); + SelectionRecord record = rr.createSelection("N"); + EXPECT_EQ(record.selectionType, SelectionRecord::STOICHIOMETRY); + vector ids = rr.getStoichiometryIds(); + EXPECT_NE(std::find(ids.begin(), ids.end(), "N"), ids.end()); +} + +TEST_F(SBMLFeatures, variable_stoich_ar_selection_type) { + RoadRunner rr((SBMLFeaturesDir / "stoich_ar.xml").string()); + SelectionRecord record = rr.createSelection("N"); + EXPECT_EQ(record.selectionType, SelectionRecord::STOICHIOMETRY); + vector ids = rr.getStoichiometryIds(); + EXPECT_NE(std::find(ids.begin(), ids.end(), "N"), ids.end()); +} + +TEST_F(SBMLFeatures, variable_stoich_ar_const_selection_type) { + RoadRunner rr((SBMLFeaturesDir / "stoich_ar_const.xml").string()); + SelectionRecord record = rr.createSelection("N"); + EXPECT_EQ(record.selectionType, SelectionRecord::STOICHIOMETRY); + vector ids = rr.getStoichiometryIds(); + EXPECT_NE(std::find(ids.begin(), ids.end(), "N"), ids.end()); +} + +TEST_F(SBMLFeatures, variable_stoich_event_selection_type) { + RoadRunner rr((SBMLFeaturesDir / "stoich_event.xml").string()); + SelectionRecord record = rr.createSelection("N"); + EXPECT_EQ(record.selectionType, SelectionRecord::STOICHIOMETRY); + vector ids = rr.getStoichiometryIds(); + EXPECT_NE(std::find(ids.begin(), ids.end(), "N"), ids.end()); +} diff --git a/wrappers/Python/roadrunner/roadrunner.i b/wrappers/Python/roadrunner/roadrunner.i index 9a940cd0f2..db246acce7 100644 --- a/wrappers/Python/roadrunner/roadrunner.i +++ b/wrappers/Python/roadrunner/roadrunner.i @@ -1912,7 +1912,8 @@ namespace std { class ostream{}; } SelectionRecord.FLOATING | SelectionRecord.BOUNDARY | SelectionRecord.COMPARTMENT | - SelectionRecord.GLOBAL_PARAMETER) + SelectionRecord.GLOBAL_PARAMETER | + SelectionRecord.STOICHIOMETRY) def resetParameter(self): """ Reset parameters to CURRENT init(X) values.