Fixes for several bugs found by @Daniel-James-Cairns-Biology

source/CVODEIntegrator.cpp

@@ -318,22 +318,22 @@ namespace rr {
     void CVODEIntegrator::setIndividualTolerance(std::string sid, double value) {
 
         // the tolerance std::vector that will be stored
-        // [0, numIndFloatingSpecies) stores tolerances for independent floating species
-        // [numIndFloatingSpecies, numIndFloatingSpecies+numRateRule) stores tolerances for variables that have rate rule
+        // Match ExecutableModel::getStateVector(): rate rules first, followed by
+        // independent floating species.
         std::vector<double> v = getAbsoluteToleranceVector();
 
         int speciesIndex = mModel->getFloatingSpeciesIndex(sid);
         std::ptrdiff_t index;
         if (speciesIndex > -1 && speciesIndex < mModel->getNumIndFloatingSpecies()) {
             // sid is an independent floating species
-            v[speciesIndex] = value;
+            v[mModel->getNumRateRules() + speciesIndex] = value;
         } else {
             // sid might has a rate rule
             std::vector<std::string> symbols = mModel->getRateRuleSymbols();
             std::vector<std::string>::iterator it = std::find(symbols.begin(), symbols.end(), sid);
             if (it != symbols.end()) {
                 // found it
-                index = mModel->getNumIndFloatingSpecies() + std::distance(symbols.begin(), it);
+                index = std::distance(symbols.begin(), it);
                 v[index] = value;
             } else {
                 throw std::invalid_argument("CVODEIntegrator::setIndividualTolerance failed, given sid " + sid +
@@ -856,11 +856,11 @@ namespace rr {
             if (initSpecies[s] == 0) {
                 int comp = mModel->getCompartmentIndexForFloatingSpecies(s);
                 if (volumes[comp] != 0.0) {
-                    amount_tolerances[s] = amount_tolerances[s] * abs(volumes[comp]);
+                    amount_tolerances[rrs + s] = amount_tolerances[rrs + s] * abs(volumes[comp]);
                 }
             }
             else {
-                amount_tolerances[s] = amount_tolerances[s] * abs(initSpecies[s]);
+                amount_tolerances[rrs + s] = amount_tolerances[rrs + s] * abs(initSpecies[s]);
             }
         }
 
@@ -873,13 +873,13 @@ namespace rr {
                     // the symbol defined by the rate rule is a species
                     int comp = mModel->getCompartmentIndexForFloatingSpecies(speciesIndex);
                     if (volumes[comp] != 0.0) {
-                        amount_tolerances[species + rr] = amount_tolerances[species + rr] * abs(volumes[comp]);
+                        amount_tolerances[rr] = amount_tolerances[rr] * abs(volumes[comp]);
                     }
                 }
                 //Otherwise just leave amount_tolerances as it is.
             }
             else {
-                amount_tolerances[species + rr] = amount_tolerances[species + rr] * abs(initRRs[rr]);
+                amount_tolerances[rr] = amount_tolerances[rr] * abs(initRRs[rr]);
             }
         }
 

source/EulerIntegrator.h

@@ -61,43 +61,48 @@ namespace rr {
         */
         EulerIntegrator(ExecutableModel *m)
                 : Integrator(m),
-                  eventStatus(std::vector<unsigned char>(m->getNumEvents(), false)),
-                  previousEventStatus(std::vector<unsigned char>(m->getNumEvents(), false)) {
+                  rateBuffer(nullptr),
+                  stateBufferBegin(nullptr),
+                  stateBufferEnd(nullptr),
+                  stateVectorSize(0) {
             EulerIntegrator::resetSettings();
 
-            mModel = m;
             exampleParameter1 = 3.14;
             exampleParameter2 = "hello";
             rrLog(Logger::LOG_WARNING) << "Euler integrator is inaccurate";
-            //std::cerr << "Number of event triggers: " << m->getEventTriggers(0, 0, 0) << std::endl;
-
-            if (mModel) {
-                // calling the getStateVector with a NULL argument returns
-                // the size of teh state std::vector.
-                stateVectorSize = mModel->getStateVector(NULL);
-                rateBuffer = new double[stateVectorSize];
-                stateBufferBegin = new double[stateVectorSize];
-                stateBufferEnd = new double[stateVectorSize];
-            } else {
-                rateBuffer = NULL;
-                stateBufferBegin = NULL;
-                stateBufferEnd = NULL;
-            }
+            syncWithModel(m);
         }
 
         /**
         * delete any memory we allocated
         */
         ~EulerIntegrator() override {
+            delete[] rateBuffer;
+            delete[] stateBufferBegin;
+            delete[] stateBufferEnd;
+        };
+
+        void syncWithModel(ExecutableModel *m) override {
+            delete[] rateBuffer;
+            delete[] stateBufferBegin;
+            delete[] stateBufferEnd;
+            rateBuffer = nullptr;
+            stateBufferBegin = nullptr;
+            stateBufferEnd = nullptr;
+            stateVectorSize = 0;
+
+            mModel = m;
+            eventStatus.clear();
+            previousEventStatus.clear();
             if (mModel) {
-                delete[] rateBuffer;
-                delete[] stateBufferBegin;
-                delete[] stateBufferEnd;
-                rateBuffer = nullptr;
-                stateBufferBegin = nullptr;
-                stateBufferEnd = nullptr;
+                stateVectorSize = mModel->getStateVector(nullptr);
+                rateBuffer = new double[stateVectorSize];
+                stateBufferBegin = new double[stateVectorSize];
+                stateBufferEnd = new double[stateVectorSize];
+                eventStatus.assign(mModel->getNumEvents(), false);
+                previousEventStatus.assign(mModel->getNumEvents(), false);
             }
-        };
+        }
 
         /**
         * integrates the model from t0 to t0 + hstep
@@ -362,12 +367,14 @@ namespace rr {
         * two buffers to store the state std::vector rate, and
         * new state std::vector
         */
-        double *rateBuffer, *stateBufferBegin, *stateBufferEnd;
+        double *rateBuffer = nullptr;
+        double *stateBufferBegin = nullptr;
+        double *stateBufferEnd = nullptr;
 
         /**
         * size of state std::vector
         */
-        int stateVectorSize;
+        int stateVectorSize = 0;
 
         std::vector<unsigned char> eventStatus;
         std::vector<unsigned char> previousEventStatus;

source/Matrix.h

@@ -7,6 +7,7 @@
 
 #include "rr-libstruct/lsMatrix.h"
 #include <algorithm>
+#include <cmath>
 
 namespace rr {
 
@@ -209,7 +210,7 @@ namespace rr {
         bool equals = true;
         for (int i = 0; i < numRows(); i++) {
             for (int j = 0; j < numCols(); j++) {
-                if ((this->operator()(i, j) - other(i, j)) > tolerance) {
+                if (!(std::abs(this->operator()(i, j) - other(i, j)) <= tolerance)) {
                     equals = false;
                     break;
                 }

source/Matrix3D.h

@@ -318,7 +318,7 @@ namespace rr {
             }
             bool equal = true;
             for (int i = 0; i < numZ(); i++) {
-                if ((index_[i] - other.index_[i]) > tol) {
+                if (!(std::abs(index_[i] - other.index_[i]) <= tol)) {
                     equal = false;
                     break;
                 }

source/c/rrCompiledModelGenerator.cpp

@@ -242,7 +242,9 @@ bool CompiledModelGenerator::expressionContainsSymbol(const std::string& express
           return false;
       }
       ASTNode *ast = SBML_parseFormula(expression.c_str());
-      return expressionContainsSymbol(ast, symbol);
+      bool contains = expressionContainsSymbol(ast, symbol);
+      delete ast;
+      return contains;
 }
 
 const Symbol* CompiledModelGenerator::getSpecies(const std::string& id)

source/llvm/Jit.h

@@ -102,8 +102,8 @@ namespace rrllvm {
     using rr_minFnTy = FnPtr_d2;
 
     // for a sparse matrix used in llvm world
-    using csr_matrix_set_nz_FnTy = rr::csr_matrix *(*)(int, int, double);
-    using csr_matrix_get_nz_FnTy = rr::csr_matrix *(*)(int, int);
+    using csr_matrix_set_nz_FnTy = bool (*)(rr::csr_matrix *, unsigned, unsigned, double);
+    using csr_matrix_get_nz_FnTy = double (*)(const rr::csr_matrix *, unsigned, unsigned);
 
     // function signatures for distrib
     using DistribFnTy_d1 = double (*)(Random *, double);

source/rrIniFile.cpp

@@ -162,7 +162,9 @@ namespace rr
 				else if (Line.find_first_of('[') == 0) // Found a section
 				{
 					Line.erase(0, 1);
-					Line.erase(Line.find_last_of(']'), 1);
+					std::string::size_type closingBracket = Line.find_last_of(']');
+					if (closingBracket != std::string::npos)
+						Line.erase(closingBracket, 1);
 					pSection = GetSection(Line, true);
 					rrLog(lDebug3) << "Located section: " + pSection->mName;
 					Comment = std::string("");
@@ -296,7 +298,9 @@ namespace rr
 					if (Line.find_first_of('[') == 0) // Found a section
 					{
 						Line.erase(0, 1);
-						Line.erase(Line.find_last_of(']'), 1);
+						std::string::size_type closingBracket = Line.find_last_of(']');
+						if (closingBracket != std::string::npos)
+							Line.erase(closingBracket, 1);
 
 						if (theSection == Line)
 						{

source/rrRoadRunner.cpp

@@ -74,6 +74,23 @@ namespace rr {
 
     typedef std::vector<std::string> string_vector;
 
+    class ConfigValueGuard {
+    public:
+        explicit ConfigValueGuard(Config::Keys key)
+            : key_(key), value_(Config::getValue(key)) {}
+
+        ~ConfigValueGuard() {
+            Config::setValue(key_, value_);
+        }
+
+        ConfigValueGuard(const ConfigValueGuard&) = delete;
+        ConfigValueGuard& operator=(const ConfigValueGuard&) = delete;
+
+    private:
+        Config::Keys key_;
+        Setting value_;
+    };
+
 
 // we can write a single function to pick the std::string lists out
 // of the model instead of duplicating it 6 times with
@@ -3234,7 +3251,7 @@ namespace rr {
         check_model();
 
         get_self();
-        std::int32_t savedJacobianMode = Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs<std::int32_t>();
+        ConfigValueGuard restoreJacobianMode(Config::ROADRUNNER_JACOBIAN_MODE);
         Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, Config::ROADRUNNER_JACOBIAN_MODE_CONCENTRATIONS); 
 
         if (self.model->getNumReactions() == 0 && self.model->getNumRateRules() > 0) {
@@ -3337,9 +3354,6 @@ namespace rr {
                 }
 
             }
-            // Put back User selected JACOBIAN_MODE:
-            Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, savedJacobianMode);
-
             // get the row/column ids, independent floating species
             std::list<std::string> list;
             self.model->getIds(SelectionRecord::FLOATING_AMOUNT, list);
@@ -3362,7 +3376,7 @@ namespace rr {
             h = self.roadRunnerOptions.jacobianStepSize;
         }
 
-        std::int32_t savedJacobianMode = Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs<std::int32_t>();
+        ConfigValueGuard restoreJacobianMode(Config::ROADRUNNER_JACOBIAN_MODE);
 
         // For our purposes here, we want all independent floating species
         // plus all floating species that have rate rules.
@@ -3457,8 +3471,6 @@ namespace rr {
                 jac(ri, ci) = origVal / compVol;
             }
 
-            // Put back User selected JACOBIAN_MODE:
-            Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, savedJacobianMode);
         }
         return jac;
     }

source/rrSparse.cpp

@@ -128,7 +128,7 @@ csr_matrix* csr_matrix_new(unsigned m, unsigned n,
 
 bool csr_matrix_set_nz(csr_matrix* mat, unsigned row, unsigned col, double val)
 {
-    if (mat && row <= mat->m && col <= mat->n)
+    if (mat && row < mat->m && col < mat->n)
     {
         for (unsigned k = mat->rowptr[row]; k < mat->rowptr[row + 1]; k++)
         {
@@ -144,7 +144,7 @@ bool csr_matrix_set_nz(csr_matrix* mat, unsigned row, unsigned col, double val)
 
 double csr_matrix_get_nz(const csr_matrix* mat, unsigned row, unsigned col)
 {
-    if (mat && row <= mat->m && col <= mat->n)
+    if (mat && row < mat->m && col < mat->n)
     {
         for (unsigned k = mat->rowptr[row]; k < mat->rowptr[row + 1]; k++)
         {
@@ -325,4 +325,3 @@ std::ostream& operator <<(std::ostream& os, const csr_matrix* mat)
 
 }
 
-

test/cxx_api_tests/CMakeLists.txt

@@ -2,6 +2,9 @@ add_test_executable(test_cxx_api_SettingTests test_targets SettingTests.cpp)
 add_test_executable(test_cxx_api_RoadRunnerAPITestsMCJit test_targets RoadRunnerAPITests.h RoadRunnerAPITestsWithMCJit.cpp ${SharedTestFiles})
 add_test_executable(test_cxx_api_RoadRunnerAPITestsLLJit test_targets RoadRunnerAPITests.h RoadRunnerAPITestsWithLLJit.cpp ${SharedTestFiles})
 add_test_executable(test_cxx_api_GillespieTests test_targets GillespieTests.cpp ${SharedTestFiles})
+add_test_executable(test_cxx_api_EulerIntegratorTests test_targets EulerIntegratorTests.cpp ${SharedTestFiles})
+add_test_executable(test_cxx_api_SparseMatrixTests test_targets SparseMatrixTests.cpp ${SharedTestFiles})
+add_test_executable(test_cxx_api_IniFileTests test_targets IniFileTests.cpp ${SharedTestFiles})
 add_test_executable(test_cxx_api_LoggerTests test_targets LoggerTests.cpp ${SharedTestFiles})
 add_test_executable(test_cxx_api_BasicDictionaryTests test_targets BasicDictionaryTests.cpp ${SharedTestFiles})
 add_test_executable(test_cxx_api_SelectionRecordTests test_targets SelectionRecordTests.cpp ${SharedTestFiles})