Skip to content

Commit

Permalink
Merge pull request #912 from sys-bio/issue-910-fix-jacobians
Browse files Browse the repository at this point in the history 
Issue #910: Fix jacobians
  • Loading branch information
@luciansmith
luciansmith committed Nov 10, 2021
2 parents 83f063e + aa4dcbb commit 000f05b
Show file tree
Hide file tree
Showing 7 changed files with 860 additions and 32 deletions.
6 changes: 3 additions & 3 deletions azure-pipelines.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
68 changes: 40 additions & 28 deletions source/rrRoadRunner.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::int32_t>() ==
Config::ROADRUNNER_JACOBIAN_MODE_AMOUNTS) {
Expand All @@ -3724,70 +3727,79 @@ namespace rr {
}

double value;
double originalConc = 0;
double origCurrentVal = 0;
double result = std::numeric_limits<double>::quiet_NaN();

// note setting init values auotmatically sets the current values to the
// init values

// this causes a reset, so need to save the current amounts to set them back
// as init conditions.
std::vector<double> conc(self.model->getNumFloatingSpecies());
std::vector<double> currentVals(self.model->getNumFloatingSpecies());
(self.model.get()->*getValuePtr)(currentVals.size(), 0, &currentVals[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<double> 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<double> currentVolumes(self.model->getNumCompartments());
(self.model.get()->*getCurrentVolumesPtr)(currentVolumes.size(), 0, &currentVolumes[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<double> initConc(self.model->getNumFloatingSpecies());
(self.model.get()->*getInitValuePtr)(initConc.size(), 0, &initConc[0]);
std::vector<double> 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, &currentVolumes[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, &currentVals[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);
Expand All @@ -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, &currentVals[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, &currentVals[0]);

return result;
}
Expand Down
2 changes: 1 addition & 1 deletion test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}")

Expand Down
99 changes: 99 additions & 0 deletions test/model_analysis/model_analysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down

0 comments on commit 000f05b

Please sign in to comment.