Skip to content

Commit

Permalink
Update tests for new Jacobian calcs
Browse files Browse the repository at this point in the history
  • Loading branch information
@NSR-Physiome
NSR-Physiome committed Jan 13, 2025
1 parent 9f747bc commit 9091c58
Show file tree
Hide file tree
Showing 3 changed files with 324 additions and 14 deletions.
236 changes: 236 additions & 0 deletions test/model_analysis/model_analysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -451,7 +451,243 @@ TEST_F(ModelAnalysisTests, SameJacobians3) {
EXPECT_NEAR(origVext, rr.getValue("Vext"), 0.0001);
EXPECT_NEAR(origapap, rr.getValue("apap"), 0.0001);
}
TEST_F(ModelAnalysisTests, jacobian_multiComp_diffVols_JACOB_MODE_FALSE) {
// Rxn over two compartments of differing volumes.
std::int32_t savedJacobianMode = Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs<std::int32_t>();
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, FALSE);
rr::RoadRunner rr((modelAnalysisModelsDir / "jacobian_multi_comps_2sp_3rxn.xml").string());

ls::DoubleMatrix ues = rr.getUnscaledElasticityMatrix();
ls::DoubleMatrix jacobian = rr.getFullJacobian();

// Expected unscaled elasticity matrix:
// S1, S2
// J0 [[ 0, 0],
// J1 [0.333333, 0],
// J2 [0, 0.961538]]

// Expected full jacobian matrix:
// S1 S2
// S1 [ -0.333333333, 0 ]
// S2 [ 0.192307692, -0.961538462 ]

int rows = 3;
int cols = 2;
ASSERT_EQ(ues.CSize(), cols);
ASSERT_EQ(ues.RSize(), rows);

double** ues_expect = new double* [rows];
for (int i = 0; i < rows; ++i) {
ues_expect[i] = new double[cols]();
}
ues_expect[0][0] = 0;
ues_expect[0][1] = 0;
ues_expect[1][0] = 0.33333;
ues_expect[1][1] = 0;
ues_expect[2][0] = 0;
ues_expect[2][1] = 0.961538;

ls::DoubleMatrix* ues_expected = new ls::DoubleMatrix(3, 2);
ues_expected->initializeFrom2DMatrix(ues_expect, 3, 2);

for (unsigned int c = 0; c < ues.CSize(); c++)
{
for (unsigned int r = 0; r < ues.RSize(); r++)
{
EXPECT_NEAR(ues.Element(r, c), ues_expected->Element(r, c), 0.0001);
}
}

double** jac_expect = new double* [2]; // 2x2
for (int i = 0; i < 2; ++i) {
jac_expect[i] = new double[cols]();
}
jac_expect[0][0] = -0.333333333;
jac_expect[0][1] = 0;
jac_expect[1][0] = 0.192307692;
jac_expect[1][1] = -0.961538462;

ls::DoubleMatrix* jac_expected = new ls::DoubleMatrix(2, 2);
jac_expected->initializeFrom2DMatrix(jac_expect, 2, 2);

ASSERT_EQ(jacobian.CSize(), 2);
ASSERT_EQ(jacobian.RSize(), 2);
for (unsigned int c = 0; c < jacobian.CSize(); c++)
{
for (unsigned int r = 0; r < jacobian.RSize(); r++)
{
// std::cout <<"Jacobian row, col: " << r << ", " << c << ": " << jacobian.Element(r, c) << std::endl;
// std::cout << "Expected Jacobian: " << jac_expected->Element(r, c) << std::endl;
EXPECT_NEAR(jacobian.Element(r, c), jac_expected->Element(r, c), 0.0001);
}
}
// Put back default JACOBIAN_MODE:
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, savedJacobianMode);
}

TEST_F(ModelAnalysisTests, jacobian_multiComp_diffVols_JACOB_MODE_TRUE) {
// Rxn over two compartments of differing volumes with JACOBIAN_MODE flag set to TRUE.
std::int32_t savedJacobianMode = Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs<std::int32_t>();
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, TRUE);
rr::RoadRunner rr((modelAnalysisModelsDir / "jacobian_multi_comps_2sp_3rxn.xml").string());

ls::DoubleMatrix ues = rr.getUnscaledElasticityMatrix();
ls::DoubleMatrix jacobian = rr.getFullJacobian();

// Expected unscaled elasticity matrix:
// S1, S2
// J0 [[ 0, 0],
// J1 [0.5, 0],
// J2 [0, 2.5]]

// Expected full jacobian matrix:
// S1 S2
// S1 [ -0.333333333, 0 ]
// S2 [ 0.192307692, -0.961538462 ]

int rows = 3;
int cols = 2;
ASSERT_EQ(ues.CSize(), cols);
ASSERT_EQ(ues.RSize(), rows);

double** ues_expect = new double* [rows];
for (int i = 0; i < rows; ++i) {
ues_expect[i] = new double[cols]();
}
ues_expect[0][0] = 0;
ues_expect[0][1] = 0;
ues_expect[1][0] = 0.5;
ues_expect[1][1] = 0;
ues_expect[2][0] = 0;
ues_expect[2][1] = 2.5;

ls::DoubleMatrix* ues_expected = new ls::DoubleMatrix(3, 2);
ues_expected->initializeFrom2DMatrix(ues_expect, 3, 2);

for (unsigned int c = 0; c < ues.CSize(); c++)
{
for (unsigned int r = 0; r < ues.RSize(); r++)
{
EXPECT_NEAR(ues.Element(r, c), ues_expected->Element(r, c), 0.0001);
}
}

double** jac_expect = new double* [2]; // 2x2
for (int i = 0; i < 2; ++i) {
jac_expect[i] = new double[cols]();
}
jac_expect[0][0] = -0.333333333;
jac_expect[0][1] = 0;
jac_expect[1][0] = 0.192307692;
jac_expect[1][1] = -0.961538462;

ls::DoubleMatrix* jac_expected = new ls::DoubleMatrix(2, 2);
jac_expected->initializeFrom2DMatrix(jac_expect, 2, 2);

ASSERT_EQ(jacobian.CSize(), 2);
ASSERT_EQ(jacobian.RSize(), 2);
for (unsigned int c = 0; c < jacobian.CSize(); c++)
{
for (unsigned int r = 0; r < jacobian.RSize(); r++)
{
// std::cout <<"Jacobian row, col: " << r << ", " << c << ": " << jacobian.Element(r, c) << std::endl;
// std::cout << "Expected Jacobian: " << jac_expected->Element(r, c) << std::endl;
EXPECT_NEAR(jacobian.Element(r, c), jac_expected->Element(r, c), 0.0001);
}
}
// Put back default selected JACOBIAN_MODE:
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, savedJacobianMode);
}

TEST_F(ModelAnalysisTests, jacobian_multiComp_diffVols_CompsEqOne) {
// Rxn over two compartments of volume one.
std::int32_t savedJacobianMode = Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs<std::int32_t>();
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, FALSE);
rr::RoadRunner rr((modelAnalysisModelsDir / "jacobian_multi_comps_2sp_3rxn.xml").string());
ASSERT_EQ(rr.getNumberOfCompartments(), 2);
for (int i = 0; i < rr.getNumberOfCompartments(); i++) {
rr.setCompartmentByIndex(i, 1);
}

ls::DoubleMatrix jacobian = rr.getFullJacobian();
// Expected full jacobian matrix:
// S1 S2
// S1 [ -0.5, 0 ]
// S2 [ 0.5, -2.5 ]

int rows = 2;
int cols = 2;
double** jac_expect = new double* [rows]; // 2x2
for (int i = 0; i < rows; ++i) {
jac_expect[i] = new double[cols]();
}
jac_expect[0][0] = -0.5;
jac_expect[0][1] = 0;
jac_expect[1][0] = 0.5;
jac_expect[1][1] = -2.5;

ls::DoubleMatrix* jac_expected = new ls::DoubleMatrix(rows, cols);
jac_expected->initializeFrom2DMatrix(jac_expect, rows, cols);

ASSERT_EQ(jacobian.CSize(), cols);
ASSERT_EQ(jacobian.RSize(), rows);
for (unsigned int c = 0; c < jacobian.CSize(); c++)
{
for (unsigned int r = 0; r < jacobian.RSize(); r++)
{
// std::cout <<"Jacobian row, col: " << r << ", " << c << ": " << jacobian.Element(r, c) << std::endl;
// std::cout << "Expected Jacobian: " << jac_expected->Element(r, c) << std::endl;
EXPECT_NEAR(jacobian.Element(r, c), jac_expected->Element(r, c), 0.0001);
}
}
// Put back default JACOBIAN_MODE:
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, savedJacobianMode);
}

TEST_F(ModelAnalysisTests, jacobian_multiComp_diffVols_CompsEqSame) {
// Rxn over two compartments of volume two.
std::int32_t savedJacobianMode = Config::getValue(Config::ROADRUNNER_JACOBIAN_MODE).getAs<std::int32_t>();
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, FALSE);
rr::RoadRunner rr((modelAnalysisModelsDir / "jacobian_multi_comps_2sp_3rxn.xml").string());
ASSERT_EQ(rr.getNumberOfCompartments(), 2);
for (int i = 0; i < rr.getNumberOfCompartments(); i++) {
rr.setCompartmentByIndex(i, 2);
}

ls::DoubleMatrix jacobian = rr.getFullJacobian();
// Expected full jacobian matrix:
// S1 S2
// S1 [ -0.25, 0 ]
// S2 [ 0.25, -1.25]

int rows = 2;
int cols = 2;
double** jac_expect = new double* [rows]; // 2x2
for (int i = 0; i < rows; ++i) {
jac_expect[i] = new double[cols]();
}
jac_expect[0][0] = -0.25;
jac_expect[0][1] = 0;
jac_expect[1][0] = 0.25;
jac_expect[1][1] = -1.25;

ls::DoubleMatrix* jac_expected = new ls::DoubleMatrix(rows, cols);
jac_expected->initializeFrom2DMatrix(jac_expect, rows, cols);

ASSERT_EQ(jacobian.CSize(), cols);
ASSERT_EQ(jacobian.RSize(), rows);
for (unsigned int c = 0; c < jacobian.CSize(); c++)
{
for (unsigned int r = 0; r < jacobian.RSize(); r++)
{
// std::cout <<"Jacobian row, col: " << r << ", " << c << ": " << jacobian.Element(r, c) << std::endl;
// std::cout << "Expected Jacobian: " << jac_expected->Element(r, c) << std::endl;
EXPECT_NEAR(jacobian.Element(r, c), jac_expected->Element(r, c), 0.0001);
}
}
// Put back default JACOBIAN_MODE:
Config::setValue(Config::ROADRUNNER_JACOBIAN_MODE, savedJacobianMode);
}

TEST_F(ModelAnalysisTests, SimulateCVODEFromNegativeStartGeneral) {
//Event: at S1 > 500: S1 = 300;
Expand Down
74 changes: 74 additions & 0 deletions test/models/ModelAnalysis/jacobian_multi_comps_2sp_3rxn.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v2.13.2 with libSBML version 5.19.5. -->
<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1">
<model metaid="multiCompart_2sp_3rxn" id="multiCompart_2sp_3rxn">
<listOfCompartments>
<compartment id="c1" spatialDimensions="3" size="1.5" constant="true"/>
<compartment id="c2" spatialDimensions="3" size="2.6" constant="true"/>
</listOfCompartments>
<listOfSpecies>
<species id="X0" compartment="c1" initialConcentration="10" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
<species id="S1" compartment="c1" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species id="S2" compartment="c2" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species id="X1" compartment="c2" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
</listOfSpecies>
<listOfParameters>
<parameter id="k1" value="0.4" constant="true"/>
<parameter id="k2" value="0.5" constant="true"/>
<parameter id="k3" value="2.5" constant="true"/>
</listOfParameters>
<listOfReactions>
<reaction id="J0" reversible="true" fast="false">
<listOfReactants>
<speciesReference species="X0" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
<speciesReference species="S1" stoichiometry="1" constant="true"/>
</listOfProducts>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> k1 </ci>
<ci> X0 </ci>
</apply>
</math>
</kineticLaw>
</reaction>
<reaction id="J1" reversible="true" fast="false">
<listOfReactants>
<speciesReference species="S1" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
<speciesReference species="S2" stoichiometry="1" constant="true"/>
</listOfProducts>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> k2 </ci>
<ci> S1 </ci>
</apply>
</math>
</kineticLaw>
</reaction>
<reaction id="J2" reversible="true" fast="false">
<listOfReactants>
<speciesReference species="S2" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
<speciesReference species="X1" stoichiometry="1" constant="true"/>
</listOfProducts>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> k3 </ci>
<ci> S2 </ci>
</apply>
</math>
</kineticLaw>
</reaction>
</listOfReactions>
</model>
</sbml>
28 changes: 14 additions & 14 deletions test/rrtest_files/jacobian_1.rrtest
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,13 +158,13 @@

[Full Jacobian]
-2.15 0.27 0.09
1.1 -1.07 0.09
0.4 0.4 -0.64
0.55 -0.535 0.045
0.2 0.2 -0.32

[Amount Jacobian]
-2.15 0.135 0.045
1.1 -0.535 0.045
0.4 0.2 -0.32
-2.15 0.27 0.09
0.55 -0.535 0.045
0.2 0.2 -0.32

[Steady State Fluxes]
J1 = 0.427394807520143
Expand All @@ -181,24 +181,24 @@ J4 = 0.213697403760072
"[S1] [S2] [S3]"

[Individual Eigenvalues]
S1= -2.37948 0
S2= -1.02054 0
S3= -0.459983 0
S1= -2.24281 0
S2= -0.237138 0
S3= -0.525048 0

[Eigenvalue Matrix]
-2.37948 0
-1.02054 0
-0.459983 0
-2.24281 0
-0.237138 0
-0.525048 0

[Individual Amount Eigenvalues]
S1= -2.242810 0
S2= -0.525048 0
S3= -0.237138 0
S2= -0.237138 0
S3= -0.525048 0

[Eigenvalue Amount Matrix]
-2.242810 0
-0.525048 0
-0.237138 0
-0.525048 0

[Unscaled Elasticity Matrix]
-0.25 0 0
Expand Down

0 comments on commit 9091c58

Please sign in to comment.