Merge pull request #1158 from sys-bio/fix-1135

Fix 1135
sys-bio · Nov 2, 2023 · 5ac0b78 · 5ac0b78
2 parents 8f77756 + 5f78f22
commit 5ac0b78
Show file tree

Hide file tree

Showing 10 changed files with 94 additions and 111 deletions.
diff --git a/source/GillespieIntegrator.cpp b/source/GillespieIntegrator.cpp
@@ -61,23 +61,13 @@ namespace rr
 
         assert(floatingSpeciesStart >= 0);
 
-        // get rows and columns
-        mModel->getStoichiometryMatrix(&stoichRows, &stoichCols, nullptr);
-        stoichData = new double[stoichRows * stoichCols];
-
-        // fill stoichData
-        mModel->getStoichiometryMatrix(&stoichRows, &stoichCols, &stoichData);
-
         setEngineSeed(getValue("seed").get<std::uint64_t>());
 	}
 
 	GillespieIntegrator::GillespieIntegrator(ExecutableModel* m)
 	    :   Integrator(m),
 			timeScale(1.0),
 			stoichScale(1.0),
-			stoichRows(0),
-			stoichCols(0),
-			stoichData(nullptr),
 			reactionRates(nullptr),
 			reactionRatesBuffer(nullptr),
 			stateVector(nullptr),
@@ -98,12 +88,10 @@ namespace rr
             delete[] reactionRatesBuffer;
             delete[] stateVector;
             delete[] stateVectorRate;
-            delete[] stoichData;
             reactionRates = nullptr;
             reactionRatesBuffer = nullptr;
             stateVector = nullptr;
             stateVectorRate = nullptr;
-            stoichData = nullptr;
         }
 	}
 
@@ -114,12 +102,10 @@ namespace rr
         delete[] reactionRatesBuffer;
         delete[] stateVector;
         delete[] stateVectorRate;
-        delete[] stoichData;
         reactionRates = nullptr;
         reactionRatesBuffer = nullptr;
         stateVector = nullptr;
         stateVectorRate = nullptr;
-        stoichData = nullptr;
 
         mModel = m;
         mModel->reset();
@@ -130,10 +116,6 @@ namespace rr
         timeScale = 1.;
         stoichScale = 1.;
 
-        stoichRows = 0;
-        stoichCols = 0;
-        stoichData = nullptr;
-
         initializeFromModel();
     }
 

diff --git a/source/GillespieIntegrator.h b/source/GillespieIntegrator.h
@@ -133,11 +133,6 @@ namespace rr
         int stateVectorSize;
         double* stateVector;
         double* stateVectorRate;
-        // m rows x n cols
-        // offset = row*NUMCOLS + column
-        int stoichRows;
-        int stoichCols;
-        double* stoichData;
         std::vector<unsigned char> eventStatus;
         std::vector<unsigned char> previousEventStatus;
 
@@ -146,11 +141,10 @@ namespace rr
 
         double urand();
         void setEngineSeed(std::uint64_t seed);
-        std::uint64_t getSeed() const;
 
         inline double getStoich(uint species, uint reaction)
         {
-            return stoichData[species * stoichCols + reaction];
+            return mModel->getStoichiometry(species, reaction);
         }
 
         /**

diff --git a/source/llvm/LLVMExecutableModel.cpp b/source/llvm/LLVMExecutableModel.cpp
@@ -2278,46 +2278,6 @@ double LLVMExecutableModel::getStoichiometry(int speciesIndex, int reactionIndex
     return isnan(result) ? 0 : result;
 }
 
-int LLVMExecutableModel::getStoichiometryMatrix(int* pRows, int* pCols,
-        double** pData)
-{
-    // asking for matrix size
-    if (pRows && pCols && pData == 0)
-    {
-        *pRows = modelData->stoichiometry->m;
-        *pCols = modelData->stoichiometry->n;
-        return modelData->stoichiometry->m * modelData->stoichiometry->n;
-    }
-
-    // allocate data
-    if (pRows && pCols && pData && *pData == 0)
-    {
-        double* data = (double*)malloc(modelData->stoichiometry->m *
-                modelData->stoichiometry->n * sizeof(double));
-        *pRows = modelData->stoichiometry->m;
-        *pCols = modelData->stoichiometry->n;
-        *pData = data;
-        csr_matrix_fill_dense(modelData->stoichiometry, data);
-
-        return modelData->stoichiometry->m * modelData->stoichiometry->n;
-    }
-
-    // use user data
-    if (pRows && *pRows == modelData->stoichiometry->m &&
-            pCols && *pCols == modelData->stoichiometry->n && pData && *pData)
-    {
-        double* data = *pData;
-        csr_matrix_fill_dense(modelData->stoichiometry, data);
-
-        return modelData->stoichiometry->m * modelData->stoichiometry->n;
-    }
-
-
-    throw_llvm_exception("invalid args");
-}
-
-
-
 /******************************* Events Section *******************************/
 #if (1) /**********************************************************************/
 /******************************************************************************/

diff --git a/source/llvm/LLVMExecutableModel.h b/source/llvm/LLVMExecutableModel.h
@@ -348,19 +348,6 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel
 
     virtual double getStoichiometry(int speciesIndex, int reactionIndex);
 
-    /**
-     * allocate a block of memory and copy the stochiometric values into it,
-     * and return it.
-     *
-     * The caller is responsible for freeing the memory that is referenced by data.
-     *
-     * @param[out] rows will hold the number of rows in the matrix.
-     * @param[out] cols will hold the number of columns in the matrix.
-     * @param[out] data a pointer which will hold a newly allocated memory block.
-     */
-    virtual int getStoichiometryMatrix(int* rows, int* cols, double** data);
-
-
     /******************************* Initial Conditions Section *******************/
     #if (1) /**********************************************************************/
     /******************************************************************************/

diff --git a/source/rrExecutableModel.h b/source/rrExecutableModel.h
@@ -597,18 +597,6 @@ namespace rr {
         /***********************************************************************/
         /******************************************************************************/
 
-        /**
-         * allocate a block of memory and copy the stochiometric values into it,
-         * and return it.
-         *
-         * The caller is responsible for freeing the memory that is referenced by data.
-         *
-         * @param[out] rows will hold the number of rows in the matrix.
-         * @param[out] cols will hold the number of columns in the matrix.
-         * @param[out] data a pointer which will hold a newly allocated memory block.
-         */
-        virtual int getStoichiometryMatrix(int *rows, int *cols, double **data) = 0;
-
         /**
          * Get the current stiochiometry value for the given species / reaction.
          *

diff --git a/test/cxx_api_tests/GillespieTests.cpp b/test/cxx_api_tests/GillespieTests.cpp
@@ -9,6 +9,7 @@
 
 #include "GillespieIntegrator.h"
 #include "rrConfig.h"
+#include "RoadRunnerTest.h"
 using namespace rr;
 
 /**
@@ -18,13 +19,33 @@ using namespace rr;
  * that need a new test to fix a bug.
  */
 
-class GillespieTests : public ::testing::Test{
+class GillespieTests : public RoadRunnerTest{
 
 public:
     OpenLinearFlux openLinearFlux;
     GillespieTests() = default;
 };
 
+TEST_F(GillespieTests, OddStoichiometries) {
+    RoadRunner rr1((rrTestModelsDir_ / "SBMLFeatures" / "S1_S2_S2.xml").string());
+    RoadRunner rr2((rrTestModelsDir_ / "SBMLFeatures" / "S1_2S2.xml").string());
+    rr1.setIntegrator("gillespie");
+    rr2.setIntegrator("gillespie");
+    rr1.getIntegrator()->setValue("seed", 1234);
+    rr2.getIntegrator()->setValue("seed", 1234);
+    rr1.getFullStoichiometryMatrix();
+    const ls::DoubleMatrix* results1 = rr1.simulate(0, 5, 20);
+    const ls::DoubleMatrix* results2 = rr2.simulate(0, 5, 20);
+
+    ASSERT_EQ(results1->RSize(), results2->RSize());
+    ASSERT_EQ(results1->CSize(), results2->CSize());
+    for (unsigned int row = 0; row < results1->RSize(); row++) {
+        for (unsigned int col = 0; col < results1->CSize(); col++) {
+            EXPECT_NEAR(results1->Element(row, col), results2->Element(row, col), 0.00001);
+        }
+    }
+}
+
 TEST_F(GillespieTests, Seed){
     RoadRunner rr(openLinearFlux.str());
     rr.setIntegrator("gillespie");

diff --git a/test/mockups/MockExecutableModel.h b/test/mockups/MockExecutableModel.h
@@ -67,7 +67,6 @@ class MockExecutableModel : public rr::ExecutableModel{
     MOCK_METHOD(int, getSupportedIdTypes, (), (override));
     MOCK_METHOD(double, getValue, (const std::string &id), (override));
     MOCK_METHOD(void, setValue, (const std::string &id, double value), (override));
-    MOCK_METHOD(int, getStoichiometryMatrix, (int * rows, int * cols, double * *data), (override));
     MOCK_METHOD(double, getStoichiometry, (int speciesIndex, int reactionIndex), (override));
     MOCK_METHOD(int, getNumConservedMoieties, (), (override));
     MOCK_METHOD(int, getConservedMoietyIndex, (const std::string &eid), (override));

diff --git a/test/models/SBMLFeatures/S1_2S2.xml b/test/models/SBMLFeatures/S1_2S2.xml
@@ -0,0 +1,35 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!-- Created by libAntimony version v2.14.0 with libSBML version 5.20.2. -->
+<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
+  <model metaid="sep_S2" id="sep_S2">
+    <listOfCompartments>
+      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="S1" compartment="default_compartment" initialConcentration="10" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="S2" compartment="default_compartment" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="k1" value="0.3" constant="true"/>
+    </listOfParameters>
+    <listOfReactions>
+      <reaction id="_J0" reversible="true">
+        <listOfReactants>
+          <speciesReference species="S1" stoichiometry="1" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="S2" stoichiometry="2" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> k1 </ci>
+              <ci> S1 </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>
diff --git a/test/models/SBMLFeatures/S1_S2_S2.xml b/test/models/SBMLFeatures/S1_S2_S2.xml
@@ -0,0 +1,36 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!-- Created by libAntimony version v2.14.0 with libSBML version 5.20.2. -->
+<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
+  <model metaid="sep_S2" id="sep_S2">
+    <listOfCompartments>
+      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="S1" compartment="default_compartment" initialConcentration="10" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="S2" compartment="default_compartment" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="k1" value="0.3" constant="true"/>
+    </listOfParameters>
+    <listOfReactions>
+      <reaction id="_J0" reversible="true">
+        <listOfReactants>
+          <speciesReference species="S1" stoichiometry="1" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="S2" stoichiometry="1" constant="true"/>
+          <speciesReference species="S2" stoichiometry="1" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> k1 </ci>
+              <ci> S1 </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>
diff --git a/wrappers/Python/roadrunner/roadrunner.i b/wrappers/Python/roadrunner/roadrunner.i
@@ -2645,25 +2645,6 @@ namespace std { class ostream{}; }
     }
 
 
-    PyObject* getCurrentStoichiometryMatrix() {
-        int rows = 0;
-        int cols = 0;
-        double* data = 0;
-
-        $self->getStoichiometryMatrix(&rows, &cols, &data);
-
-        int nd = 2;
-        npy_intp dims[2] = {rows, cols};
-
-        PyObject *pArray = PyArray_New(&PyArray_Type, nd, dims, NPY_DOUBLE, NULL, data, 0,
-                NPY_ARRAY_CARRAY | NPY_ARRAY_OWNDATA, NULL);
-
-//        VERIFY_PYARRAY(pArray);
-
-        return pArray;
-    }
-
-
     /**
      * get values.
      */