27 #include "sbml/SBMLDocument.h"
28 #include "sbml/Model.h"
29 #include "sbml/FunctionDefinition.h"
30 #include "sbml/Rule.h"
31 #include "sbml/InitialAssignment.h"
32 #include "sbml/Reaction.h"
33 #include "sbml/KineticLaw.h"
34 #include "sbml/SpeciesReference.h"
35 #include "sbml/Event.h"
36 #include "sbml/Delay.h"
37 #include "sbml/Trigger.h"
38 #include "sbml/EventAssignment.h"
39 #include "sbml/StoichiometryMath.h"
40 #include "sbml/math/ASTNode.h"
41 #include "sbml/math/MathML.h"
58 mNumExceededFunctions(0),
59 mNumFailedFunctions(0),
63 mNumCOPASIFunctions(0),
64 mNumExceededCOPASIFunctions(0),
65 mNumFailedCOPASIFunctions(0),
67 mNumKineticFunctions(0),
68 mNumMassActionsKinetics(0),
69 mNumConstantFluxKinetics(0),
70 mNumMappedKineticExpressions(0),
71 mNumUnmappedKineticExpressions(0),
72 mDifferentNormalform(0),
108 while (it5 != endit5)
128 std::vector<std::string>::const_iterator it = filenames.begin(), endit = filenames.end();
132 std::cout <<
"Processing " << *it << std::endl;
138 std::cout <<
"number of COPASI function definitions: " <<
mNumCOPASIFunctions << std::endl;
145 std::cout <<
"the following " <<
mUnreadableFiles.size() <<
" files could not be read: " << std::endl;
150 std::cout << *it << std::endl;
155 std::cout <<
mNumFiles <<
" files have been processed." << std::endl;
158 std::cout <<
"number of failed function definitions: " <<
mNumFailedFunctions << std::endl;
159 std::cout <<
"number of expressions: " <<
mNumExpressions << std::endl;
160 std::cout <<
"number of exceeded expressions: " <<
mNumExceeded << std::endl;
161 std::cout <<
"number of failed expressions: " <<
mNumFailed << std::endl;
167 while (it2 != endit2)
172 while (it3 != endit3)
176 std::cout <<
"The functions \"" << it2->first <<
"\" and \"" << it3->first <<
"\" in the COPASI database are equal." << std::endl;
186 unsigned int numClassified = 0;
187 unsigned int numDubious = 0;
189 while (it5 != endit5)
194 while (it3 != endit3)
217 std::cout <<
"Number of function definitons that could be classified: " << numClassified << std::endl;
218 std::cout <<
"Number of function definitons that were classified incorrectly: " << numDubious << std::endl;
219 std::cout <<
"Number of function definitons that could not be classified: " <<
mNormalizedFunctionDefinitions.size() - numClassified - numDubious << std::endl;
221 std::cout <<
"Number of kinetic expressions that could be mapped to a function definition: " <<
mNumMappedKineticExpressions << std::endl;
223 std::cout <<
"List of the number of expressions mapped to a certain function definition: " << std::endl;
226 while (it6 != endit6)
228 std::cout << it6->first <<
" : " << it6->second << std::endl;
232 std::cout <<
"There are " <<
mSBOMap.size() <<
" different SBO Terms." << std::endl;
233 std::map<int, std::vector<CNormalFraction*> >::iterator sboIt =
mSBOMap.begin(), sboEndit =
mSBOMap.end();
235 while (sboIt != sboEndit)
237 std::cout <<
"There are " << sboIt->second.size() <<
" expressions for SBO term " << sboIt->first <<
"." << std::endl;
241 std::cout <<
"Number of kinetic expressions with sbo terms: " <<
mNumSBO << std::endl;
242 std::cout <<
"Number of kinetic expressions with sbo terms that could not be normalized to the same normalform: " <<
mDifferentNormalform << std::endl;
243 std::cout <<
"The expressions that could not be mapped are divided into " <<
mUnknownCategories.size() <<
" different expressions." << std::endl;
245 unsigned int moreThanOne = 0;
246 unsigned int moreThanFive = 0;
247 unsigned int moreThanTen = 0;
248 unsigned int moreThanTwenty = 0;
249 unsigned int numMoreThanOne = 0;
250 unsigned int numMoreThanFive = 0;
251 unsigned int numMoreThanTen = 0;
252 unsigned int numMoreThanTwenty = 0;
254 while (catIt != catEndit)
256 if (catIt->second > 20)
259 numMoreThanTwenty += catIt->second;
260 std::cout <<
"Expression with more than 20 instances: " << catIt->first->toString() << std::endl;
262 else if (catIt->second > 10)
265 numMoreThanTen += catIt->second;
267 else if (catIt->second > 5)
270 numMoreThanFive += catIt->second;
272 else if (catIt->second > 1)
275 numMoreThanOne += catIt->second;
281 std::cout << moreThanTwenty <<
" uncategorized expressions have more than 20 instances (" << numMoreThanTwenty <<
" total)" << std::endl;
282 std::cout << moreThanTen <<
" uncategorized expressions have more than 10 instances (" << numMoreThanTen <<
" total)" << std::endl;
283 std::cout << moreThanFive <<
" uncategorized expressions have more than 5 instances (" << numMoreThanFive <<
" total)" << std::endl;
284 std::cout << moreThanOne <<
" uncategorized expressions have more than 1 instance (" << numMoreThanOne <<
" total)" << std::endl;
285 std::cout <<
mUnknownCategories.size() - (moreThanTwenty + moreThanTen + moreThanFive + moreThanOne) <<
" uncategorized expression are unique." << std::endl;
287 unsigned int num = 0;
289 while (timeIt != timeEndit && (*timeIt) > 600)
297 std::cout <<
"Number of expression taking more than 10 minutes: " << num << std::endl;
301 while (timeIt != timeEndit && (*timeIt) > 60)
309 std::cout <<
"Number of expression taking more than 1 minute: " << num << std::endl;
312 while (timeIt != timeEndit && (*timeIt) > 10)
320 std::cout <<
"Number of expression taking more than 10 seconds: " << num << std::endl;
323 while (timeIt != timeEndit && (*timeIt) > 1)
331 std::cout <<
"Number of expression taking more than 1 second: " << num << std::endl;
334 while (timeIt != timeEndit && (*timeIt) > 0.1)
342 std::cout <<
"Number of expression taking more than 1/10 seconds: " << num << std::endl;
345 while (timeIt != timeEndit && (*timeIt) > 0.01)
353 std::cout <<
"Number of expression taking more than 1/100 seconds: " << num << std::endl;
356 while (timeIt != timeEndit)
364 std::cout <<
"Number of expression taking less than 1/100 seconds: " << num << std::endl;
392 if (pDocument != NULL)
394 const Model* pModel = pDocument->getModel();
416 const InitialAssignment* pInitialAssignment = NULL;
417 const ASTNode* pMath = NULL;
418 ASTNode* pNewMath = NULL;
421 unsigned int i, iMax = pModel->getListOfInitialAssignments()->size();
423 for (i = 0; i < iMax; ++i)
425 pInitialAssignment = pModel->getInitialAssignment(i);
426 assert(pInitialAssignment != NULL);
427 pMath = pInitialAssignment->getMath();
432 assert(pNewMath != NULL);
439 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
443 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
447 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
454 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
459 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
465 assert(pFraction != NULL);
471 const Rule* pRule = NULL;
472 iMax = pModel->getListOfRules()->size();
474 for (i = 0; i < iMax; ++i)
476 pRule = pModel->getRule(i);
477 assert(pRule != NULL);
478 pMath = pRule->getMath();
483 assert(pNewMath != NULL);
489 assert(pFraction != NULL);
492 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
496 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
500 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
507 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
512 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
522 const Reaction* pReaction = NULL;
523 iMax = pModel->getListOfReactions()->size();
525 for (i = 0; i < iMax; ++i)
527 pReaction = pModel->getReaction(i);
528 assert(pReaction != NULL);
529 const KineticLaw* pLaw = pReaction->getKineticLaw();
530 int sboTerm = pLaw->getSBOTerm();
534 pMath = pLaw->getMath();
541 assert(result ==
true);
542 assert(pNewMath != NULL);
549 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
553 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
557 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
561 assert(pFraction != NULL);
564 std::string
id = pReaction->getId();
579 assert(pCOPASIReaction != NULL);
581 if (dynamic_cast<const CMassAction*>(pCOPASIReaction->
getFunction()) == NULL)
636 while (catIt != catEndit)
675 std::vector<CNormalFraction*> v;
676 v.push_back(pFraction);
683 std::cout <<
"Some expressions for SBO term " << sboTerm <<
" can not be normalized to the same normal form." << std::endl;
687 mSBOMap[sboTerm].push_back(pFraction);
695 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
700 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
709 const SpeciesReference* pSpeciesReference = NULL;
712 unsigned j, jMax = pReaction->getListOfReactants()->size();
714 for (j = 0; j < jMax; ++j)
716 pSpeciesReference = pReaction->getReactant(j);
717 assert(pSpeciesReference != NULL);
719 if (pSpeciesReference->isSetStoichiometryMath())
721 const StoichiometryMath* pSMath = pSpeciesReference->getStoichiometryMath();
722 assert(pSMath != NULL);
723 pMath = pSMath->getMath();
728 assert(pNewMath != NULL);
735 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
739 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
743 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
747 assert(pFraction != NULL);
752 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
757 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
768 jMax = pReaction->getListOfProducts()->size();
770 for (j = 0; j < jMax; ++j)
772 pSpeciesReference = pReaction->getProduct(j);
773 assert(pSpeciesReference != NULL);
775 if (pSpeciesReference->isSetStoichiometryMath())
777 const StoichiometryMath* pSMath = pSpeciesReference->getStoichiometryMath();
778 assert(pSMath != NULL);
779 pMath = pSMath->getMath();
784 assert(pNewMath != NULL);
791 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
795 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
799 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
803 assert(pFraction != NULL);
808 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
813 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
825 const Event* pEvent = NULL;
826 const Trigger* pTrigger = NULL;
827 const Delay* pDelay = NULL;
828 iMax = pModel->getListOfEvents()->size();
830 for (i = 0; i < iMax; ++i)
832 pEvent = pModel->getEvent(i);
833 assert(pEvent != NULL);
835 pTrigger = pEvent->getTrigger();
836 assert(pTrigger != NULL);
837 pMath = pTrigger->getMath();
842 assert(pNewMath != NULL);
849 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
853 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
857 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
861 assert(pFraction != NULL);
866 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
871 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
880 if (pEvent->isSetDelay())
882 pDelay = pEvent->getDelay();
883 assert(pDelay != NULL);
884 pMath = pDelay->getMath();
889 assert(pNewMath != NULL);
896 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
900 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
904 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
908 assert(pFraction != NULL);
913 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
918 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
928 unsigned int j, jMax = pEvent->getListOfEventAssignments()->size();
929 const EventAssignment* pEventAssignment = NULL;
931 for (j = 0; j < jMax; ++j)
933 pEventAssignment = pEvent->getEventAssignment(j);
934 assert(pEventAssignment != NULL);
935 pMath = pEventAssignment->getMath();
940 assert(pNewMath != NULL);
947 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
951 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
955 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
959 assert(pFraction != NULL);
964 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewMath) <<
"\"." << std::endl;
969 std::cerr <<
"expression \"" << writeMathMLToString(pNewMath) <<
"\" could not be normalized." << std::endl;
985 const FunctionDefinition* pFunDef = NULL;
986 const ASTNode* pRoot = NULL;
987 ASTNode* pNewRoot = NULL;
989 unsigned int i, iMax = pModel->getListOfFunctionDefinitions()->size();
991 for (i = 0; i < iMax; ++i)
993 pFunDef = pModel->getFunctionDefinition(i);
994 pRoot = pFunDef->getMath();
997 if (pRoot != NULL && pRoot->getNumChildren() > 0)
1000 const ASTNode* pMath = pRoot->getChild(pRoot->getNumChildren() - 1);
1001 assert(pMath != NULL);
1003 assert(pNewRoot != NULL);
1010 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
1014 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
1018 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
1022 assert(pFraction != NULL);
1027 std::cerr <<
"recursion limit exceeded for expression \"" << writeMathMLToString(pNewRoot) <<
"\"." << std::endl;
1032 std::cerr <<
"expression \"" << writeMathMLToString(pNewRoot) <<
"\" could not be normalized." << std::endl;
1251 assert(pFunctionDB != NULL);
1253 unsigned int i = 0, iMax = loadedFunctions.
size();
1258 assert(pTree != NULL);
1260 if (dynamic_cast<const CFunction*>(pTree) != NULL && dynamic_cast<const CMassAction*>(pTree) == NULL)
1266 assert(pExpanded != NULL);
1270 double timeDiff =
mTV2.tv_sec -
mTV1.tv_sec;
1274 timeDiff -= (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
1278 timeDiff += (
mTV2.tv_usec -
mTV1.tv_usec) / 1e6;
1283 assert(pFraction != NULL);
1289 std::cerr <<
"recursion limit exceeded for expression:" << std::endl;
1295 std::cerr <<
"expression could not be normalized: " << std::endl;
1308 std::vector<std::pair<ASTNode*, unsigned int> > nodeStack;
1309 std::map<std::string, std::string> replacementMap;
1312 unsigned int compartment_index = 1;
1313 unsigned int species_index = 1;
1314 unsigned int reaction_index = 1;
1315 unsigned int parameter_index = 1;
1316 std::ostringstream os;
1318 ASTNode* pCurrent = pNode;
1320 nodeStack.push_back(std::pair<ASTNode*, unsigned int>(NULL, 0));
1321 std::map<std::string, std::string>::const_iterator pos;
1323 while (pCurrent != NULL)
1327 if (pCurrent == nodeStack.back().first)
1329 ++nodeStack.back().second;
1331 if (nodeStack.back().second == pCurrent->getNumChildren())
1334 nodeStack.erase(--nodeStack.end());
1335 pCurrent = nodeStack.back().first;
1340 pCurrent = pCurrent->getChild(nodeStack.back().second);
1345 iMax = pCurrent->getNumChildren();
1351 nodeStack.push_back(std::pair<ASTNode*, unsigned int>(pCurrent, 0));
1352 pCurrent = pCurrent->getChild(0);
1356 if (pCurrent->isName())
1363 std::string
id = pCurrent->getName();
1365 pos = replacementMap.find(
id);
1367 if (pos == replacementMap.end())
1370 if ((pReaction->getKineticLaw() != NULL && pReaction->getKineticLaw()->getParameter(
id) != NULL) || pModel->getParameter(
id) != NULL)
1373 os <<
"K_" << parameter_index;
1374 replacementMap.insert(std::pair<std::string, std::string>(
id, os.str()));
1375 pCurrent->setName(os.str().c_str());
1379 else if (pModel->getSpecies(
id) != NULL)
1382 os <<
"S_" << species_index;
1383 replacementMap.insert(std::pair<std::string, std::string>(
id, os.str()));
1384 pCurrent->setName(os.str().c_str());
1387 else if (pModel->getCompartment(
id) != NULL)
1390 os <<
"C_" << compartment_index;
1391 replacementMap.insert(std::pair<std::string, std::string>(
id, os.str()));
1392 pCurrent->setName(os.str().c_str());
1393 ++compartment_index;
1395 else if (pModel->getReaction(
id) != NULL)
1398 os <<
"R_" << reaction_index;
1399 replacementMap.insert(std::pair<std::string, std::string>(
id, os.str()));
1400 pCurrent->setName(os.str().c_str());
1411 pCurrent->setName(pos->second.c_str());
1416 pCurrent = nodeStack.back().first;
1431 std::vector<std::string> filenames;
1435 filenames.push_back(argv[i]);
1440 test.
run(filenames);
1444 std::cerr <<
"Usage: stresstest SBMLFILE1 [SBMLFILE2 ...]" << std::endl;
void normalizeFunctionDB()
std::set< double > mProcessTimes
CCopasiVectorN< CFunction > & loadedFunctions()
std::map< int, std::vector< CNormalFraction * > > mSBOMap
SBMLDocument * getCurrentSBMLDocument()
static CNormalFraction * normAndSimplifyReptdly(const CEvaluationTree *tree0, unsigned int depth=0)
static bool normalize_names(ASTNode *pNode, const Reaction *pReaction, const Model *pModel)
const std::string & getObjectName() const
unsigned int mNumExceededCOPASIFunctions
unsigned int mNumExceeded
virtual size_t size() const
unsigned int mNumFailedCOPASIFunctions
unsigned int mNumExpressions
CCopasiDataModel * mpDataModel
unsigned int mNumExceededFunctions
void normalizeMath(const std::string &filename)
int main(int argc, char **argv)
bool importSBML(const std::string &fileName, CProcessReport *pImportHandler=NULL, const bool &deleteOldData=true)
void normalizeAndSimplifyExpressions(const Model *pModel)
bool are_equal(const CNormalFraction *pLHS, const CNormalFraction *pRHS)
unsigned int mDifferentNormalform
unsigned int mNumUnmappedKineticExpressions
const CFunction * getFunction() const
CNormalFraction * create_simplified_normalform(const ASTNode *pSource)
unsigned int mNumCOPASIFunctions
std::map< std::string, unsigned int > mExpressionMappings
unsigned int mNumMappedKineticExpressions
void run(const std::vector< std::string > &filenames)
std::multimap< std::string, CNormalFraction * > mNormalizedCOPASIFunctionDefinitions
unsigned int mNumKineticFunctions
static CFunctionDB * getFunctionList()
static CCopasiDataModel * addDatamodel()
ASTNode * create_expression(const ASTNode *pSource, const ListOfFunctionDefinitions *pFunctions)
unsigned int mNumFunctionDefinitions
unsigned int mNumFailedFunctions
void normalizeAndSimplifyFunctionDefinitions(const Model *pModel)
CCopasiVectorNS< CReaction > & getReactions()
std::vector< std::string > mUnreadableFiles
static void init(int argc, char *argv[], const bool &withGui=false)
std::vector< std::pair< std::string, CNormalFraction * > > mNormalizedFunctionDefinitions
ASTNode * expand_function_calls(const ASTNode *pNode, const ListOfFunctionDefinitions *pFunctions)
void printRecursively(std::ostream &os, int indent=0) const
std::vector< CNormalFraction * > mNormalizedExpressions
CEvaluationNode * getRoot()
std::map< CNormalFraction *, int > mUnknownCategories