Skip to content

Commit 245c8a0

Browse files
skeatingskeating
committed
Improves SBML rate rule conversion
Addresses issues in SBML rate rule conversion by refining stoichiometry handling and variable replacement. The changes ensure correct reaction stoichiometry by considering the `useStoichiometryFromMath` option. It also prevents replacement of CSymbolFunction to avoid unexpected behavior.
1 parent 3e84bf5 commit 245c8a0

File tree

7 files changed

+134
-74
lines changed

7 files changed

+134
-74
lines changed

‎src/sbml/conversion/SBMLConverter.cpp‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -330,14 +330,15 @@ ASTNode* SBMLConverter::replaceAssignedVariablesWithMath(ASTNode* original)
330330
if (ar != NULL && ar->isSetMath() == true)
331331
{
332332
ASTNode* arMath = ar->getMath()->deepCopy();
333+
if (arMath->isCSymbolFunction() || arMath->getType() == AST_NAME_TIME)
334+
{
335+
// we cannot replace a CSymbolFunction as it may have different values in different contexts
336+
delete arMath;
337+
continue;
338+
}
333339
ASTNode* variable = new ASTNode(AST_NAME);
334340
variable->setName(ar->getVariable().c_str());
335-
cout << "assignment rule " << i << ": " << SBML_formulaToL3String(arMath) << " variable "
336-
<< SBML_formulaToL3String(variable) << " original " << SBML_formulaToL3String(newMath) << endl;
337-
338341
analyser.replaceExpressionInNodeWithNode(newMath, variable, arMath);
339-
cout << "afterwards assignment rule " << i << ": " << SBML_formulaToL3String(arMath) << " variable "
340-
<< SBML_formulaToL3String(variable) << " original " << SBML_formulaToL3String(newMath) << endl;
341342
}
342343
}
343344
return newMath;

‎src/sbml/conversion/SBMLRateRuleConverter.cpp‎

Lines changed: 65 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -235,7 +235,7 @@ SBMLRateRuleConverter::getDefaultProperties() const
235235
{
236236
prop.addOption("inferReactions", true,
237237
"Infer reactions from rateRules in the model");
238-
prop.addOption("useStoichiometryFromMath", false,
238+
prop.addOption("useStoichiometryFromMath", true,
239239
"If a number appears in the math use it as the stoichiometry");
240240

241241
init = true;
@@ -278,7 +278,7 @@ SBMLRateRuleConverter::setDocument(SBMLDocument* doc)
278278
{
279279
if (SBMLConverter::setDocument(doc) == LIBSBML_OPERATION_SUCCESS)
280280
{
281-
if (mDocument != NULL)
281+
if (mDocument != NULL && mDocument->getModel() != NULL)
282282
{
283283
mOriginalModel = mDocument->getModel()->clone();
284284
return LIBSBML_OPERATION_SUCCESS;
@@ -468,11 +468,11 @@ void SBMLRateRuleConverter::populateTerms()
468468
// Fages algo 3.6 Step 2
469469
createTerms(node);
470470
}
471-
for (unsigned int n = 0; n < mTerms.size(); n++)
472-
{
473-
ASTNode* node = mTerms.at(n);
474-
cout << "Term " << n << ": " << SBML_formulaToL3String(node) << endl;
475-
}
471+
//for (unsigned int n = 0; n < mTerms.size(); n++)
472+
//{
473+
// ASTNode* node = mTerms.at(n);
474+
// cout << "Term " << n << ": " << SBML_formulaToL3String(node) << endl;
475+
//}
476476
//print_vectors(mCoefficients);
477477
}
478478

@@ -561,18 +561,15 @@ SBMLRateRuleConverter::determineCoefficient(ASTNode* ode, unsigned int termN, do
561561
// first child should be a number
562562
// take it out of the term
563563
// it will be used as a coefficient
564-
//if (ode_node->getType() == AST_REAL)
565-
//{
566-
// coeff = ode_node->getValue();
567-
// found = true;
568-
//}
564+
// unless we do not want it used as stoichiometry
565+
569566
if (ode_node->getType() == AST_TIMES && ode_node->getNumChildren() > 0)
570567
{
571568
if (ode_node->getChild(0)->isNumber())
572569
{
573570
coeff = ode_node->getChild(0)->getValue();
574571
ode_node->removeChild(0, true);
575-
// we don't want to be loved with the times node if it has only one child
572+
// we don't want to be left with the times node if it has only one child
576573
if (ode_node->getNumChildren() == 1)
577574
{
578575
ASTNode* child = ode_node->getChild(0)->deepCopy();
@@ -951,10 +948,10 @@ SBMLRateRuleConverter::populateInitialODEinfo()
951948
}
952949
}
953950

954-
for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++)
955-
{
956-
cout << mODEs[odeIndex].first << ": " << SBML_formulaToL3String(mODEs[odeIndex].second) << endl;
957-
}
951+
//for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++)
952+
//{
953+
// cout << mODEs[odeIndex].first << ": " << SBML_formulaToL3String(mODEs[odeIndex].second) << endl;
954+
//}
958955
}
959956

960957

@@ -1019,7 +1016,7 @@ SBMLRateRuleConverter::populateReactionCoefficients()
10191016

10201017
bool SBMLRateRuleConverter::useStoichiometryFromMath()
10211018
{
1022-
bool value = false;
1019+
bool value = true;
10231020
if (getProperties() == NULL || getProperties()->hasOption("useStoichiometryFromMath") == false)
10241021
{
10251022
return value;
@@ -1149,10 +1146,53 @@ SBMLRateRuleConverter::dealWithSpecies()
11491146
}
11501147
}
11511148

1152-
void
1153-
SBMLRateRuleConverter::dealWithStoichiometry()
1149+
bool
1150+
SBMLRateRuleConverter::needToAdjustStoichiometryAndMath(unsigned int odeNumber)
11541151
{
1155-
// TO DO
1152+
bool adjust = false;
1153+
if (mODEs.at(odeNumber).second->getType() == AST_TIMES &&
1154+
mODEs.at(odeNumber).second->getNumChildren() > 0 &&
1155+
mODEs.at(odeNumber).second->getChild(0)->isNumber() &&
1156+
(mODEs.at(odeNumber).second->getChild(0)->getValue() > 1.0 ||
1157+
mODEs.at(odeNumber).second->getChild(0)->getValue() < -1.0)
1158+
)
1159+
{
1160+
adjust = true;
1161+
}
1162+
return adjust;
1163+
}
1164+
1165+
double
1166+
SBMLRateRuleConverter::dealWithStoichiometry(double stoichiometry, ASTNode& math, unsigned int odeNumber)
1167+
{
1168+
if (useStoichiometryFromMath())
1169+
{
1170+
return stoichiometry;
1171+
}
1172+
else
1173+
{
1174+
if (stoichiometry > 1.0 && needToAdjustStoichiometryAndMath(odeNumber) )
1175+
{
1176+
1177+
ASTNode* newMath = mODEs.at(odeNumber).second->deepCopy();
1178+
if (newMath->getType() == AST_TIMES && newMath->getNumChildren() > 0 &&
1179+
newMath->getChild(0)->isNumber())
1180+
{
1181+
double value = newMath->getChild(0)->getValue();
1182+
if (value < 0.0)
1183+
{
1184+
// we don't want a value that is less than zero
1185+
ASTNode* child = newMath->getChild(0)->deepCopy();
1186+
child->setValue(-1.0 * value);
1187+
newMath->replaceChild(0, child);
1188+
}
1189+
}
1190+
math = *newMath;
1191+
delete newMath;
1192+
stoichiometry = 1.0;
1193+
}
1194+
}
1195+
return stoichiometry;
11561196
}
11571197

11581198
void
@@ -1162,6 +1202,7 @@ SBMLRateRuleConverter::createReactions()
11621202
char number[4];
11631203
for (setCoeffIt it = mCoefficients.begin(); it != mCoefficients.end(); ++it)
11641204
{
1205+
ASTNode* math = (it->first)->deepCopy();
11651206
Reaction *r = mDocument->getModel()->createReaction();
11661207
r->setReversible(false);
11671208
r->setFast(false);
@@ -1175,7 +1216,7 @@ SBMLRateRuleConverter::createReactions()
11751216
double stoichiometry = 1.0;
11761217
if (mReactants[i][j] > 0)
11771218
{
1178-
stoichiometry = mReactants[i][j];
1219+
stoichiometry = dealWithStoichiometry(mReactants[i][j], *math, j);
11791220
SpeciesReference *sr = r->createReactant();
11801221
sr->setSpecies(mODEs[j].first);
11811222
sr->setStoichiometry(stoichiometry);
@@ -1184,7 +1225,7 @@ SBMLRateRuleConverter::createReactions()
11841225
}
11851226
if (mProducts[i][j] > 0)
11861227
{
1187-
stoichiometry = mProducts[i][j];
1228+
stoichiometry = dealWithStoichiometry(mProducts[i][j], *math, j);
11881229
SpeciesReference *sr = r->createProduct();
11891230
sr->setSpecies(mODEs[j].first);
11901231
sr->setStoichiometry(stoichiometry);
@@ -1201,7 +1242,7 @@ SBMLRateRuleConverter::createReactions()
12011242
if (itemAdded && !r->isSetKineticLaw())
12021243
{
12031244
KineticLaw *kl = r->createKineticLaw();
1204-
kl->setMath(it->first);
1245+
kl->setMath(math);
12051246
}
12061247

12071248
// check whetherkinetic law uses a species not listed as p/r/m

‎src/sbml/conversion/SBMLRateRuleConverter.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -263,10 +263,11 @@ class LIBSBML_EXTERN SBMLRateRuleConverter : public SBMLConverter
263263
void dealWithSpecies();
264264
void createReactions();
265265
void removeRules();
266-
void dealWithStoichiometry();
266+
double dealWithStoichiometry(double stoichiometry, ASTNode& math, unsigned int odeNumber);
267+
bool needToAdjustStoichiometryAndMath(unsigned int odeNumber);
267268

269+
// member variables populated during analysis;
268270

269-
// member variables populated during analysis
270271
pairODEs mODEs;
271272
std::vector<ASTNode*> mTerms;
272273
setCoeff mCoefficients;

‎src/sbml/conversion/test/TestRoundtripConverter.cpp‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -110,7 +110,7 @@ bool test_rule_to_reaction_to_rule(const std::string& raterule_file,
110110
rule_rn_converter->setDocument(d);
111111
if (rule_rn_converter->convert() != LIBSBML_OPERATION_SUCCESS)
112112
{
113-
cout << "rule_reaction_rule converter rule->reaction failed" << endl;
113+
//cout << "rule_reaction_rule converter rule->reaction failed" << endl;
114114
delete d_rn;
115115
delete d_rule;
116116
delete d_rule1;
@@ -134,7 +134,7 @@ bool test_rule_to_reaction_to_rule(const std::string& raterule_file,
134134
rn_rule_converter->setDocument(d);
135135
if (rn_rule_converter->convert() != LIBSBML_OPERATION_SUCCESS)
136136
{
137-
cout << "rule_reaction_rule converter reaction->rule failed" << endl;
137+
//cout << "rule_reaction_rule converter reaction->rule failed" << endl;
138138
delete d_rn;
139139
delete d_rule;
140140
delete d_rule1;
@@ -172,7 +172,7 @@ bool test_reaction_to_rule_to_reaction(const std::string& raterule_file,
172172
rn_rule_converter->setDocument(d);
173173
if (rn_rule_converter->convert() != LIBSBML_OPERATION_SUCCESS)
174174
{
175-
cout << "reaction_rule_reaction converter reaction->rule failed" << endl;
175+
//cout << "reaction_rule_reaction converter reaction->rule failed" << endl;
176176
delete d_rn;
177177
delete d_rn1;
178178
delete d_rule;
@@ -196,7 +196,7 @@ bool test_reaction_to_rule_to_reaction(const std::string& raterule_file,
196196
rule_rn_converter->setDocument(d);
197197
if (rule_rn_converter->convert() != LIBSBML_OPERATION_SUCCESS)
198198
{
199-
cout << "reaction_rule_reaction converter rule->reaction failed" << endl;
199+
//cout << "reaction_rule_reaction converter rule->reaction failed" << endl;
200200
delete d_rn;
201201
delete d_rule;
202202
delete d_rn1;
@@ -825,7 +825,7 @@ create_suite_TestSBMLRoundtripConverter(void)
825825
tcase_add_test(tcase, test_roundtrip_01_reverse);
826826
tcase_add_test(tcase, test_roundtrip_02);
827827
tcase_add_test(tcase, test_roundtrip_02_reverse);
828-
//tcase_add_test(tcase, test_roundtrip_03); // rule to reaction problem with stoichiometry
828+
tcase_add_test(tcase, test_roundtrip_03); // rule to reaction problem with stoichiometry
829829
tcase_add_test(tcase, test_roundtrip_03_reverse);
830830
tcase_add_test(tcase, test_roundtrip_04);
831831
tcase_add_test(tcase, test_roundtrip_04_reverse);

‎src/sbml/conversion/test/TestRunner.c‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -121,11 +121,13 @@ main (void)
121121
int num_failed;
122122

123123
setTestDataDirectory();
124-
SRunner *runner = srunner_create(create_suite_TestSBMLReactionConverter());
124+
////SRunner* runner = srunner_create(create_suite_TestSBMLRateRuleConverter());
125+
// // SRunner* runner = srunner_create(create_suite_TestSBMLReactionConverter());
126+
// SRunner* runner = srunner_create(create_suite_TestSBMLRateRuleConverter());
127+
// srunner_add_suite(runner, create_suite_TestSBMLReactionConverter());
128+
// srunner_add_suite(runner, create_suite_TestSBMLRoundtripConverter());
125129

126-
//srunner_add_suite(runner, create_suite_TestSBMLRateRuleConverter());
127-
128-
/*SRunner *runner = srunner_create( create_suite_TestConversionOption() );
130+
SRunner *runner = srunner_create( create_suite_TestConversionOption() );
129131
srunner_add_suite( runner, create_suite_TestSBMLRuleConverter () );
130132
srunner_add_suite( runner, create_suite_TestConversionProperties () );
131133
srunner_add_suite( runner, create_suite_TestSBMLConverterRegistry () );
@@ -136,10 +138,10 @@ main (void)
136138
srunner_add_suite( runner, create_suite_TestStripPackageConverter () );
137139
srunner_add_suite( runner, create_suite_TestLevelVersionConverter () );
138140
srunner_add_suite( runner, create_suite_TestRateOfConverter () );
141+
srunner_add_suite(runner, create_suite_TestExpressionAnalyser());
139142
srunner_add_suite( runner, create_suite_TestSBMLRateRuleConverter () );
140143
srunner_add_suite( runner, create_suite_TestSBMLRoundtripConverter () );
141144
srunner_add_suite( runner, create_suite_TestSBMLReactionConverter () );
142-
srunner_add_suite( runner, create_suite_TestExpressionAnalyser () );*/
143145

144146
/* srunner_set_fork_status(runner, CK_NOFORK); */
145147

0 commit comments

Comments
 (0)