C++ Species::getBoundaryCondition方法代码示例

bool variesIn(string id, Model* model,  const map<string, vector<double> >& results)
{
  if (variesIn(id, results)) return true;
  Rule* rule = model->getRule(id);
  if (rule != NULL) {
    if (rule->getTypeCode() == SBML_RATE_RULE) return true;
    return variesIn(rule->getMath(), model, results);
  }
  if (!getConstant(id, model, results)) {
    //The variable might be set by an algebraic rule
    for (unsigned long r=0; r<model->getNumRules(); r++) {
      if (model->getRule(r)->getTypeCode() == SBML_ALGEBRAIC_RULE) {
        if (appearsIn(id, model->getRule(r)->getMath()) &&
            variesBesides(id, model->getRule(r)->getMath(), model, results)) return true;
      }
    }
  }
  for (unsigned long e=0; e<model->getNumEvents(); e++) {
    Event* event = model->getEvent(e);
    for (unsigned long ea=0; ea<event->getNumEventAssignments(); ea++) {
      if (event->getEventAssignment(ea)->getVariable() == id) return true;
    }
  }
  Species* species = model->getSpecies(id);
  if (species==NULL) return false;
  if (species->isSetBoundaryCondition() && species->getBoundaryCondition()==true) return false;
  for (unsigned long r=0; r<model->getNumReactions(); r++) {
    Reaction* rxn = model->getReaction(r);
    if (variesIn(id, rxn->getListOfReactants(), results)) return true;
    if (variesIn(id, rxn->getListOfProducts(), results)) return true;
  }
  return false;
}

void checkSpecies(Model* model, set<string>& components, set<string>& tests,  const map<string, vector<double> >& results, int type)
{
  //Must call this after 'checkCompartments' because we look in 'tests' for 'NonUnityCompartment'.
  if (model->getNumSpecies() > 0) {
    components.insert("Species");
    if (type==0) {
      tests.insert("Amount||Concentration");
    }
    else if (type==2) {
      tests.insert("Amount");
    }
	set<string> compartments;
    for (unsigned int s=0; s<model->getNumSpecies(); s++) {
      Species* species = model->getSpecies(s);
      if (species->isSetBoundaryCondition() && species->getBoundaryCondition()) {
        tests.insert("BoundaryCondition");
      }
      if (species->getConstant()) {
        tests.insert("ConstantSpecies");
      }
      if (species->isSetConversionFactor()) {
        tests.insert("ConversionFactors");
      }
      if (species->isSetHasOnlySubstanceUnits() && species->getHasOnlySubstanceUnits()) {
        tests.insert("HasOnlySubstanceUnits");
      }
      if (!species->isSetInitialAmount() && !species->isSetInitialConcentration()) {
        tests.insert("InitialValueReassigned");
      }
      else if (species->isSetId() && initialOverriddenIn(species->getId(), model, results, tests)) {
        tests.insert("InitialValueReassigned");
      }
      if (species->isSetCompartment()) {
        compartments.insert(species->getCompartment());
      }
    }
    if (tests.find("MultiCompartment") != tests.end() && compartments.size()==1 && model->getNumSpecies() > 1) {
      cerr << "Error:  multiple compartments discovered, but all species are in a single compartment." << endl;
      tests.insert("ERRORMultiCompartment");
    }
  }
}

void
Model::assignRequiredValues()
{
  // when converting to L3 some attributes which have default values in L1/L2
  // but are required in L3 are not present or set
  unsigned int i, n;

  if (getNumUnitDefinitions() > 0)
  {
    for (i = 0; i < getNumUnitDefinitions(); i++)
    {
      for (n = 0; n < getUnitDefinition(i)->getNumUnits(); n++)
      {
        Unit *u = getUnitDefinition(i)->getUnit(n);
        if (!u->isSetExponent())
          u->setExponent(1.0);
        if (!u->isSetScale())
          u->setScale(0);
        if (!u->isSetMultiplier())
          u->setMultiplier(1.0);
      }
    }
  }
  
  if (getNumCompartments() > 0)
  {
    for (i = 0; i < getNumCompartments(); i++)
    {
      Compartment *c = getCompartment(i);
      c->setConstant(c->getConstant());
    }
  }
  if (getNumSpecies() > 0)
  {
    for (i = 0; i < getNumSpecies(); i++)
    {
      Species * s = getSpecies(i);
      s->setBoundaryCondition(s->getBoundaryCondition());
      s->setHasOnlySubstanceUnits(s->getHasOnlySubstanceUnits());
      s->setConstant(s->getConstant());
    }
  }
  if (getNumParameters() > 0)
  {
    for (i = 0; i < getNumParameters(); i++)
    {
      Parameter * p = getParameter(i);
      p->setConstant(p->getConstant());
    }
  }
  if (getNumReactions() > 0)
  {
    for (i = 0; i < getNumReactions(); i++)
    {
      Reaction * r = getReaction(i);
      r->setFast(r->getFast());
      r->setReversible(r->getReversible());
      if (r->getNumReactants() > 0)
      {
        for (n = 0; n < r->getNumReactants(); n++)
        {
          SpeciesReference *sr = r->getReactant(n);
          if (sr->isSetStoichiometryMath())
          {
            sr->setConstant(false);
          }
          else
          {
            sr->setConstant(true);
          }
        }
      }
      if (r->getNumProducts() > 0)
      {
        for (n = 0; n < r->getNumProducts(); n++)
        {
          SpeciesReference *sr = r->getProduct(n);
          if (sr->isSetStoichiometryMath())
          {
            sr->setConstant(false);
          }
          else
          {
            sr->setConstant(true);
          }
        }
      }
    }
  }
  if (getNumEvents() > 0)
  {
    for (i = 0; i < getNumEvents(); i++)
    {
      Event * e = getEvent(i);
      e->setUseValuesFromTriggerTime(e->getUseValuesFromTriggerTime());

      if (e->isSetTrigger())
      {
        Trigger *t = e->getTrigger();
        t->setPersistent(true);
//.........这里部分代码省略.........

/* create MOLECULE */
map< string,Id > SbmlReader::createMolecule( map< string,Id > &idMap )
{	
	map<string,Id>molMap;
	static const Cinfo* moleculeCinfo = initMoleculeCinfo();
	static const Finfo* modeFinfo = moleculeCinfo->findFinfo( "mode" );
	static const Finfo* nInitFinfo = moleculeCinfo->findFinfo( "nInit" );	
	static const Cinfo* kincomptCinfo = initKinComptCinfo();
	static const Finfo* dimensionFinfo = kincomptCinfo->findFinfo( "numDimensions" );
	static const Finfo* sizeFinfo = kincomptCinfo->findFinfo( "size" );
	int num_species = model_->getNumSpecies();
	//cout << "num species: " << num_species << endl;
	for ( int sindex = 0; sindex < num_species; sindex++ )
	{
		Species* s = model_->getSpecies(sindex);
		if (!s){
			//cout << "species " << sindex << " is nul" << endl;
			continue;
		}
		std::string compt = "";		
		if ( s->isSetCompartment() )		
			compt = s->getCompartment();
		if (compt.length()< 1){
			cout << "compt is empty for species "<< sindex << endl;
			continue;
		}
		string id = s->getId();
		//cout<<"species is :"<<id<<endl;
		if (id.length() < 1){
			continue;
		}
		/*std::string name = "";
		if ( s->isSetName() ){
			name = s->getName();
		}*/ 
		Id comptEl = idMap[compt];
		molecule_ = Neutral::create( "Molecule",id,comptEl,Id::scratchId() );//create Molecule
		molMap[id] = comptEl; 
		elmtMap_[id] = Eref( molecule_ );
		//printNotes(s);
		UnitDefinition * ud = s->getDerivedUnitDefinition();
		assert(ud != NULL);
		/*string udef = UnitDefinition::printUnits(ud);
		//cout << "species unit :" << udef << endl;*/
		double initvalue =0.0;
		if ( s->isSetInitialConcentration() )
			initvalue = s->getInitialConcentration();
		else if ( s->isSetInitialAmount() )
			initvalue = s->getInitialAmount() ;
		else {
			cout << "Invalid SBML: Either initialConcentration or initialAmount must be set." << endl;
			return molMap;
		}
		double transvalue = transformUnits(1,ud);
		//cout << "initAmount: " << initvalue << endl;
		transvalue *= initvalue;
		//cout << "transvalue: " << transvalue << endl;
		unsigned int dimension;
        	get< unsigned int >( comptEl.eref(), dimensionFinfo,dimension );
		bool initconc = s->isSetInitialConcentration();
		if ( dimension > 0 &&  initconc) { 
			double size;
			get< double > (comptEl.eref(),sizeFinfo,size); 		
			transvalue *= size;	
		}
		::set< double >( molecule_, nInitFinfo, transvalue ); //initialAmount 	
		bool cons=s->getConstant(); 
		bool bcondition = s->getBoundaryCondition();
		if (cons){
			::set< int >( molecule_,modeFinfo,4 );
			//cout << "set constant 4" << endl;
		}else if (bcondition) {
			::set< int >( molecule_,modeFinfo,1 );
			//cout << "set constant 1" << endl;
		}else {
			::set< int >( molecule_,modeFinfo,0 );
			//cout << "set constant 0" << endl;	
		}
	}
	return molMap;
}

/* create MOLECULE  */
const SbmlReader::sbmlStr_mooseId SbmlReader::createMolecule( map< string,Id > &comptSidMIdMap) {
    Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
    map< string, Id >molSidcmptMIdMap;
    double transvalue = 0.0;
    int num_species = model_->getNumSpecies();

    if (num_species == 0) {
        baseId = Id();
        errorFlag_ = true;
        return molSidcmptMIdMap;
    }

    for ( int sindex = 0; sindex < num_species; sindex++ ) {
        Species* spe = model_->getSpecies(sindex);
        if (!spe) {
            continue;
        }
        std::string compt = "";
        if ( spe->isSetCompartment() ) {
            compt = spe->getCompartment();
        }
        if (compt.length()< 1) {
            //cout << "compt is empty for species "<< sindex << endl;
            continue;
        }
        string id = spe->getId();
        if (id.length() < 1) {
            continue;
        }
        std::string name = "";
        if ( spe->isSetName() ) {
            name = spe->getName();
            name = nameString(name);
        }
        if (name.empty())
            name = id;

        double initvalue =0.0;
        if ( spe->isSetInitialConcentration() )
            initvalue = spe->getInitialConcentration();
        else if ( spe->isSetInitialAmount() )
            initvalue = spe->getInitialAmount() ;
        else {
            unsigned int nr = model_->getNumRules();
            bool found = false;
            for ( unsigned int r = 0; r < nr; r++ ) {
                Rule * rule = model_->getRule(r);
                bool assignRule = rule->isAssignment();
                if ( assignRule ) {
                    string rule_variable = rule->getVariable();
                    if (rule_variable.compare(id) == 0) {
                        found = true;
                        break;
                    }
                }
            }
            if (found == false) {
                cout << "Invalid SBML: Either initialConcentration or initialAmount must be set or it should be found in assignmentRule but non happening for " << spe->getName() <<endl;
                return molSidcmptMIdMap;
            }
        }
        Id comptEl = comptSidMIdMap[compt];
        Id meshEntry = Neutral::child( comptEl.eref(), "mesh" );
        bool constant = spe->getConstant();
        bool boundaryCondition = spe->getBoundaryCondition();
        if (boundaryCondition == true)
            cout << "Pools having BoundaryCondition true " << name <<endl;
        Id pool;
        //If constant is true then its equivalent to BuffPool in moose
        if (constant == true)
            //if( (boundaryCondition == true) && (constant==false))
            pool = shell->doCreate("BufPool",comptEl,name,1);
        else
            pool = shell->doCreate("Pool", comptEl, name ,1);
        molSidcmptMIdMap[id] = comptEl;
        //Map to Molecule SBML id to Moose Id
        molSidMIdMap_[id] = pool;

        //shell->doAddMsg( "OneToOne",pool, "mesh", meshEntry, "mesh" );
        bool bcondition = spe->getBoundaryCondition();
        if ( constant == true && bcondition == false)
            cout <<"The species "<< name << " should not appear in reactant or product as per sbml Rules"<< endl;

        unsigned int spatialDimen =Field< unsigned int >::get( comptEl, "numDimensions");

        UnitDefinition * ud = spe->getDerivedUnitDefinition();
        assert(ud != NULL);
        bool hasonlySubUnit = spe->getHasOnlySubstanceUnits();

        //double v = Field< double >::get( comptEl.path(), "volume" );
        transvalue = transformUnits(1,ud,"substance",hasonlySubUnit);
        if (hasonlySubUnit) {
            // In Moose, no. of molecules (nInit) and unit is "item"
            if (spatialDimen > 0 && spe->isSetInitialAmount() ) {
                //transvalue *= initvalue;
                initvalue *=transvalue;
                Field < double> :: set( pool, "nInit", initvalue);
            }
        } else {
//.........这里部分代码省略.........