From 70069de2af3590102e6e0205822cd455c962a101 Mon Sep 17 00:00:00 2001 From: Lucian Smith Date: Wed, 19 Aug 2026 16:38:29 -0700 Subject: [PATCH 1/2] Fixes for several bugs found by @Daniel-James-Cairns-Biology Presumably discovered by Claude, since one of the fixes is in completely dead code. Used the submitted patches, and added tests myself. --- source/CVODEIntegrator.cpp | 16 ++-- source/EulerIntegrator.h | 59 ++++++------ source/Matrix.h | 3 +- source/Matrix3D.h | 2 +- source/c/rrCompiledModelGenerator.cpp | 4 +- source/llvm/Jit.h | 4 +- source/rrIniFile.cpp | 8 +- source/rrRoadRunner.cpp | 26 ++++-- source/rrSparse.cpp | 5 +- test/cxx_api_tests/CMakeLists.txt | 3 + test/cxx_api_tests/JacobianTests.cpp | 55 +++++++++++ test/cxx_api_tests/Matrix3DTests.cpp | 45 +++++++++ test/cxx_api_tests/MatrixTests.cpp | 29 ++++++ test/llvm_tests/CMakeLists.txt | 1 + .../CVODEIntegratorTests/CvodeUnitTest.cpp | 92 +++++++++++++++++++ 15 files changed, 301 insertions(+), 51 deletions(-) diff --git a/source/CVODEIntegrator.cpp b/source/CVODEIntegrator.cpp index 2b37835b30..a9e33b115a 100644 --- a/source/CVODEIntegrator.cpp +++ b/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 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 symbols = mModel->getRateRuleSymbols(); std::vector::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]); } } diff --git a/source/EulerIntegrator.h b/source/EulerIntegrator.h index f88e3cb4d9..e8651c43fb 100644 --- a/source/EulerIntegrator.h +++ b/source/EulerIntegrator.h @@ -61,43 +61,48 @@ namespace rr { */ EulerIntegrator(ExecutableModel *m) : Integrator(m), - eventStatus(std::vector(m->getNumEvents(), false)), - previousEventStatus(std::vector(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 eventStatus; std::vector previousEventStatus; diff --git a/source/Matrix.h b/source/Matrix.h index 4a1b4a512c..ac1b4b3da9 100644 --- a/source/Matrix.h +++ b/source/Matrix.h @@ -7,6 +7,7 @@ #include "rr-libstruct/lsMatrix.h" #include +#include 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; } diff --git a/source/Matrix3D.h b/source/Matrix3D.h index d5bc85cb97..2894511efd 100644 --- a/source/Matrix3D.h +++ b/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; } diff --git a/source/c/rrCompiledModelGenerator.cpp b/source/c/rrCompiledModelGenerator.cpp index 8916fab81d..162f8a5824 100644 --- a/source/c/rrCompiledModelGenerator.cpp +++ b/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) diff --git a/source/llvm/Jit.h b/source/llvm/Jit.h index f41e421435..055bf2efe2 100644 --- a/source/llvm/Jit.h +++ b/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); diff --git a/source/rrIniFile.cpp b/source/rrIniFile.cpp index 86b366d8e2..21638c01e8 100644 --- a/source/rrIniFile.cpp +++ b/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) { diff --git a/source/rrRoadRunner.cpp b/source/rrRoadRunner.cpp index 5655dbb2c8..807d42dce4 100644 --- a/source/rrRoadRunner.cpp +++ b/source/rrRoadRunner.cpp @@ -74,6 +74,23 @@ namespace rr { typedef std::vector 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(); + 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 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(); + 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; } diff --git a/source/rrSparse.cpp b/source/rrSparse.cpp index 35d58e8ae9..bb2fabfa70 100644 --- a/source/rrSparse.cpp +++ b/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) } - diff --git a/test/cxx_api_tests/CMakeLists.txt b/test/cxx_api_tests/CMakeLists.txt index c1bc1d17cc..bdbc085f70 100644 --- a/test/cxx_api_tests/CMakeLists.txt +++ b/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}) diff --git a/test/cxx_api_tests/JacobianTests.cpp b/test/cxx_api_tests/JacobianTests.cpp index fc18870bad..9075880d58 100644 --- a/test/cxx_api_tests/JacobianTests.cpp +++ b/test/cxx_api_tests/JacobianTests.cpp @@ -10,6 +10,30 @@ using namespace rr; +// A model with no reactions, governed entirely by a rate rule (dS0/dt = -S0). +// getFullJacobian() takes an early-return path for this case. +static const std::string PureRateRuleNoReactionsSBML = R"( + + + + + + + + S0 + + +)"; + +// A model with no species, no reactions, and no rate rules. Both +// getFullJacobian() (the nr == 0 branch) and getReducedJacobian() (zero +// included species) take early/degenerate paths for this case. +static const std::string NoSpeciesNoReactionsSBML = R"( + + + +)"; + class JacobianTests : public RoadRunnerTest { public: @@ -190,4 +214,35 @@ TEST_F(JacobianTests, SimpleFluxReducedConc) { checkJacobianReducedConc("SimpleFlux", 1e-4); } +/** + * Regression tests for the ROADRUNNER_JACOBIAN_MODE leak: getFullJacobian() + * and getReducedJacobian() force CONCENTRATIONS mode and are supposed to + * restore the caller's setting before returning, on every path. Each of + * these hits an early/degenerate return that historically skipped the + * restore, leaving the config stuck at CONCENTRATIONS. + */ +TEST_F(JacobianTests, FullJacobianRestoresConfigOnPureRateRuleEarlyReturn) { + RoadRunner rr(PureRateRuleNoReactionsSBML); + Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS); + rr.getFullJacobian(); + ASSERT_EQ((std::int32_t) Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS, + Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs()); +} + +TEST_F(JacobianTests, FullJacobianRestoresConfigWhenNoReactionsOrRateRules) { + RoadRunner rr(NoSpeciesNoReactionsSBML); + Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS); + rr.getFullJacobian(); + ASSERT_EQ((std::int32_t) Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS, + Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs()); +} + +TEST_F(JacobianTests, ReducedJacobianRestoresConfigWhenNoIncludedSpecies) { + RoadRunner rr(NoSpeciesNoReactionsSBML); + Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS); + rr.getReducedJacobian(); + ASSERT_EQ((std::int32_t) Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS, + Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs()); +} + diff --git a/test/cxx_api_tests/Matrix3DTests.cpp b/test/cxx_api_tests/Matrix3DTests.cpp index ac9df2afc9..2cb54950c7 100644 --- a/test/cxx_api_tests/Matrix3DTests.cpp +++ b/test/cxx_api_tests/Matrix3DTests.cpp @@ -4,6 +4,7 @@ #include "gtest/gtest.h" #include "Matrix3D.h" +#include using namespace rr; @@ -417,6 +418,50 @@ TEST_F(Matrix3DTests, AlmostEqualsWhenFalse) { ASSERT_FALSE(first.almostEquals(second, 1e-4)); } +TEST_F(Matrix3DTests, AlmostEqualsWhenFalseReverseDirection) { + // almostEquals must be symmetric. The one-sided check + // (index_[i] - other.index_[i]) > tol only caught the case where the + // callee's index was smaller than the argument's. + Matrix3D small( + {0.0, 1.0}, + { + {{0.0}}, + {{0.0}} + } + ); + Matrix3D large( + {0.0, 1000.0}, + { + {{0.0}}, + {{0.0}} + } + ); + ASSERT_FALSE(small.almostEquals(large, 1e-4)); + ASSERT_FALSE(large.almostEquals(small, 1e-4)); +} + +TEST_F(Matrix3DTests, AlmostEqualsWithNaNIndexDifferenceIsFalse) { + // NaN comparisons are always false, so the old check treated a NaN + // difference between index values as "not greater than tolerance", + // silently reporting equality. + double nan = std::numeric_limits::quiet_NaN(); + Matrix3D first( + {nan, 1.0}, + { + {{0.0}}, + {{0.0}} + } + ); + Matrix3D second( + {2.0, 1.0}, + { + {{0.0}}, + {{0.0}} + } + ); + ASSERT_FALSE(first.almostEquals(second, 1e-4)); +} + TEST_F(Matrix3DTests, CheckRowNames) { Matrix3D matrix3D( {0.0, 1.0}, diff --git a/test/cxx_api_tests/MatrixTests.cpp b/test/cxx_api_tests/MatrixTests.cpp index 170abe58f2..3e2339c5d3 100644 --- a/test/cxx_api_tests/MatrixTests.cpp +++ b/test/cxx_api_tests/MatrixTests.cpp @@ -4,6 +4,7 @@ #include "gtest/gtest.h" #include "Matrix.h" +#include using namespace rr; @@ -62,6 +63,34 @@ TEST_F(MatrixTests, TestAlmostEqualsWhenFalse) { ASSERT_FALSE(first.almostEquals(second, 1e-5)); } +TEST_F(MatrixTests, TestAlmostEqualsWhenFalseReverseDirection) { + // almostEquals must be symmetric: if a.almostEquals(b) is false because + // b is bigger, b.almostEquals(a) must also be false. The one-sided check + // (lhs - rhs) > tolerance only caught the first ordering. + Matrix small({ + {0.0, 0.0}, + }); + Matrix large({ + {1000.0, 0.0}, + }); + ASSERT_FALSE(small.almostEquals(large, 1e-4)); + ASSERT_FALSE(large.almostEquals(small, 1e-4)); +} + +TEST_F(MatrixTests, TestAlmostEqualsWithNaNDifferenceIsFalse) { + // NaN comparisons are always false, so the old check + // (lhs - rhs) > tolerance treated any NaN difference as "not greater + // than tolerance", silently reporting equality. + double nan = std::numeric_limits::quiet_NaN(); + Matrix first({ + {nan, 1.0}, + }); + Matrix second({ + {2.0, 1.0}, + }); + ASSERT_FALSE(first.almostEquals(second, 1e-4)); +} + TEST_F(MatrixTests, TestAlmostEqualsWhenFalseButAcceptableTolerance) { Matrix first({ {1.12346, 2.1234, 3.1234}, diff --git a/test/llvm_tests/CMakeLists.txt b/test/llvm_tests/CMakeLists.txt index 5a9b1be499..cc12ed1a48 100644 --- a/test/llvm_tests/CMakeLists.txt +++ b/test/llvm_tests/CMakeLists.txt @@ -3,6 +3,7 @@ add_test_executable(test_llvm_MCJitTests test_targets MCJitTests.cpp JitTests.h) add_test_executable(test_llvm_LLJitTests test_targets LLJitTests.cpp JitTests.h) add_test_executable(test_llvm_MCJitMapFunctionsToJitSymbolsTests test_targets FunctionMapping/MCJitMapFunctionsToJitSymbolsTests.cpp) add_test_executable(test_llvm_LLJitMapFunctionsToJitSymbolsTests test_targets FunctionMapping/LLJitMapFunctionsToJitSymbolsTests.cpp) +add_test_executable(test_llvm_CsrMatrixJitAliasSignatureTests test_targets FunctionMapping/CsrMatrixJitAliasSignatureTests.cpp) add_test_executable(test_llvm_RandomTests test_targets RandomTests.cpp) add_test_executable( test_llvm_EvalInitialConditionsTests test_targets diff --git a/test/sundials-tests/CVODEIntegratorTests/CvodeUnitTest.cpp b/test/sundials-tests/CVODEIntegratorTests/CvodeUnitTest.cpp index 6d30de027f..853553119d 100644 --- a/test/sundials-tests/CVODEIntegratorTests/CvodeUnitTest.cpp +++ b/test/sundials-tests/CVODEIntegratorTests/CvodeUnitTest.cpp @@ -226,6 +226,98 @@ TEST_F(CVODEIntegratorUnitTests, restart) { cvodeIntegrator.restart(0); } +/** + * Regression tests for the tolerance-vector packing order. + * + * ExecutableModel::getStateVector() packs the state as + * [rate-rule variables, independent floating species]. The tolerance + * vector must follow the same layout, but setIndividualTolerance() and + * getAbsoluteToleranceVector() historically assumed the reverse order + * ([species, rate rules]), silently mapping tolerances to the wrong + * state variable whenever a model had both. + */ +TEST_F(CVODEIntegratorUnitTests, SetIndividualToleranceRateRuleComesBeforeSpecies) { + // 2 independent species (S0, S1) and 1 rate-rule variable (R0): state + // vector length 3, overriding the fixture's default of 2. + EXPECT_CALL(mockExecutableModel, getStateVector).WillRepeatedly(Return(3)); + EXPECT_CALL(mockExecutableModel, getNumIndFloatingSpecies).WillRepeatedly(Return(2)); + EXPECT_CALL(mockExecutableModel, getNumRateRules).WillRepeatedly(Return(1)); + EXPECT_CALL(mockExecutableModel, getRateRuleSymbols) + .WillRepeatedly(Return(std::vector{"R0"})); + EXPECT_CALL(mockExecutableModel, getFloatingSpeciesIndex("R0")).WillRepeatedly(Return(-1)); + EXPECT_CALL(mockExecutableModel, getNumCompartments).WillRepeatedly(Return(1)); + + CVODEIntegrator cvodeIntegrator(&mockExecutableModel); + cvodeIntegrator.setValue("absolute_tolerance", 1e-6); + cvodeIntegrator.setIndividualTolerance("R0", 0.5); + + std::vector tol = cvodeIntegrator.getValue("absolute_tolerance").get>(); + ASSERT_EQ(3, tol.size()); + // Rate rules occupy [0, numRateRules); R0 is the only one, so its + // tolerance must land at index 0, not at index 2 (numIndFloatingSpecies). + ASSERT_NEAR(0.5, tol[0], 1e-9); + ASSERT_NEAR(1e-6, tol[1], 1e-9); + ASSERT_NEAR(1e-6, tol[2], 1e-9); +} + +TEST_F(CVODEIntegratorUnitTests, SetIndividualToleranceForSpeciesIsOffsetByRateRuleCount) { + // 2 independent species (S0, S1) and 1 rate-rule variable: state + // vector length 3, overriding the fixture's default of 2. + EXPECT_CALL(mockExecutableModel, getStateVector).WillRepeatedly(Return(3)); + EXPECT_CALL(mockExecutableModel, getNumIndFloatingSpecies).WillRepeatedly(Return(2)); + EXPECT_CALL(mockExecutableModel, getNumRateRules).WillRepeatedly(Return(1)); + // getAbsoluteToleranceVector() scales rate-rule tolerances unconditionally, + // even though this test's tolerance update itself never touches a rate + // rule. With getRateRuleValues() left at its default of 0, that scaling + // takes the branch that also calls getFloatingSpeciesIndex("R0") to see + // whether the rate-rule variable is itself a species -- so that needs a + // matching expectation too, or gmock treats it as an unexpected call to + // a method it only has an "S1" expectation for. + EXPECT_CALL(mockExecutableModel, getRateRuleSymbols) + .WillRepeatedly(Return(std::vector{"R0"})); + EXPECT_CALL(mockExecutableModel, getFloatingSpeciesIndex("R0")).WillRepeatedly(Return(-1)); + EXPECT_CALL(mockExecutableModel, getFloatingSpeciesIndex("S1")).WillRepeatedly(Return(1)); + EXPECT_CALL(mockExecutableModel, getNumCompartments).WillRepeatedly(Return(1)); + + CVODEIntegrator cvodeIntegrator(&mockExecutableModel); + cvodeIntegrator.setValue("absolute_tolerance", 1e-6); + cvodeIntegrator.setIndividualTolerance("S1", 0.25); + + std::vector tol = cvodeIntegrator.getValue("absolute_tolerance").get>(); + ASSERT_EQ(3, tol.size()); + // S1 is independent-species index 1. With 1 rate rule packed ahead of + // the species block, its slot is numRateRules + 1 == 2, not 1. + ASSERT_NEAR(0.25, tol[2], 1e-9); + ASSERT_NEAR(1e-6, tol[0], 1e-9); + ASSERT_NEAR(1e-6, tol[1], 1e-9); +} + +TEST_F(CVODEIntegratorUnitTests, AbsoluteToleranceVectorPacksRateRulesBeforeSpecies) { + // 1 independent species (S0, amount 5) and 1 rate-rule variable + // (R0, value 3). Both amounts are non-zero, so getAbsoluteToleranceVector + // scales each slot by abs(value) directly, without needing compartment + // volumes. + EXPECT_CALL(mockExecutableModel, getNumIndFloatingSpecies).WillRepeatedly(Return(1)); + EXPECT_CALL(mockExecutableModel, getNumRateRules).WillRepeatedly(Return(1)); + EXPECT_CALL(mockExecutableModel, getNumCompartments).WillRepeatedly(Return(1)); + EXPECT_CALL(mockExecutableModel, getRateRuleSymbols) + .WillRepeatedly(Return(std::vector{"R0"})); + EXPECT_CALL(mockExecutableModel, getFloatingSpeciesAmounts) + .WillRepeatedly(DoAll(SetArgPointee<2>(5.0), Return(0))); + EXPECT_CALL(mockExecutableModel, getRateRuleValues) + .WillRepeatedly(SetArgPointee<0>(3.0)); + + CVODEIntegrator cvodeIntegrator(&mockExecutableModel); + cvodeIntegrator.setValue("absolute_tolerance", 1e-6); + std::vector tol = cvodeIntegrator.getAbsoluteToleranceVector(); + + ASSERT_EQ(2, tol.size()); + // Packed order is [rate rules, species]: R0 (value 3) at index 0, + // S0 (amount 5) at index 1. + ASSERT_NEAR(1e-6 * 3.0, tol[0], 1e-12); + ASSERT_NEAR(1e-6 * 5.0, tol[1], 1e-12); +} + From 76aa2d2e1f5ad5a830e0fc26b3af7e5ee1b1e6ce Mon Sep 17 00:00:00 2001 From: Lucian Smith Date: Wed, 19 Aug 2026 17:17:08 -0700 Subject: [PATCH 2/2] Actually add new tests. --- test/cxx_api_tests/EulerIntegratorTests.cpp | 108 ++++++++++++++++++ test/cxx_api_tests/IniFileTests.cpp | 60 ++++++++++ test/cxx_api_tests/SparseMatrixTests.cpp | 54 +++++++++ .../CsrMatrixJitAliasSignatureTests.cpp | 61 ++++++++++ 4 files changed, 283 insertions(+) create mode 100644 test/cxx_api_tests/EulerIntegratorTests.cpp create mode 100644 test/cxx_api_tests/IniFileTests.cpp create mode 100644 test/cxx_api_tests/SparseMatrixTests.cpp create mode 100644 test/llvm_tests/FunctionMapping/CsrMatrixJitAliasSignatureTests.cpp diff --git a/test/cxx_api_tests/EulerIntegratorTests.cpp b/test/cxx_api_tests/EulerIntegratorTests.cpp new file mode 100644 index 0000000000..eb7eb770e5 --- /dev/null +++ b/test/cxx_api_tests/EulerIntegratorTests.cpp @@ -0,0 +1,108 @@ +#include "gtest/gtest.h" +#include "rrRoadRunner.h" +#include "RoadRunnerTest.h" +#include "EulerIntegrator.h" + +using namespace rr; + +// One species, no events. Small enough that the Euler integrator's +// internal buffers end up sized for exactly one state variable and zero +// events. +static const std::string OneSpeciesNoEventSBML = R"( + + + + + + + + + + kS0 + + +)"; + +// Five species and one event. Reloading into this model after the one +// above should force the Euler integrator to rebuild its state/rate +// buffers (now sized for 5 states) and its event-status vector (now sized +// for 1 event), and drop its pointer to the model RoadRunner::load() just +// replaced. +static const std::string FiveSpeciesWithEventSBML = R"( + + + + + + + + + + + + + + + kA + + + + + + A0.9 + + + + 99 + + + + +)"; + +class EulerIntegratorTests : public RoadRunnerTest { +public: + EulerIntegratorTests() = default; +}; + +/** + * Regression test for EulerIntegrator inheriting Integrator's empty + * syncWithModel(). RoadRunner::load() deletes the old ExecutableModel and + * calls syncWithModel() on every already-constructed integrator so they can + * pick up the new one and resize their internal state. Euler didn't + * override it, so after a reload it kept a dangling pointer to the freed + * model plus rate/state buffers and an event-status vector sized for the + * model that no longer exists. + * + * Under a plain (non-sanitized) build this may silently produce wrong + * results rather than crash outright -- the reliable signal is a + * heap-use-after-free / heap-buffer-overflow under ASan or valgrind. The + * functional assertion below (the event never fires) is a bug that doesn't + * depend on a sanitizer to observe, but running this test under ASan is + * the more decisive check for the underlying memory-safety fix. + */ +TEST_F(EulerIntegratorTests, ReloadResizesBuffersAndEventsWithoutStaleModel) { + RoadRunner rr(OneSpeciesNoEventSBML); + rr.setIntegrator("euler"); + + // Run once so Euler actually allocates its buffers/event vector + // against the 1-species/0-event model. + rr.simulate(0, 1, 10); + + rr.load(FiveSpeciesWithEventSBML); + ASSERT_EQ("euler", rr.getIntegrator()->getName()) + << "reload should not have reset the selected integrator"; + + ASSERT_NO_THROW({ + rr.simulate(0, 20, 200); + }); + + // A decays below 0.9 well within this window, so the event should have + // fired and set E = 99. If Euler's event-status vector never grew past + // its stale size of zero (carried over from the event-free model), the + // trigger can never be observed and E stays at its initial value of 1. + double finalE = rr.getValue("E"); + EXPECT_NEAR(99.0, finalE, 1e-6) + << "event did not fire after reload; Euler's event bookkeeping " + "likely still reflects the pre-reload model."; +} diff --git a/test/cxx_api_tests/IniFileTests.cpp b/test/cxx_api_tests/IniFileTests.cpp new file mode 100644 index 0000000000..7ae4e269ab --- /dev/null +++ b/test/cxx_api_tests/IniFileTests.cpp @@ -0,0 +1,60 @@ +#include "gtest/gtest.h" +#include "rrIniFile.h" +#include +#include + +using namespace rr; + +class IniFileTests : public ::testing::Test { +public: + std::filesystem::path tempFile; + + IniFileTests() { + tempFile = std::filesystem::temp_directory_path() / "roadrunner_malformed_ini_test.ini"; + } + + ~IniFileTests() override { + std::error_code ec; + std::filesystem::remove(tempFile, ec); + } + + void writeFile(const std::string &contents) { + std::ofstream out(tempFile, std::ios::out | std::ios::trunc); + out << contents; + out.close(); + } +}; + +/** + * Regression tests: Load() and LoadSection() both unconditionally erase() + * at find_last_of(']') for any line starting with '['. For a section + * header with no closing bracket, find_last_of returns npos, and + * erase(npos, 1) throws std::out_of_range -- turning a tolerantly-parsed + * malformed file into a crash. + */ +TEST_F(IniFileTests, LoadDoesNotThrowOnSectionMissingClosingBracket) { + writeFile( + "[missing-close\n" + "key1=value1\n" + "[goodsection]\n" + "key2=value2\n" + ); + + IniFile ini; + ASSERT_NO_THROW(ini.Load(tempFile.string())); + + // Parsing should have recovered and continued past the malformed line. + EXPECT_TRUE(ini.SectionExists("goodsection")); +} + +TEST_F(IniFileTests, LoadSectionDoesNotThrowOnSectionMissingClosingBracket) { + writeFile( + "[missing-close\n" + "key1=value1\n" + "[goodsection]\n" + "key2=value2\n" + ); + + IniFile ini(tempFile.string()); + ASSERT_NO_THROW(ini.LoadSection("goodsection")); +} diff --git a/test/cxx_api_tests/SparseMatrixTests.cpp b/test/cxx_api_tests/SparseMatrixTests.cpp new file mode 100644 index 0000000000..9e58547c22 --- /dev/null +++ b/test/cxx_api_tests/SparseMatrixTests.cpp @@ -0,0 +1,54 @@ +#include "gtest/gtest.h" +#include "rrSparse.h" +#include + +using namespace rr; + +class SparseMatrixTests : public ::testing::Test { +public: + SparseMatrixTests() = default; +}; + +TEST_F(SparseMatrixTests, GetNzValidIndexReturnsStoredValue) { + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {42.0}); + ASSERT_NEAR(42.0, csr_matrix_get_nz(mat, 0, 0), 1e-9); + csr_matrix_delete(mat); +} + +/** + * Regression tests for a strict-bounds bug: csr_matrix_get_nz/set_nz + * accepted row == m or col == n -- one past the last valid index -- because + * they checked row <= m / col <= n instead of strict less-than. rowptr is + * only allocated with m + 1 entries, so row == m reads rowptr[m + 1], one + * past the allocation: an ASan-detectable heap-buffer-overflow on a 1x1 + * matrix. Under a non-sanitized build this may not crash and could return + * whatever garbage follows the allocation instead of NaN/false, so these + * are most decisive when run under ASan or valgrind. + */ +TEST_F(SparseMatrixTests, GetNzRowEqualToDimensionIsRejected) { + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {42.0}); + double result = csr_matrix_get_nz(mat, 1, 0); // row == m + ASSERT_TRUE(std::isnan(result)); + csr_matrix_delete(mat); +} + +TEST_F(SparseMatrixTests, GetNzColEqualToDimensionIsRejected) { + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {42.0}); + double result = csr_matrix_get_nz(mat, 0, 1); // col == n + ASSERT_TRUE(std::isnan(result)); + csr_matrix_delete(mat); +} + +TEST_F(SparseMatrixTests, SetNzRowEqualToDimensionIsRejected) { + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {0.0}); + bool wrote = csr_matrix_set_nz(mat, 1, 0, 7.0); // row == m + ASSERT_FALSE(wrote); + csr_matrix_delete(mat); +} + +TEST_F(SparseMatrixTests, SetNzColEqualToDimensionIsRejected) { + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {0.0}); + bool wrote = csr_matrix_set_nz(mat, 0, 1, 7.0); // col == n + ASSERT_FALSE(wrote); + csr_matrix_delete(mat); +} diff --git a/test/llvm_tests/FunctionMapping/CsrMatrixJitAliasSignatureTests.cpp b/test/llvm_tests/FunctionMapping/CsrMatrixJitAliasSignatureTests.cpp new file mode 100644 index 0000000000..5d24309c2a --- /dev/null +++ b/test/llvm_tests/FunctionMapping/CsrMatrixJitAliasSignatureTests.cpp @@ -0,0 +1,61 @@ +// +// Regression tests for the JIT sparse-matrix function-pointer aliases. +// +// csr_matrix_get_nz_FnTy / csr_matrix_set_nz_FnTy (declared in Jit.h) exist +// so that code can safely call through a JIT-looked-up function address +// for csr_matrix_get_nz/csr_matrix_set_nz using these aliases. Before the +// fix the aliases didn't match the real functions at all -- wrong return +// type, wrong parameter count and types -- so calling through them was +// undefined behavior. The existing tests in MapFunctionsToJitSymbolsTests.h +// only checked that the looked-up address was non-null, which can't catch +// a signature mismatch. +// +// This file intentionally calls through the aliases with the real +// arguments (a csr_matrix*, plus row/col indices). Pre-fix, this does not +// even compile -- the alias only accepts two plain ints and returns a +// csr_matrix* -- so a compile failure here is the expected pre-patch +// signal, not a sign these tests are broken. Post-fix it compiles and the +// calls behave correctly. +// + +#include "gtest/gtest.h" +#include "llvm/LLJit.h" +#include "llvm/MCJit.h" +#include "Jit.h" +#include "rrSparse.h" +#include "rrRoadRunnerOptions.h" + +using namespace rr; +using namespace rrllvm; + +TEST(CsrMatrixJitAliasSignatureTests, LLJitGetAndSetNzInvokeCorrectly) { + LLJit llJit(LoadSBMLOptions().modelGeneratorOpt); + + csr_matrix_get_nz_FnTy getNz = + (csr_matrix_get_nz_FnTy) llJit.lookupFunctionAddress("csr_matrix_get_nz"); + csr_matrix_set_nz_FnTy setNz = + (csr_matrix_set_nz_FnTy) llJit.lookupFunctionAddress("csr_matrix_set_nz"); + ASSERT_FALSE(getNz == nullptr); + ASSERT_FALSE(setNz == nullptr); + + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {0.0}); + ASSERT_TRUE(setNz(mat, 0, 0, 7.5)); + ASSERT_NEAR(7.5, getNz(mat, 0, 0), 1e-9); + csr_matrix_delete(mat); +} + +TEST(CsrMatrixJitAliasSignatureTests, MCJitGetAndSetNzInvokeCorrectly) { + MCJit mcJit(LoadSBMLOptions().modelGeneratorOpt); + + csr_matrix_get_nz_FnTy getNz = + (csr_matrix_get_nz_FnTy) mcJit.lookupFunctionAddress("csr_matrix_get_nz"); + csr_matrix_set_nz_FnTy setNz = + (csr_matrix_set_nz_FnTy) mcJit.lookupFunctionAddress("csr_matrix_set_nz"); + ASSERT_FALSE(getNz == nullptr); + ASSERT_FALSE(setNz == nullptr); + + csr_matrix *mat = csr_matrix_new(1, 1, {0}, {0}, {0.0}); + ASSERT_TRUE(setNz(mat, 0, 0, 7.5)); + ASSERT_NEAR(7.5, getNz(mat, 0, 0), 1e-9); + csr_matrix_delete(mat); +}