diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 09fe6ca00f..fb9d82cf1c 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -957,9 +957,9 @@ stages: echo "sudo yum install -y ncurses-devel" sudo yum install -y ncurses-devel - # && cmake -DLLVM_INSTALL_PREFIX=/llvm-6.x/install-llvm-6.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/libroadrunner-deps/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$roadrunnerInstall37 -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=/Miniconda3/envs/py37 -DSWIG_EXECUTABLE=$swigExec .. - echo "cmake command: cmake -DLLVM_INSTALL_PREFIX=/llvm-6.x/install-llvm-6.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/libroadrunner-deps/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$(INSTALL_DIRECTORY) -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=$(PythonRoot) -DSWIG_EXECUTABLE=$(SwigExecutable) .." - cmake -DLLVM_INSTALL_PREFIX=/llvm-6.x/install-llvm-6.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/libroadrunner-deps/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$(INSTALL_DIRECTORY) -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=$(PythonRoot) -DSWIG_EXECUTABLE=$(SwigExecutable) .. + # && cmake -DLLVM_INSTALL_PREFIX=/llvm-6.x/install-llvm-6.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/libroadrunner-deps/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$roadrunnerInstall37 -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=/Miniconda3/envs/py37 -DSWIG_EXECUTABLE=$swigExec -DCMAKE_BUILD_TYPE=Release .. + echo "cmake command: cmake -DLLVM_INSTALL_PREFIX=/llvm-6.x/install-llvm-6.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/libroadrunner-deps/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$(INSTALL_DIRECTORY) -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=$(PythonRoot) -DSWIG_EXECUTABLE=$(SwigExecutable) -DCMAKE_BUILD_TYPE=Release .." + cmake -DLLVM_INSTALL_PREFIX=/llvm-6.x/install-llvm-6.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/libroadrunner-deps/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$(INSTALL_DIRECTORY) -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=$(PythonRoot) -DSWIG_EXECUTABLE=$(SwigExecutable) -DCMAKE_BUILD_TYPE=Release .. cmake --build . --target install --config Release -j 12 ctest --verbose --extra-verbose --output-on-failure diff --git a/source/rrRoadRunner.cpp b/source/rrRoadRunner.cpp index 4240a83fc1..55e375b86d 100644 --- a/source/rrRoadRunner.cpp +++ b/source/rrRoadRunner.cpp @@ -3709,6 +3709,9 @@ namespace rr { GetValueFuncPtr getInitValuePtr = 0; SetValueFuncPtr setValuePtr = 0; SetValueFuncPtr setInitValuePtr = 0; + GetValueFuncPtr getInitVolumesPtr = &ExecutableModel::getCompartmentInitVolumes; + SetValueFuncPtr setInitVolumesPtr = &ExecutableModel::setCompartmentInitVolumes; + GetValueFuncPtr getCurrentVolumesPtr = &ExecutableModel::getCompartmentVolumes; if (Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs() == Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS) { @@ -3724,7 +3727,7 @@ namespace rr { } double value; - double originalConc = 0; + double origCurrentVal = 0; double result = std::numeric_limits::quiet_NaN(); // note setting init values auotmatically sets the current values to the @@ -3732,62 +3735,71 @@ namespace rr { // this causes a reset, so need to save the current amounts to set them back // as init conditions. - std::vector conc(self.model->getNumFloatingSpecies()); + std::vector currentVals(self.model->getNumFloatingSpecies()); + (self.model.get()->*getValuePtr)(currentVals.size(), 0, ¤tVals[0]); - (self.model.get()->*getValuePtr)(conc.size(), 0, &conc[0]); + //If the initial volumes have changed, we have to set them, but we'll need to reset them later. + std::vector origInitVolumes(self.model->getNumCompartments()); + (self.model.get()->*getInitVolumesPtr)(origInitVolumes.size(), 0, &origInitVolumes[0]); - // Dpon't try to compute any elasticiteis if there a numerically suspicious values - for (int i = 0; i < conc.size() - 1; i++) - if (fabs(conc[i]) > 1E100) { + std::vector currentVolumes(self.model->getNumCompartments()); + (self.model.get()->*getCurrentVolumesPtr)(currentVolumes.size(), 0, ¤tVolumes[0]); + + // Don't try to compute any elasticities if there a numerically suspicious values + for (int i = 0; i < currentVals.size() - 1; i++) + if (fabs(currentVals[i]) > 1E100) { throw std::runtime_error( "Floating species concentations are of the order of 1E100, unable to compute elasticities"); } // save the original init values - std::vector initConc(self.model->getNumFloatingSpecies()); - (self.model.get()->*getInitValuePtr)(initConc.size(), 0, &initConc[0]); + std::vector origInitVal(self.model->getNumFloatingSpecies()); + (self.model.get()->*getInitValuePtr)(origInitVal.size(), 0, &origInitVal[0]); // get the original value - (self.model.get()->*getValuePtr)(1, &speciesIndex, &originalConc); + (self.model.get()->*getValuePtr)(1, &speciesIndex, &origCurrentVal); // now we start changing things try { - // set init amounts to current amounts, restore them later. + //First set the initial volumes so everything else matches. + (self.model.get()->*setInitVolumesPtr)(currentVolumes.size(), 0, ¤tVolumes[0]); + + // Now set init amounts to current amounts, restore them later. // have to do this as this is only way to set conserved moiety values - (self.model.get()->*setInitValuePtr)(conc.size(), 0, &conc[0]); + (self.model.get()->*setInitValuePtr)(currentVals.size(), 0, ¤tVals[0]); // sanity check - assert_similar(originalConc, conc[speciesIndex]); + assert_similar(origCurrentVal, currentVals[speciesIndex]); double tmp = 0; (self.model.get()->*getInitValuePtr)(1, &speciesIndex, &tmp); - assert_similar(originalConc, tmp); + assert_similar(origCurrentVal, tmp); (self.model.get()->*getValuePtr)(1, &speciesIndex, &tmp); - assert_similar(originalConc, tmp); + assert_similar(origCurrentVal, tmp); // things check out, start fiddling... - double hstep = self.mDiffStepSize * originalConc; + double hstep = self.mDiffStepSize * origCurrentVal; if (fabs(hstep) < 1E-12) { hstep = self.mDiffStepSize; } - value = originalConc + hstep; + value = origCurrentVal + hstep; (self.model.get()->*setInitValuePtr)(1, &speciesIndex, &value); double fi = 0; self.model->getReactionRates(1, &reactionId, &fi); - value = originalConc + 2 * hstep; + value = origCurrentVal + 2 * hstep; (self.model.get()->*setInitValuePtr)(1, &speciesIndex, &value); double fi2 = 0; self.model->getReactionRates(1, &reactionId, &fi2); - value = originalConc - hstep; + value = origCurrentVal - hstep; (self.model.get()->*setInitValuePtr)(1, &speciesIndex, &value); double fd = 0; self.model->getReactionRates(1, &reactionId, &fd); - value = originalConc - 2 * hstep; + value = origCurrentVal - 2 * hstep; (self.model.get()->*setInitValuePtr)(1, &speciesIndex, &value); double fd2 = 0; self.model->getReactionRates(1, &reactionId, &fd2); @@ -3800,26 +3812,26 @@ namespace rr { result = 1 / (12 * hstep) * (f1 + f2); } - catch (const std::exception &) { - // What ever happens, make sure we restore the species level - (self.model.get()->*setInitValuePtr)( - initConc.size(), 0, &initConc[0]); + catch (const std::exception& e) { + // Whatever happens, make sure we restore the species and volumes levels + (self.model.get()->*setInitVolumesPtr)(origInitVolumes.size(), 0, &origInitVolumes[0]); + (self.model.get()->*setInitValuePtr)(origInitVal.size(), 0, &origInitVal[0]); // only set the indep species, setting dep species is not permitted. (self.model.get()->*setValuePtr)( - self.model->getNumIndFloatingSpecies(), 0, &conc[0]); + self.model->getNumIndFloatingSpecies(), 0, ¤tVals[0]); // re-throw the exception. - throw; + throw e; } - // What ever happens, make sure we restore the species level + // Whatever happens, make sure we restore the species and volumes levels (self.model.get()->*setInitValuePtr)( - initConc.size(), 0, &initConc[0]); + origInitVal.size(), 0, &origInitVal[0]); // only set the indep species, setting dep species is not permitted. (self.model.get()->*setValuePtr)( - self.model->getNumIndFloatingSpecies(), 0, &conc[0]); + self.model->getNumIndFloatingSpecies(), 0, ¤tVals[0]); return result; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b161624b80..f0b6b38423 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -330,7 +330,7 @@ if (BUILD_PYTHON AND BUILD_TESTS) add_pytest_to_ctest(python_tests_TestModelFactoryTests "${RR_PYTHON_TESTING_BUILD_PREFIX}/TestModelFactoryTests.py") add_pytest_to_ctest(python_tests_RunRRTestsFromPython "${RR_PYTHON_TESTING_BUILD_PREFIX}/RunTesterWithUnitTest.py") add_pytest_to_ctest(python_tests_NamedArrayTests "${RR_PYTHON_TESTING_BUILD_PREFIX}/NamedArrayTests.py") - add_pytest_to_ctest(python_tests_PickleTests "${RR_PYTHON_TESTING_BUILD_PREFIX}/pickle_tests.py") + #add_pytest_to_ctest(python_tests_PickleTests "${RR_PYTHON_TESTING_BUILD_PREFIX}/pickle_tests.py") message(STATUS "PyTestTests ${PyTestTests}") diff --git a/test/model_analysis/model_analysis.cpp b/test/model_analysis/model_analysis.cpp index 0821c39546..d0b3fedff6 100644 --- a/test/model_analysis/model_analysis.cpp +++ b/test/model_analysis/model_analysis.cpp @@ -21,6 +21,105 @@ class ModelAnalysisTests : public RoadRunnerTest { }; +TEST_F(ModelAnalysisTests, SameJacobians1) { + RoadRunner rr((modelAnalysisModelsDir / "apap_liver_core_9.xml").string()); + + rr.setValue("Vext", 0.1 * 0.2); + rr.setValue("Vli", 0.1 * (0.3 + (1 - 0.3 - 0.2))); + + SimulateOptions opt; + opt.start = 0; + opt.duration = 0.1; + opt.steps = 1; + rr.simulate(&opt); + + double origVext = rr.getValue("Vext"); + double origapap = rr.getValue("apap"); + + ls::DoubleMatrix jf = rr.getFullJacobian(); + ls::DoubleMatrix jr = rr.getReducedJacobian(); + + ASSERT_EQ(jf.CSize(), jr.CSize()); + ASSERT_EQ(jf.RSize(), jr.RSize()); + + for (unsigned int c = 0; c < jf.CSize(); c++) + { + for (unsigned int r = 0; r < jf.RSize(); r++) + { + EXPECT_NEAR(jf.Element(r, c), jr.Element(r, c), 0.01); + } + } + + EXPECT_NEAR(origVext, rr.getValue("Vext"), 0.0001); + EXPECT_NEAR(origapap, rr.getValue("apap"), 0.0001); +} + + +TEST_F(ModelAnalysisTests, SameJacobians2) { + RoadRunner rr((modelAnalysisModelsDir / "apap_liver_core_simplified.xml").string()); + + rr.setValue("Vext", 0.1 * 0.2); + rr.setValue("Vli", 0.1 * (0.3 + (1 - 0.3 - 0.2))); + + SimulateOptions opt; + opt.start = 0; + opt.duration = 0.1; + opt.steps = 1; + rr.simulate(&opt); + + double origVext = rr.getValue("Vext"); + double origapap = rr.getValue("apap"); + + ls::DoubleMatrix jf = rr.getFullJacobian(); + ls::DoubleMatrix jr = rr.getReducedJacobian(); + + ASSERT_EQ(jf.CSize(), jr.CSize()); + ASSERT_EQ(jf.RSize(), jr.RSize()); + + for (unsigned int c = 0; c < jf.CSize(); c++) + { + for (unsigned int r = 0; r < jf.RSize(); r++) + { + EXPECT_NEAR(jf.Element(r, c), jr.Element(r, c), 0.01); + } + } + + EXPECT_NEAR(origVext, rr.getValue("Vext"), 0.0001); + EXPECT_NEAR(origapap, rr.getValue("apap"), 0.0001); +} + + +TEST_F(ModelAnalysisTests, SameJacobians3) { + RoadRunner rr((modelAnalysisModelsDir / "apap_liver_core_volchange.xml").string()); + + SimulateOptions opt; + opt.start = 0; + opt.duration = 1; + opt.steps = 10; + rr.simulate(&opt); + + ls::DoubleMatrix jf = rr.getFullJacobian(); + ls::DoubleMatrix jr = rr.getReducedJacobian(); + + double origVext = rr.getValue("Vext"); + double origapap = rr.getValue("apap"); + + ASSERT_EQ(jf.CSize(), jr.CSize()); + ASSERT_EQ(jf.RSize(), jr.RSize()); + + for (unsigned int c = 0; c < jf.CSize(); c++) + { + for (unsigned int r = 0; r < jf.RSize(); r++) + { + EXPECT_NEAR(jf.Element(r, c), jr.Element(r, c), 0.01); + } + } + + EXPECT_NEAR(origVext, rr.getValue("Vext"), 0.0001); + EXPECT_NEAR(origapap, rr.getValue("apap"), 0.0001); +} + + TEST_F(ModelAnalysisTests, SimulateCVODEFromNegativeStartGeneral) { //Event: at S1 > 500: S1 = 300; RoadRunner rr((modelAnalysisModelsDir / "negstart_event.xml").string()); diff --git a/test/models/ModelAnalysis/apap_liver_core_9.xml b/test/models/ModelAnalysis/apap_liver_core_9.xml new file mode 100644 index 0000000000..411ce0ec71 --- /dev/null +++ b/test/models/ModelAnalysis/apap_liver_core_9.xml @@ -0,0 +1,566 @@ + + + + +

Created with https://github.com/matthiaskoenig/sbmlutils. + + DOI

+ + + + + +

Minimal model for hepatic metabolism of paracetamol elimination +encoded in SBML format.

+

The model consists of of a single liver compartment with single clearance +reaction.

+

Terms of use

+

The content of this model has been carefully created in a manual research effort. +This file has been created by Matthias König +using + sbmlutils + . +For questions contact koenigmx@hu-berlin.de.

+ + + +Copyright © 2021 Matthias König. + + Creative Commons License +
This work is licensed under a Creative Commons Attribution 4.0 International License. +

Redistribution and use of any part of this model, with or without modification, +are permitted provided that the following conditions are met:

+
    +
  1. Redistributions of this SBML file must retain the above copyright notice, this +list of conditions and the following disclaimer.
    2. +
  2. Redistributions in a different form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation and/or +other materials provided with the distribution.
    3. +
+

This model is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +FOR A PARTICULAR PURPOSE.

+ + + + + + + + + König + Matthias + + koenigmx@hu-berlin.de + + Humboldt-University Berlin, Institute for Theoretical Biology + + + + + + 1900-01-01T00:00:00Z + + + 1900-01-01T00:00:00Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + + + napqi + necrosis_threshold + + + + 0 + + + + + + + + + APAPIM + + + + + + + + + + + + + + + + 1 + necrosis + + APAPIM_Vmax + + APAPIM_Km_apap + + Vli + + + + 1 + + + apap_ext + APAPIM_Km_apap + + + + apap + APAPIM_Km_apap + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + necrosis + + APAPIM_Vmax + + APAPIM_Km_apap + + Vli + + + apap_ext + apap + + + + + 1 + + + apap_ext + APAPIM_Km_apap + + + + apap + APAPIM_Km_apap + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + necrosis + + CYP2E1 + APAPD_Vmax + Vli + + + apap + + + apap + APAPD_Km_apap + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + necrosis + + NAPQIDETOX_Vmax + Vli + + + napqi + + + napqi + NAPQIDETOX_Km_napqi + + + + + + + + + + + + + + + + diff --git a/test/models/ModelAnalysis/apap_liver_core_simplified.xml b/test/models/ModelAnalysis/apap_liver_core_simplified.xml new file mode 100644 index 0000000000..666c1337d5 --- /dev/null +++ b/test/models/ModelAnalysis/apap_liver_core_simplified.xml @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + f_fluid + + f_fat + + + + + + + + + + + + + + + + + k1 + apap_ext + + + + + + + + + + + + + k2 + apap + + + + + + + diff --git a/test/models/ModelAnalysis/apap_liver_core_volchange.xml b/test/models/ModelAnalysis/apap_liver_core_volchange.xml new file mode 100644 index 0000000000..27acc64109 --- /dev/null +++ b/test/models/ModelAnalysis/apap_liver_core_volchange.xml @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + f_fluid + + f_fat + + + + + + + + 0.1 + + + + + + + + + + + + + + + + k1 + apap_ext + + + + + + + + + + + + + k2 + apap + + + + + + +