add getExtendedStoichiometryMatrix function, fix #509

sys-bio · Dec 21, 2018 · df5535a · df5535a
1 parent 138098b
commit df5535a
Show file tree

Hide file tree

Showing 2 changed files with 126 additions and 0 deletions.
diff --git a/source/rrRoadRunner.cpp b/source/rrRoadRunner.cpp
@@ -2220,6 +2220,130 @@ DoubleMatrix RoadRunner::getFullStoichiometryMatrix()
     return m;
 }
 
+
+DoubleMatrix RoadRunner::getExtendedStoichiometryMatrix()
+{
+    check_model();
+    get_self();
+    ls::LibStructural *ls = getLibStruct();
+
+    if (self.loadOpt.getConservedMoietyConversion()) {
+        // pointer to mat owned by ls
+        DoubleMatrix m = *ls->getReorderedStoichiometryMatrix();
+        ls->getReorderedStoichiometryMatrixLabels(
+                m.getRowNames(), m.getColNames());
+        return m;
+    }
+
+    // pointer to owned matrix
+    DoubleMatrix *mptr = ls->getStoichiometryMatrix();
+    if (!mptr)
+    {
+        throw CoreException("Error: Stoichiometry matrix does not exist for this model");
+    }
+    DoubleMatrix m = *mptr;
+    ls->getStoichiometryMatrixLabels(m.getRowNames(), m.getColNames());
+
+    libsbml::SBMLReader reader;
+    libsbml::SBMLDocument *doc = reader.readSBMLFromString(self.mCurrentSBML);
+    libsbml::Model *model = doc->getModel();
+    typedef cxx11_ns::unordered_map<int,int> RxnIdxToRowMap;
+    typedef cxx11_ns::unordered_map<int,libsbml::Reaction*> RxnIdxToRxnMap;
+    typedef cxx11_ns::unordered_map<libsbml::Species*,int> SpcToRowMap;
+    RxnIdxToRowMap missing_rct_row_map;
+    RxnIdxToRowMap missing_prd_row_map;
+    RxnIdxToRxnMap rxn_map;
+    SpcToRowMap boundary_spc_row_map;
+
+    int n_rows = m.numRows();
+
+    for (int i=0; i<m.getColNames().size(); ++i) {
+      libsbml::Reaction* r = model->getReaction(m.getColNames().at(i));
+      assert(r);
+      rxn_map[i] = r;
+      if (!r->getNumReactants()) {
+        missing_rct_row_map[i] = n_rows++;
+      } else {
+        for (int j=0; j<r->getNumReactants(); ++j) {
+          libsbml::SpeciesReference* ref = r->getReactant(j);
+          assert(ref);
+          libsbml::Species* s = model->getSpecies(ref->getSpecies());
+          assert(s);
+          if (s->getBoundaryCondition() && boundary_spc_row_map.find(s) == boundary_spc_row_map.end())
+            boundary_spc_row_map[s] = n_rows++;
+        }
+      }
+      if (!r->getNumProducts()) {
+        missing_prd_row_map[i] = n_rows++;
+      } else {
+        for (int j=0; j<r->getNumProducts(); ++j) {
+          libsbml::SpeciesReference* ref = r->getProduct(j);
+          assert(ref);
+          libsbml::Species* s = model->getSpecies(ref->getSpecies());
+          assert(s);
+          if (s->getBoundaryCondition() && boundary_spc_row_map.find(s) == boundary_spc_row_map.end())
+            boundary_spc_row_map[s] = n_rows++;
+        }
+      }
+    }
+
+    DoubleMatrix extended_matrix(n_rows, m.numCols());
+    extended_matrix.setRowNames(m.getRowNames());
+    extended_matrix.setColNames(m.getColNames());
+    extended_matrix.getRowNames().resize(n_rows);
+    for (int i=0; i<m.numRows(); ++i) {
+      for (int j=0; j<m.numCols(); ++j) {
+        extended_matrix(i,j) = m(i,j);
+      }
+    }
+    for (int i=m.numRows(); i<n_rows; ++i) {
+      for (int j=0; j<m.numCols(); ++j) {
+        extended_matrix(i,j) = 0;
+      }
+    }
+    for (RxnIdxToRowMap::const_iterator i=missing_rct_row_map.begin(); i!=missing_rct_row_map.end(); ++i){
+      extended_matrix(i->second,i->first) = -1;
+      RxnIdxToRxnMap::const_iterator r = rxn_map.find(i->first);
+      if (r != rxn_map.end())
+        extended_matrix.getRowNames().at(i->second) = r->second->getId() + "_source";
+    }
+    for (RxnIdxToRowMap::const_iterator i=missing_prd_row_map.begin(); i!=missing_prd_row_map.end(); ++i){
+      extended_matrix(i->second,i->first) = 1;
+      RxnIdxToRxnMap::const_iterator r = rxn_map.find(i->first);
+      if (r != rxn_map.end())
+        extended_matrix.getRowNames().at(i->second) = r->second->getId() + "_sink";
+    }
+    for (SpcToRowMap::const_iterator i=boundary_spc_row_map.begin(); i!=boundary_spc_row_map.end(); ++i) {
+      for (int j=0; j<m.getColNames().size(); ++j) {
+        libsbml::Reaction* r = model->getReaction(m.getColNames().at(j));
+        assert(r);
+        for (int k=0; k<r->getNumReactants(); ++k) {
+          libsbml::SpeciesReference* ref = r->getReactant(k);
+          assert(ref);
+          libsbml::Species* s = model->getSpecies(ref->getSpecies());
+          assert(s);
+          if (s == i->first) {
+            extended_matrix(i->second,j) = -1;
+            extended_matrix.getRowNames().at(i->second) = s->getId();
+          }
+        }
+        for (int k=0; k<r->getNumProducts(); ++k) {
+          libsbml::SpeciesReference* ref = r->getProduct(k);
+          assert(ref);
+          libsbml::Species* s = model->getSpecies(ref->getSpecies());
+          assert(s);
+          if (s == i->first) {
+            extended_matrix(i->second,j) = 1;
+            extended_matrix.getRowNames().at(i->second) = s->getId();
+          }
+        }
+      }
+    }
+
+    return extended_matrix;
+}
+
+
 DoubleMatrix RoadRunner::getL0Matrix()
 {
     check_model();

diff --git a/source/rrRoadRunner.h b/source/rrRoadRunner.h
@@ -543,6 +543,8 @@ class RR_DECLSPEC RoadRunner
      */
     ls::DoubleMatrix getFullStoichiometryMatrix();
 
+    ls::DoubleMatrix getExtendedStoichiometryMatrix();
+
 
     ls::DoubleMatrix getL0Matrix();