C++ Shell::doAddMsg方法代码示例

/* Enzymatic Reaction  */
void SbmlReader::setupEnzymaticReaction( const EnzymeInfo & einfo,string enzname, const map< string, Id > &molSidcmptMIdMap,string name1) {
    string enzPool = einfo.enzyme;

    Id comptRef = molSidcmptMIdMap.find(enzPool)->second; //gives compartment of sp
    Id meshEntry = Neutral::child( comptRef.eref(), "mesh" );
    Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );

    //Creating enz pool to enzyme site
    Id enzPoolId = molSidMIdMap_.find(enzPool)->second;
    Id enzyme_ = shell->doCreate("Enz", enzPoolId, name1, 1);
    //shell->doAddMsg( "Single", meshEntry, "remeshReacs", enzyme_, "remesh");

    Id complex = einfo.complex;
    //Moving enzyme site under enzyme
    shell->doMove(complex,enzyme_);
    shell->doAddMsg("OneToAll",enzyme_,"cplx",complex,"reac");

    shell->doAddMsg("OneToOne",enzyme_,"enz",enzPoolId,"reac");

    vector< Id >::const_iterator sub_itr;
    for ( sub_itr = einfo.substrates.begin(); sub_itr != einfo.substrates.end(); sub_itr++ ) {
        Id S = (*sub_itr);
        Id b = shell->doAddMsg( "OneToOne", enzyme_, "sub" ,S , "reac" );
    }

    vector< Id >::const_iterator prd_itr;
    for ( prd_itr = einfo.products.begin(); prd_itr != einfo.products.end(); prd_itr++ ) {
        Id P = (*prd_itr);
        shell->doAddMsg ("OneToOne",enzyme_,"prd", P,"reac");
    }
    // populate k3,k2,k1 in this order only.
    Field < double > :: set( enzyme_, "k3", einfo.k3 );
    Field < double > :: set( enzyme_, "k2", einfo.k2 );
    Field < double > :: set( enzyme_, "k1", einfo.k1 );
}

void SbmlReader::getRules() {
    unsigned int nr = model_->getNumRules();
    //if (nr > 0)
    //  cout << "\n ##### Need to populate funcpool and sumtotal which is pending due to equations \n";
    Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
    for ( unsigned int r = 0; r < nr; r++ ) {

        Rule * rule = model_->getRule(r);
        bool assignRule = rule->isAssignment();
        if ( assignRule ) {
            string rule_variable = rule->getVariable();
            map< string,Id >::iterator v_iter;
            map< string,Id >::iterator m_iter;
            v_iter = molSidMIdMap_.find( rule_variable );
            if (v_iter != molSidMIdMap_.end()) {
                Id rVariable = molSidMIdMap_.find(rule_variable)->second;
                string rstring =molSidMIdMap_.find(rule_variable)->first;
                Id sumId = shell->doCreate( "SumFunc", rVariable, "func", 1 );
                rVariable.element()->zombieSwap( FuncPool::initCinfo() );
                ObjId ret = shell->doAddMsg( "single",
                                             ObjId( sumId, 0 ), "output",
                                             ObjId( rVariable, 0 ), "input" );
                assert( ret != ObjId() );
                const ASTNode * ast = rule->getMath();
                vector< string > ruleMembers;
                ruleMembers.clear();
                printMembers( ast,ruleMembers );
                for ( unsigned int rm = 0; rm < ruleMembers.size(); rm++ ) {
                    m_iter = molSidMIdMap_.find( ruleMembers[rm] );
                    if ( m_iter != molSidMIdMap_.end() ) {
                        Id rMember = molSidMIdMap_.find(ruleMembers[rm])->second;
                        ObjId ret = shell->doAddMsg( "single",
                                                     ObjId( rMember, 0 ), "nOut",
                                                     ObjId( sumId, 0 ), "input" );
                        string test = molSidMIdMap_.find(ruleMembers[rm])->first;
                    } else {
                        cerr << "SbmlReader::getRules: Assignment rule member is not a species" << endl;
                        // In assignment rule there are constants instead of molecule which is yet to deal in moose.
                        errorFlag_ = true;
                    }
                }
            }
        }
        bool rateRule = rule->isRate();
        if ( rateRule ) {
            cout << "warning : for now Rate Rule is not handled " << endl;
            errorFlag_ = true;
        }
        bool  algebRule = rule->isAlgebraic();
        if ( algebRule ) {
            cout << "warning: for now Algebraic Rule is not handled" << endl;
            errorFlag_ = true;
        }
    }
}

void testMMenzProcess()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	//////////////////////////////////////////////////////////////////////
	// This set is the test kinetic calculation using MathFunc
	//////////////////////////////////////////////////////////////////////
	Id nid = shell->doCreate( "Neutral", Id(), "n", 1 );
	//////////////////////////////////////////////////////////////////////
	// This set is the reference kinetic calculation using MMEnz
	//////////////////////////////////////////////////////////////////////
	Id pid = shell->doCreate( "Pool", nid, "p", 1 ); // substrate
	Id qid = shell->doCreate( "Pool", nid, "q", 1 );	// enz mol
	Id rid = shell->doCreate( "Pool", nid, "r", 1 ); // product
	Id mmid = shell->doCreate( "MMenz", nid, "mm", 1 ); // mmenz

	Id tabid2 = shell->doCreate( "Table", nid, "tab2", 1 ); //output plot

	Field< double >::set( mmid, "Km", 1.0 );
	Field< double >::set( mmid, "kcat", 1.0 );
	Field< double >::set( pid, "nInit", 1.0 );
	Field< double >::set( qid, "nInit", 1.0 );
	Field< double >::set( rid, "nInit", 0.0 );

	shell->doAddMsg( "Single", ObjId( mmid ), "sub", ObjId( pid ), "reac" );
	shell->doAddMsg( "Single", ObjId( mmid ), "prd", ObjId( rid ), "reac" );
	shell->doAddMsg( "Single", ObjId( qid ), "nOut", ObjId( mmid ), "enzDest" );
	shell->doAddMsg( "Single", ObjId( pid ), "nOut", ObjId( tabid2 ), "input" );
	shell->doSetClock( 0, 0.01 );
	shell->doSetClock( 1, 0.01 );
	shell->doUseClock( "/n/mm,/n/tab2", "process", 0 );
	shell->doUseClock( "/n/#[ISA=Pool]", "process", 1 );
	
	//////////////////////////////////////////////////////////////////////
	// Now run models and compare outputs
	//////////////////////////////////////////////////////////////////////

	shell->doReinit();
	shell->doStart( 10 );

	vector< double > vec = Field< vector< double > >::get( tabid2, "vec" );
	assert( vec.size() == 1001 );
	for ( unsigned int i = 0; i < vec.size(); ++i ) {
		double t = 0.01 * i;
		double et = estT( vec[i] );
		assert( doubleApprox( t, et ) );
	}

	shell->doDelete( nid );
	cout << "." << flush;
}

void SbmlReader::addSubPrd(Reaction * reac,Id reaction_,string type) {
    map< string,double > rctMap;
    map< string,double >::iterator rctMap_iter;
    double rctcount=0.0;
    Shell * shell = reinterpret_cast< Shell* >( Id().eref().data() );
    rctMap.clear();
    unsigned int nosubprd;
    const SpeciesReference* rct;
    if (type == "sub") {
        nosubprd = reac->getNumReactants();
    } else
        nosubprd = reac->getNumProducts();
    for ( unsigned int rt=0; rt<nosubprd; rt++ ) {
        if (type == "sub")
            rct = reac->getReactant(rt);
        else
            rct = reac->getProduct(rt);
        std:: string sp = rct->getSpecies();
        rctMap_iter = rctMap.find(sp);
        if ( rctMap_iter != rctMap.end() )
            rctcount = rctMap_iter->second;
        else
            rctcount = 0.0;
        rctcount += rct->getStoichiometry();
        rctMap[sp] = rctcount;
        if (type =="sub")
            noOfsub_ +=rctcount;
        for ( int i=0; (int)i<rct->getStoichiometry(); i++ )
            shell->doAddMsg( "OneToOne", reaction_, type ,molSidMIdMap_[sp] , "reac" );
    }
}

bool ReadSwc::build( Id parent, 
				double lambda, double RM, double RA, double CM )
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	vector< Id > compts( segs_.size() );
	for ( unsigned int i = 0; i < branches_.size(); ++i ) {
		SwcBranch& br = branches_[i];
		for ( unsigned int j = 0; j < br.segs_.size(); ++j ) {
			Id compt;
			SwcSegment& seg = segs_[ br.segs_[j] -1 ];
			unsigned int paIndex = seg.parent();
			if ( paIndex == ~0U ) { // soma
				compt = makeCompt( parent, seg, seg, RM, RA, CM, i, j );
			} else {
				SwcSegment& pa = segs_[ paIndex - 1 ];
				compt = makeCompt( parent, seg, pa, RM, RA, CM, i, j );
				assert( compt != Id() );
				assert( compts[ paIndex -1 ] != Id() );
				shell->doAddMsg( "Single", 
					compts[paIndex-1], "axial", compt, "raxial" );
			}
			assert( compt != Id() );
			compts[ seg.myIndex() -1 ] = compt;
		}
	}
	return true;
}

/*    set up Michalies Menten reaction  */
void SbmlReader::setupMMEnzymeReaction( Reaction * reac,string rid,string rname,const map< string, Id > &molSidcmptMIdMap ) {
    string::size_type loc = rid.find( "_MM_Reaction_" );
    if( loc != string::npos ) {
        int strlen = rid.length();
        rid.erase( loc,strlen-loc );
    }
    unsigned int num_modifr = reac->getNumModifiers();
    for ( unsigned int m = 0; m < num_modifr; m++ ) {
        const ModifierSpeciesReference* modifr=reac->getModifier( m );
        std::string sp = modifr->getSpecies();
        Id enzyme_;
        Id E = molSidMIdMap_.find(sp)->second;

        Id comptRef = molSidcmptMIdMap.find(sp)->second; //gives compartment of sp
        Id meshEntry = Neutral::child( comptRef.eref(), "mesh" );
        Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
        enzyme_ = shell->doCreate("MMenz",E,rname,1);
        //shell->doAddMsg( "Single", meshEntry, "remeshReacs", enzyme_, "remesh");
        shell->doAddMsg("Single",E,"nOut",enzyme_,"enzDest");

        KineticLaw * klaw=reac->getKineticLaw();
        vector< double > rate;
        rate.clear();
        getKLaw( klaw,true,rate );
        if ( errorFlag_ )
            return;
        else if ( !errorFlag_ ) {
            for ( unsigned int rt = 0; rt < reac->getNumReactants(); rt++ ) {
                const SpeciesReference* rct = reac->getReactant( rt );
                sp=rct->getSpecies();
                Id S = molSidMIdMap_.find(sp)->second;
                shell->doAddMsg( "OneToOne", enzyme_, "sub" ,S , "reac" );
            }
            for ( unsigned int pt = 0; pt < reac->getNumProducts(); pt++ ) {
                const SpeciesReference* pdt = reac->getProduct(pt);
                sp = pdt->getSpecies();
                Id P = molSidMIdMap_.find(sp)->second;
                shell->doAddMsg( "OneToOne", enzyme_, "prd" ,P, "reac" );
            }
            Field < double > :: set( enzyme_, "kcat", rate[0] );
            Field < double > :: set( enzyme_, "numKm", rate[1] );
        }

    }
}

void testReacVolumeScaling()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id comptId = shell->doCreate( "CubeMesh", Id(), "cube", 1 );
	Id meshId( comptId.value() + 1 );
	Id subId = shell->doCreate( "Pool", comptId, "sub", 1 );
	Id prdId = shell->doCreate( "Pool", comptId, "prd", 1 );
	Id reacId = shell->doCreate( "Reac", comptId, "reac", 1 );

	double vol1 = 1e-15;

	ObjId mid = shell->doAddMsg( "OneToOne", 
		subId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );
	mid = shell->doAddMsg( "OneToOne", 
		prdId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );

	vector< double > coords( 9, 10.0e-6 );
	coords[0] = coords[1] = coords[2] = 0;

	Field< vector< double > >::set( comptId, "coords", coords );

	double volume = Field< double >::get( comptId, "volume" );
	assert( doubleEq( volume, vol1 ) );

	ObjId ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	ret = shell->doAddMsg( "Single", reacId, "prd", prdId, "reac" );
	assert( ret != ObjId() );

	Field< double >::set( reacId, "Kf", 2 );
	Field< double >::set( reacId, "Kb", 3 );
	double x = Field< double >::get( reacId, "kf" );
	assert( doubleEq( x, 2 ) );
	x = Field< double >::get( reacId, "kb" );
	assert( doubleEq( x, 3 ) );
	
	ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	double conv = 1.0 / ( NA * vol1 );
	x = Field< double >::get( reacId, "kf" );
	assert( doubleEq( x, 2 * conv ) );
	x = Field< double >::get( reacId, "kb" );
	assert( doubleEq( x, 3 ) );

	ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	ret = shell->doAddMsg( "Single", reacId, "prd", prdId, "reac" );
	assert( ret != ObjId() );
	x = Field< double >::get( reacId, "kf" );
	assert( doubleEq( x, 2 * conv * conv ) );
	x = Field< double >::get( reacId, "kb" );
	assert( doubleEq( x, 3 * conv ) );

	shell->doDelete( comptId );
	cout << "." << flush;
}

// See what Element::getNeighbors does with 2 sub <----> prd.
void testTwoReacGetNeighbors()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id comptId = shell->doCreate( "CubeMesh", Id(), "cube", 1 );
	Id meshId( comptId.value() + 1 );
	Id subId = shell->doCreate( "Pool", comptId, "sub", 1 );
	Id prdId = shell->doCreate( "Pool", comptId, "prd", 1 );
	Id reacId = shell->doCreate( "Reac", comptId, "reac", 1 );

	ObjId mid = shell->doAddMsg( "OneToOne", 
		subId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );
	mid = shell->doAddMsg( "OneToOne", 
		prdId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );

	ObjId ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );

	ret = shell->doAddMsg( "Single", reacId, "prd", prdId, "reac" );
	assert( ret != ObjId() );

	vector< Id > pools;
	unsigned int num = reacId.element()->getNeighbors( pools, 
		Reac::initCinfo()->findFinfo( "toSub" ) );
	assert( num == 2 );
	assert( pools[0] == subId );
	assert( pools[1] == subId );

	pools.clear();
	num = reacId.element()->getNeighbors( pools, 
		Reac::initCinfo()->findFinfo( "sub" ) );
	assert( num == 2 );
	assert( pools[0] == subId );
	assert( pools[1] == subId );

	shell->doDelete( comptId );
	cout << "." << flush;
}

void testPoolVolumeScaling()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id comptId = shell->doCreate( "CylMesh", Id(), "cyl", 1 );
	Id meshId( comptId.value() + 1 );
	Id poolId = shell->doCreate( "Pool", comptId, "pool", 1 );

	ObjId mid = shell->doAddMsg( "OneToOne", 
		ObjId( poolId, 0 ), "requestVolume",
		ObjId( meshId, 0 ), "get_volume" );

	assert( mid != ObjId() );

	vector< double > coords( 9, 0.0 );
	double x1 = 100e-6;
	double r0 = 10e-6;
	double r1 = 5e-6;
	double lambda = x1;
	coords[3] = x1;
	coords[6] = r0;
	coords[7] = r1;
	coords[8] = lambda;

	Field< vector< double > >::set( comptId, "coords", coords );

	double volume = Field< double >::get( poolId, "volume" );
	assert( doubleEq( volume, PI * x1 * (r0+r1) * (r0+r1) / 4.0 ) );

	Field< double >::set( poolId, "n", 400 );
	double volscale = 1 / ( NA * volume );
	double conc = Field< double >::get( poolId, "conc" );
	assert( doubleEq( conc, 400 * volscale ) );
	Field< double >::set( poolId, "conc", 500 * volscale );
	double n = Field< double >::get( poolId, "n" );
	assert( doubleEq( n, 500 ) );

	Field< double >::set( poolId, "nInit", 650 );
	double concInit = Field< double >::get( poolId, "concInit" );
	assert( doubleEq( concInit, 650 * volscale ) );
	Field< double >::set( poolId, "concInit", 10 * volscale );
	n = Field< double >::get( poolId, "nInit" );
	assert( doubleEq( n, 10 ) );

	shell->doDelete( comptId );
	cout << "." << flush;
}

void ReadCell::addChannelMessage( Id chan )
{
	/*
	* Get child objects of type Mstring, named addmsg1, 2, etc.
	* These define extra messages to be assembled at setup.
	* Similar to what was done with GENESIS.
	*/
	vector< Id > kids;
	Neutral::children( chan.eref(), kids );
	Shell *shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id cwe = shell->getCwe();
	shell->setCwe( chan );
	for ( vector< Id >::iterator i = kids.begin(); i != kids.end(); ++i )
	{
		// Ignore kid if its name does not begin with "addmsg"..
		const string& name = i->element()->getName();
		if ( name.find( "addmsg", 0 ) != 0 )
			continue;

		string s = Field< string >::get( *i, "value" );
		vector< string > token;
		tokenize( s, " 	", token );
		assert( token.size() == 4 );
		ObjId src = shell->doFind( token[0] );
		ObjId dest = shell->doFind( token[2] );

		// I would like to assert, or warn here, but there are legitimate 
		// cases where not all possible messages are actually available 
		// to set up. So I just bail.
		if ( src.bad() || dest.bad()) {
#ifndef NDEBUG
				/*
			cout << "ReadCell::addChannelMessage( " << chan.path() << 
				"): " << name << " " << s << 
				": Bad src " << src << " or dest " << dest << endl;
				*/
#endif
			continue; 
		}
		ObjId mid = 
			shell->doAddMsg( "single", src, token[1], dest, token[3] );
		assert( !mid.bad());
	}
	shell->setCwe( cwe );
}

void testCompartmentProcess()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	unsigned int size = 100;
	double Rm = 1.0;
	double Ra = 0.01;
	double Cm = 1.0;
	double dt = 0.01;
	double runtime = 10;
	double lambda = sqrt( Rm / Ra );

	Id cid = shell->doCreate( "Compartment", Id(), "compt", size );
	assert( Id::isValid(cid));
	assert( cid.eref().element()->numData() == size );

	bool ret = Field< double >::setRepeat( cid, "initVm", 0.0 );
	assert( ret );
	Field< double >::setRepeat( cid, "inject", 0 );
	// Only apply current injection in first compartment
	Field< double >::set( ObjId( cid, 0 ), "inject", 1.0 ); 
	Field< double >::setRepeat( cid, "Rm", Rm );
	Field< double >::setRepeat( cid, "Ra", Ra );
	Field< double >::setRepeat( cid, "Cm", Cm );
	Field< double >::setRepeat( cid, "Em", 0 );
	Field< double >::setRepeat( cid, "Vm", 0 );

	// The diagonal message has a default stride of 1, so it connects
	// successive compartments.
	// Note that the src and dest elements here are identical, so we cannot
	// use a shared message. The messaging system will get confused about
	// direction to send data. So we split up the shared message that we
	// might have used, below, into two individual messages.
	// MsgId mid = shell->doAddMsg( "Diagonal", ObjId( cid ), "raxial", ObjId( cid ), "axial" );
	ObjId mid = shell->doAddMsg( "Diagonal", ObjId( cid ), "axialOut", ObjId( cid ), "handleAxial" );
	assert( !mid.bad());
	// mid = shell->doAddMsg( "Diagonal", ObjId( cid ), "handleRaxial", ObjId( cid ), "raxialOut" );
	mid = shell->doAddMsg( "Diagonal", ObjId( cid ), "raxialOut", ObjId( cid ), "handleRaxial" );
	assert( !mid.bad() );
	// ObjId managerId = Msg::getMsg( mid )->manager().objId();
	// Make the raxial data go from high to lower index compartments.
	Field< int >::set( mid, "stride", -1 );

#ifdef DO_SPATIAL_TESTS
	shell->doSetClock( 0, dt );
	shell->doSetClock( 1, dt );
	// Ensure that the inter_compt msgs go between nodes once every dt.
	shell->doSetClock( 9, dt ); 

	shell->doUseClock( "/compt", "init", 0 );
	shell->doUseClock( "/compt", "process", 1 );

	shell->doReinit();
	shell->doStart( runtime );

	double Vmax = Field< double >::get( ObjId( cid, 0 ), "Vm" );

	double delta = 0.0;
	// We measure only the first 50 compartments as later we 
	// run into end effects because it is not an infinite cable
	for ( unsigned int i = 0; i < size; i++ ) {
		double Vm = Field< double >::get( ObjId( cid, i ), "Vm" );
		double x = Vmax * exp( - static_cast< double >( i ) / lambda );
		delta += ( Vm - x ) * ( Vm - x );
	 	// cout << i << " (x, Vm) = ( " << x << ", " << Vm << " )\n";
	}
	assert( delta < 1.0e-5 );
#endif // DO_SPATIAL_TESTS
	shell->doDelete( cid );
	cout << "." << flush;
}

//.........这里部分代码省略.........
                                    // if plots exist then will be placing at "/data"
                                    Id graphs;
                                    //Id dataId;
                                    Id dataIdTest;
                                    if (parentId2 == Id())
                                        graphs = s->doCreate( "Neutral", parentId2, "data", 1);
                                    else
                                        // need to check how to put / while coming from gui as the path is /model/modelName??? 27 jun 2014
                                        findModelParent ( Id(), modelName, dataIdTest, modelName ) ;
                                        string test = "/data";
                                        Id tgraphs(test);
                                        graphs=tgraphs;
                                        //graphs = s->doCreate("Neutral",parentId,"data",1);
                                        //Id dataId;
                                        //if (dataId == Id())
                                        //    cout << "Id " << dataId;
                                        //    graphs = s->doCreate( "Neutral",dataId, "data", 1);
                                        assert( graphs != Id() );
                                    plotValue = (grandChildNode.getChild(0).toXMLString()).c_str();
                                    istringstream pltVal(plotValue);
                                    string pltClean;
                                    while (getline(pltVal,pltClean, ';')) {
                                        pltClean.erase( remove( pltClean.begin(), pltClean.end(), ' ' ), pltClean.end() );
                                        //string plotPath = location+pltClean;
                                        string plotPath = base_.path()+pltClean;
                                        Id plotSId(plotPath);
                                        unsigned pos = pltClean.find('/');
                                        if (pos != std::string::npos)
                                            pltClean = pltClean.substr(pos+1,pltClean.length());
                                        replace(pltClean.begin(),pltClean.end(),'/','_');
                                        string plotName =  pltClean + ".conc";
                                        Id pltPath(graphs.path());
                                        Id tab = s->doCreate( "Table", pltPath, plotName, 1 );
                                        if (tab != Id())
                                            s->doAddMsg("Single",tab,"requestOut",plotSId,"getConc");
                                    }//while
                                    /* passing /model and /data         */
                                    string comptPath =base_.path()+"/##";
                                    s->doUseClock(comptPath,"process",4);

                                    string tablePath = graphs.path()+"/##[TYPE=Table]";
                                    s->doUseClock( tablePath, "process",8 );
                                }//plots
                                /*else
                                  nodeValue = atof ((grandChildNode.getChild(0).toXMLString()).c_str());

                                  if (nodeName == "runTime")
                                  sm->setRunTime(nodeValue);
                                  else if (nodeName == "simdt")
                                  sm->setSimDt(nodeValue);
                                  else if(nodeName == "plotdt")
                                  sm->setPlotDt(nodeValue);
                                  */

                            } //grandChild
                            else
                                cout << "Warning: expected exactly ONE child of " << nodeName << " but none found "<<endl;
                        } //gchild
                    } //moose and modelAnnotation
                }
            }//annotation Node
            else {
                //4 for simdt and 8 for plotdt
                s->doUseClock(base_.path()+"/##","process",4);
                s->doUseClock(+"/data/##[TYPE=Table]","process",8);
                s->doSetClock(4,0.1);
                s->doSetClock(8,0.1);

            }
            vector< ObjId > compts;
            string comptpath = base_.path()+"/##[ISA=ChemCompt]";
            wildcardFind( comptpath, compts );
            vector< ObjId >::iterator i = compts.begin();
            string comptName = nameString(Field<string> :: get(ObjId(*i),"name"));
            string simpath = base_.path() + "/##";
            s->doUseClock( simpath, "process", 4 );

            //wildcardFind( plotpath, plots );
            //Id pathexist(base_.path()+"/kinetics");
            /*
               if (solverClass.empty())
               {
               if( pathexist != Id())
               sm->build(base_.eref(),&q,"rk5");
               else
               sm->buildForSBML(base_.eref(),&q,"rk5");
               }
               else
               { if(pathexist != Id())
               sm->build(base_.eref(),&q,solverClass);
               else
               sm->buildForSBML(base_.eref(),&q,solverClass);
               }
               */
            return base_;
        }

    } else
        return baseId;
}

void testHSolveUtils( )
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	bool success;
	
	Id n = shell->doCreate( "Neutral", Id(), "n" );
	
	/**
	 *  First we test the functions which return the compartments linked to a
	 *  given compartment: adjacent(), and children().
	 *  
	 *  A small tree is created for this:
	 *  
	 *               c0
	 *                L c1
	 *                   L c2
	 *                   L c3
	 *                   L c4
	 *                   L c5
	 *  
	 *  (c0 is the parent of c1. c1 is the parent of c2, c3, c4, c5.)
	 */
	Id c[ 6 ];
	c[ 0 ] = shell->doCreate( "Compartment", n, "c0" );
	c[ 1 ] = shell->doCreate( "Compartment", n, "c1" );
	c[ 2 ] = shell->doCreate( "Compartment", n, "c2" );
	c[ 3 ] = shell->doCreate( "Compartment", n, "c3" );
	c[ 4 ] = shell->doCreate( "Compartment", n, "c4" );
	c[ 5 ] = shell->doCreate( "Compartment", n, "c5" );
	
	MsgId mid;
	mid = shell->doAddMsg( "Single", c[ 0 ], "axial", c[ 1 ], "raxial" );
	ASSERT( mid != Msg::bad, "Linking compartments" );
	mid = shell->doAddMsg( "Single", c[ 1 ], "axial", c[ 2 ], "raxial" );
	ASSERT( mid != Msg::bad, "Linking compartments" );
	mid = shell->doAddMsg( "Single", c[ 1 ], "axial", c[ 3 ], "raxial" );
	ASSERT( mid != Msg::bad, "Linking compartments" );
	mid = shell->doAddMsg( "Single", c[ 1 ], "axial", c[ 4 ], "raxial" );
	ASSERT( mid != Msg::bad, "Linking compartments" );
	mid = shell->doAddMsg( "Single", c[ 1 ], "axial", c[ 5 ], "raxial" );
	ASSERT( mid != Msg::bad, "Linking compartments" );
	
	vector< Id > found;
	unsigned int nFound;
	
	/* 
	 * Testing version 1 of HSolveUtils::adjacent.
	 * It finds all neighbours of given compartment.
	 */
	// Neighbours of c0
	nFound = HSolveUtils::adjacent( c[ 0 ], found );
	ASSERT( nFound == found.size(), "Finding adjacent compartments" );
	// c1 is adjacent
	ASSERT( nFound == 1, "Finding adjacent compartments" );
	ASSERT( found[ 0 ] == c[ 1 ], "Finding adjacent compartments" );
	
	// Neighbours of c1
	found.clear();
	nFound = HSolveUtils::adjacent( c[ 1 ], found );
	ASSERT( nFound == 5, "Finding adjacent compartments" );
	// c0 is adjacent
	success =
		find( found.begin(), found.end(), c[ 0 ] ) != found.end();
	ASSERT( success, "Finding adjacent compartments" );
	// c2 - c5 are adjacent
	for ( int i = 2; i < 6; i++ ) {
		success =
			find( found.begin(), found.end(), c[ i ] ) != found.end();
		ASSERT( success, "Finding adjacent compartments" );
	}
	
	// Neighbours of c2
	found.clear();
	nFound = HSolveUtils::adjacent( c[ 2 ], found );
	// c1 is adjacent
	ASSERT( nFound == 1, "Finding adjacent compartments" );
	ASSERT( found[ 0 ] == c[ 1 ], "Finding adjacent compartments" );
	
	/*
	 * Testing version 2 of HSolveUtils::adjacent.
	 * It finds all but one neighbours of given compartment.
	 * The the second argument to 'adjacent' is the one that is excluded.
	 */
	// Neighbours of c1 (excluding c0)
	found.clear();
	nFound = HSolveUtils::adjacent( c[ 1 ], c[ 0 ], found );
	ASSERT( nFound == 4, "Finding adjacent compartments" );
	// c2 - c5 are adjacent
	for ( int i = 2; i < 6; i++ ) {
		success =
			find( found.begin(), found.end(), c[ i ] ) != found.end();
		ASSERT( success, "Finding adjacent compartments" );
	}
	
	// Neighbours of c1 (excluding c2)
	found.clear();
	nFound = HSolveUtils::adjacent( c[ 1 ], c[ 2 ], found );
	ASSERT( nFound == 4, "Finding adjacent compartments" );
	// c0 is adjacent
	success =
//.........这里部分代码省略.........

/**
 * Tab controlled by table
 * A + Tab <===> B
 * A + B -sumtot--> tot1
 * 2B <===> C
 * 
 * C ---e1Pool ---> D
 * D ---e2Pool ----> E
 *
 * All these are created on /kinetics, a cube compartment of vol 1e-18 m^3
 *
 */
Id makeReacTest()
{
	double simDt = 0.1;
	double plotDt = 0.1;
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );

	Id pools[10];
	unsigned int i = 0;
	// Make the objects.
	Id kin = s->doCreate( "CubeMesh", Id(), "kinetics", 1 );
	Id tab = s->doCreate( "StimulusTable", kin, "tab", 1 );
	Id T = pools[i++] = s->doCreate( "BufPool", kin, "T", 1 );
	Id A = pools[i++] = s->doCreate( "Pool", kin, "A", 1 );
	Id B = pools[i++] = s->doCreate( "Pool", kin, "B", 1 );
	Id C = pools[i++] = s->doCreate( "Pool", kin, "C", 1 );
	Id D = pools[i++] = s->doCreate( "Pool", kin, "D", 1 );
	Id E = pools[i++] = s->doCreate( "Pool", kin, "E", 1 );
	Id tot1 = pools[i++] = s->doCreate( "BufPool", kin, "tot1", 1 );
	Id sum = s->doCreate( "Function", tot1, "func", 1 ); // Silly that it has to have this name. 
	Id sumInput( sum.value() + 1 );
	Id e1Pool = s->doCreate( "Pool", kin, "e1Pool", 1 );
	Id e2Pool = s->doCreate( "Pool", kin, "e2Pool", 1 );
	Id e1 = s->doCreate( "Enz", e1Pool, "e1", 1 );
	Id cplx = s->doCreate( "Pool", e1, "cplx", 1 );
	Id e2 = s->doCreate( "MMenz", e2Pool, "e2", 1 );
	Id r1 = s->doCreate( "Reac", kin, "r1", 1 );
	Id r2 = s->doCreate( "Reac", kin, "r2", 1 );
	Id plots = s->doCreate( "Table2", kin, "plots", 7 );

	// Connect them up
	s->doAddMsg( "Single", tab, "output", T, "setN" );

	s->doAddMsg( "Single", r1, "sub", T, "reac" );
	s->doAddMsg( "Single", r1, "sub", A, "reac" );
	s->doAddMsg( "Single", r1, "prd", B, "reac" );

	Field< unsigned int >::set( sum, "numVars", 2 );
	s->doAddMsg( "Single", A, "nOut", ObjId( sumInput, 0, 0 ), "input" );
	s->doAddMsg( "Single", B, "nOut", ObjId( sumInput, 0, 1 ), "input" );
	s->doAddMsg( "Single", sum, "valueOut", tot1, "setN" );

	s->doAddMsg( "Single", r2, "sub", B, "reac" );
	s->doAddMsg( "Single", r2, "sub", B, "reac" );
	s->doAddMsg( "Single", r2, "prd", C, "reac" );

	s->doAddMsg( "Single", e1, "sub", C, "reac" );
	s->doAddMsg( "Single", e1, "enz", e1Pool, "reac" );
	s->doAddMsg( "Single", e1, "cplx", cplx, "reac" );
	s->doAddMsg( "Single", e1, "prd", D, "reac" );

	s->doAddMsg( "Single", e2, "sub", D, "reac" );
	s->doAddMsg( "Single", e2Pool, "nOut", e2, "enzDest" );
	s->doAddMsg( "Single", e2, "prd", E, "reac" );

	// Set parameters.
	Field< double >::set( A, "concInit", 2 );
	Field< double >::set( e1Pool, "concInit", 1 );
	Field< double >::set( e2Pool, "concInit", 1 );
	Field< string >::set( sum, "expr", "x0+x1" );
	Field< double >::set( r1, "Kf", 0.2 );
	Field< double >::set( r1, "Kb", 0.1 );
	Field< double >::set( r2, "Kf", 0.1 );
	Field< double >::set( r2, "Kb", 0.0 );
	Field< double >::set( e1, "Km", 5 );
	Field< double >::set( e1, "kcat", 1 );
	Field< double >::set( e1, "ratio", 4 );
	Field< double >::set( e2, "Km", 5 );
	Field< double >::set( e2, "kcat", 1 );
	vector< double > stim( 100, 0.0 );
	double vol = Field< double >::get( kin, "volume" );
	for ( unsigned int i = 0; i< 100; ++i ) {
		stim[i] = vol * NA * (1.0 + sin( i * 2.0 * PI / 100.0 ) );
	}
	Field< vector< double > >::set( tab, "vector", stim );
	Field< double >::set( tab, "stepSize", 0.0 );
	Field< double >::set( tab, "stopTime", 10.0 );
	Field< double >::set( tab, "loopTime", 10.0 );
	Field< bool >::set( tab, "doLoop", true );
	

	// Connect outputs
	for ( unsigned int i = 0; i < 7; ++i )
		s->doAddMsg( "Single", ObjId( plots,i), 
						"requestOut", pools[i], "getConc" );

	// Schedule it.
	for ( unsigned int i = 11; i < 18; ++i )
		s->doSetClock( i, simDt );
//.........这里部分代码省略.........

//.........这里部分代码省略.........

        int count = -1;
        for ( unsigned int a = 0; a < arraySize; a++ )
            if ( array[ a ] == -1 )
                count++;
            else
                tree[ count ].children.push_back( array[ a ] );

        //////////////////////////////////////////
        // Create cell inside moose; setup solver.
        //////////////////////////////////////////
        Id n = shell->doCreate( "Neutral", Id(), "n", 1 );

        vector< Id > c( nCompt );
        for ( i = 0; i < nCompt; i++ )
        {
            ostringstream name;
            name << "c" << i;
            c[ i ] = shell->doCreate( "Compartment", n, name.str() , 1);

            Field< double >::set( c[ i ], "Ra", tree[ i ].Ra );
            Field< double >::set( c[ i ], "Rm", tree[ i ].Rm );
            Field< double >::set( c[ i ], "Cm", tree[ i ].Cm );
            Field< double >::set( c[ i ], "Em", Em[ i ] );
            Field< double >::set( c[ i ], "initVm", V[ i ] );
            Field< double >::set( c[ i ], "Vm", V[ i ] );
        }

        for ( i = 0; i < nCompt; i++ )
        {
            vector< unsigned int >& child = tree[ i ].children;
            for ( j = 0; j < ( int )( child.size() ); j++ )
            {
                ObjId mid = shell->doAddMsg(
                                "Single", c[ i ], "axial", c[ child[ j ] ], "raxial" );
                ASSERT( ! mid.bad(), "Creating test model" );
            }
        }

        HP.setup( c[ 0 ], dt );

        /*
         * Here we check if the cell was read in correctly by the solver.
         * This test only checks if all the created compartments were read in.
         * It doesn't check if they have been assigned hines' indices correctly.
         */
        vector< Id >& hc = HP.compartmentId_;
        ASSERT( ( int )( hc.size() ) == nCompt, "Tree traversal" );
        for ( i = 0; i < nCompt; i++ )
            ASSERT(
                find( hc.begin(), hc.end(), c[ i ] ) != hc.end(), "Tree traversal"
            );

        //////////////////////////////////////////
        // Setup local matrix
        //////////////////////////////////////////

        /*
         * First we need to ensure that the hines' indices for the local model
         * and those inside the solver match. If the numbering is different,
         * then the matrices will not agree.
         *
         * In the following, we find out the indices assigned by the solver,
         * and impose them on the local data structures.
         */