Get variable stoichiometries fully working

.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" ]

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}")
 

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

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);

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));

source/llvm/Jit.cpp

@@ -135,6 +135,8 @@ namespace rrllvm {
                 "getPiecewiseTrigger"));
         modelResources->evalVolatileStoichPtr = reinterpret_cast<EvalVolatileStoichCodeGen::FunctionPtr>(lookupFunctionAddress(
                 "evalVolatileStoich"));
+        modelResources->resetRateRuleStoichPtr = reinterpret_cast<ResetRateRuleStoichCodeGen::FunctionPtr>(lookupFunctionAddress(
+                "resetRateRuleStoich"));
         modelResources->evalConversionFactorPtr = reinterpret_cast<EvalConversionFactorCodeGen::FunctionPtr>(lookupFunctionAddress(
                 "evalConversionFactor"));
         if (options & LoadSBMLOptions::READ_ONLY) {

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<std::string> &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<LLVMModelDataSymbols::SpeciesReferenceInfo>::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<LLVMModelDataSymbols::SpeciesReferenceInfo>::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<LLVMModelDataSymbols::SpeciesReferenceInfo> stoichiometryIndx = symbols->getStoichiometryList();
     std::list<LLVMModelDataSymbols::SpeciesReferenceInfo>::const_iterator stoichiometry = stoichiometryIndx.begin();
     for (int i = 0; i < index; i++)

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.

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);

source/llvm/LLVMModelDataSymbols.h

@@ -710,6 +710,13 @@ class LLVMModelDataSymbols
 
     std::vector<SpeciesReferenceType> 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.

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();
+                }
               }
             }
 

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;