C++ KineticLaw类代码示例

/* convert from L1 to L3 */
void 
Model::convertL3ToL1 ()
{
  //
  // Level 3 allows a model to be specified without a Compartment.  However
  // this is not valid in Level 1.  Thus if a L3 model has no Compartment
  // one must be included 
  //
  if (getNumCompartments() == 0)
  {
    createCompartment()->setId(ASSIGNED_COMPARTMENT);

  }
  dealWithModelUnits();
  
  dealWithAssigningL1Stoichiometry(*this, false);
  for (unsigned int i = 0; i < getNumReactions(); i++)
  {
    Reaction *r = getReaction(i);
    if (r->isSetKineticLaw())
    {
      KineticLaw *kl = r->getKineticLaw();
      for (unsigned int j = 0; j < kl->getNumLocalParameters(); j++)
      {
        Parameter *lp = new Parameter(getLevel(), getVersion());
        (*lp) = *(kl->getLocalParameter(j));
        kl->addParameter(lp);
      }
    }
  }
}

END_TEST


START_TEST ( test_Reaction_parent_NULL )
{
    SBMLDocument *d = new SBMLDocument();
    Model *m = d->createModel();
    Reaction *c = m->createReaction();
    SpeciesReference *sr = c->createReactant();
    KineticLaw *kl = c->createKineticLaw();

    fail_unless(c->getAncestorOfType(SBML_MODEL) == m);
    fail_unless (c->getSBMLDocument() == d);
    fail_unless(sr->getAncestorOfType(SBML_REACTION) == c);
    fail_unless(kl->getAncestorOfType(SBML_REACTION) == c);

    Reaction *c1 = c->clone();
    delete d;

    fail_unless(c1->getAncestorOfType(SBML_MODEL) == NULL);
    fail_unless(c1->getParentSBMLObject() == NULL);
    fail_unless (c1->getSBMLDocument() == NULL);

    SpeciesReference *sr1 = c1->getReactant(0);
    fail_unless(sr1->getAncestorOfType(SBML_MODEL) == NULL);
    fail_unless(sr1->getAncestorOfType(SBML_REACTION) == c1);
    fail_unless (sr1->getSBMLDocument() == NULL);

    fail_unless(c1->getKineticLaw()->getAncestorOfType(SBML_MODEL) == NULL);
    fail_unless(c1->getKineticLaw()->getAncestorOfType(SBML_REACTION) == c1);
    fail_unless (c1->getKineticLaw()->getSBMLDocument() == NULL);


    delete c1;
}

void checkReactions(Model* model, set<string>& components, set<string>& tests,  const map<string, vector<double> >& results, int type)
{
  if (model->getNumReactions() > 0) {
    components.insert("Reaction");
    for (unsigned int r=0; r<model->getNumReactions(); r++) {
      Reaction* rxn = model->getReaction(r);
      if (rxn->isSetFast() && rxn->getFast()) {
        tests.insert("FastReaction");
      }
      if (rxn->isSetReversible() && rxn->getReversible()) {
        if (type!=1) {
          tests.insert("ReversibleReaction [?]");
        }
      }
      ListOfSpeciesReferences* reactants = rxn->getListOfReactants();
      checkSpeciesRefs(model, reactants, components, tests, results);
      ListOfSpeciesReferences* products = rxn->getListOfProducts();
      checkSpeciesRefs(model, products, components, tests, results);
      if (rxn->isSetKineticLaw()) {
        KineticLaw* kl = rxn->getKineticLaw();
        if (kl->getNumParameters() > 0) {
          tests.insert("LocalParameters");
        }
        if (kl->isSetMath() == false) {
          tests.insert("NoMathML");
        }
      }
    }
  }
}

END_TEST

START_TEST ( test_KineticLaw )
{
  KineticLaw* kl = new KineticLaw(2, 4);
  
  fail_unless (kl->hasRequiredAttributes());

  delete kl;
}

/* adds species referred to in a KineticLaw to the ListOfModifiers
 * this will only be applicable when up converting an L1 model
 */
void 
Model::addModifiers ()
{
  //
  // Level 2/3 has a listOfModifiers associated with a Reaction
  // which are not listed in a L1 Model.
  // For each symbol in the Reaction's KineticLaw,
  // that symbol is a modifier iff:
  //
  //   1. It is defined as a Species in the Model
  //   2. It is not a Reactant or Product in this Reaction.
  //
  // Thus modifiers must be added where appropriate.
  //
  const char *id;

  unsigned int size;
  unsigned int n, l;

  const ASTNode *node;
  List          *names;
  KineticLaw* kl;

  for (n = 0; n < getNumReactions(); n++)
  {
    kl = getReaction(n)->getKineticLaw();

    if (kl == NULL || kl->isSetMath() == false) continue;
   
    node  = kl->getMath();
    names = node->getListOfNodes((ASTNodePredicate) ASTNode_isName);
    size  = names->getSize();

    for (l = 0; l < size; l++)
    {
      node = (ASTNode *) names->get(l);
      id   = node->getName();

      // 1. It is an AST_NAME (not AST_NAME_TIME), and
      if (node->getType() != AST_NAME) continue;

      // 2. It refers to a Species in this Model, and
      if (id == NULL || getSpecies(id) == NULL) continue;

      // 3. It is not a Reactant, Product, or (already) a Modifier
      if (getReaction(n)->getReactant(id) != NULL) continue;
      if (getReaction(n)->getProduct (id) != NULL) continue;
      if (getReaction(n)->getModifier(id) != NULL) continue;

      getReaction(n)->createModifier()->setSpecies(id);
    }

    delete names;
  }
}

END_TEST

START_TEST ( test_KineticLaw )
{
  KineticLaw* kl = new KineticLaw(2, 4);
   
  fail_unless (!(kl->hasRequiredElements()));

  kl->setMath(SBML_parseFormula("kl"));

  fail_unless (kl->hasRequiredElements());

  delete kl;
}

void 
convertPow(SBMLDocument* doc, 
           bool shouldChangePow, 
           bool inlineCompartmentSizes)
{

  Model* model = SBMLDocument_getModel(doc);
  if (model == NULL)
  {
    return;
  }

  std::map<string, double> compartmentValueMap;
  if (inlineCompartmentSizes)
  {
    for(unsigned int i = 0; i < model->getNumCompartments(); ++i)
    {
      Compartment* c = model->getCompartment(i);
      compartmentValueMap[c->getId()] = c->getSize();
    }
  }

  for(unsigned int i = 0; i < model->getNumReactions(); ++i)
  {
    Reaction* r  = model->getReaction(i);
    KineticLaw* kl = r->getKineticLaw();
    const char* strKineticFormula;
    if (kl == NULL)
    {
      strKineticFormula = "";
    }
    else
    {
      strKineticFormula = KineticLaw_getFormula(kl);
      if (strKineticFormula == NULL)
      {
        continue;
      }
    }

    ASTNode* ast_Node	= SBML_parseFormula (strKineticFormula);
    if (ast_Node != NULL)
    {
      changePow(ast_Node, compartmentValueMap, shouldChangePow);
      kl->setMath (ast_Node);
    }
    delete ast_Node;
  }
}

void createKineticLawForReaction(Reaction* reaction)
{
  if (reaction == NULL)
    return;
  reaction->unsetKineticLaw();
  KineticLaw *law = reaction->getKineticLaw();
  if (law == NULL)
  {
    law = reaction->createKineticLaw();
    LocalParameter* fluxValue = law->createLocalParameter();
    fluxValue->initDefaults();
    fluxValue->setId("FLUX_VALUE");
    fluxValue->setValue(0);
    fluxValue->setUnits("dimensionless");
    ASTNode* astn = SBML_parseFormula("FLUX_VALUE");
    law->setMath(astn);
    delete astn;
  }

  LocalParameter* LB = law->getLocalParameter("LOWER_BOUND");
  if (LB == NULL)
  {
    LB = law->createLocalParameter();
    LB->initDefaults();
    LB->setId("LOWER_BOUND");
    LB->setUnits("dimensionless");
    LB->setValue(-std::numeric_limits<double>::infinity());
  }

  LocalParameter* UB = law->getLocalParameter("UPPER_BOUND");
  if (UB == NULL)
  {
    UB = law->createLocalParameter();
    UB->initDefaults();
    UB->setId("UPPER_BOUND");
    UB->setUnits("dimensionless");
    LB->setValue(std::numeric_limits<double>::infinity());
  }

  LocalParameter* param = law->getLocalParameter("OBJECTIVE_COEFFICIENT");
  if (param == NULL)
  {
    param = law->createLocalParameter();
    param->initDefaults();
    param->setId("OBJECTIVE_COEFFICIENT");
    param->setUnits("dimensionless");
    param->setValue(0);
  }

}

void updateKineticLawFromBound(Reaction* reaction, FluxBound* current)
{
  if (reaction == NULL || current == NULL)
    return;
  const string operation = current -> getOperation();

  KineticLaw *law = reaction->getKineticLaw();
  LocalParameter* LB = law->getLocalParameter("LOWER_BOUND");
  LocalParameter* UB = law->getLocalParameter("UPPER_BOUND");

  if (operation == "less" || operation == "lessEqual" || operation == "equal")
  {
    UB->setValue(current->getValue());
  }
  if (operation == "greater" || operation == "greaterEqual" || operation == "equal")
  {
    LB->setValue(current->getValue());
  }
  
}

END_TEST


START_TEST (test_WriteL3SBML_KineticLaw_ListOfParameters)
{
  const char* expected = 
    "<kineticLaw>\n"
    "  <listOfLocalParameters>\n"
    "    <localParameter id=\"n\" value=\"1.2\"/>\n"
    "  </listOfLocalParameters>\n"
    "</kineticLaw>";

  KineticLaw *kl = D->createModel()->createReaction()->createKineticLaw();

  LocalParameter *p = kl->createLocalParameter();
  p->setId("n");
  p->setValue(1.2);

  char* sbml = kl->toSBML();
  fail_unless( equals(expected, sbml) );
  safe_free(sbml);
}

void setObjectiveCoefficient(FbcModelPlugin* plugin, Model* model)
{
  if (plugin == NULL || model == NULL)
    return;

  Objective* obj = plugin->getActiveObjective();
  if (obj == NULL)
    return;

  for (unsigned int i = 0; i < obj->getNumFluxObjectives(); ++i)
  {
    FluxObjective* fluxObj = obj->getFluxObjective(i);
    if (fluxObj == NULL)
      continue;
    Reaction* reaction = model->getReaction(fluxObj->getReaction());
    if (reaction == NULL)
      continue;
    KineticLaw* law = reaction->getKineticLaw();
    if (law == NULL)
      continue;
    LocalParameter* param = law->getLocalParameter("OBJECTIVE_COEFFICIENT");
    param->setValue(fluxObj->getCoefficient());
  }
}

/* convert from L1 to L3 */
void 
Model::convertL3ToL2 (bool strict)
{
  dealWithModelUnits();

  dealWithStoichiometry();

  dealWithEvents(strict);

  for (unsigned int i = 0; i < getNumReactions(); i++)
  {
    Reaction *r = getReaction(i);
    if (r->isSetKineticLaw())
    {
      KineticLaw *kl = r->getKineticLaw();
      for (unsigned int j = 0; j < kl->getNumLocalParameters(); j++)
      {
        Parameter *lp = new Parameter(getLevel(), getVersion());
        (*lp) = *(kl->getLocalParameter(j));
        kl->addParameter(lp);
      }
    }
  }
}

int Submodel::convertTimeAndExtentWith(const ASTNode* tcf, const ASTNode* xcf, const ASTNode* klmod)
{
  if (tcf==NULL && xcf==NULL) return LIBSBML_OPERATION_SUCCESS;
  Model* model = getInstantiation();
  if (model==NULL) {
    //getInstantiation sets its own error messages.
    return LIBSBML_OPERATION_FAILED;
  }
  ASTNode tcftimes(AST_TIMES);
  ASTNode tcfdiv(AST_DIVIDE);
  if (tcf != NULL) {
    tcftimes.addChild(tcf->deepCopy());
    tcfdiv.addChild(tcf->deepCopy());
  }
  ASTNode rxndivide(AST_DIVIDE);
  if (klmod != NULL) {
    ASTNode rxnref(AST_NAME);
    rxndivide.addChild(rxnref.deepCopy());
    rxndivide.addChild(klmod->deepCopy());
  }
  List* allElements = model->getAllElements();
  for (ListIterator iter = allElements->begin(); iter != allElements->end(); ++iter)
  {
    SBase* element = static_cast<SBase*>(*iter);
    assert(element != NULL);
    ASTNode* ast1 = NULL;
    ASTNode* ast2 = NULL;
    Constraint* constraint = NULL;
    Delay* delay = NULL;
    EventAssignment* ea = NULL;
    InitialAssignment* ia = NULL;
    KineticLaw* kl = NULL;
    Priority* priority = NULL;
    RateRule* rrule = NULL;
    Rule* rule = NULL;
    Submodel* submodel = NULL;
    Trigger* trigger = NULL;
    string cf = "";
    //Reaction math will be converted below, in the bits with the kinetic law.  But because of that, we need to handle references *to* the reaction:  even if it has no kinetic law, the units have changed, and this needs to be reflected by the flattening routine.
    if (rxndivide.getNumChildren() != 0 && element->getTypeCode()==SBML_REACTION && element->isSetId()) {
      rxndivide.getChild(0)->setName(element->getId().c_str());
      for (ListIterator iter = allElements->begin(); iter != allElements->end(); ++iter)
      {
        SBase* subelement = static_cast<SBase*>(*iter);
        subelement->replaceSIDWithFunction(element->getId(), &rxndivide);
      }
    }

    //Submodels need their timeConversionFactor and extentConversionFactor attributes converted.  We're moving top-down, so all we need to do here is fix the conversion factor attributes themselves, pointing them to new parameters if need be.
    if ((tcf !=NULL || xcf != NULL) && element->getTypeCode()==SBML_COMP_SUBMODEL) {
      submodel = static_cast<Submodel*>(element);
      if (tcf != NULL) {
        if (submodel->isSetTimeConversionFactor()) {
          createNewConversionFactor(cf, tcf, submodel->getTimeConversionFactor(), model);
          submodel->setTimeConversionFactor(cf);
        }
        else {
          submodel->setTimeConversionFactor(tcf->getName());
        }
      }
      if (xcf != NULL) {
        if (submodel->isSetExtentConversionFactor()) {
          createNewConversionFactor(cf, xcf, submodel->getExtentConversionFactor(), model);
          submodel->setExtentConversionFactor(cf);
        }
        else {
          submodel->setExtentConversionFactor(xcf->getName());
        }
      }
    }
    if (tcf==NULL) {
      if (klmod !=NULL && element->getTypeCode()==SBML_KINETIC_LAW) {
        kl = static_cast<KineticLaw*>(element);
        if (kl->isSetMath()) {
          ast1 = new ASTNode(AST_TIMES);
          ast1->addChild(klmod->deepCopy());
          ast1->addChild(kl->getMath()->deepCopy());
          kl->setMath(ast1);
          delete ast1;
        }
      }
    }
    else {
      // All math 'time' and 'delay' csymbols must still be converted.
      // Also, several constructs are modified directly.
      switch(element->getTypeCode()) {
        //This would be a WHOLE LOT SIMPLER if there was a 'hasMath' class in libsbml.  But even so, it would have to
        // handle the kinetic laws, rate rules, and delays separately.
      case SBML_KINETIC_LAW:
        //Kinetic laws are multiplied by 'klmod'.
        kl = static_cast<KineticLaw*>(element);
        ast1 = kl->getMath()->deepCopy();
        convertCSymbols(ast1, &tcfdiv, &tcftimes);
        if (klmod !=NULL) {
          kl = static_cast<KineticLaw*>(element);
          if (kl->isSetMath()) {
            ast2 = new ASTNode(AST_TIMES);
            ast2->addChild(klmod->deepCopy());
            ast2->addChild(ast1);
            kl->setMath(ast2);
//.........这里部分代码省略.........

int
main (int argc, char* argv[])
{
  if (argc != 2)
  {
    cout << endl << "Usage: printNotes filename" << endl << endl;
    return 1;
  }

  unsigned int i,j;
  const char* filename   = argv[1];
  SBMLDocument* document;
  SBMLReader reader;

  document = reader.readSBML(filename);

  unsigned int errors = document->getNumErrors();

  cout << endl;
  cout << "filename: " << filename << endl;
  cout << endl;

  if(errors > 0)
  {
    document->printErrors(cerr);
    delete document;

    return errors;
  }

  /* Model */

  Model* m = document->getModel();
  printNotes(m);

  for(i=0; i < m->getNumReactions(); i++)
  {
    Reaction* re = m->getReaction(i);
    printNotes(re);

    /* SpeciesReference (Reacatant) */

    for(j=0; j < re->getNumReactants(); j++)
    {
      SpeciesReference* rt = re->getReactant(j);
      if (rt->isSetNotes()) cout << "   ";
      printNotes(rt, (rt->isSetSpecies() ? rt->getSpecies() : std::string("")) );
    }

    /* SpeciesReference (Product) */

    for(j=0; j < re->getNumProducts(); j++)
    {
      SpeciesReference* rt = re->getProduct(j);
      if (rt->isSetNotes()) cout << "   ";
      printNotes(rt, (rt->isSetSpecies() ? rt->getSpecies() : std::string("")) );
    }

    /* ModifierSpeciesReference (Modifier) */

    for(j=0; j < re->getNumModifiers(); j++)
    {
      ModifierSpeciesReference* md = re->getModifier(j);
      if (md->isSetNotes()) cout << "   ";
      printNotes(md, (md->isSetSpecies() ? md->getSpecies() : std::string("")) );
    }

    /* Kineticlaw */

    if(re->isSetKineticLaw())
    {
      KineticLaw* kl = re->getKineticLaw();
      if (kl->isSetNotes()) cout << "   ";
      printNotes(kl);

      /* Parameter */

      for(j=0; j < kl->getNumParameters(); j++)
      {
        Parameter* pa = kl->getParameter(j);
        if (pa->isSetNotes()) cout << "      ";
        printNotes(pa);
      }
    }

  }

  /* Species */

  for(i=0; i < m->getNumSpecies(); i++)
  {
    Species* sp = m->getSpecies(i);
    printNotes(sp);
  }

  /* Compartment */

  for(i=0; i < m->getNumCompartments(); i++)
  {
    Compartment* sp = m->getCompartment(i);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                      assert(pSBase != NULL);

                      // since L2V2 species references can have ids
                      if (pSBase->isSetId())
                        {
                          id = pSBase->getId();

                          if (ids.find(id) == ids.end())
                            {
                              ids.insert(std::pair<const std::string, const SBase*>(id, pSBase));
                            }
                          else
                            {
                              CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 68, id.c_str());
                            }
                        }

                      if (pSBase->isSetMetaId())
                        {
                          id = pSBase->getMetaId();

                          if (metaIds.find(id) == metaIds.end())
                            {
                              metaIds.insert(std::pair<const std::string, const SBase*>(id, pSBase));
                            }
                          else
                            {
                              CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 67, id.c_str());
                            }
                        }
                    }
                }

              KineticLaw* pKLaw = pReaction->getKineticLaw();

              if (pKLaw != NULL)
                {
                  if (pKLaw->isSetMetaId())
                    {
                      id = pKLaw->getMetaId();

                      if (metaIds.find(id) == metaIds.end())
                        {
                          metaIds.insert(std::pair<const std::string, const SBase*>(id, pKLaw));
                        }
                      else
                        {
                          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 67, id.c_str());
                        }
                    }

                  pSBase = pKLaw->getListOfParameters();

                  if (pSBase != NULL)
                    {
                      if (pSBase->isSetMetaId())
                        {
                          id = pSBase->getMetaId();

                          if (metaIds.find(id) == metaIds.end())
                            {
                              metaIds.insert(std::pair<const std::string, const SBase*>(id, pSBase));
                            }
                          else
                            {
                              CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 67, id.c_str());