diff --git a/src/sbml/SBMLErrorTable.h b/src/sbml/SBMLErrorTable.h index 4f07775286..53d2a1d6e4 100644 --- a/src/sbml/SBMLErrorTable.h +++ b/src/sbml/SBMLErrorTable.h @@ -11047,8 +11047,7 @@ static const sbmlErrorTableEntry errorTable[] = "LibSBML expected to read the annotation into a ModelHistory " "object. Unfortunately, some attributes were not present or correct " "and the resulting ModelHistory object will not correctly " - "produce the annotation. This functionality will be improved in " - "later versions of libSBML. ", + "produce the annotation. ", {"", "", "L2V2 Section 6.3", @@ -11076,8 +11075,7 @@ static const sbmlErrorTableEntry errorTable[] = "LibSBML expected to read the annotation into a ModelHistory " "object. Unfortunately, some attributes were not present or correct " "and the resulting ModelHistory object is NULL. Thus it will fail to " - "produce the annotation. This functionality will be improved in " - "later versions of libSBML. ", + "produce the annotation. ", {"", "", "L2V2 Section 6.3", diff --git a/src/sbml/SBase.cpp b/src/sbml/SBase.cpp index 4e1ada2ac5..3ba470d240 100644 --- a/src/sbml/SBase.cpp +++ b/src/sbml/SBase.cpp @@ -4576,6 +4576,7 @@ SBase::read (XMLInputStream& stream) else if ( next.isStart() ) { const std::string nextName = next.getName(); + const std::string nextURI = next.getURI(); #if 0 cout << "[DEBUG] SBase::read " << nextName << " uri " << stream.peek().getURI() << endl; @@ -4620,7 +4621,7 @@ SBase::read (XMLInputStream& stream) || readAnnotation(stream) || readNotes(stream) )) { - logUnknownElement(nextName, getLevel(), getVersion()); + logUnknownElement(nextName, getLevel(), getVersion(), nextURI); stream.skipPastEnd( stream.next() ); } } @@ -5343,11 +5344,17 @@ SBase::logUnknownAttribute( const string& attribute, void SBase::logUnknownElement( const string& element, const unsigned int level, - const unsigned int version ) + const unsigned int version, + const string& URI) { bool logged = false; ostringstream msg; + // if we have an unknown element that is not in an SBML name space + // this is perfectly allowed XML + if (getPackageName() == "core" && SBMLNamespaces::getSBMLNamespaceURI(level, version) != URI) + return; + if (level > 2 && getTypeCode() == SBML_LIST_OF) { int tc = static_cast(this)->getItemTypeCode(); diff --git a/src/sbml/SBase.h b/src/sbml/SBase.h index 5ed23d9dcb..df17ac3e5c 100644 --- a/src/sbml/SBase.h +++ b/src/sbml/SBase.h @@ -3532,7 +3532,8 @@ newModel.addSpecies(s1); */ void logUnknownElement( const std::string& element, const unsigned int level, - const unsigned int version ); + const unsigned int version, + const std::string& URI = NULL); /** diff --git a/src/sbml/conversion/ExpressionAnalyser.cpp b/src/sbml/conversion/ExpressionAnalyser.cpp index 66cf981944..d181ec2d61 100644 --- a/src/sbml/conversion/ExpressionAnalyser.cpp +++ b/src/sbml/conversion/ExpressionAnalyser.cpp @@ -1,57 +1,57 @@ - -/** - * @file ExpressionAnalyser.cpp - * @brief Implementation of ExpressionAnalyser - * @author Sarah Keating - * @author Alessandro Felder - * - * - */ - - -#include - - -#include -#include -#include -#include -#include -#include -#include - -#ifdef __cplusplus - -using namespace std; - -LIBSBML_CPP_NAMESPACE_BEGIN - - + +/** + * @file ExpressionAnalyser.cpp + * @brief Implementation of ExpressionAnalyser + * @author Sarah Keating + * @author Alessandro Felder + * + * + */ + + +#include + + +#include +#include +#include +#include +#include +#include +#include + +#ifdef __cplusplus + +using namespace std; + +LIBSBML_CPP_NAMESPACE_BEGIN + + ExpressionAnalyser::ExpressionAnalyser() { } @@ -61,739 +61,759 @@ ExpressionAnalyser::ExpressionAnalyser(Model * m, pairODEs odes) { mModel = m; mODEs = odes; - SBMLTransforms::mapComponentValues(mModel); - mModel->populateAllElementIdList(); - mNewVarName = "newVar"; - mNewVarCount = 1; + SBMLTransforms::mapComponentValues(mModel); + mModel->populateAllElementIdList(); + mNewVarName = "newVar"; + mNewVarCount = 1; +} + +ExpressionAnalyser::ExpressionAnalyser(const ExpressionAnalyser& orig) : + mModel( orig.mModel) +{ +} + +/* +* Assignment operator for SBMLLevelVersionConverter. +*/ +ExpressionAnalyser& +ExpressionAnalyser::operator=(const ExpressionAnalyser& rhs) +{ + if (&rhs != this) + { + mModel = rhs.mModel; + } + + return *this; +} + +ExpressionAnalyser* +ExpressionAnalyser::clone() const +{ + return new ExpressionAnalyser(*this); +} + +/* + * Destroy this object. + */ +ExpressionAnalyser::~ExpressionAnalyser () +{ + for (std::vector >::iterator it = mODEs.begin(); it != mODEs.end(); ++it) + { + if (it->second != NULL) + { + delete it->second; + it->second = NULL; + } + } + mODEs.clear(); + SBMLTransforms::clearComponentValues(mModel); +} + +/* +* Set ode pairs +*/ +int +ExpressionAnalyser::setODEPairs(std::vector< std::pair< std::string, ASTNode*> > odes) +{ + mODEs = odes; + return LIBSBML_OPERATION_SUCCESS; +} + + +/* +* Set ode model +*/ +int +ExpressionAnalyser::setModel(Model* model) +{ + SBMLTransforms::clearComponentValues(mModel); + mModel = model; + SBMLTransforms::mapComponentValues(model); + return LIBSBML_OPERATION_SUCCESS; +} + +/* +* Check whether the expression has a parent expression which may already have been analysed +* in which case we do not need to re analyse the child expression +* e.g. if we have k-x-y do not need to analyse k-x +*/ +bool +ExpressionAnalyser::hasExpressionAlreadyRecorded(SubstitutionValues_t* value) +{ + for (unsigned int i = mExpressions.size(); i > 0; i--) + { + SubstitutionValues_t* exp = mExpressions.at(i - 1); + + // if we already have parent eg K+v-x-y dont record children eg K+v-x + std::pair parent = getParentNode(value->current, exp->current); + if (parent.first != NULL) + { + return true; + } + switch (value->type) + { + case TYPE_K_MINUS_X_MINUS_Y: + if (value->k_value == exp->k_value && + value->x_value == exp->x_value && + value->y_value == exp->y_value && + value->dxdt_expression == exp->dxdt_expression && + value->dydt_expression == exp->dydt_expression && + value->type == exp->type) + { + return true; + } + break; + case TYPE_K_PLUS_V_MINUS_X_MINUS_Y: + if (value->k_value == exp->k_value && + value->x_value == exp->x_value && + value->y_value == exp->y_value && + value->dxdt_expression == exp->dxdt_expression && + value->dydt_expression == exp->dydt_expression && + value->v_expression == exp->v_expression && + value->type == exp->type) + { + return true; + } + break; + case TYPE_K_MINUS_X_PLUS_W_MINUS_Y: + if (value->k_value == exp->k_value && + value->x_value == exp->x_value && + value->y_value == exp->y_value && + value->dxdt_expression == exp->dxdt_expression && + value->dydt_expression == exp->dydt_expression && + value->w_expression == exp->w_expression && + value->type == exp->type) + { + return true; + } + break; + case TYPE_K_MINUS_X: + if (value->k_value == exp->k_value && + value->x_value == exp->x_value && + value->dxdt_expression == exp->dxdt_expression && + value->type == exp->type) + { + return true; + } + break; + case TYPE_K_PLUS_V_MINUS_X: + if (value->k_value == exp->k_value && + value->x_value == exp->x_value && + value->dxdt_expression == exp->dxdt_expression && + value->v_expression == exp->v_expression && + value->type == exp->type) + { + return true; + } + break; + case TYPE_MINUS_X_PLUS_Y: + if (value->x_value == exp->x_value && + value->y_value == exp->y_value && + value->dxdt_expression == exp->dxdt_expression && + value->dydt_expression == exp->dydt_expression && + value->type == exp->type) + { + return true; + } + break; + default: + break; + } + } + return false; +} + + +bool +ExpressionAnalyser::analyseNode(ASTNode* node, SubstitutionValues_t *value) +{ + unsigned int numChildren = node->getNumChildren(); + ASTNodeType_t type = node->getType(); + ASTNode* rightChild = node->getRightChild(); + ASTNode* leftChild = node->getLeftChild(); + //node->printMath(); + switch (type) + { + case AST_PLUS: + // -x+y node binary; plus; left child type minus; rightchild var/const + // + + // - y + // x + if (numChildren != 2 || rightChild->getType() != AST_NAME + || leftChild->getType() != AST_MINUS + || leftChild->getNumChildren() != 1) + return false; + + // if we get to this point, the only thing left to check is + // whether the ->left->right grandchild (the x in -x+y) is a variable species. + if (isVariableSpeciesOrParameter(leftChild->getChild(0))) + { + value->x_value = leftChild->getChild(0)->getName(); + value->y_value = rightChild->getName(); + value->dydt_expression = getODEFor(rightChild->getName()); + value->dxdt_expression = getODEFor(leftChild->getChild(0)->getName()); + value->type = TYPE_MINUS_X_PLUS_Y; + value->current = node; + return true; + } + break; + + case AST_MINUS: + // - - + // k x - y + // k x + // k-x or k-x-y node binary; right child (x,y) is variable + // - - - + // + x - y + y + // k v + x - w + // k v k x + // k+v-x; right child x var; left child plus node with left child k constant + // k+v-x-y; rightchild var y; left child minus == k+v-x + // k-x+w-y; rightchild var y; left child plus with left child == k-x + if (numChildren != 2 || !isVariableSpeciesOrParameter(rightChild)) + return false; + // if left child is numerical constant or a parameter and right child variable, it IS k-x + if (isNumericalConstantOrConstantParameter(leftChild, isNumber) + && isVariableSpeciesOrParameter(rightChild)) + { + + if (isNumber) + { + value->k_value = "number"; + value->k_real_value = leftChild->getValue(); + } + value->x_value = rightChild->getName(); + value->dxdt_expression = getODEFor(rightChild->getName()); + value->type = TYPE_K_MINUS_X; + value->current = node; + return true; + } + // left child + with it's left child const we have k+v-x + // left child + with it's left child k-x we have k-x+w-y + else if (leftChild->getType() == AST_PLUS) + { + if (isNumericalConstantOrConstantParameter(leftChild->getChild(0), isNumber)) + { + value->k_value = leftChild->getChild(0)->getName(); + value->x_value = rightChild->getName(); + value->dxdt_expression = getODEFor(rightChild->getName()); + value->v_expression = leftChild->getChild(1); + value->type = TYPE_K_PLUS_V_MINUS_X; + value->current = node; + return true; + } + else if (analyseNode(leftChild->getChild(0), value) && value->type == TYPE_K_MINUS_X) + { + value->y_value = rightChild->getName(); + value->dydt_expression = getODEFor(rightChild->getName()); + value->w_expression = leftChild->getChild(1); + value->type = TYPE_K_MINUS_X_PLUS_W_MINUS_Y; + value->current = node; + return true; + } + } + else if (leftChild->getType() == AST_MINUS + && isVariableSpeciesOrParameter(leftChild->getRightChild())) + { + // if left child is k+v-x we have k+v-x-y + // or if left child is k-x we have k-x-y + if (analyseNode(leftChild, value)) + { + if (value->type == TYPE_K_PLUS_V_MINUS_X) + { + value->type = TYPE_K_PLUS_V_MINUS_X_MINUS_Y; + value->y_value = rightChild->getName(); + value->dydt_expression = getODEFor(value->y_value); + value->current = node; + return true; + } + else if (value->type == TYPE_K_MINUS_X) + { + value->y_value = rightChild->getName(); + value->dydt_expression = getODEFor(rightChild->getName()); + value->type = TYPE_K_MINUS_X_MINUS_Y; + value->current = node; + return true; + } + } + return false; + } + break; + default: + return false; + } + return false; +} + +/* +* Return the ODE for the given variable +* or an ASTNode representing zero if there is no time derivative +*/ +ASTNode* +ExpressionAnalyser::getODEFor(std::string name) +{ + for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++) + { + std::pair ode = mODEs.at(odeIndex); + if (name == ode.first) + { + return ode.second; + } + } + ASTNode* zero = new ASTNode(AST_REAL); + zero->setValue(0.0); + return zero->deepCopy(); +} +void +ExpressionAnalyser::analyse(bool minusXPlusYOnly) +{ + for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++) + { + std::pair ode = mODEs.at(odeIndex); + ASTNode* odeRHS = ode.second; + odeRHS->reduceToBinary(); + List* operators = odeRHS->getListOfNodes((ASTNodePredicate)ASTNode_isOperator); + ListIterator it = operators->begin(); + + while (it != operators->end()) + { + ASTNode* currentNode = (ASTNode*)*it; + cout << SBML_formulaToL3String(currentNode) << endl; + if (minusXPlusYOnly && currentNode->getType() != AST_PLUS) + { + it++; + continue; + } + SubstitutionValues_t* value = new SubstitutionValues_t; + value->type = TYPE_UNKNOWN; + value->k_real_value = util_NaN(); + value->dxdt_expression = NULL; + value->dydt_expression = NULL; + value->v_expression = NULL; + value->w_expression = NULL; + if (analyseNode(currentNode, value)) + { + value->odeIndex = odeIndex; + if (!hasExpressionAlreadyRecorded(value)) + { + mExpressions.push_back(value); + } + } + it++; + } + } +} + +void +ExpressionAnalyser::detectHiddenSpecies(List * hiddenSpecies) +{ + // find -x+y and replace with y-x + analyse(true); + reorderMinusXPlusYIteratively(); + mExpressions.clear(); + + // find cases of k-x/k-x-y/k+v-x/k+v-x-y/k-x+w-y + analyse(); + + for (unsigned int i = 0; i < mExpressions.size(); i++) + { + SubstitutionValues_t *exp = mExpressions.at(i); + ASTNode* currentNode = exp->current; + for (unsigned int j = 0; j < mODEs.size(); j++) + { + std::pair ode = mODEs.at(j); + ASTNode* odeRHS = ode.second; + int index = parameterAlreadyCreated(exp); + if (index >= 0) + { + exp->z_value = mExpressions.at(index)->z_value; + replaceExpressionWithNewParameter(odeRHS, exp); + } + else + { + std::string zName = getUniqueNewParameterName(); + exp->z_value = zName; + replaceExpressionWithNewParameter(odeRHS, exp); + } +// cout << "ode in main: " << SBML_formulaToL3String(odeRHS) << endl; + } + } + addParametersAndRateRules(hiddenSpecies); } -ExpressionAnalyser::ExpressionAnalyser(const ExpressionAnalyser& orig) : - mModel( orig.mModel) -{ -} - -/* -* Assignment operator for SBMLLevelVersionConverter. -*/ -ExpressionAnalyser& -ExpressionAnalyser::operator=(const ExpressionAnalyser& rhs) -{ - if (&rhs != this) - { - mModel = rhs.mModel; - } - - return *this; -} - -ExpressionAnalyser* -ExpressionAnalyser::clone() const -{ - return new ExpressionAnalyser(*this); -} - -/* - * Destroy this object. - */ -ExpressionAnalyser::~ExpressionAnalyser () -{ - for (std::vector >::iterator it = mODEs.begin(); it != mODEs.end(); ++it) - { - if (it->second != NULL) - { - delete it->second; - it->second = NULL; - } - } - mODEs.clear(); - SBMLTransforms::clearComponentValues(mModel); -} - -/* -* Set ode pairs -*/ -int -ExpressionAnalyser::setODEPairs(std::vector< std::pair< std::string, ASTNode*> > odes) -{ - mODEs = odes; - return LIBSBML_OPERATION_SUCCESS; -} - - -/* -* Set ode model -*/ -int -ExpressionAnalyser::setModel(Model* model) -{ - SBMLTransforms::clearComponentValues(mModel); - mModel = model; - SBMLTransforms::mapComponentValues(model); - return LIBSBML_OPERATION_SUCCESS; -} - -/* -* Check whether the expression has a parent expression which may already have been analysed -* in which case we do not need to re analyse the child expression -* e.g. if we have k-x-y do not need to analyse k-x -*/ -bool -ExpressionAnalyser::hasExpressionAlreadyRecorded(SubstitutionValues_t* value) -{ - for (unsigned int i = mExpressions.size(); i > 0; i--) - { - SubstitutionValues_t* exp = mExpressions.at(i - 1); - - // if we already have parent eg K+v-x-y dont record children eg K+v-x - std::pair parent = getParentNode(value->current, exp->current); - if (parent.first != NULL) - { - return true; - } - switch (value->type) - { - case TYPE_K_MINUS_X_MINUS_Y: - if (value->k_value == exp->k_value && - value->x_value == exp->x_value && - value->y_value == exp->y_value && - value->dxdt_expression == exp->dxdt_expression && - value->dydt_expression == exp->dydt_expression && - value->type == exp->type) - { - return true; - } - break; - case TYPE_K_PLUS_V_MINUS_X_MINUS_Y: - if (value->k_value == exp->k_value && - value->x_value == exp->x_value && - value->y_value == exp->y_value && - value->dxdt_expression == exp->dxdt_expression && - value->dydt_expression == exp->dydt_expression && - value->v_expression == exp->v_expression && - value->type == exp->type) - { - return true; - } - break; - case TYPE_K_MINUS_X_PLUS_W_MINUS_Y: - if (value->k_value == exp->k_value && - value->x_value == exp->x_value && - value->y_value == exp->y_value && - value->dxdt_expression == exp->dxdt_expression && - value->dydt_expression == exp->dydt_expression && - value->w_expression == exp->w_expression && - value->type == exp->type) - { - return true; - } - break; - case TYPE_K_MINUS_X: - if (value->k_value == exp->k_value && - value->x_value == exp->x_value && - value->dxdt_expression == exp->dxdt_expression && - value->type == exp->type) - { - return true; - } - break; - case TYPE_K_PLUS_V_MINUS_X: - if (value->k_value == exp->k_value && - value->x_value == exp->x_value && - value->dxdt_expression == exp->dxdt_expression && - value->v_expression == exp->v_expression && - value->type == exp->type) - { - return true; - } - break; - case TYPE_MINUS_X_PLUS_Y: - if (value->x_value == exp->x_value && - value->y_value == exp->y_value && - value->dxdt_expression == exp->dxdt_expression && - value->dydt_expression == exp->dydt_expression && - value->type == exp->type) - { - return true; - } - break; - default: - break; - } - } - return false; -} - - -bool -ExpressionAnalyser::analyseNode(ASTNode* node, SubstitutionValues_t *value) -{ - unsigned int numChildren = node->getNumChildren(); - ASTNodeType_t type = node->getType(); - ASTNode* rightChild = node->getRightChild(); - ASTNode* leftChild = node->getLeftChild(); - //node->printMath(); - switch (type) - { - case AST_PLUS: - // -x+y node binary; plus; left child type minus; rightchild var/const - // + - // - y - // x - if (numChildren != 2 || rightChild->getType() != AST_NAME - || leftChild->getType() != AST_MINUS - || leftChild->getNumChildren() != 1) - return false; - - // if we get to this point, the only thing left to check is - // whether the ->left->right grandchild (the x in -x+y) is a variable species. - if (isVariableSpeciesOrParameter(leftChild->getChild(0))) - { - value->x_value = leftChild->getChild(0)->getName(); - value->y_value = rightChild->getName(); - value->dydt_expression = getODEFor(rightChild->getName()); - value->dxdt_expression = getODEFor(leftChild->getChild(0)->getName()); - value->type = TYPE_MINUS_X_PLUS_Y; - value->current = node; - return true; - } - break; - - case AST_MINUS: - // - - - // k x - y - // k x - // k-x or k-x-y node binary; right child (x,y) is variable - // - - - - // + x - y + y - // k v + x - w - // k v k x - // k+v-x; right child x var; left child plus node with left child k constant - // k+v-x-y; rightchild var y; left child minus == k+v-x - // k-x+w-y; rightchild var y; left child plus with left child == k-x - if (numChildren != 2 || !isVariableSpeciesOrParameter(rightChild)) - return false; - // if left child is numerical constant or a parameter and right child variable, it IS k-x - if (isNumericalConstantOrConstantParameter(leftChild) - && isVariableSpeciesOrParameter(rightChild)) - { - value->k_value = leftChild->getName(); - value->x_value = rightChild->getName(); - value->dxdt_expression = getODEFor(rightChild->getName()); - value->type = TYPE_K_MINUS_X; - value->current = node; - return true; - } - // left child + with it's left child const we have k+v-x - // left child + with it's left child k-x we have k-x+w-y - else if (leftChild->getType() == AST_PLUS) - { - if (isNumericalConstantOrConstantParameter(leftChild->getChild(0))) - { - value->k_value = leftChild->getChild(0)->getName(); - value->x_value = rightChild->getName(); - value->dxdt_expression = getODEFor(rightChild->getName()); - value->v_expression = leftChild->getChild(1); - value->type = TYPE_K_PLUS_V_MINUS_X; - value->current = node; - return true; - } - else if (analyseNode(leftChild->getChild(0), value) && value->type == TYPE_K_MINUS_X) - { - value->y_value = rightChild->getName(); - value->dydt_expression = getODEFor(rightChild->getName()); - value->w_expression = leftChild->getChild(1); - value->type = TYPE_K_MINUS_X_PLUS_W_MINUS_Y; - value->current = node; - return true; - } - } - else if (leftChild->getType() == AST_MINUS - && isVariableSpeciesOrParameter(leftChild->getRightChild())) - { - // if left child is k+v-x we have k+v-x-y - // or if left child is k-x we have k-x-y - if (analyseNode(leftChild, value)) - { - if (value->type == TYPE_K_PLUS_V_MINUS_X) - { - value->type = TYPE_K_PLUS_V_MINUS_X_MINUS_Y; - value->y_value = rightChild->getName(); - value->dydt_expression = getODEFor(value->y_value); - value->current = node; - return true; - } - else if (value->type == TYPE_K_MINUS_X) - { - value->y_value = rightChild->getName(); - value->dydt_expression = getODEFor(rightChild->getName()); - value->type = TYPE_K_MINUS_X_MINUS_Y; - value->current = node; - return true; - } - } - return false; - } - break; - default: - return false; - } - return false; -} - -/* -* Return the ODE for the given variable -* or an ASTNode representing zero if there is no time derivative -*/ -ASTNode* -ExpressionAnalyser::getODEFor(std::string name) -{ - for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++) - { - std::pair ode = mODEs.at(odeIndex); - if (name == ode.first) - { - return ode.second; - } - } - ASTNode* zero = new ASTNode(AST_REAL); - zero->setValue(0.0); - return zero->deepCopy(); -} -void -ExpressionAnalyser::analyse(bool minusXPlusYOnly) -{ - for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++) - { - std::pair ode = mODEs.at(odeIndex); - ASTNode* odeRHS = ode.second; - odeRHS->reduceToBinary(); - List* operators = odeRHS->getListOfNodes((ASTNodePredicate)ASTNode_isOperator); - ListIterator it = operators->begin(); - - while (it != operators->end()) - { - ASTNode* currentNode = (ASTNode*)*it; - if (minusXPlusYOnly && currentNode->getType() != AST_PLUS) - { - it++; - continue; - } - SubstitutionValues_t* value = new SubstitutionValues_t; - value->type = TYPE_UNKNOWN; - value->dxdt_expression = NULL; - value->dydt_expression = NULL; - value->v_expression = NULL; - value->w_expression = NULL; - if (analyseNode(currentNode, value)) - { - value->odeIndex = odeIndex; - if (!hasExpressionAlreadyRecorded(value)) - { - mExpressions.push_back(value); - } - } - it++; - } - } -} - -void -ExpressionAnalyser::detectHiddenSpecies(List * hiddenSpecies) -{ - // find -x+y and replace with y-x - analyse(true); - reorderMinusXPlusYIteratively(); - mExpressions.clear(); - - // find cases of k-x/k-x-y/k+v-x/k+v-x-y/k-x+w-y - analyse(); - - for (unsigned int i = 0; i < mExpressions.size(); i++) - { - SubstitutionValues_t *exp = mExpressions.at(i); - ASTNode* currentNode = exp->current; - for (unsigned int j = 0; j < mODEs.size(); j++) - { - std::pair ode = mODEs.at(j); - ASTNode* odeRHS = ode.second; - int index = parameterAlreadyCreated(exp); - if (index >= 0) - { - exp->z_value = mExpressions.at(index)->z_value; - replaceExpressionWithNewParameter(odeRHS, exp); - } - else - { - std::string zName = getUniqueNewParameterName(); - exp->z_value = zName; - replaceExpressionWithNewParameter(odeRHS, exp); - } -// cout << "ode in main: " << SBML_formulaToL3String(odeRHS) << endl; - } - } - addParametersAndRateRules(hiddenSpecies); -} - -/* -* Replace a child node within a node with the given replacement mode -* -* param node ASTNode * parent node containing node to be replaced -* param replaced ASTNode * node to be replaced if found in parent node -* param replacement -*/ -void -ExpressionAnalyser::replaceExpressionInNodeWithNode(ASTNode* node, ASTNode* replaced, ASTNode* replacement) -{ - if (node == NULL) - { - return; - } - //cout << "node: " << SBML_formulaToL3String(node) << endl; - //cout << "with: " << SBML_formulaToL3String(replaced) << endl; - //cout << "by: " << SBML_formulaToL3String(replacement) << endl; - // we might be replcing the whole node - if (node == replaced) - { - replaced = node->deepCopy(); - (*node) = *replacement; - } - else - { - std::paircurrentParentAndIndex = make_pair((ASTNode*)NULL, (int)(NAN)); - ASTNode* currentParent; - int index; - do - { - currentParentAndIndex = getParentNode(replaced, node); - currentParent = currentParentAndIndex.first; - index = currentParentAndIndex.second; - if (currentParent != NULL) - { - currentParent->replaceChild(index, replacement->deepCopy(), false); - // intentionally, don't delete replacement as it's now owned by currentParent! - } - } while (currentParent != NULL); - } -} - -void -ExpressionAnalyser::replaceExpressionInNodeWithVar(ASTNode* node, ASTNode* replaced, std::string var) -{ - ASTNode* z = new ASTNode(AST_NAME); - z->setName(var.c_str()); - replaceExpressionInNodeWithNode(node, replaced, z); -} - -std::string -ExpressionAnalyser::getUniqueNewParameterName() -{ - char number[4]; - sprintf(number, "%u", mNewVarCount); - - std::string name = mNewVarName + string(number); - IdList ids = mModel->getAllElementIdList(); - while (ids.contains(name)) - { - mNewVarCount++; - sprintf(number, "%u", mNewVarCount); - name = mNewVarName + string(number); - } - return name; -} - - -void -ExpressionAnalyser::addParametersAndRateRules(List* hiddenSpecies) -{ - for (unsigned int i = 0; i < mExpressions.size(); i++) - { - SubstitutionValues_t *exp = mExpressions.at(i); - if (mModel->getParameter(exp->z_value) == NULL) - { - // create expression for z - ASTNode* kx = new ASTNode(AST_MINUS); - ASTNode* k = new ASTNode(AST_NAME); - k->setName(exp->k_value.c_str()); - ASTNode* x = new ASTNode(AST_NAME); - x->setName(exp->x_value.c_str()); - kx->addChild(k); - kx->addChild(x); - - ASTNode* zNode = new ASTNode(AST_MINUS); - - - // add raterule defining dz/dt - ASTNode* dxdt = exp->dxdt_expression->deepCopy(); - RateRule* raterule = mModel->createRateRule(); - raterule->setVariable(exp->z_value); - ASTNode* math = new ASTNode(AST_TIMES); - ASTNode* minus1 = new ASTNode(AST_REAL); - minus1->setValue(-1.0); - - ASTNode* dydt = NULL; - ASTNode* plus = NULL; - ASTNode* y = NULL; - switch (exp->type) - { - case TYPE_K_MINUS_X: - case TYPE_K_PLUS_V_MINUS_X: - // dz/dt = -dx/dt - math->addChild(minus1); - math->addChild(dxdt); - - // z = k - x - (*zNode) = *kx; - - break; - case TYPE_K_MINUS_X_MINUS_Y: - case TYPE_K_PLUS_V_MINUS_X_MINUS_Y: - case TYPE_K_MINUS_X_PLUS_W_MINUS_Y: - // dz/dt = - (dx/dt + dy/dt) - dydt = exp->dydt_expression->deepCopy(); - plus = new ASTNode(AST_PLUS); - plus->addChild(dxdt); - plus->addChild(dydt); - math->addChild(minus1); - math->addChild(plus); - - // z = k-x-y - y = new ASTNode(AST_NAME); - y->setName(exp->y_value.c_str()); - zNode->addChild(kx); - zNode->addChild(y); - - break; - default: - break; - } - raterule->setMath(math); - - // introduce z - Parameter* zParam = mModel->createParameter(); - zParam->setId(exp->z_value); - zParam->setConstant(false); - zParam->setValue(SBMLTransforms::evaluateASTNode(zNode, mModel)); - hiddenSpecies->add(zParam); - - delete zNode; - delete math; //its children dxdt and minus1 deleted as part of this. - } - } -} - - -void -ExpressionAnalyser::replaceExpressionWithNewParameter(ASTNode* ode, SubstitutionValues_t* exp) -{ - if (exp->type == TYPE_K_MINUS_X || exp->type == TYPE_K_MINUS_X_MINUS_Y) - { - replaceExpressionInNodeWithVar(ode, exp->current, exp->z_value); - //cout << "ode in new param var: " << SBML_formulaToL3String(ode) << endl; - for (unsigned int i = 0; i < mExpressions.size(); i++) - { - SubstitutionValues_t *thisexp = mExpressions.at(i); - - if (thisexp->dxdt_expression != NULL) - { - replaceExpressionInNodeWithVar(thisexp->dxdt_expression, exp->current, exp->z_value); - } - if (thisexp->dydt_expression != NULL) - { - replaceExpressionInNodeWithVar(thisexp->dydt_expression, exp->current, exp->z_value); - } - } - - } - if (exp->type == TYPE_K_PLUS_V_MINUS_X || exp->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y) - { - ASTNode* replacement = new ASTNode(AST_PLUS); - ASTNode* z = new ASTNode(AST_NAME); - z->setName(exp->z_value.c_str()); - ASTNode *v = exp->v_expression->deepCopy(); - replacement->addChild(z); - replacement->addChild(v); - replaceExpressionInNodeWithNode(ode, exp->current, replacement); - //cout << "ode in new param node: " << SBML_formulaToL3String(ode) << endl; - for (unsigned int i = 0; i < mExpressions.size(); i++) - { - SubstitutionValues_t *thisexp = mExpressions.at(i); - - if (thisexp->dxdt_expression != NULL) - { - //cout << "dxdt_b4: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; - replaceExpressionInNodeWithNode(thisexp->dxdt_expression, exp->current, replacement); - //cout << "dxdt: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; - } - if (thisexp->dydt_expression != NULL) - { - //cout << "dydt_b4: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; - replaceExpressionInNodeWithNode(thisexp->dydt_expression, exp->current, replacement); - //cout << "dydt: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; - } - } - } - if (exp->type == TYPE_K_MINUS_X_PLUS_W_MINUS_Y) - { - ASTNode* replacement = new ASTNode(AST_PLUS); - ASTNode* z = new ASTNode(AST_NAME); - z->setName(exp->z_value.c_str()); - ASTNode *v = exp->w_expression->deepCopy(); - replacement->addChild(z); - replacement->addChild(v); - //cout << "ode in new param node: " << SBML_formulaToL3String(ode) << endl; - //cout << "current in new param node: " << SBML_formulaToL3String(exp->current) << endl; - //cout << "replace in new param node: " << SBML_formulaToL3String(replacement) << endl; - replaceExpressionInNodeWithNode(ode, exp->current, replacement); - //cout << "ode in new param node: " << SBML_formulaToL3String(ode) << endl; - for (unsigned int i = 0; i < mExpressions.size(); i++) - { - SubstitutionValues_t *thisexp = mExpressions.at(i); - - if (thisexp->dxdt_expression != NULL) - { - //cout << "dxdt_b4: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; - replaceExpressionInNodeWithNode(thisexp->dxdt_expression, exp->current, replacement); - //cout << "dxdt: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; - } - if (thisexp->dydt_expression != NULL) - { - //cout << "dydt_b4: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; - replaceExpressionInNodeWithNode(thisexp->dydt_expression, exp->current, replacement); - //cout << "dydt: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; - } - } - } -} - - -bool -ExpressionAnalyser::addHiddenVariablesForKMinusX(List* hiddenSpecies) -{ - return true; -} - -/* -* check that the variable names in these two expressions match -*/ -bool -ExpressionAnalyser::matchesVariables(SubstitutionValues_t* exp, SubstitutionValues_t* exp1) -{ - if (exp->k_value == exp1->k_value && - exp->x_value == exp1->x_value && - exp->y_value == exp1->y_value) - { - return true; - } - return false; -} - - -/* -* Have we already created a parameter for this expression -* if so, return name -*/ -int -ExpressionAnalyser::parameterAlreadyCreated(SubstitutionValues_t* exp) -{ - int match = -1; - unsigned int i = 0; - while (!match && i < mExpressions.size()) - { - SubstitutionValues_t* thisexp = mExpressions.at(i); - if (thisexp->z_value.empty() || !matchesVariables(thisexp, exp)) - { - i++; - continue; - } - if (thisexp->type == exp->type) - { - match = i; - } - else if (thisexp->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y && exp->type == TYPE_K_MINUS_X_MINUS_Y) - { - match = i; - } - i++; - } - return match; -} - -/* - * Check whether for node is a name node representing species or a non constant parameter -*/ -bool ExpressionAnalyser::isVariableSpeciesOrParameter(ASTNode* node) -{ - if (!node->isName()) // some nodes, like * operators, don't seem to have a name in the first place - return false; - Species* species = mModel->getSpecies(node->getName()); - Parameter* parameter = mModel->getParameter(node->getName()); // some species in rate rules may be defined as variable parameters - bool isVariableSpeciesOrParameter = (species != NULL && !species->getConstant()); - bool isVariableParameter = (parameter!=NULL && !parameter->getConstant()); - return isVariableSpeciesOrParameter || isVariableParameter; -} - -/* -* Check whether for node is a name node representing a constant parameter or a numerical node -*/ -bool ExpressionAnalyser::isNumericalConstantOrConstantParameter(ASTNode* node) -{ - if (!node->isName()) // some nodes, like * operators, don't seem to have a name in the first place - return false; - Parameter* parameter = mModel->getParameter(node->getName()); - bool isConstantParameter = (parameter != NULL) && (parameter->getConstant()); - bool isNumericalConstant = node->isNumber() && node->isConstant(); - return isNumericalConstant || isConstantParameter; -} - -/* -* Reorder any instance of - x + y with y - x in the set of ODEs. -* Fages Algorithm 3.1 Step 1 -*/ -void ExpressionAnalyser::reorderMinusXPlusYIteratively() -{ - for (unsigned int i = 0; i < mExpressions.size(); i++) - { - SubstitutionValues_t* exp = mExpressions.at(i); - if (exp->type != TYPE_MINUS_X_PLUS_Y) - continue; - ASTNode* ode = (mODEs.at(exp->odeIndex)).second; - ASTNode* replacement = new ASTNode(AST_MINUS); - ASTNode* y = new ASTNode(AST_NAME); - y->setName((exp->y_value).c_str()); - ASTNode* x = new ASTNode(AST_NAME); - x->setName((exp->x_value).c_str()); - replacement->addChild(y); - replacement->addChild(x); - replaceExpressionInNodeWithNode(ode, exp->current, replacement); - } -} - -std::pair ExpressionAnalyser::getParentNode(const ASTNode* child, const ASTNode* root) -{ - //cout << "root " << SBML_formulaToL3String(root) << endl; - //cout << "child " << SBML_formulaToL3String(child) << endl; - for (unsigned int i = 0; i < root->getNumChildren(); i++) - { - if (root->getChild(i)->exactlyEqual(*(child))) - { - return std::pair(const_cast(root), i); - } - } - for (unsigned int i = 0; i < root->getNumChildren(); i++) - { - std::pair parent = getParentNode(child, root->getChild(i)); - if (parent.first != NULL) - { - return parent; - } - } - return std::pair(NULL, (int)(NAN)); -} - -/** @endcond */ - -LIBSBML_CPP_NAMESPACE_END - -#endif /* __cplusplus */ - - +/* +* Replace a child node within a node with the given replacement mode +* +* param node ASTNode * parent node containing node to be replaced +* param replaced ASTNode * node to be replaced if found in parent node +* param replacement +*/ +void +ExpressionAnalyser::replaceExpressionInNodeWithNode(ASTNode* node, ASTNode* replaced, ASTNode* replacement) +{ + if (node == NULL) + { + return; + } + //cout << "node: " << SBML_formulaToL3String(node) << endl; + //cout << "with: " << SBML_formulaToL3String(replaced) << endl; + //cout << "by: " << SBML_formulaToL3String(replacement) << endl; + // we might be replcing the whole node + if (node == replaced) + { + replaced = node->deepCopy(); + (*node) = *replacement; + } + else + { + std::paircurrentParentAndIndex = make_pair((ASTNode*)NULL, (int)(NAN)); + ASTNode* currentParent; + int index; + do + { + currentParentAndIndex = getParentNode(replaced, node); + currentParent = currentParentAndIndex.first; + index = currentParentAndIndex.second; + if (currentParent != NULL) + { + currentParent->replaceChild(index, replacement->deepCopy(), false); + // intentionally, don't delete replacement as it's now owned by currentParent! + } + } while (currentParent != NULL); + } +} + +void +ExpressionAnalyser::replaceExpressionInNodeWithVar(ASTNode* node, ASTNode* replaced, std::string var) +{ + ASTNode* z = new ASTNode(AST_NAME); + z->setName(var.c_str()); + replaceExpressionInNodeWithNode(node, replaced, z); +} + +std::string +ExpressionAnalyser::getUniqueNewParameterName() +{ + char number[4]; + sprintf(number, "%u", mNewVarCount); + + std::string name = mNewVarName + string(number); + IdList ids = mModel->getAllElementIdList(); + while (ids.contains(name)) + { + mNewVarCount++; + sprintf(number, "%u", mNewVarCount); + name = mNewVarName + string(number); + } + return name; +} + + +void +ExpressionAnalyser::addParametersAndRateRules(List* hiddenSpecies) +{ + for (unsigned int i = 0; i < mExpressions.size(); i++) + { + SubstitutionValues_t *exp = mExpressions.at(i); + if (mModel->getParameter(exp->z_value) == NULL) + { + // create expression for z + ASTNode* kx = new ASTNode(AST_MINUS); + ASTNode* k = new ASTNode(AST_NAME); + k->setName(exp->k_value.c_str()); + ASTNode* x = new ASTNode(AST_NAME); + x->setName(exp->x_value.c_str()); + kx->addChild(k); + kx->addChild(x); + + ASTNode* zNode = new ASTNode(AST_MINUS); + + + // add raterule defining dz/dt + ASTNode* dxdt = exp->dxdt_expression->deepCopy(); + RateRule* raterule = mModel->createRateRule(); + raterule->setVariable(exp->z_value); + ASTNode* math = new ASTNode(AST_TIMES); + ASTNode* minus1 = new ASTNode(AST_REAL); + minus1->setValue(-1.0); + + ASTNode* dydt = NULL; + ASTNode* plus = NULL; + ASTNode* y = NULL; + switch (exp->type) + { + case TYPE_K_MINUS_X: + case TYPE_K_PLUS_V_MINUS_X: + // dz/dt = -dx/dt + math->addChild(minus1); + math->addChild(dxdt); + + // z = k - x + (*zNode) = *kx; + + break; + case TYPE_K_MINUS_X_MINUS_Y: + case TYPE_K_PLUS_V_MINUS_X_MINUS_Y: + case TYPE_K_MINUS_X_PLUS_W_MINUS_Y: + // dz/dt = - (dx/dt + dy/dt) + dydt = exp->dydt_expression->deepCopy(); + plus = new ASTNode(AST_PLUS); + plus->addChild(dxdt); + plus->addChild(dydt); + math->addChild(minus1); + math->addChild(plus); + + // z = k-x-y + y = new ASTNode(AST_NAME); + y->setName(exp->y_value.c_str()); + zNode->addChild(kx); + zNode->addChild(y); + + break; + default: + break; + } + raterule->setMath(math); + + // introduce z + Parameter* zParam = mModel->createParameter(); + zParam->setId(exp->z_value); + zParam->setConstant(false); + zParam->setValue(SBMLTransforms::evaluateASTNode(zNode, mModel)); + hiddenSpecies->add(zParam); + + delete zNode; + delete math; //its children dxdt and minus1 deleted as part of this. + } + } +} + + +void +ExpressionAnalyser::replaceExpressionWithNewParameter(ASTNode* ode, SubstitutionValues_t* exp) +{ + if (exp->type == TYPE_K_MINUS_X || exp->type == TYPE_K_MINUS_X_MINUS_Y) + { + replaceExpressionInNodeWithVar(ode, exp->current, exp->z_value); + //cout << "ode in new param var: " << SBML_formulaToL3String(ode) << endl; + for (unsigned int i = 0; i < mExpressions.size(); i++) + { + SubstitutionValues_t *thisexp = mExpressions.at(i); + + if (thisexp->dxdt_expression != NULL) + { + replaceExpressionInNodeWithVar(thisexp->dxdt_expression, exp->current, exp->z_value); + } + if (thisexp->dydt_expression != NULL) + { + replaceExpressionInNodeWithVar(thisexp->dydt_expression, exp->current, exp->z_value); + } + } + + } + if (exp->type == TYPE_K_PLUS_V_MINUS_X || exp->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y) + { + ASTNode* replacement = new ASTNode(AST_PLUS); + ASTNode* z = new ASTNode(AST_NAME); + z->setName(exp->z_value.c_str()); + ASTNode *v = exp->v_expression->deepCopy(); + replacement->addChild(z); + replacement->addChild(v); + replaceExpressionInNodeWithNode(ode, exp->current, replacement); + //cout << "ode in new param node: " << SBML_formulaToL3String(ode) << endl; + for (unsigned int i = 0; i < mExpressions.size(); i++) + { + SubstitutionValues_t *thisexp = mExpressions.at(i); + + if (thisexp->dxdt_expression != NULL) + { + //cout << "dxdt_b4: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; + replaceExpressionInNodeWithNode(thisexp->dxdt_expression, exp->current, replacement); + //cout << "dxdt: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; + } + if (thisexp->dydt_expression != NULL) + { + //cout << "dydt_b4: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; + replaceExpressionInNodeWithNode(thisexp->dydt_expression, exp->current, replacement); + //cout << "dydt: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; + } + } + } + if (exp->type == TYPE_K_MINUS_X_PLUS_W_MINUS_Y) + { + ASTNode* replacement = new ASTNode(AST_PLUS); + ASTNode* z = new ASTNode(AST_NAME); + z->setName(exp->z_value.c_str()); + ASTNode *v = exp->w_expression->deepCopy(); + replacement->addChild(z); + replacement->addChild(v); + //cout << "ode in new param node: " << SBML_formulaToL3String(ode) << endl; + //cout << "current in new param node: " << SBML_formulaToL3String(exp->current) << endl; + //cout << "replace in new param node: " << SBML_formulaToL3String(replacement) << endl; + replaceExpressionInNodeWithNode(ode, exp->current, replacement); + //cout << "ode in new param node: " << SBML_formulaToL3String(ode) << endl; + for (unsigned int i = 0; i < mExpressions.size(); i++) + { + SubstitutionValues_t *thisexp = mExpressions.at(i); + + if (thisexp->dxdt_expression != NULL) + { + //cout << "dxdt_b4: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; + replaceExpressionInNodeWithNode(thisexp->dxdt_expression, exp->current, replacement); + //cout << "dxdt: " << SBML_formulaToL3String(thisexp->dxdt_expression) << endl; + } + if (thisexp->dydt_expression != NULL) + { + //cout << "dydt_b4: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; + replaceExpressionInNodeWithNode(thisexp->dydt_expression, exp->current, replacement); + //cout << "dydt: " << SBML_formulaToL3String(thisexp->dydt_expression) << endl; + } + } + } +} + + +bool +ExpressionAnalyser::addHiddenVariablesForKMinusX(List* hiddenSpecies) +{ + return true; +} + +/* +* check that the variable names in these two expressions match +*/ +bool +ExpressionAnalyser::matchesVariables(SubstitutionValues_t* exp, SubstitutionValues_t* exp1) +{ + if (exp->k_value == exp1->k_value && + exp->x_value == exp1->x_value && + exp->y_value == exp1->y_value) + { + return true; + } + return false; +} + + +/* +* Have we already created a parameter for this expression +* if so, return name +*/ +int +ExpressionAnalyser::parameterAlreadyCreated(SubstitutionValues_t* exp) +{ + int match = -1; + unsigned int i = 0; + while (!match && i < mExpressions.size()) + { + SubstitutionValues_t* thisexp = mExpressions.at(i); + if (thisexp->z_value.empty() || !matchesVariables(thisexp, exp)) + { + i++; + continue; + } + if (thisexp->type == exp->type) + { + match = i; + } + else if (thisexp->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y && exp->type == TYPE_K_MINUS_X_MINUS_Y) + { + match = i; + } + i++; + } + return match; +} + +/* + * Check whether for node is a name node representing species or a non constant parameter +*/ +bool ExpressionAnalyser::isVariableSpeciesOrParameter(ASTNode* node) +{ + if (!node->isName()) // some nodes, like * operators, don't seem to have a name in the first place + return false; + Species* species = mModel->getSpecies(node->getName()); + Parameter* parameter = mModel->getParameter(node->getName()); // some species in rate rules may be defined as variable parameters + bool isVariableSpeciesOrParameter = (species != NULL && !species->getConstant()); + bool isVariableParameter = (parameter!=NULL && !parameter->getConstant()); + return isVariableSpeciesOrParameter || isVariableParameter; +} + +/* +* Check whether for node is a name node representing a constant parameter or a numerical node +*/ +bool ExpressionAnalyser::isNumericalConstantOrConstantParameter(ASTNode* node, bool& isNumber) +{ + bool isConstantParameter = false; + isNumber = false; + + if (node->isName()) // some nodes, like * operators, don't seem to have a name in the first place + { + Parameter* parameter = mModel->getParameter(node->getName()); + isConstantParameter = (parameter != NULL) && (parameter->getConstant()); + } + bool isNumericalConstant = node->isNumber() || node->isConstant(); + + if (isConstantParameter) + return true; + else if (isNumericalConstant) + { + isNumber = true; + return true; + } + else + return false; +} + +/* +* Reorder any instance of - x + y with y - x in the set of ODEs. +* Fages Algorithm 3.1 Step 1 +*/ +void ExpressionAnalyser::reorderMinusXPlusYIteratively() +{ + for (unsigned int i = 0; i < mExpressions.size(); i++) + { + SubstitutionValues_t* exp = mExpressions.at(i); + if (exp->type != TYPE_MINUS_X_PLUS_Y) + continue; + ASTNode* ode = (mODEs.at(exp->odeIndex)).second; + ASTNode* replacement = new ASTNode(AST_MINUS); + ASTNode* y = new ASTNode(AST_NAME); + y->setName((exp->y_value).c_str()); + ASTNode* x = new ASTNode(AST_NAME); + x->setName((exp->x_value).c_str()); + replacement->addChild(y); + replacement->addChild(x); + replaceExpressionInNodeWithNode(ode, exp->current, replacement); + } +} + +std::pair ExpressionAnalyser::getParentNode(const ASTNode* child, const ASTNode* root) +{ + //cout << "root " << SBML_formulaToL3String(root) << endl; + //cout << "child " << SBML_formulaToL3String(child) << endl; + for (unsigned int i = 0; i < root->getNumChildren(); i++) + { + if (root->getChild(i)->exactlyEqual(*(child))) + { + return std::pair(const_cast(root), i); + } + } + for (unsigned int i = 0; i < root->getNumChildren(); i++) + { + std::pair parent = getParentNode(child, root->getChild(i)); + if (parent.first != NULL) + { + return parent; + } + } + return std::pair(NULL, (int)(NAN)); +} + +/** @endcond */ + +LIBSBML_CPP_NAMESPACE_END + +#endif /* __cplusplus */ + + diff --git a/src/sbml/conversion/ExpressionAnalyser.h b/src/sbml/conversion/ExpressionAnalyser.h index 3b2f30f266..28bd531370 100644 --- a/src/sbml/conversion/ExpressionAnalyser.h +++ b/src/sbml/conversion/ExpressionAnalyser.h @@ -79,6 +79,7 @@ typedef enum |*/ struct SubstitutionValues_t { std::string k_value; + double k_real_value; std::string x_value; std::string y_value; ASTNode * dxdt_expression; @@ -223,7 +224,7 @@ class LIBSBML_EXTERN ExpressionAnalyser * @param node the node to check * @return true if the node is a constant number/parameter */ - bool isNumericalConstantOrConstantParameter(ASTNode* node); + bool isNumericalConstantOrConstantParameter(ASTNode* node, bool& isNumber); /* * Have we already created a parameter for this expression diff --git a/src/sbml/conversion/SBMLRateRuleConverter.cpp b/src/sbml/conversion/SBMLRateRuleConverter.cpp index ed94d4089f..33c43b1232 100644 --- a/src/sbml/conversion/SBMLRateRuleConverter.cpp +++ b/src/sbml/conversion/SBMLRateRuleConverter.cpp @@ -196,6 +196,74 @@ SBMLRateRuleConverter::matchesProperties(const ConversionProperties &props) cons } +void +SBMLRateRuleConverter::catchAnomalies() +{ + // if a rate rule is actually assigned to a variable that is itself + // the subject of an assignment rule we want to replace that math + // before we do any analysis + AssignmentRule* assignment; + Model* model = mDocument->getModel(); + for (unsigned int i = 0; i < mDocument->getModel()->getNumRules(); ++i) + { + Rule* rule = model->getRule(i); + if (rule->getTypeCode() != SBML_RATE_RULE) + { + continue; + } + ASTNode* math = const_cast(rule->getMath()); + if (math->getType() == AST_NAME) + { + assignment = model->getAssignmentRuleByVariable(math->getName()); + if (assignment == NULL) + { + continue; + } + + const ASTNode* assignment_math = assignment->getMath(); + math->setType(assignment_math->getType()); + for (unsigned int n = 0; n < assignment_math->getNumChildren(); ++n) + { + math->addChild(assignment_math->getChild(n)->deepCopy()); + } + } + + } + + + //if (this->getType() == AST_NAME) + //{ + // // here the node is a variable + // std::string name = this->getName(); + // const SBase* parent = this->getParentSBMLObject(); + // const Model* model = static_cast(parent->getAncestorOfType(SBML_MODEL)); + // const AssignmentRule* rule = model->getAssignmentRuleByVariable(name); + + // // and the variable is assigned by rule + // // we want to replace this node by the rule + // if (rule != NULL) + // { + // const ASTNode* math = rule->getMath(); + // // if this is a name we just want to replace + // ASTNodeType_t type = math->getType(); + // if (type == AST_NAME) + // + // { + // this->setName(math->getName()); + // } + // else + // { + // this->setType(math->getType()); + // for (unsigned int i = 0; i < math->getNumChildren(); ++i) + // { + // this->addChild(math->getChild(i)->deepCopy()); + // } + // } + + // } + //} +} + int SBMLRateRuleConverter::convert() { @@ -205,6 +273,7 @@ SBMLRateRuleConverter::convert() { return returnValue; } + catchAnomalies(); // Fages algo 3.6 Steps 1-2 populateODEinfo(); @@ -376,7 +445,12 @@ SBMLRateRuleConverter::determineCoefficient(ASTNode* ode, unsigned int termN, do // first child should be a number // take it out of the term // it will be used as a coefficient - if (term1->getType() == AST_TIMES && term1->getNumChildren() > 0) + if (term1->getType() == AST_REAL) + { + coeff = term1->getValue(); + found = true; + } + else if (term1->getType() == AST_TIMES && term1->getNumChildren() > 0) { if (term1->getChild(0)->isNumber()) { @@ -567,7 +641,7 @@ SBMLRateRuleConverter::isPositive(const ASTNode* node, bool& posDeriv) } else { - posDeriv = false; + if (mDerivSign == NEGATIVE_DERIVATIVE) posDeriv = false; } signDetermined = true; } @@ -664,8 +738,13 @@ SBMLRateRuleConverter::addToTerms(ASTNode* node) } else if (term->isNumber()) { - delete term; - return; + // however if the term is a constant we don't want to delete + if (term->getNumChildren() != 0) + { + delete term; + return; + } + term->setValue(1.0); } if (mTerms.size() == 0) diff --git a/src/sbml/conversion/SBMLRateRuleConverter.h b/src/sbml/conversion/SBMLRateRuleConverter.h index 2253976fa8..f2904bbf31 100644 --- a/src/sbml/conversion/SBMLRateRuleConverter.h +++ b/src/sbml/conversion/SBMLRateRuleConverter.h @@ -62,6 +62,7 @@ #include #include #include +#include typedef enum { @@ -74,6 +75,7 @@ typedef enum LIBSBML_CPP_NAMESPACE_BEGIN + typedef std::vector< std::pair< std::string, ASTNode*> > pairODEs; class LIBSBML_EXTERN SBMLRateRuleConverter : public SBMLConverter @@ -244,6 +246,8 @@ class LIBSBML_EXTERN SBMLRateRuleConverter : public SBMLConverter void createReactions(); void removeRules(); + void catchAnomalies(); + // member variables populated during analysis pairODEs mODEs; diff --git a/src/sbml/conversion/test/TestRunner.c b/src/sbml/conversion/test/TestRunner.c index 5c6f8ee80f..7e6a99a19f 100644 --- a/src/sbml/conversion/test/TestRunner.c +++ b/src/sbml/conversion/test/TestRunner.c @@ -118,9 +118,9 @@ main (void) int num_failed; setTestDataDirectory(); -// SRunner *runner = srunner_create(create_suite_TestUnitsConverterL2()); + SRunner *runner = srunner_create(create_suite_TestSBMLRateRuleConverter()); - SRunner *runner = srunner_create( create_suite_TestConversionOption() ); + /*SRunner *runner = srunner_create( create_suite_TestConversionOption() ); srunner_add_suite( runner, create_suite_TestSBMLRuleConverter () ); srunner_add_suite( runner, create_suite_TestConversionProperties () ); srunner_add_suite( runner, create_suite_TestSBMLConverterRegistry () ); @@ -131,7 +131,7 @@ main (void) srunner_add_suite( runner, create_suite_TestStripPackageConverter () ); srunner_add_suite( runner, create_suite_TestLevelVersionConverter () ); srunner_add_suite( runner, create_suite_TestRateOfConverter () ); - srunner_add_suite( runner, create_suite_TestSBMLRateRuleConverter () ); + srunner_add_suite( runner, create_suite_TestSBMLRateRuleConverter () );*/ /* srunner_set_fork_status(runner, CK_NOFORK); */ diff --git a/src/sbml/math/ASTNode.cpp b/src/sbml/math/ASTNode.cpp index 9bf7ec7e8f..497ec9a228 100644 --- a/src/sbml/math/ASTNode.cpp +++ b/src/sbml/math/ASTNode.cpp @@ -4086,6 +4086,7 @@ LIBSBML_EXTERN void ASTNode::refactor() { + // we need to look at whether we have got a single variable that is assigned by an assignment rule refactorNumbers(); encompassUnaryMinus(); createNonBinaryTree(); diff --git a/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-04-10102.xml b/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-06.xml similarity index 100% rename from src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-04-10102.xml rename to src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-06.xml diff --git a/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-06-10102.xml b/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-07.xml similarity index 100% rename from src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-06-10102.xml rename to src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-07.xml diff --git a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-27-20517.xml b/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-08.xml similarity index 80% rename from src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-27-20517.xml rename to src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-08.xml index 69370fd243..f60c1bb550 100644 --- a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-27-20517.xml +++ b/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-08.xml @@ -1,6 +1,6 @@ +xmlns:appear="http://www.sbml.org/sbml/level3/version1/appear" appear:required="true"> diff --git a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-29-20517.xml b/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-09.xml similarity index 86% rename from src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-29-20517.xml rename to src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-09.xml index ff0a4a7137..cce61f113a 100644 --- a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-29-20517.xml +++ b/src/sbml/validator/test/test-data/libsbml-constraints/99107-fail-01-09.xml @@ -1,6 +1,6 @@ diff --git a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-28-20517.xml b/src/sbml/validator/test/test-data/libsbml-constraints/99108-fail-01-06.xml similarity index 86% rename from src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-28-20517.xml rename to src/sbml/validator/test/test-data/libsbml-constraints/99108-fail-01-06.xml index bdc090f245..c53076cd87 100644 --- a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-28-20517.xml +++ b/src/sbml/validator/test/test-data/libsbml-constraints/99108-fail-01-06.xml @@ -1,6 +1,6 @@ diff --git a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-30-20517.xml b/src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-07.xml similarity index 87% rename from src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-30-20517.xml rename to src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-07.xml index 1877ba8743..43726846a6 100644 --- a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-30-20517.xml +++ b/src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-07.xml @@ -1,6 +1,6 @@ +xmlns:appear="http://www.sbother"> diff --git a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-32-20517.xml b/src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-08.xml similarity index 88% rename from src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-32-20517.xml rename to src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-08.xml index 2da0266d66..d8544069c2 100644 --- a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-32-20517.xml +++ b/src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-08.xml @@ -1,6 +1,6 @@ diff --git a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-31-20517.xml b/src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-09.xml similarity index 88% rename from src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-31-20517.xml rename to src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-09.xml index 7815afde46..5c83a7d1de 100644 --- a/src/sbml/validator/test/test-data/sbml-xml-constraints/10102-fail-01-31-20517.xml +++ b/src/sbml/validator/test/test-data/sbml-general-consistency-constraints/20517-fail-01-09.xml @@ -1,6 +1,6 @@