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

void testRunKsolve()
{
	double simDt = 0.1;
	// double plotDt = 0.1;
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	Id kin = makeReacTest();
	Id ksolve = s->doCreate( "Ksolve", kin, "ksolve", 1 );
	Id stoich = s->doCreate( "Stoich", ksolve, "stoich", 1 );
	Field< Id >::set( stoich, "compartment", kin );
	Field< Id >::set( stoich, "ksolve", ksolve );
	Field< string >::set( stoich, "path", "/kinetics/##" );
	s->doUseClock( "/kinetics/ksolve", "process", 4 ); 
	s->doSetClock( 4, simDt );

	s->doReinit();
	s->doStart( 20.0 );
	Id plots( "/kinetics/plots" );
	for ( unsigned int i = 0; i < 7; ++i ) {
		stringstream ss;
		ss << "plot." << i;
		SetGet2< string, string >::set( ObjId( plots, i ), "xplot", 
						"tsr2.plot", ss.str() );
	}
	s->doDelete( kin );
	cout << "." << flush;
}

int main( int argc, char** argv )
{
	bool doUnitTests = 0;
	bool doRegressionTests = 0;
	unsigned int benchmark = 0;
	// This reorders the OpFunc to Fid mapping to ensure it is node and
	// compiler independent.
	Id shellId = init( argc, argv, doUnitTests, doRegressionTests, benchmark );
	// Note that the main loop remains the parser loop, though it may
	// spawn a lot of other stuff.
	Element* shelle = shellId.element();
	Shell* s = reinterpret_cast< Shell* >( shelle->data( 0 ) );
	if ( doUnitTests )
		nonMpiTests( s ); // These tests do not need the process loop.

	if ( Shell::myNode() == 0 ) {
		if ( Shell::numNodes() > 1 ) {
			// Use the last clock for the postmaster, so that it is called
			// after everything else has been processed and all messages
			// are ready to send out.
			s->doUseClock( "/postmaster", "process", 9 );
			s->doSetClock( 9, 1.0 ); // Use a sensible default.
		}
#ifdef DO_UNIT_TESTS
		if ( doUnitTests ) {
			mpiTests();
			processTests( s );
		}
		// if ( doRegressionTests ) regressionTests();
#endif
		// These are outside unit tests because they happen in optimized
		// mode, using a command-line argument. As soon as they are done
		// the system quits, in order to estimate timing.
		if ( benchmark != 0 ) {
			mooseBenchmarks( benchmark );
			s->doQuit();
		} else {
			// Here we set off a little event loop to poll user input. 
			// It deals with the doQuit call too.
                        if(! quitFlag)
                            Shell::launchParser(); 
		}
	} else {
		PostMaster* p = reinterpret_cast< PostMaster* >( ObjId( 3 ).data());
		while ( Shell::keepLooping() ) {
			p->clearPending();
		}
	}
	Msg::clearAllMsgs();
	Id::clearAllElements();
#ifdef USE_MPI
	MPI_Finalize();
#endif
	return 0;
}

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 testTaperingCylDiffn()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 25e-6;
	double r0 = 2e-6;
	double r1 = 1e-6;
	double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
	double runtime = 10.0;
	double dt = 0.1; // 0.2 is the highest dt for which the error is in bounds
	double diffConst = 1.0e-12; 
	// Should set explicitly, currently during creation of DiffPoolVec
	//double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
	Field< double >::set( cyl, "r0", r0 );
	Field< double >::set( cyl, "r1", r1 );
	Field< double >::set( cyl, "x0", 0 );
	Field< double >::set( cyl, "x1", len );
	Field< double >::set( cyl, "diffLength", diffLength );
	unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
	assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
	Id pool = s->doCreate( "Pool", cyl, "pool", 1 );
	Field< double >::set( pool, "diffConst", diffConst );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", cyl );
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	// Next: build by setting the path of the dsolve.
	Field< string >::set( dsolve, "path", "/model/cyl/pool" );
	// Then find a way to test it.
	assert( pool.element()->numData() == ndc );
	Field< double >::set( ObjId( pool, 0 ), "nInit", 1.0 );

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

	double myTot = 0.0;
	vector< double > poolVec;
   	Field< double >::getVec( pool, "n", poolVec );
	for ( unsigned int i = 0; i < poolVec.size(); ++i ) {
		myTot += poolVec[i];
	} 
	assert( doubleEq( myTot, 1.0 ) );

	s->doDelete( model );
	cout << "." << flush;
}

/**
 * The readcell function implements the old GENESIS cellreader
 * functionality. Although it is really a parser operation, I
 * put it here in Kinetics because the cell format is independent
 * of parser and is likely to remain a legacy for a while.
 */
Id ReadKkit::read(
	const string& filename, 
	const string& modelname,
	Id pa, const string& methodArg )
{
	string method = methodArg;
	ifstream fin( filename.c_str() );
	if (!fin){
		cerr << "ReadKkit::read: could not open file " << filename << endl;
		return Id();
    }

	if ( method.substr(0, 4) == "old_" ) {
		moveOntoCompartment_ = false;
		method = method.substr( 4 );
	}

	Shell* s = reinterpret_cast< Shell* >( ObjId().data() );
	Id mgr = makeStandardElements( pa, modelname );
	assert( mgr != Id() );

	baseId_ = mgr;
	basePath_ = mgr.path();
	enzCplxMols_.resize( 0 );

	innerRead( fin );

	assignPoolCompartments();
	assignReacCompartments();
	assignEnzCompartments();
	assignMMenzCompartments();

	convertParametersToConcUnits();

	s->doSetClock( 8, plotdt_ );
	
	string plotpath = basePath_ + "/graphs/##[TYPE=Table]," +
			basePath_ + "/moregraphs/##[TYPE=Table]";
	s->doUseClock( plotpath, "process", 8 );

	setMethod( s, mgr, simdt_, plotdt_, method );

	s->doReinit();
	return mgr;
}

void testRunGsolve()
{
	double simDt = 0.1;
	// double plotDt = 0.1;
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	Id kin = makeReacTest();
	double volume = 1e-21;
	Field< double >::set( kin, "volume", volume );
	Field< double >::set( ObjId( "/kinetics/A" ), "concInit", 2 );
	Field< double >::set( ObjId( "/kinetics/e1Pool" ), "concInit", 1 );
	Field< double >::set( ObjId( "/kinetics/e2Pool" ), "concInit", 1 );
	Id e1( "/kinetics/e1Pool/e1" );
	Field< double >::set( e1, "Km", 5 );
	Field< double >::set( e1, "kcat", 1 );
	vector< double > stim( 100, 0.0 );
	for ( unsigned int i = 0; i< 100; ++i ) {
		stim[i] = volume * NA * (1.0 + sin( i * 2.0 * PI / 100.0 ) );
	}
	Field< vector< double > >::set( ObjId( "/kinetics/tab" ), "vector", stim );


	Id gsolve = s->doCreate( "Gsolve", kin, "gsolve", 1 );
	Id stoich = s->doCreate( "Stoich", gsolve, "stoich", 1 );
	Field< Id >::set( stoich, "compartment", kin );
	Field< Id >::set( stoich, "ksolve", gsolve );
	
	Field< string >::set( stoich, "path", "/kinetics/##" );
	s->doUseClock( "/kinetics/gsolve", "process", 4 ); 
	s->doSetClock( 4, simDt );

	s->doReinit();
	s->doStart( 20.0 );
	Id plots( "/kinetics/plots" );
	for ( unsigned int i = 0; i < 7; ++i ) {
		stringstream ss;
		ss << "plot." << i;
		SetGet2< string, string >::set( ObjId( plots, i ), "xplot", 
						"tsr3.plot", ss.str() );
	}
	s->doDelete( kin );
	cout << "." << flush;
}

/**
 * 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 );
//.........这里部分代码省略.........

void testCalcJunction()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	// Make a neuron with same-size dend and spine. PSD is tiny.
	// Put a, b, c in dend, b, c, d in spine, c, d, f in psd. No reacs.
	// See settling of all concs by diffusion, pairwise.
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id dend = s->doCreate( "Compartment", model, "dend", 1 );
	Id neck = s->doCreate( "Compartment", model, "spine_neck", 1 );
	Id head = s->doCreate( "Compartment", model, "spine_head", 1 );
	Field< double >::set( dend, "x", 10e-6 );
	Field< double >::set( dend, "diameter", 2e-6 );
	Field< double >::set( dend, "length", 10e-6 );
	Field< double >::set( neck, "x0", 9e-6 );
	Field< double >::set( neck, "x", 9e-6 );
	Field< double >::set( neck, "y", 1e-6 );
	Field< double >::set( neck, "diameter", 0.5e-6 );
	Field< double >::set( neck, "length", 1.0e-6 );
	Field< double >::set( head, "x0", 9e-6 );
	Field< double >::set( head, "x", 9e-6 );
	Field< double >::set( head, "y0", 1e-6 );
	Field< double >::set( head, "y", 11e-6 );
	Field< double >::set( head, "diameter", 2e-6 );
	Field< double >::set( head, "length", 10e-6 );
	s->doAddMsg( "Single", ObjId( dend ), "raxial", ObjId( neck ), "axial");
	s->doAddMsg( "Single", ObjId( neck ), "raxial", ObjId( head ), "axial");

	Id nm = s->doCreate( "NeuroMesh", model, "nm", 1 );
	Field< double >::set( nm, "diffLength", 10e-6 );
	Field< bool >::set( nm, "separateSpines", true );
	Id sm = s->doCreate( "SpineMesh", model, "sm", 1 );
	Id pm = s->doCreate( "PsdMesh", model, "pm", 1 );
	ObjId mid = s->doAddMsg( "Single", ObjId( nm ), "spineListOut", ObjId( sm ), "spineList" );
	assert( !mid.bad() );
	mid = s->doAddMsg( "Single", ObjId( nm ), "psdListOut", ObjId( pm ), "psdList" );
	Field< Id >::set( nm, "cell", model );

	vector< Id > pools( 9 );
	static string names[] = {"a", "b", "c", "b", "c", "d", "c", "d", "e" };
	static Id parents[] = {nm, nm, nm, sm, sm, sm, pm, pm, pm};
	for ( unsigned int i = 0; i < 9; ++i ) {
		pools[i] = s->doCreate( "Pool", parents[i], names[i], 1 );
		assert( pools[i] != Id() );
		Field< double >::set( pools[i], "concInit", 1.0 + 1.0 * i );
		Field< double >::set( pools[i], "diffConst", 1e-11 );
		if ( i < 6 ) {
			double vol = Field< double >::get( pools[i], "volume" );
			assert( doubleEq( vol, 10e-6 * 1e-12 * PI ) );
		}
	}
	Id dendsolve = s->doCreate( "Dsolve", model, "dendsolve", 1 );
	Id spinesolve = s->doCreate( "Dsolve", model, "spinesolve", 1 );
	Id psdsolve = s->doCreate( "Dsolve", model, "psdsolve", 1 );
	Field< Id >::set( dendsolve, "compartment", nm );
	Field< Id >::set( spinesolve, "compartment", sm );
	Field< Id >::set( psdsolve, "compartment", pm );
	Field< string >::set( dendsolve, "path", "/model/nm/#" );
	Field< string >::set( spinesolve, "path", "/model/sm/#" );
	Field< string >::set( psdsolve, "path", "/model/pm/#" );
	assert( Field< unsigned int >::get( dendsolve, "numAllVoxels" ) == 1 );
	assert( Field< unsigned int >::get( spinesolve, "numAllVoxels" ) == 1 );
	assert( Field< unsigned int >::get( psdsolve, "numAllVoxels" ) == 1 );
	assert( Field< unsigned int >::get( dendsolve, "numPools" ) == 3 );
	assert( Field< unsigned int >::get( spinesolve, "numPools" ) == 3 );
	assert( Field< unsigned int >::get( psdsolve, "numPools" ) == 3 );
	SetGet2< Id, Id >::set( dendsolve, "buildNeuroMeshJunctions", 
					spinesolve, psdsolve );
	s->doSetClock( 0, 0.01 );
	s->doUseClock( "/model/#solve", "process", 0 );
	s->doReinit();
	s->doStart( 100 );

	for ( unsigned int i = 0; i < 9; ++i ) {
		double c = Field< double >::get( pools[i], "conc" );
		double n = Field< double >::get( pools[i], "n" );
		double v = Field< double >::get( pools[i], "volume" );
		cout << pools[i].path() << ": " << c << ", " << n << ", " <<
				n / v << ", " <<
				v << endl;
	}
	s->doDelete( model );
	cout << "." << flush;
}

void testCylDiffnWithStoich()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 25e-6;
	double r0 = 1e-6;
	double r1 = 1e-6;
	double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
	double runtime = 10.0;
	double dt0 = 0.1; // Used for diffusion. 0.2 is the highest dt for which the error is in bounds
	double dt1 = 1; // Used for chem.
	double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
	Field< double >::set( cyl, "r0", r0 );
	Field< double >::set( cyl, "r1", r1 );
	Field< double >::set( cyl, "x0", 0 );
	Field< double >::set( cyl, "x1", len );
	Field< double >::set( cyl, "diffLength", diffLength );
	unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
	assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
	Id pool1 = s->doCreate( "Pool", cyl, "pool1", 1 );
	Id pool2 = s->doCreate( "Pool", cyl, "pool2", 1 );
	Field< double >::set( pool1, "diffConst", diffConst );
	Field< double >::set( pool2, "diffConst", diffConst/2 );

	Id stoich = s->doCreate( "Stoich", model, "stoich", 1 );
	Id ksolve = s->doCreate( "Ksolve", model, "ksolve", 1 );
	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( stoich, "compartment", cyl );
	Field< Id >::set( stoich, "ksolve", ksolve );
	Field< Id >::set( stoich, "dsolve", dsolve );
	Field< string >::set( stoich, "path", "/model/cyl/#" );
	assert( pool1.element()->numData() == ndc );

	// Then find a way to test it.
	vector< double > poolVec;
	Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
	Field< double >::set( ObjId( pool2, 0 ), "nInit", 1.0 );
   	Field< double >::getVec( pool1, "nInit", poolVec );
	assert( poolVec.size() == ndc );
	assert( doubleEq( poolVec[0], 1.0 ) );
	assert( doubleEq( poolVec[1], 0.0 ) );

	vector< double > nvec = 
		LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( nvec.size() == ndc );

	// Next: build by doing reinit
	s->doUseClock( "/model/dsolve", "process", 0 );
	s->doUseClock( "/model/ksolve", "process", 1 );
	s->doSetClock( 0, dt0 );
	s->doSetClock( 1, dt1 );
	s->doReinit();
	s->doStart( runtime );

	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
   	Field< double >::getVec( pool1, "n", poolVec );
	assert( nvec.size() == poolVec.size() );
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		assert( doubleEq( nvec[i], poolVec[i] ) );
	/*
	cout << endl;
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		cout << nvec[i] << "	";
	cout << endl;
	*/

	double dx = diffLength;
	double err = 0.0;
	double analyticTot = 0.0;
	double myTot = 0.0;
	for ( unsigned int i = 0; i < nvec.size(); ++i ) {
		double x = i * dx + dx * 0.5;
		// This part is the solution as a func of x,t.
		double y = dx *  // This part represents the init n of 1 in dx
			( 1.0 / sqrt( PI * diffConst * runtime ) ) * 
			exp( -x * x / ( 4 * diffConst * runtime ) ); 
		err += ( y - nvec[i] ) * ( y - nvec[i] );
		//cout << i << "	" << x << "	" << y << "	" << conc[i] << endl;
		analyticTot += y;
		myTot += nvec[i];
	} 
	assert( doubleEq( myTot, 1.0 ) );
	// cout << "analyticTot= " << analyticTot << ", myTot= " << myTot << endl;
	assert( err < 1.0e-5 );


	s->doDelete( model );
	cout << "." << flush;
}

void testCellDiffn()
{
	Id makeCompt( Id parentCompt, Id parentObj,
		string name, double len, double dia, double theta );
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 40e-6;
	double dia = 10e-6;
	double diffLength = 1e-6;
	double dt = 1.0e-1;
	double runtime = 100.0;
	double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id soma = makeCompt( Id(), model, "soma", dia, dia, 90 );
	Id dend = makeCompt( soma, model, "dend", len, 3e-6, 0 );
	Id branch1 = makeCompt( dend, model, "branch1", len, 2e-6, 45.0 );
	Id branch2 = makeCompt( dend, model, "branch2", len, 2e-6, -45.0 );
	Id twig1 = makeCompt( branch1, model, "twig1", len, 1.5e-6, 90.0 );
	Id twig2 = makeCompt( branch1, model, "twig2", len, 1.5e-6, 0.0 );

	Id nm = s->doCreate( "NeuroMesh", model, "neuromesh", 1 );
	Field< double >::set( nm, "diffLength", diffLength );
	Field< string >::set( nm, "geometryPolicy", "cylinder" );
	Field< Id >::set( nm, "cell", model );
	unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
	assert( ns == 6 );
	unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
	assert( ndc == 210  );
	Id pool1 = s->doCreate( "Pool", nm, "pool1", 1 );
	Field< double >::set( pool1, "diffConst", diffConst );
	Id pool2 = s->doCreate( "Pool", nm, "pool2", 1 );
	Field< double >::set( pool2, "diffConst", diffConst );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", nm );
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	// Next: build by setting path
	Field< string >::set( dsolve, "path", "/model/neuromesh/pool#" );

	vector< double > nvec = 
		LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( nvec.size() == ndc );
	assert( pool1.element()->numData() == ndc );
	Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
	Field< double >::set( ObjId( pool2, ndc - 1 ), "nInit", 2.0 );

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

	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	vector< double > pool1Vec;
	Field< double >::getVec( pool1, "n", pool1Vec );
	assert( pool1Vec.size() == ndc );

	vector< double > pool2Vec;
	Field< double >::getVec( pool2, "n", pool2Vec );
	assert( pool2Vec.size() == ndc );
	double myTot1 = 0;
	double myTot2 = 0;
	for ( unsigned int i = 0; i < nvec.size(); ++i ) {
		assert( doubleEq( pool1Vec[i], nvec[i] ) );
		myTot1 += nvec[i];
		myTot2 += pool2Vec[i];
	}
	assert( doubleEq( myTot1, 1.0 ) );
	assert( doubleEq( myTot2, 2.0 ) );

	/*
	cout << endl;
	cout << "Big cell: " << endl;
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		cout << nvec[i] << ", " << pool2Vec[i] << endl;
	cout << endl;
	*/


	s->doDelete( model );
	cout << "." << flush;
}

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;
}

/**
 * In all cases we set up the same amount of data transfer by the msgs, that
 * is, equivalent to a fully recurrently connected network.
 * Used in regressionTests/benchmarkTests.cpp
 */
void benchmarkMsg( unsigned int n, string msgType )
{
	Eref sheller = Id().eref();
	Shell* shell = reinterpret_cast< Shell* >( sheller.data() );
	vector< double > init( n );
	for ( unsigned int i = 0; i < n; ++i )
		init[i] = (i + 1) * 1e6;

	Id a1 = shell->doCreate( "Arith", Id(), "a1", n );

	if ( msgType == "Single" ) {
		for ( unsigned int i = 0; i < n; ++i ) {
			for ( unsigned int j = 0; j < n; ++j ) {
				ObjId m1 = shell->doAddMsg( "Single", 
					ObjId( a1, i ), "output", ObjId( a1, j ), "arg3" );
				assert( !m1.bad() );
			}
		}
	} else if ( msgType == "OneToAll" ) {
		for ( unsigned int i = 0; i < n; ++i ) {
			ObjId m1 = shell->doAddMsg( "OneToAll", 
				ObjId( a1, i ), "output", ObjId( a1, 0 ), "arg3" );
			assert( !m1.bad() );
		}
	} else if ( msgType == "OneToOne" ) {
		for ( unsigned int i = 0; i < n; ++i ) { // just repeat it n times
			ObjId m1 = shell->doAddMsg( "OneToOne", 
				ObjId( a1, 0 ), "output", ObjId( a1, 0 ), "arg3" );
			assert( !m1.bad() );
		}
	} else if ( msgType == "Diagonal" ) {
		for ( unsigned int i = 0; i < 2 * n; ++i ) { // Set up all offsets
			ObjId m1 = shell->doAddMsg( "Diagonal", 
				ObjId( a1, 0 ), "output", ObjId( a1, 0 ), "arg3" );
			Field< int >::set( m1, "stride", n - i );
		}
	} else if ( msgType == "Sparse" ) {
		ObjId m1 = shell->doAddMsg( "Sparse", 
			ObjId( a1, 0 ), "output", ObjId( a1, 0 ), "arg3" );
	
		SetGet2< double, long >::set( m1, 
			"setRandomConnectivity", 1.0, 1234 );
	} 

	shell->doUseClock( "/a1", "proc", 0 );
	for ( unsigned int i = 0; i < 10; ++i )
		shell->doSetClock( i, 0 );
	shell->doSetClock( 0, 1 );
	shell->doReinit();
	SetGet1< double >::setVec( a1, "arg1", init );
	shell->doStart( 100 );
	for ( unsigned int i = 0; i < n; ++i )
		init[i] = 0; // be sure we don't retain old info.
	init.clear();
	Field< double >::getVec( a1, "outputValue", init );
	cout << endl;
	for ( unsigned int i = 0; i < n; ++i ) {
		cout << i << " " << init[i] << "	";
		if ( i % 5 == 4 )
			cout << endl;
	}

	shell->doDelete( a1 );
}

//.........这里部分代码省略.........
                                    // 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;
}