Merge pull request #1145 from sys-bio/adding_get_set_stoichiometry_va…

…lues_feature

The feature to get and set parametric stoichiometries is added

source/CVODEIntegrator.h

@@ -243,7 +243,7 @@ namespace rr {
          * 
          * The absolute tolerance vector is either set directly by the user as a vector, or is generated from the single absolute tolerance value (either default or set by the user) multiplied by the initial value of every element in the state vector (independent floating species, and elements controlled by rate rules).  If that initial value is zero, the corresponding element of the absolute tolerance vector is the single value multiplied by the compartment volume for species, or by one for all other values.
          */
-        virtual std::vector<double> getAbsoluteToleranceVector();
+        std::vector<double> getAbsoluteToleranceVector() override;
 
 
     private:

source/llvm/EvalInitialConditionsCodeGen.cpp

@@ -189,7 +189,7 @@ void EvalInitialConditionsCodeGen::codeGenStoichiometry(
     Value *stoich = builder.CreateLoad(stoichEP, "stoichiometry");
 
     std::list<LLVMModelDataSymbols::SpeciesReferenceInfo> stoichEntries =
-            dataSymbols.getStoichiometryIndx();
+            dataSymbols.getStoichiometryList();
 
     for (std::list<LLVMModelDataSymbols::SpeciesReferenceInfo>::iterator i =
             stoichEntries.begin(); i != stoichEntries.end(); i++)

source/llvm/LLVMExecutableModel.cpp

@@ -678,6 +678,30 @@ std::string LLVMExecutableModel::getReactionId(size_t id)
     }
 }
 
+int LLVMExecutableModel::getStoichiometryIndex(const std::string& id)
+{
+    return symbols->getStoichiometryIndex(id);
+}
+
+int LLVMExecutableModel::getStoichiometryIndex(const std::string& speciesId, const std::string& reactionId)
+{
+    return symbols->getStoichiometryIndex(speciesId, reactionId);
+}
+
+std::string LLVMExecutableModel::getStoichiometryId(size_t id)
+{
+    std::vector<std::string> ids = symbols->getStoichiometryIds();
+    if (id < ids.size())
+    {
+        return ids[id];
+    }
+    else
+    {
+        throw_llvm_exception("index out of range");
+        return "";
+    }
+}
+
 void LLVMExecutableModel::evalInitialConditions(uint32_t flags)
 {
     evalInitialConditionsPtr(modelData, flags);
@@ -1514,6 +1538,10 @@ const rr::SelectionRecord& LLVMExecutableModel::getSelection(const std::string&
                 sel.selectionType = SelectionRecord::EVENT;
                 sel.index = index;
                 break;
+            case LLVMModelDataSymbols::STOICHIOMETRY:
+                sel.selectionType = SelectionRecord::STOICHIOMETRY;
+                sel.index = index;
+                break;
             default:
                 throw LLVMException("No sbml element exists for symbol '" + str + "'");
                 break;
@@ -1596,6 +1624,9 @@ const rr::SelectionRecord& LLVMExecutableModel::getSelection(const std::string&
                 break;
             }
             break;
+        case SelectionRecord::STOICHIOMETRY:
+            sel.index = getStoichiometryIndex(sel.p1, sel.p2);
+            break;
 
         default:
             rrLog(Logger::LOG_ERROR) << "A new SelectionRecord should not have this value: "
@@ -1702,6 +1733,9 @@ void LLVMExecutableModel::setValue(const std::string& id, double value)
     case SelectionRecord::INITIAL_GLOBAL_PARAMETER:
         setGlobalParameterInitValues(1, &index, &value);
         break;
+    case SelectionRecord::STOICHIOMETRY:
+        setStoichiometries(1, &index, &value);
+        break;
     default:
         throw LLVMException("Invalid selection '" + sel.to_string() + "' for setting value");
         break;
@@ -2270,7 +2304,74 @@ 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, bool strict)
+{
+    if (len == 1) {
+        int index = *indx;
+        double value = *values;
+        return setStoichiometry(index, value);
+    }
+
+    return -1;
+}
+
+int LLVMExecutableModel::setStoichiometry(int index, double 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 speciesIndex, int reactionIndex, double value)
+{
+    double result = csr_matrix_set_nz(modelData->stoichiometry, speciesIndex, reactionIndex, value);
+    return isnan(result) ? 0 : result;
+}
+
+double LLVMExecutableModel::getStoichiometry(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 * getStoichiometry(stoichiometry->row, stoichiometry->column);
+    else if (stoichiometry->type == LLVMModelDataSymbols::SpeciesReferenceType::Product)
+        return getStoichiometry(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::getStoichiometry(int speciesIndex, int reactionIndex)
 {

source/llvm/LLVMExecutableModel.h

@@ -328,6 +328,9 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel
     virtual std::string getCompartmentId(size_t);
     virtual int getReactionIndex(const std::string&);
     virtual std::string getReactionId(size_t);
+    virtual int getStoichiometryIndex(const std::string&);
+    virtual int getStoichiometryIndex(const std::string& speciesId, const std::string& reactionId);
+    virtual std::string getStoichiometryId(size_t);
 
     virtual void print(std::ostream &stream);
 
@@ -338,13 +341,23 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel
     virtual int setConservedMoietyValues(size_t len, int const *indx,
             const double *values);
 
-
     virtual int setCompartmentVolumes(size_t len, int const* indx,
         const double* values);
 
     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, bool strict);
+
+    virtual int setStoichiometry(int index, double value);
+
+    virtual int setStoichiometry(int speciesIndex, int reactionIndex, double value);
+
+    virtual double getStoichiometry(int index);
 
     virtual double getStoichiometry(int speciesIndex, int reactionIndex);
 

source/llvm/LLVMModelDataSymbols.cpp

@@ -266,6 +266,15 @@ LLVMModelDataSymbols::SymbolIndexType LLVMModelDataSymbols::getSymbolIndex(
         result = i->second;
         return EVENT;
     }
+    else
+    {
+        for (int stoichIndex = 0 ; stoichIndex < stoichIds.size(); ++stoichIndex) {
+            if (stoichIds.at(stoichIndex) == name) {
+                result = stoichIndex;
+                return STOICHIOMETRY;
+            }
+        }
+    }
 
     result = -1;
     return INVALID_SYMBOL;
@@ -405,13 +414,51 @@ size_t LLVMModelDataSymbols::getReactionSize() const
     return reactionsMap.size();
 }
 
+int LLVMModelDataSymbols::getStoichiometryIndex(const std::string& id) const
+{
+    for (int i = 0; i < stoichIds.size(); ++i)
+    {
+        if (stoichIds.at(i) == id)
+            return i;
+    }
+
+    return -1;
+}
+
+int LLVMModelDataSymbols::getStoichiometryIndex(const std::string& speciesId, const std::string& reactionId) const
+{
+    std::list<LLVMModelDataSymbols::SpeciesReferenceInfo> stoichiometryIndx = getStoichiometryList();
+    int speciesIndex = getFloatingSpeciesIndex(speciesId);
+    int reactionIndex = getReactionIndex(reactionId);
+    std::list<LLVMModelDataSymbols::SpeciesReferenceInfo>::const_iterator stoichiometry = stoichiometryIndx.begin();
+    int stoichiometryIndex = 0;
+    while (stoichiometry != stoichiometryIndx.end()) {
+        if (stoichiometry->row == speciesIndex && stoichiometry->column == reactionIndex)
+            return stoichiometryIndex;
+        ++stoichiometry;
+        ++stoichiometryIndex;
+    }
+
+    return -1;
+}
+
+std::vector<std::string> LLVMModelDataSymbols::getStoichiometryIds() const
+{
+    return stoichIds;
+}
+
+size_t LLVMModelDataSymbols::getStoichiometrySize() const
+{
+    return getStoichiometryList().size();
+}
+
 size_t  LLVMModelDataSymbols::getFloatingSpeciesSize() const
 {
     return floatingSpeciesMap.size();
 }
 
 std::list<LLVMModelDataSymbols::SpeciesReferenceInfo>
-    LLVMModelDataSymbols::getStoichiometryIndx() const
+    LLVMModelDataSymbols::getStoichiometryList() const
 {
     std::list<SpeciesReferenceInfo> result;
 
@@ -1436,6 +1483,10 @@ void LLVMModelDataSymbols::initReactions(const libsbml::Model* model)
                     stoichIds.push_back(p->isSetId() ? p->getId() : "");
                     stoichTypes.push_back(Product);
 
+                    // in case this species is both a product and reactant, can look up
+                    // index of the just added Reactant
+                    speciesMap[speciesIdx] = stoichTypes.size() - 1;
+
                     if (p->isSetId() && p->getId().length() > 0)
                     {
                         if (namedSpeciesReferenceInfo.find(p->getId())
@@ -1883,6 +1934,13 @@ int LLVMModelDataSymbols::getConservedMoietyIndex(
     return -1;
 }
 
+bool LLVMModelDataSymbols::isConservedMoietyAnalysis() const {
+    if (conservedMoietyGlobalParameter.size())
+        return true;
+
+    return false;
+}
+
 void LLVMModelDataSymbols::saveState(std::ostream& out) const
 {
 	rr::saveBinary(out, conservedMoietySpeciesSet);

source/llvm/LLVMModelDataSymbols.h

@@ -151,6 +151,7 @@ class LLVMModelDataSymbols
         GLOBAL_PARAMETER,
         REACTION,
         EVENT,
+        STOICHIOMETRY,
         INVALID_SYMBOL
     };
 
@@ -244,6 +245,11 @@ class LLVMModelDataSymbols
     std::vector<std::string> getReactionIds() const;
     size_t getReactionSize() const;
 
+    int getStoichiometryIndex(std::string const&) const;
+    int getStoichiometryIndex(const std::string&, const std::string&) const;
+    std::vector<std::string> getStoichiometryIds() const;
+    size_t getStoichiometrySize() const;
+
 
     std::vector<std::string> getGlobalParameterIds() const;
 
@@ -310,7 +316,7 @@ class LLVMModelDataSymbols
      * in the list of pairs, first is the row (species) index,
      * and second is the column (reaction) index.
      */
-    std::list<SpeciesReferenceInfo> getStoichiometryIndx() const;
+    std::list<SpeciesReferenceInfo> getStoichiometryList() const;
 
     /**
      * initialize and allocate the buffers (including the stoich matrix)
@@ -422,6 +428,11 @@ class LLVMModelDataSymbols
      */
     int getConservedMoietyIndex(const std::string& name) const;
 
+    /**
+     * check if the conserved moiety is turned on for this model
+     */
+    bool isConservedMoietyAnalysis() const;
+
 private:
 
     /**