Fixes for named stoichiometries

source/llvm/LLVMExecutableModel.cpp

@@ -1246,7 +1246,7 @@ void LLVMExecutableModel::getIds(int types, std::list<std::string> &ids)
                 if (getStoichiometry(s, r) != 0)
                 {
                     string rid = getReactionId(r);
-                    ids.push_back("stoich(" + sid + ", " + rid + ")");
+                    ids.push_back(symbols->getStoichiometryIdFor(sid, rid));
                 }
             }
         }
@@ -1603,6 +1603,10 @@ const rr::SelectionRecord& LLVMExecutableModel::getSelection(const std::string&
                 sel.selectionType = SelectionRecord::INITIAL_GLOBAL_PARAMETER;
                 sel.index = index;
                 break;
+            case LLVMModelDataSymbols::STOICHIOMETRY:
+                sel.selectionType = SelectionRecord::INITIAL_STOICHIOMETRY;
+                sel.index = index;
+                break;
             default:
                 std::string msg = "Invalid Id for initial value: '" + str + "'";
                 throw LLVMException(msg);
@@ -1744,7 +1748,8 @@ void LLVMExecutableModel::setValue(const std::string& id, double value)
         setGlobalParameterInitValues(1, &index, &value);
         break;
     case SelectionRecord::STOICHIOMETRY:
-        setStoichiometries(1, &index, &value);
+    case SelectionRecord::INITIAL_STOICHIOMETRY:
+        setStoichiometry(index, value);
         break;
     default:
         throw LLVMException("Invalid selection '" + sel.to_string() + "' for setting value");
@@ -2319,25 +2324,73 @@ int LLVMExecutableModel::setCompartmentVolumes(size_t len, const int* indx,
     return result;
 }
 
-int LLVMExecutableModel::setStoichiometries(size_t len, const int* indx,
-                                               const double* values)
-{
-    return setStoichiometries(len, indx, values, true);
-}
+//int LLVMExecutableModel::setStoichiometries(size_t len, const int* indx,
+//                                               const double* values)
+//{
+//    return setStoichiometries(len, indx, values, true);
+//}
+
+//int LLVMExecutableModel::setStoichiometries(size_t len, const int* indx,
+//                                               const double* values, bool strict)
+//{
+//    if (len == 1) {
+//        int index = *indx;
+//        double value = *values;
+//        return setStoichiometry(index, value);
+//    }
+//
+//    return -1;
+//}
 
-int LLVMExecutableModel::setStoichiometries(size_t len, const int* indx,
-                                               const double* values, bool strict)
+int LLVMExecutableModel::setStoichiometry(int index, double value)
 {
-    if (len == 1) {
-        int index = *indx;
-        double value = *values;
-        return setStoichiometry(index, value);
-    }
+    if (std::signbit(value))
+        throw LLVMException("Invalid stoichiometry value");
+
+    if (symbols->isConservedMoietyAnalysis())
+        throw LLVMException("Unable to set stoichiometries when conserved moieties are on");
+
+    std::list <LLVMModelDataSymbols::SpeciesReferenceInfo> stoichiometryIndx = symbols->getStoichiometryList();
+    std::list<LLVMModelDataSymbols::SpeciesReferenceInfo>::const_iterator stoichiometry = stoichiometryIndx.begin();
+    for (int i = 0; i < index; i++)
+        ++stoichiometry;
+    if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Product)
+        return setStoichiometry(stoichiometry->row, stoichiometry->column, value);
+    else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant)
+        return setStoichiometry(stoichiometry->row, stoichiometry->column, -1 * value);
+    else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::MultiReactantProduct)
+        throw LLVMException("Cannot set stoichiometry for a MultiReactantProduct");
+    else
+        throw LLVMException("Cannot set stoichiometry for a Modifier");
 
     return -1;
 }
 
-int LLVMExecutableModel::setStoichiometry(int index, double value)
+int LLVMExecutableModel::setStoichiometry(int speciesIndex, int reactionIndex, double value)
+{
+    double result = csr_matrix_set_nz(modelData->stoichiometry, speciesIndex, reactionIndex, value);
+    return isnan(result) ? 0 : result;
+}
+
+//int LLVMExecutableModel::setInitStoichiometries(size_t len, const int* indx,
+//                                               const double* values)
+//{
+//    return setInitStoichiometries(len, indx, values, true);
+//}
+
+//int LLVMExecutableModel::setInitStoichiometries(size_t len, const int* indx,
+//                                               const double* values, bool strict)
+//{
+//    if (len == 1) {
+//        int index = *indx;
+//        double value = *values;
+//        return setInitStoichiometry(index, value);
+//    }
+//
+//    return -1;
+//}
+
+int LLVMExecutableModel::setInitStoichiometry(int index, double value)
 {
     if (std::signbit(value))
         throw LLVMException("Invalid stoichiometry value");
@@ -2350,9 +2403,9 @@ int LLVMExecutableModel::setStoichiometry(int index, double value)
     for (int i = 0; i < index; i++)
         ++stoichiometry;
     if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Product)
-        return setStoichiometry(stoichiometry->row, stoichiometry->column, value);
+        return setInitStoichiometry(stoichiometry->row, stoichiometry->column, value);
     else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant)
-        return setStoichiometry(stoichiometry->row, stoichiometry->column, -1 * value);
+        return setInitStoichiometry(stoichiometry->row, stoichiometry->column, -1 * value);
     else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::MultiReactantProduct)
         throw LLVMException("Cannot set stoichiometry for a MultiReactantProduct");
     else
@@ -2361,8 +2414,10 @@ int LLVMExecutableModel::setStoichiometry(int index, double value)
     return -1;
 }
 
-int LLVMExecutableModel::setStoichiometry(int speciesIndex, int reactionIndex, double value)
+int LLVMExecutableModel::setInitStoichiometry(int speciesIndex, int reactionIndex, double value)
 {
+    //For now, we don't store 'initStoichiometry' separately.  Essentially, every stoichiometry set is to the initial value.
+    //double result = csr_matrix_set_nz(modelData->initStoichiometry, speciesIndex, reactionIndex, value);
     double result = csr_matrix_set_nz(modelData->stoichiometry, speciesIndex, reactionIndex, value);
     return isnan(result) ? 0 : result;
 }
@@ -2394,6 +2449,35 @@ double LLVMExecutableModel::getStoichiometry(int speciesIndex, int reactionIndex
     return isnan(result) ? 0 : result;
 }
 
+double LLVMExecutableModel::getInitStoichiometry(int index)
+{
+    if (symbols->isConservedMoietyAnalysis())
+        throw LLVMException("Unable to get stoichiometries when conserved moieties are on");
+
+    if (index < 0)
+        throw LLVMException("The stoichiometry index is not valid");
+    std::list<LLVMModelDataSymbols::SpeciesReferenceInfo> stoichiometryIndx = symbols->getStoichiometryList();
+    std::list<LLVMModelDataSymbols::SpeciesReferenceInfo>::const_iterator stoichiometry = stoichiometryIndx.begin();
+    for (int i = 0; i < index; i++)
+        ++stoichiometry;
+    if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Reactant)
+        return -1 * getInitStoichiometry(stoichiometry->row, stoichiometry->column);
+    else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Product)
+        return getInitStoichiometry(stoichiometry->row, stoichiometry->column);
+    else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::MultiReactantProduct)
+        throw LLVMException("Cannot return stoichiometry for a MultiReactantProduct");
+    else
+        throw LLVMException("Cannot return stoichiometry for a Modifier");
+}
+
+double LLVMExecutableModel::getInitStoichiometry(int speciesIndex, int reactionIndex)
+{
+    //For now, we don't store 'initStoichiometry' separately.  Essentially, every stoichiometry set is to the initial value.
+    //double result = csr_matrix_get_nz(modelData->initStoichiometry, speciesIndex, reactionIndex);
+    double result = csr_matrix_get_nz(modelData->stoichiometry, speciesIndex, reactionIndex);
+    return isnan(result) ? 0 : result;
+}
+
 /******************************* Events Section *******************************/
 #if (1) /**********************************************************************/
 /******************************************************************************/

source/llvm/LLVMExecutableModel.h

@@ -348,20 +348,34 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel
     virtual int setCompartmentVolumes(size_t len, int const* indx,
         const double* values, bool strict);
 
-    virtual int setStoichiometries(size_t len, int const* indx,
-                                      const double* values);
+    //virtual int setStoichiometries(size_t len, int const* indx,
+    //                                  const double* values);
 
-    virtual int setStoichiometries(size_t len, int const* indx,
-                                      const double* values, bool strict);
+    //virtual int setStoichiometries(size_t len, int const* indx,
+    //                                  const double* values, bool strict);
 
     virtual int setStoichiometry(int index, double value);
 
     virtual int setStoichiometry(int speciesIndex, int reactionIndex, double value);
 
+    //virtual int setInitStoichiometries(size_t len, int const* indx,
+    //  const double* values);
+
+    //virtual int setInitStoichiometries(size_t len, int const* indx,
+    //  const double* values, bool strict);
+
+    virtual int setInitStoichiometry(int index, double value);
+
+    virtual int setInitStoichiometry(int speciesIndex, int reactionIndex, double value);
+
     virtual double getStoichiometry(int index);
 
     virtual double getStoichiometry(int speciesIndex, int reactionIndex);
 
+    virtual double getInitStoichiometry(int index);
+
+    virtual double getInitStoichiometry(int speciesIndex, int reactionIndex);
+
     /******************************* Initial Conditions Section *******************/
     #if (1) /**********************************************************************/
     /******************************************************************************/

source/llvm/LLVMModelData.cpp

@@ -125,6 +125,7 @@ void LLVMModelData_save(LLVMModelData *data, std::ostream& out)
 	rr::saveBinary(out, data->stateVectorSize);
 	//Save the stoichiometry matrix
 	rr::csr_matrix_dump_binary(data->stoichiometry, out);
+  //rr::csr_matrix_dump_binary(data->initStoichiometry, out);
 
 	//We do not need to save random because LLVMExecutableModel will make a new one if it is null
 
@@ -205,6 +206,7 @@ LLVMModelData* LLVMModelData_from_save(std::istream& in)
 
 	//Load the stoichiometry matrix
 	data->stoichiometry = rr::csr_matrix_new_from_binary(in);
+  //data->initStoichiometry = rr::csr_matrix_new_from_binary(in);
 
 	//Alias pointer offsets
 

source/llvm/LLVMModelData.h

@@ -223,25 +223,30 @@ struct LLVMModelData
      */
     double*                             floatingSpeciesAmountsAlias;      // 31
 
+    /*
+     * initial stoichiometry matrix
+     */
+    //rr::csr_matrix*                     initStoichiometry;                // 32
+
     /**
 	 * binary data layout:
      *
-     * compartmentVolumes                [numIndCompartmentVolumes]       // 32
-     * initCompartmentVolumes            [numInitCompartmentVolumes]      // 33
-     * initFloatingSpeciesAmounts        [numInitFloatingSpecies]         // 34
-     * boundarySpeciesAmounts            [numIndBoundarySpecies]          // 35
-     * initBoundarySpeciesAmounts        [numInitBoundarySpecies]         // 36
-     * globalParameters                  [numIndGlobalParameters]         // 37
-     * initGlobalParameters              [numInitGlobalParameters]        // 38
-     * reactionRates                     [numReactions]                   // 39
+     * compartmentVolumes                [numIndCompartmentVolumes]
+     * initCompartmentVolumes            [numInitCompartmentVolumes]
+     * initFloatingSpeciesAmounts        [numInitFloatingSpecies]
+     * boundarySpeciesAmounts            [numIndBoundarySpecies]
+     * initBoundarySpeciesAmounts        [numInitBoundarySpecies]
+     * globalParameters                  [numIndGlobalParameters]
+     * initGlobalParameters              [numInitGlobalParameters]
+     * reactionRates                     [numReactions]
      *
-     * rateRuleValues                    [numRateRules]                   // 40
-     * floatingSpeciesAmounts            [numIndFloatingSpecies]          // 41
+     * rateRuleValues                    [numRateRules]
+     * floatingSpeciesAmounts            [numIndFloatingSpecies]
      */
 
 	 /**
 	 * This dynamic-sized array will be allocated while this ModelData is allocated.
-	 * Ten array in the permanent data section is stored contiously in this chunck.
+	 * Ten array in the permanent data section is stored contiously in this chunk.
 	 * Size of each array is defined by ten unsigned integer above.
 	 * Values can be accessed using ten alias pointers defined above.
 	 * 

source/llvm/LLVMModelDataSymbols.cpp

@@ -402,6 +402,19 @@ int LLVMModelDataSymbols::getStoichiometryIndex(const std::string& speciesId, co
     return -1;
 }
 
+std::string LLVMModelDataSymbols::getStoichiometryIdFor(const std::string& speciesId, const std::string& reactionId) const
+{
+  int index = getStoichiometryIndex(speciesId, reactionId);
+  if (index == -1) {
+    return "";
+  }
+  string stoichid = stoichIds[index];
+  if (stoichid == "") {
+    return "stoich(" + speciesId + ", " + reactionId + ")";
+  }
+  return stoichid;
+}
+
 std::vector<std::string> LLVMModelDataSymbols::getStoichiometryIds() const
 {
     return stoichIds;
@@ -1328,6 +1341,7 @@ typedef cxx11_ns::unordered_map<uint, ssize_type> UIntUMap;
 void LLVMModelDataSymbols::initReactions(const libsbml::Model* model)
 {
     // get the reactions
+    std::vector<std::string> namedStoichiometryIds;
     const ListOfReactions *reactions = model->getListOfReactions();
     for (size_t i = 0; i < reactions->size(); i++)
     {
@@ -1357,6 +1371,7 @@ void LLVMModelDataSymbols::initReactions(const libsbml::Model* model)
                 SpeciesReferenceInfo info =
                 { static_cast<uint>(speciesIdx), static_cast<uint>(i), Reactant, r->getId() };
                 namedSpeciesReferenceInfo[r->getId()] = info;
+                namedStoichiometryIds.push_back(r->getId());
               }
               else
               {
@@ -1419,6 +1434,7 @@ void LLVMModelDataSymbols::initReactions(const libsbml::Model* model)
                 SpeciesReferenceInfo info =
                 { static_cast<uint>(speciesIdx), static_cast<uint>(i), Product, p->getId() };
                 namedSpeciesReferenceInfo[p->getId()] = info;
+                namedStoichiometryIds.push_back(p->getId());
               }
               else
               {
@@ -1462,6 +1478,22 @@ 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);
+      }
+    }
+
 }
 
 

source/llvm/LLVMModelDataSymbols.h

@@ -248,6 +248,7 @@ class LLVMModelDataSymbols
 
     int getStoichiometryIndex(std::string const&) const;
     int getStoichiometryIndex(const std::string&, const std::string&) const;
+    std::string getStoichiometryIdFor(const std::string&, const std::string&) const;
     std::vector<std::string> getStoichiometryIds() const;
     size_t getStoichiometrySize() const;
 

source/llvm/LLVMModelGenerator.cpp

@@ -286,10 +286,10 @@ namespace rrllvm {
     {
         SharedModelResourcesPtr modelResources = std::make_shared<ModelResources>();
 
-        char* docSBML = doc->toSBML();
 
         ModelGeneratorContext modelGeneratorContext(doc, options, JitFactory::makeJitEngine(options));
 
+        char* docSBML = doc->toSBML();
         std::string sbmlMD5 = getSBMLMD5(std::string((const char*)docSBML), options);
         if (modelResources->sbmlMD5.empty()) {
             modelResources->sbmlMD5 = sbmlMD5;
@@ -381,6 +381,23 @@ namespace rrllvm {
                     }
                 }
             }
+
+            //Copy any changed stoichiometries.  Make sure we know if we're in 'conserved moieties' mode.
+            for (int r = 0; r < oldModel->getNumReactions(); r++) {
+              string rID = oldModel->getReactionId(r);
+              int new_r = newModel->getReactionIndex(rID);
+              if (new_r < 0) continue;
+              int s_max = oldModel->getNumIndFloatingSpecies();
+              for (int s = 0; s < s_max; s++) {
+                string sID = oldModel->getFloatingSpeciesId(s);
+                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);
+              }
+            }
+
             newModel->setTime(oldModel->getTime());
         }
 

source/llvm/ModelDataIRBuilder.cpp

@@ -744,7 +744,7 @@ llvm::StructType *ModelDataIRBuilder::getStructType(llvm::Module *module)
 unsigned ModelDataIRBuilder::getModelDataSize(llvm::Module *module, const DataLayout& dl)
 {
     StructType *structType = getStructType(module);
-	uint64_t llvm_size = dl.getTypeStoreSize(structType);
+    uint64_t llvm_size = dl.getTypeStoreSize(structType);
 
     // the model data struct will NEVER be bigger than a 32 bit pointer!
     return (unsigned)llvm_size;