diff --git a/moose-core/.travis.yml b/moose-core/.travis.yml
index b5f3f688e66c00d97989013b897183b478ca3c70..a89ee791337440c7828cccb68c5424a341c5a962 100644
--- a/moose-core/.travis.yml
+++ b/moose-core/.travis.yml
@@ -15,6 +15,10 @@ notifications:
on_success: change
on_failure: always
+addons:
+ apt:
+ update: true
+
before_script:
- echo "OSX related"
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then nvm get head || true; fi
diff --git a/moose-core/.travis/travis_prepare_osx.sh b/moose-core/.travis/travis_prepare_osx.sh
index 65bff743617d0883dddf034c8cc7b7e17fcae747..d8aa29a1719ebd80e12a0e5b2b6920313cf1cceb 100755
--- a/moose-core/.travis/travis_prepare_osx.sh
+++ b/moose-core/.travis/travis_prepare_osx.sh
@@ -39,6 +39,6 @@ echo 'import sys; sys.path.insert(1, "/usr/local/lib/python2.7/site-packages")'
PATH=/usr/local/bin:/usr/bin:$PATH
# ensurepip
python -m ensurepip
-python -m pip install matplotlib --user
+python -m pip install matplotlib --user --upgrade
python -m pip install pyNeuroML libNeuroML --user
python -m pip install scipy --user
diff --git a/moose-core/biophysics/MarkovSolver.cpp b/moose-core/biophysics/MarkovSolver.cpp
index 1ab4adc1b6646109457d759b30b10d283aa895df..b5e3ddbba288cd3ecfd8ada9b2f40601d444dd1c 100644
--- a/moose-core/biophysics/MarkovSolver.cpp
+++ b/moose-core/biophysics/MarkovSolver.cpp
@@ -81,12 +81,12 @@ Matrix* MarkovSolver::computePadeApproximant( Matrix* Q1,
unsigned int degreeIndex )
{
Matrix *expQ;
- Matrix *U, *V, *VplusU, *VminusU, *invVminusU, *Qpower;
+ Matrix *U, *VplusU, *VminusU, *invVminusU, *Qpower;
vector< unsigned int >* swaps = new vector< unsigned int >;
unsigned int n = Q1->size();
unsigned int degree = mCandidates[degreeIndex];
double *padeCoeffs;
-
+ Matrix *V = matAlloc(n);
//Vector of Matrix pointers. Each entry is an even power of Q.
vector< Matrix* > QevenPowers;
@@ -131,7 +131,6 @@ Matrix* MarkovSolver::computePadeApproximant( Matrix* Q1,
case 7 :
case 9 :
U = matAlloc( n );
- V = matAlloc( n );
QevenPowers.push_back( Q1 );
@@ -199,6 +198,7 @@ Matrix* MarkovSolver::computePadeApproximant( Matrix* Q1,
matMatAdd( temp, Q6, 1.0, b13[6], FIRST );
matMatAdd( temp, Q4, 1.0, b13[4], FIRST );
matMatAdd( temp, Q2, 1.0, b13[2], FIRST );
+ delete( V );
V = matEyeAdd( temp, b13[0] );
delete temp;
diff --git a/moose-core/biophysics/MatrixOps.cpp b/moose-core/biophysics/MatrixOps.cpp
index a983ef748304b857df89b969236257e5258d8c89..f0f478c3540243630216e0efe34a3c5f65bd24ef 100644
--- a/moose-core/biophysics/MatrixOps.cpp
+++ b/moose-core/biophysics/MatrixOps.cpp
@@ -128,12 +128,14 @@ void matMatAdd( Matrix* A, Matrix* B, double alpha, double beta,
Matrix *C;
unsigned int n = A->size();
- if ( resIndex == FIRST )
+ if ( resIndex == FIRST ) {
C = A;
- else if ( resIndex == SECOND )
+ } else if ( resIndex == SECOND ) {
C = B;
- else
+ } else {
cerr << "matMatAdd : Invalid index supplied to store result.\n";
+ return;
+ }
for( unsigned int i = 0; i < n; ++i )
{
@@ -481,7 +483,7 @@ Matrix* matAlloc( unsigned int n )
A->resize( n );
for ( unsigned int i = 0; i < n; ++i )
- (*A)[i].resize( n );
+ (*A)[i].resize( n, 0.0 );
return A;
}
@@ -490,7 +492,7 @@ Vector* vecAlloc( unsigned int n )
{
Vector* vec = new Vector;
- vec->resize( n );
+ vec->resize( n, 0.0 );
return vec;
}
diff --git a/moose-core/biophysics/RandSpike.cpp b/moose-core/biophysics/RandSpike.cpp
index 0c2e09ffff1825cd06407a8931fd79dd41db58d7..4678970f2f7836f99f3fb6e7bdb2710e48b3a9a4 100644
--- a/moose-core/biophysics/RandSpike.cpp
+++ b/moose-core/biophysics/RandSpike.cpp
@@ -67,6 +67,14 @@ const Cinfo* RandSpike::initCinfo()
&RandSpike::setRefractT,
&RandSpike::getRefractT
);
+ static ValueFinfo< RandSpike, bool > doPeriodic( "doPeriodic",
+ "Flag: when false, do Poisson process with specified mean rate.\n"
+ "When true, fire periodically at specified rate.\n"
+ "Defaults to false. Note that refractory time overrides this: "
+ "Rate cannot exceed 1/refractT.",
+ &RandSpike::setDoPeriodic,
+ &RandSpike::getDoPeriodic
+ );
static ReadOnlyValueFinfo< RandSpike, bool > hasFired( "hasFired",
"True if RandSpike has just fired",
&RandSpike::getFired
@@ -79,6 +87,7 @@ const Cinfo* RandSpike::initCinfo()
&rate, // Value
&refractT, // Value
&absRefract, // Value
+ &doPeriodic, // Value
&hasFired, // ReadOnlyValue
};
@@ -86,7 +95,8 @@ const Cinfo* RandSpike::initCinfo()
{
"Name", "RandSpike",
"Author", "Upi Bhalla",
- "Description", "RandSpike object, generates random spikes at."
+ "Description", "RandSpike object, generates random or regular "
+ "spikes at "
"specified mean rate. Based closely on GENESIS randspike. "
};
static Dinfo< RandSpike > dinfo;
@@ -111,7 +121,8 @@ RandSpike::RandSpike()
refractT_(0.0),
lastEvent_(0.0),
threshold_(0.0),
- fired_( 0 )
+ fired_( false ),
+ doPeriodic_( false )
{;}
//////////////////////////////////////////////////////////////////
@@ -121,6 +132,10 @@ RandSpike::RandSpike()
// Value Field access function definitions.
void RandSpike::setRate( double rate )
{
+ if ( rate < 0.0 ) {
+ cout <<"Warning: RandSpike::setRate: Rate must be >= 0. Using 0.\n";
+ rate = 0.0;
+ }
rate_ = rate;
double prob = 1.0 - rate * refractT_;
if ( prob <= 0.0 ) {
@@ -149,6 +164,15 @@ bool RandSpike::getFired() const
return fired_;
}
+void RandSpike::setDoPeriodic( bool val )
+{
+ doPeriodic_ = val;
+}
+bool RandSpike::getDoPeriodic() const
+{
+ return doPeriodic_;
+}
+
//////////////////////////////////////////////////////////////////
// RandSpike::Dest function definitions.
@@ -156,16 +180,23 @@ bool RandSpike::getFired() const
void RandSpike::process( const Eref& e, ProcPtr p )
{
- if ( refractT_ > p->currTime - lastEvent_ )
+ if ( refractT_ > p->currTime - lastEvent_ || rate_ <= 0.0 )
return;
- double prob = realRate_ * p->dt;
- if ( prob >= 1.0 || prob >= mtrand() )
- {
- lastEvent_ = p->currTime;
- spikeOut()->send( e, p->currTime );
- fired_ = true;
+ fired_ = false;
+ if (doPeriodic_) {
+ if ( (p->currTime - lastEvent_) > 1.0/rate_ ) {
+ lastEvent_ = p->currTime;
+ spikeOut()->send( e, p->currTime );
+ fired_ = true;
+ }
} else {
- fired_ = false;
+ double prob = realRate_ * p->dt;
+ if ( prob >= 1.0 || prob >= mtrand() )
+ {
+ lastEvent_ = p->currTime;
+ spikeOut()->send( e, p->currTime );
+ fired_ = true;
+ }
}
}
diff --git a/moose-core/biophysics/RandSpike.h b/moose-core/biophysics/RandSpike.h
index ea0cf595d9d2097351b832092a971149628d2fb6..2a7da3de72c1ead13a743483faf5bbc0496c331a 100644
--- a/moose-core/biophysics/RandSpike.h
+++ b/moose-core/biophysics/RandSpike.h
@@ -24,6 +24,9 @@ class RandSpike
void setRefractT( double val );
double getRefractT() const;
+ void setDoPeriodic( bool val );
+ bool getDoPeriodic() const;
+
bool getFired() const;
//////////////////////////////////////////////////////////////////
@@ -42,6 +45,7 @@ class RandSpike
double lastEvent_;
double threshold_;
bool fired_;
+ bool doPeriodic_;
};
#endif // _RANDSPIKE_H
diff --git a/moose-core/builtins/Function.cpp b/moose-core/builtins/Function.cpp
index 79413c1bb0fea17bee2ebaff277423c278421456..ea0c8ae872cc043db64dd0abe7c9e905559316e8 100644
--- a/moose-core/builtins/Function.cpp
+++ b/moose-core/builtins/Function.cpp
@@ -125,9 +125,19 @@ const Cinfo * Function::initCinfo()
"When *false*, disables event-driven calculation and turns on "
"Process-driven calculations. \n"
"When *true*, enables event-driven calculation and turns off "
- "Process-driven calculations. \n",
+ "Process-driven calculations. \n"
+ "Defaults to *false*. \n",
&Function::setUseTrigger,
&Function::getUseTrigger);
+ static ValueFinfo< Function, bool > doEvalAtReinit(
+ "doEvalAtReinit",
+ "When *false*, disables function evaluation at reinit, and "
+ "just emits a value of zero to any message targets. \n"
+ "When *true*, does a function evaluation at reinit and sends "
+ "the computed value to any message targets. \n"
+ "Defaults to *false*. \n",
+ &Function::setDoEvalAtReinit,
+ &Function::getDoEvalAtReinit);
static ElementValueFinfo< Function, string > expr(
"expr",
"Mathematical expression defining the function. The underlying parser\n"
@@ -257,6 +267,8 @@ const Cinfo * Function::initCinfo()
&rate,
&derivative,
&mode,
+ &useTrigger,
+ &doEvalAtReinit,
&expr,
&numVars,
&inputs,
@@ -316,7 +328,7 @@ static const Cinfo * functionCinfo = Function::initCinfo();
Function::Function(): _t(0.0), _valid(false), _numVar(0), _lastValue(0.0),
_value(0.0), _rate(0.0), _mode(1),
- _useTrigger( false ), _stoich(0)
+ _useTrigger( false ), _doEvalAtReinit( false ), _stoich(0)
{
_parser.SetVarFactory(_functionAddVar, this);
_independent = "x0";
@@ -601,6 +613,16 @@ bool Function::getUseTrigger() const
return _useTrigger;
}
+void Function::setDoEvalAtReinit(bool doEvalAtReinit )
+{
+ _doEvalAtReinit = doEvalAtReinit;
+}
+
+bool Function::getDoEvalAtReinit() const
+{
+ return _doEvalAtReinit;
+}
+
double Function::getValue() const
{
double value = 0.0;
@@ -770,8 +792,11 @@ void Function::reinit(const Eref &e, ProcPtr p)
_valid = false;
}
_t = p->currTime;
- _value = 0.0;
- _lastValue = 0.0;
+ if (_doEvalAtReinit) {
+ _lastValue = _value = getValue();
+ } else {
+ _lastValue = _value = 0.0;
+ }
_rate = 0.0;
switch (_mode){
case 1: {
diff --git a/moose-core/builtins/Function.h b/moose-core/builtins/Function.h
index ed6c4517b65f230ebf73f43544f1d6e681d0da97..83f2b90381112ba44fc84603ba8da9af5baf69eb 100644
--- a/moose-core/builtins/Function.h
+++ b/moose-core/builtins/Function.h
@@ -92,6 +92,10 @@ class Function
void setUseTrigger(bool useTrigger);
bool getUseTrigger() const;
+ // set/get flag to do function evaluation at reinit
+ void setDoEvalAtReinit(bool doEvalAtReinit);
+ bool getDoEvalAtReinit() const;
+
void setNumVar(unsigned int num);
unsigned int getNumVar() const;
@@ -143,6 +147,7 @@ protected:
double _rate;
unsigned int _mode;
bool _useTrigger;
+ bool _doEvalAtReinit;
// this stores variables received via incoming messages, identifiers of the form x{i} are included in this
vector<Variable *> _varbuf;
// this stores variable values pulled by sending request. identifiers of the form y{i} are included in this
diff --git a/moose-core/builtins/HDF5WriterBase.cpp b/moose-core/builtins/HDF5WriterBase.cpp
index 80c47cbe0dcdb371170ed11a4b9db4b1412c3f46..747a2d8dd76bd134ef37fde2380ab5bed7bcc88e 100644
--- a/moose-core/builtins/HDF5WriterBase.cpp
+++ b/moose-core/builtins/HDF5WriterBase.cpp
@@ -84,7 +84,8 @@ hid_t require_group(hid_t file, string path)
{
vector<string> pathTokens;
moose::tokenize(path, "/", pathTokens);
- hid_t prev = file, current;
+ hid_t prev = file;
+ hid_t current = -1;
htri_t exists;
// Open the container for the event maps
for (unsigned int ii = 0; ii < pathTokens.size(); ++ii){
diff --git a/moose-core/diffusion/Dsolve.cpp b/moose-core/diffusion/Dsolve.cpp
index 61ffb1bc612c01dde5f917fd15faabd54169e0be..d2f27122fe4eb2a9777200d76468671fa4e99cb8 100644
--- a/moose-core/diffusion/Dsolve.cpp
+++ b/moose-core/diffusion/Dsolve.cpp
@@ -867,9 +867,11 @@ void Dsolve::mapXfersBetweenDsolves(
Id pool( srcSolve->pools_[i].getId() );
assert( pool != Id() );
string poolName = pool.element()->getName();
- size_t prefixLen = poolName.length() - xlen;
- if ( poolName.rfind( xferPost ) == prefixLen )
- srcMap[ poolName.substr( 0, prefixLen) ] = i;
+ if ( poolName.length() > xlen ) {
+ size_t prefixLen = poolName.length() - xlen;
+ if ( poolName.rfind( xferPost ) == prefixLen )
+ srcMap[ poolName.substr( 0, prefixLen) ] = i;
+ }
}
const Dsolve* destSolve = reinterpret_cast< const Dsolve* >(
@@ -889,7 +891,8 @@ void Dsolve::mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other)
{
Dsolve* otherSolve = reinterpret_cast< Dsolve* >(
other.eref().data() );
- vector< ConcChanInfo >& ch = channels_;
+ Dsolve* selfSolve = reinterpret_cast< Dsolve* >( self.eref().data() );
+ vector< ConcChanInfo >& ch = selfSolve->channels_;
unsigned int outIndex;
for ( unsigned int i = 0; i < ch.size(); ++i ) {
outIndex = otherSolve->convertIdToPoolIndex( ch[i].otherPool );
@@ -901,7 +904,7 @@ void Dsolve::mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other)
// Now set up the other Dsolve.
vector< ConcChanInfo >& ch2 = otherSolve->channels_;
for ( unsigned int i = 0; i < ch2.size(); ++i ) {
- outIndex = convertIdToPoolIndex( ch2[i].otherPool );
+ outIndex = selfSolve->convertIdToPoolIndex( ch2[i].otherPool );
if ( outIndex != ~0U ) {
jn.otherChannels.push_back(i);
ch2[i].otherPool = outIndex; // replace the Id with the index
@@ -933,6 +936,7 @@ void Dsolve::innerBuildMeshJunctions( Id self, Id other, bool selfIsMembraneBoun
{
DiffJunction jn; // This is based on the Spine Dsolver.
jn.otherDsolve = other.value();
+ Dsolve* dself = reinterpret_cast< Dsolve* >( self.eref().data() );
if ( selfIsMembraneBound ) {
mapChansBetweenDsolves( jn, self, other );
} else {
@@ -943,8 +947,9 @@ void Dsolve::innerBuildMeshJunctions( Id self, Id other, bool selfIsMembraneBoun
mapVoxelsBetweenMeshes( jn, self, other );
+
// printJunction( self, other, jn );
- junctions_.push_back( jn );
+ dself->junctions_.push_back( jn );
}
/////////////////////////////////////////////////////////////
diff --git a/moose-core/diffusion/Dsolve.h b/moose-core/diffusion/Dsolve.h
index d6028e62a4c7fb8bdbbe267fa87f54f8174f2c16..b8197ef4bb300255584c63346db7db02fb1591e6 100644
--- a/moose-core/diffusion/Dsolve.h
+++ b/moose-core/diffusion/Dsolve.h
@@ -96,7 +96,7 @@ class Dsolve: public ZombiePoolInterface
* the junction between any specified pair of Dsolves.
* Note that it builds the junction on the 'self' Dsolve.
*/
- void innerBuildMeshJunctions( Id self, Id other,
+ static void innerBuildMeshJunctions( Id self, Id other,
bool isMembraneBound );
/// Sets up map of matching pools for diffusion.
@@ -107,7 +107,8 @@ class Dsolve: public ZombiePoolInterface
vector< unsigned int >& srcPools,
vector< unsigned int >& destPools,
Id src, Id dest );
- void mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other);
+ static void mapChansBetweenDsolves( DiffJunction& jn,
+ Id self, Id other);
/**
* Computes flux through a junction between diffusion solvers.
diff --git a/moose-core/kinetics/Enz.cpp b/moose-core/kinetics/Enz.cpp
index 2362c0b32eb83ce2247021aa762da822b0163776..f41d109ad984d7cde9275d8fea19338147e2d1ee 100644
--- a/moose-core/kinetics/Enz.cpp
+++ b/moose-core/kinetics/Enz.cpp
@@ -136,7 +136,7 @@ void Enz::vRemesh( const Eref& e )
void Enz::vSetK1( const Eref& e, double v )
{
r1_ = k1_ = v;
- double volScale = convertConcToNumRateUsingMesh( e, subOut, 1 );
+ double volScale = convertConcToNumRateUsingMesh( e, enzOut, 1 );
Km_ = ( k2_ + k3_ ) / ( k1_ * volScale );
}
@@ -160,7 +160,11 @@ double Enz::vGetK2( const Eref& e ) const
void Enz::vSetKcat( const Eref& e, double v )
{
- double ratio = k2_ / k3_;
+ double ratio = 4.0;
+ if ( v < EPSILON )
+ v = EPSILON;
+ if (k3_ > EPSILON)
+ ratio = k2_ / k3_;
k3_ = v;
k2_ = v * ratio;
vSetKm( e, Km_ ); // Update k1_ here as well.
@@ -182,7 +186,7 @@ void Enz::vSetKm( const Eref& e, double v )
{
Km_ = v;
double volScale =
- convertConcToNumRateUsingMesh( e, subOut, 1 );
+ convertConcToNumRateUsingMesh( e, enzOut, 1 );
k1_ = ( k2_ + k3_ ) / ( v * volScale );
}
@@ -193,14 +197,14 @@ double Enz::vGetKm( const Eref& e ) const
void Enz::vSetNumKm( const Eref& e, double v )
{
- double volScale = convertConcToNumRateUsingMesh( e, subOut, 1 );
+ double volScale = convertConcToNumRateUsingMesh( e, enzOut, 1 );
k1_ = ( k2_ + k3_ ) / v;
Km_ = v / volScale;
}
double Enz::vGetNumKm( const Eref& e ) const
{
- double volScale = convertConcToNumRateUsingMesh( e, subOut, 1 );
+ double volScale = convertConcToNumRateUsingMesh( e, enzOut, 1 );
return Km_ * volScale;
}
@@ -208,7 +212,7 @@ void Enz::vSetRatio( const Eref& e, double v )
{
k2_ = v * k3_;
double volScale =
- convertConcToNumRateUsingMesh( e, subOut, 1 );
+ convertConcToNumRateUsingMesh( e, enzOut, 1 );
k1_ = ( k2_ + k3_ ) / ( Km_ * volScale );
}
@@ -224,7 +228,7 @@ void Enz::vSetConcK1( const Eref& e, double v )
cout << "Enz::vSetConcK1: Warning: value " << v << " too small\n";
return;
}
- double volScale = convertConcToNumRateUsingMesh( e, subOut, 1 );
+ double volScale = convertConcToNumRateUsingMesh( e, enzOut, 1 );
r1_ = k1_ = v * volScale;
Km_ = ( k2_ + k3_ ) / ( v );
}
diff --git a/moose-core/kinetics/MMenz.cpp b/moose-core/kinetics/MMenz.cpp
index 8b396ceb107d74c824ffda06067c97ce3d6ec28c..cd03b00fa5655e1c1a32c95ce15820fdfdd433d4 100644
--- a/moose-core/kinetics/MMenz.cpp
+++ b/moose-core/kinetics/MMenz.cpp
@@ -133,6 +133,8 @@ double MMenz::vGetNumKm( const Eref& enz ) const
void MMenz::vSetKcat( const Eref& e, double v )
{
+ if ( v < EPSILON )
+ v = EPSILON;
kcat_ = v;
}
diff --git a/moose-core/kinetics/ReadKkit.cpp b/moose-core/kinetics/ReadKkit.cpp
index d3b5bd4154067e43753f4a4e60d1453da2bb7b4d..8172ae90920e225ab8337be73316c1bd0cdca1a8 100644
--- a/moose-core/kinetics/ReadKkit.cpp
+++ b/moose-core/kinetics/ReadKkit.cpp
@@ -1,3 +1,4 @@
+
/**********************************************************************
** This program is part of 'MOOSE', the
** Messaging Object Oriented Simulation Environment.
@@ -162,15 +163,17 @@ static void positionCompt( ObjId compt, double side, bool shiftUp )
}
}
-void makeSolverOnCompt( Shell* s, const vector< ObjId >& compts,
+void makeSolverOnCompt( Shell* s,const vector< ObjId >& compts,
bool isGsolve )
{
+
if ( compts.size() > 3 ) {
cout << "Warning: ReadKkit::makeSolverOnCompt: Cannot handle " <<
compts.size() << " chemical compartments\n";
return;
}
vector< Id > stoichVec;
+ /*
if ( compts.size() == 2 ) {
double side = Field< double >::get( compts[1], "dy" );
positionCompt( compts[0], side, true );
@@ -180,7 +183,8 @@ void makeSolverOnCompt( Shell* s, const vector< ObjId >& compts,
positionCompt( compts[0], side, true );
positionCompt( compts[2], side, false );
}
-
+ */
+ /*
for ( vector< ObjId >::const_iterator
i = compts.begin(); i != compts.end(); ++i ) {
string simpath = i->path() + "/##";
@@ -189,12 +193,14 @@ void makeSolverOnCompt( Shell* s, const vector< ObjId >& compts,
ksolve = s->doCreate( "Gsolve", *i, "gsolve", 1 );
else
ksolve = s->doCreate( "Ksolve", *i, "ksolve", 1 );
+ dsolve = s->doCreate("Dsolve,")
Id stoich = s->doCreate( "Stoich", *i, "stoich", 1 );
stoichVec.push_back( stoich );
Field< Id >::set( stoich, "compartment", *i );
Field< Id >::set( stoich, "ksolve", ksolve );
Field< string >::set( stoich, "path", simpath );
}
+ */
/* Not needed now that we use xfer pools to cross compartments.
if ( stoichVec.size() == 2 ) {
SetGet1< Id >::set( stoichVec[1], "buildXreacs", stoichVec[0] );
@@ -221,18 +227,25 @@ void setMethod( Shell* s, Id mgr, double simdt, double plotdt,
assert( compt != Id() );
string simpath2 = mgr.path() + "/##[ISA=StimulusTable]," +
mgr.path() + "/##[ISA=PulseGen]";
-
string m = lower( method );
+
+ if ( m == "ksolve" || m =="gsl" || m == "gssa" || m == "gsolve" ||
+ m == "gillespie" || m == "stochastic" )
+ {
+ cout << " Warning: Default solver set is Exponential Euler. To set \'gsl\' or \'gssa\' solver use function mooseaddChemSolver(modelpath,\'solverType\')"<<"\n";
+ }
+ /*
+
if ( m == "rk4" ) {
cout << "Warning, not yet implemented. Using rk5 instead\n";
m = "rk5";
}
if ( m == "ksolve" || m == "gsl" ||
m == "rk5" || m == "rkf" || m == "rk" ) {
- makeSolverOnCompt( s, ret, false );
+ makeSolverOnCompt( s, mgr, ret, false );
} else if ( m == "gssa" || m == "gsolve" ||
m == "gillespie" || m == "stochastic" ) {
- makeSolverOnCompt( s, ret, true );
+ makeSolverOnCompt( s, mgr,ret, true );
} else if ( m == "ee" || m == "neutral" ) {
// s->doUseClock( simpath, "process", 4 );
// s->doSetClock( 4, simdt );
@@ -242,7 +255,9 @@ void setMethod( Shell* s, Id mgr, double simdt, double plotdt,
// s->doUseClock( simpath, "process", 4 );
// s->doSetClock( 4, simdt );
}
+ */
s->doUseClock( simpath2, "proc", 11 );
+ s->doSetClock( 10, simdt );
s->doSetClock( 11, simdt );
s->doSetClock( 12, simdt );
s->doSetClock( 13, simdt );
diff --git a/moose-core/kinetics/WriteKkit.cpp b/moose-core/kinetics/WriteKkit.cpp
index 0f8533c05c0df96ed788913205ee10f0c4eac097..c32b1fa024b8b7f0a03fecc6da05ff5f5855607c 100644
--- a/moose-core/kinetics/WriteKkit.cpp
+++ b/moose-core/kinetics/WriteKkit.cpp
@@ -278,10 +278,12 @@ void writePlot( ofstream& fout, Id id,
{
string path = id.path();
size_t pos = path.find( "/graphs" );
- if ( pos == string::npos )
+ if ( pos == string::npos ) {
pos = path.find( "/moregraphs" );
- if ( pos == string::npos )
+ if ( pos == string::npos ) {
return;
+ }
+ }
path = path.substr( pos );
fout << "simundump xplot " << path << " 3 524288 \\\n" <<
"\"delete_plot.w <s> <d>; edit_plot.D <w>\" " << textcolour << " 0 0 1\n";
diff --git a/moose-core/ksolve/ZombieEnz.cpp b/moose-core/ksolve/ZombieEnz.cpp
index 361ac111f64eefa17f2bfda8c5eb68c06304f747..e73c27acc30fc776603ee8ea71fb61a19499494f 100644
--- a/moose-core/ksolve/ZombieEnz.cpp
+++ b/moose-core/ksolve/ZombieEnz.cpp
@@ -23,6 +23,8 @@
#include "CplxEnzBase.h"
#include "ZombieEnz.h"
+#define EPSILON 1e-15
+
const Cinfo* ZombieEnz::initCinfo()
{
//////////////////////////////////////////////////////////////
@@ -118,8 +120,11 @@ void ZombieEnz::vSetKcat( const Eref& e, double v )
double k2 = getK2( e );
double k3 = getKcat( e );
double ratio = 4.0;
- if ( k3 > 1e-10 )
+ if ( v < EPSILON )
+ v = EPSILON;
+ if ( k3 > EPSILON ) {
ratio = k2/k3;
+ }
double Km = (k2 + k3) / concK1_;
concK1_ = v * (1.0 + ratio) / Km;
diff --git a/moose-core/ksolve/ZombieMMenz.cpp b/moose-core/ksolve/ZombieMMenz.cpp
index 95f163a0f7db6115b543bc6b05c0e06bf2d7fc8d..dc400897603a9fcf9414690b40df04234846d0e4 100644
--- a/moose-core/ksolve/ZombieMMenz.cpp
+++ b/moose-core/ksolve/ZombieMMenz.cpp
@@ -22,6 +22,8 @@
#include "EnzBase.h"
#include "ZombieMMenz.h"
+#define EPSILON 1e-15
+
const Cinfo* ZombieMMenz::initCinfo()
{
//////////////////////////////////////////////////////////////
@@ -121,6 +123,8 @@ double ZombieMMenz::vGetNumKm( const Eref& e ) const
void ZombieMMenz::vSetKcat( const Eref& e, double v )
{
+ if ( v < EPSILON )
+ v = EPSILON;
stoich_->setMMenzKcat( e, v );
}
diff --git a/moose-core/python/rdesigneur/fixXreacs.py b/moose-core/python/fixXreacs.py
similarity index 99%
rename from moose-core/python/rdesigneur/fixXreacs.py
rename to moose-core/python/fixXreacs.py
index 19379f5f189b60a02a4dcc3e17899509e5ce3e9c..de7bec638be4fcdc85b03f7598fc8758a5706479 100644
--- a/moose-core/python/rdesigneur/fixXreacs.py
+++ b/moose-core/python/fixXreacs.py
@@ -181,7 +181,7 @@ def restoreXreacs( basepath ):
#print( "Deleting {}".format( i.parent.path ) )
#print msgs
moose.delete( i.parent )
- for j in msgs:
+ for j in msgs[1:]:
if len( j ) > 0:
args = j.split( ' ' )
assert( len( args ) == 4 )
diff --git a/moose-core/python/moose/SBML/readSBML.py b/moose-core/python/moose/SBML/readSBML.py
index 06de82bc3a7317fa1d8525be13e86fb1013d3026..d9b0e90f0b531fb446123c5409c3b46e6bfd8de7 100644
--- a/moose-core/python/moose/SBML/readSBML.py
+++ b/moose-core/python/moose/SBML/readSBML.py
@@ -13,14 +13,15 @@
** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS
Created : Thu May 13 10:19:00 2016(+0530)
Version
-Last-Updated: Sat Jan 6 1:21:00 2018(+0530)
+Last-Updated: Fri May 21 11:21:00 2018(+0530)
By:HarshaRani
**********************************************************************/
2018
-Jan6 : - only if valid model exists, then printing the no of compartment,pool,reaction etc
+May 18: - cleanedup and connected cplx pool to correct parent enzyme
+Jan 6: - only if valid model exists, then printing the no of compartment,pool,reaction etc
- at reaction level a check made to see if path exist while creating a new reaction
2017
-Oct 4 : - loadpath is cleaned up
+Oct 4: - loadpath is cleaned up
Sep 13: - After EnzymaticReaction's k2 is set, explicity ratio is set to 4 to make sure it balance.
- If units are defined in the rate law for the reaction then check is made and if not in milli mole the base unit
then converted to milli unit
@@ -31,11 +32,11 @@ Sep 12: - Now mooseReadSBML return model and errorFlag
- errorFlag is set for Rules (for now piecewise is set which is not read user are warned)
- rateLaw are also calculated depending on units and number of substrates/product
-Sep 8 : -functionDefinitions is read,
+Sep 8 : - functionDefinitions is read,
- if Kf and Kb unit are not defined then checked if substance units is defined and depending on this unit Kf and Kb is calculated
-kf and kb units is not defined and if substance units is also not defined validator fails
-Aug 9 : a check made to for textColor while adding to Annotator
-Aug 8 : removed "findCompartment" function to chemConnectUtil and imported the function from the same file
+Aug 9 : - a check made to for textColor while adding to Annotator
+Aug 8 : - removed "findCompartment" function to chemConnectUtil and imported the function from the same file
TODO in
-Compartment
@@ -251,7 +252,7 @@ def checkGroup(basePath,model):
groupInfo[p.getId()] =[mem.getIdRef()]
return groupInfo
def setupEnzymaticReaction(enz, groupName, enzName,
- specInfoMap, modelAnnotaInfo):
+ specInfoMap, modelAnnotaInfo,deletcplxMol):
enzPool = (modelAnnotaInfo[groupName]["enzyme"])
enzPool = str(idBeginWith(enzPool))
enzParent = specInfoMap[enzPool]["Mpath"]
@@ -259,13 +260,14 @@ def setupEnzymaticReaction(enz, groupName, enzName,
cplx = str(idBeginWith(cplx))
complx = moose.element(specInfoMap[cplx]["Mpath"].path)
enzyme_ = moose.Enz(enzParent.path + '/' + enzName)
- moose.move(complx, enzyme_)
- moose.connect(enzyme_, "cplx", complx, "reac")
+ complx1 = moose.Pool(enzyme_.path+'/'+moose.element(complx).name)
+ specInfoMap[cplx]["Mpath"] = complx1
+ moose.connect(enzyme_, "cplx", complx1, "reac")
moose.connect(enzyme_, "enz", enzParent, "reac")
-
sublist = (modelAnnotaInfo[groupName]["substrate"])
prdlist = (modelAnnotaInfo[groupName]["product"])
-
+ deletcplxMol.append(complx.path)
+ complx = complx1
for si in range(0, len(sublist)):
sl = sublist[si]
sl = str(idBeginWith(sl))
@@ -280,7 +282,6 @@ def setupEnzymaticReaction(enz, groupName, enzName,
if (enz.isSetNotes):
pullnotes(enz, enzyme_)
-
return enzyme_, True
@@ -547,7 +548,7 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue,fun
# ToDo
# -- I need to check here if any substance/product is if ( constant == true && bcondition == false)
# cout <<"The species "<< name << " should not appear in reactant or product as per sbml Rules"<< endl;
-
+ deletecplxMol = []
errorFlag = True
reactSBMLIdMooseId = {}
msg = ""
@@ -589,7 +590,7 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue,fun
if (groupName != "" and list(
modelAnnotaInfo[groupName]["stage"])[0] == 3):
reaction_, reactionCreated = setupEnzymaticReaction(
- reac, groupName, rName, specInfoMap, modelAnnotaInfo)
+ reac, groupName, rName, specInfoMap, modelAnnotaInfo,deletecplxMol)
reaction_.k3 = modelAnnotaInfo[groupName]['k3']
reaction_.concK1 = modelAnnotaInfo[groupName]['k1']
reaction_.k2 = modelAnnotaInfo[groupName]['k2']
@@ -710,6 +711,9 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue,fun
elif reaction_.className == "MMenz":
reaction_.kcat = kfvalue
reaction_.Km = kbvalue
+ for l in deletecplxMol:
+ if moose.exists(l):
+ moose.delete(moose.element(l))
return (errorFlag, msg)
@@ -1151,6 +1155,7 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap,
poolInfo.color = v
elif k == 'Color':
poolInfo.textColor = v
+
specInfoMap[sId] = {
"Mpath": poolId,
"const": constant,
diff --git a/moose-core/python/moose/SBML/writeSBML.py b/moose-core/python/moose/SBML/writeSBML.py
index 8491b46a8aba3bfb9b2359fddce84931da0c3e28..8124611b4802e3f1b27ff2af7115e45371ae174d 100644
--- a/moose-core/python/moose/SBML/writeSBML.py
+++ b/moose-core/python/moose/SBML/writeSBML.py
@@ -13,11 +13,12 @@
** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS
Created : Friday May 27 12:19:00 2016(+0530)
Version
-Last-Updated: Sat 6 Jan 01:10:00 2018(+0530)
+Last-Updated: Mon 30 Apr 15:10:00 2018(+0530)
By: HarshaRani
**********************************************************************/
/****************************
2018
+Apr 30: indentation corrected while writting annotation for enzymecomplex
Jan 6: for Product_formation_, k3 units depends on noofSub, prd was passed which is fixed
2017
Dec 15: If model path exists is checked
@@ -143,7 +144,7 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}):
if moose.exists(key.path+'/info'):
ginfo = moose.element(key.path+'/info')
else:
- moose.Annotator(key.path+'/info')
+ ginfo = moose.Annotator(key.path+'/info')
groupCompartment = findCompartment(key)
if ginfo.color != '':
grpAnno = "<moose:GroupAnnotation>"
@@ -231,7 +232,6 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
compt = comptVec.name + "_" + \
str(comptVec.getId().value) + "_" + \
str(comptVec.getDataIndex()) + "_"
-
#Writting out S+E -> SE*
enzsetId = str(idBeginWith(cleanEnzname +
"_" + str(enz.getId().value) +
@@ -254,7 +254,6 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
enzAnno = enzAnno + "<moose:enzyme>" + \
(str(idBeginWith(convertSpecialChar(
nameList_[i])))) + "</moose:enzyme>\n"
-
#Finding Substrate, (S)
enzSub = enz.neighbors["sub"]
if not enzSub:
@@ -277,7 +276,6 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
enzAnno = enzAnno + "<moose:product>" + \
nameList_[i] + "</moose:product>\n"
foundEnzymeComplex = True
-
if foundEnzymeComplex:
# Only if S+E->SE* found, reaction is created
enzyme = cremodel_.createReaction()
@@ -318,15 +316,16 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
printParameters(kl, "k1", k1, unit)
punit = parmUnit(noofPrd - 1, cremodel_)
printParameters(kl, "k2", k2, punit)
-
if groupName != moose.element('/'):
if groupName not in groupInfo:
- groupInfo[groupName]=[enzsetId]
+ #groupInfo[groupName]=[enzsetId]
+ groupInfo[groupName]=[nameList_[0]]
else:
- groupInfo[groupName].append(enzsetId)
+ #groupInfo[groupName].append(enzsetId)
+ groupInfo[groupName].append(nameList_[0])
else:
if secplxerror:
- print ("\'"+enz.name+"\' this enzyme is not written to file because,"+ secplxerror)
+ print ("\'"+enz.parent.name+ "/"+enz.name+"\' this enzyme is not written to file because,"+ secplxerror)
#Here SE* -> E+ P
foundEnzymeEP = False
@@ -335,6 +334,7 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
"_" + str(enz.getDataIndex()) +
"_" + "Product_formation_"))
cplxeperror = ""
+ enzAnno2 = ""
#enzSubt = ""
enzOut = enz.neighbors["enzOut"]
if not enzOut:
@@ -375,16 +375,15 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
for i in range(0, len(nameList_)):
enzAnno2 = enzAnno2 + "<moose:product>" + \
nameList_[i] + "</moose:product>\n"
- enzAnno2 += "<moose:groupName>" + cleanEnzname + "_" + \
- str(enz.getId().value) + "_" + \
- str(enz.getDataIndex()) + "_" + "</moose:groupName>\n"
- enzAnno2 += "<moose:stage>2</moose:stage> \n"
- if enzannoexist:
- enzAnno2 = enzAnno2 + enzGpnCorCol
- enzAnno2 += "</moose:EnzymaticReaction>"
-
- foundEnzymeEP = True
-
+ enzAnno2 += "<moose:groupName>" + cleanEnzname + "_" + \
+ str(enz.getId().value) + "_" + \
+ str(enz.getDataIndex()) + "_" + "</moose:groupName>\n"
+ enzAnno2 += "<moose:stage>2</moose:stage> \n"
+ if enzannoexist:
+ enzAnno2 = enzAnno2 + enzGpnCorCol
+ enzAnno2 += "</moose:EnzymaticReaction>"
+
+ foundEnzymeEP = True
if foundEnzymeEP:
enzyme = cremodel_.createReaction()
enzyme.setId(enzsetIdP)
@@ -407,9 +406,11 @@ def writeEnz(modelpath, cremodel_, sceneitems,groupInfo):
printParameters(kl, "k3", k3, unit)
if groupName != moose.element('/'):
if groupName not in groupInfo:
- groupInfo[groupName]=[enzsetIdP]
+ #groupInfo[groupName]=[nameList_[0]]
+ pass
else:
- groupInfo[groupName].append(enzsetIdP)
+ #groupInfo[groupName].append(nameList_[0])
+ pass
else:
print (cplxeperror)
elif(enz.className == "MMenz" or enz.className == "ZombieMMenz"):
@@ -776,6 +777,7 @@ def writeFunc(modelpath, cremodel_):
foundFunc = False
for func in funcs:
if func:
+ #not neccessary parent is compartment can be group also
if func.parent.className == "CubeMesh" or func.parent.className == "CyclMesh":
if len(moose.element(func).neighbors["valueOut"]) > 0:
funcEle = moose.element(
@@ -785,8 +787,16 @@ def writeFunc(modelpath, cremodel_):
funcEle.name + "_" + str(funcEle.getId().value) + "_" + str(funcEle.getDataIndex()) + "_"))
expr = " "
expr = str(moose.element(func).expr)
- if not expr:
+ if expr:
foundFunc = True
+ item = func.path + '/x[0]'
+ sumtot = moose.element(item).neighbors["input"]
+ for i in range(0, len(sumtot)):
+ v = "x" + str(i)
+ if v in expr:
+ z = str(idBeginWith(str(convertSpecialChar(sumtot[i].name + "_" + str(moose.element(
+ sumtot[i]).getId().value) + "_" + str(moose.element(sumtot[i]).getDataIndex())) + "_")))
+ expr = expr.replace(v, z)
else:
foundFunc = True
fName = idBeginWith(convertSpecialChar(func.parent.name +
diff --git a/moose-core/python/moose/chemMerge/merge.py b/moose-core/python/moose/chemMerge/merge.py
index 49d7d7162b9ee6ee4c98d86b5769d6bc68967343..d3d800582fd3428c253ce0fed2716b64889b18d7 100644
--- a/moose-core/python/moose/chemMerge/merge.py
+++ b/moose-core/python/moose/chemMerge/merge.py
@@ -12,7 +12,7 @@
#** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS
#Created : Friday Dec 16 23:19:00 2016(+0530)
#Version
-#Last-Updated: Wed Oct 25 16:25:33 2017(+0530)
+#Last-Updated: Wed May 21 11:52:33 2018(+0530)
# By: Harsha
#**********************************************************************/
@@ -36,17 +36,21 @@
# --- if same Enz Name and same sub and prd then nothing is copied
# --- if same Enz Name but sub or prd is different then duplicated and copied
# -- Function are copied only if destination pool to which its suppose to connect doesn't exist with function of its own
+# which is a limitation in moose that no function can be connected to same pool
#
'''
Change log
-Oct 25: line to load SBML file was commented, which is uncommented now and also while copying MMenz had a problem which also cleaned up
-Oct 14: absolute import with mtypes just for python3
+May 21 : Instead of A and B changed to S and D (source and destination),
+ - at Pool level Neutral and its objected where copied, but again at Enz and Reac level copying was done which caused duplicates which is taken care
+ - If Neutral object name is changed for destination, then change made to source file which would be easy to copy else path issue will come
+Oct 25 : line to load SBML file was commented, which is uncommented now and also while copying MMenz had a problem which also cleaned up
+Oct 14 : absolute import with mtypes just for python3
-Oct 12: clean way of cheking the type of path provided, filepath,moose obj, moose path are taken,
+Oct 12 : clean way of cheking the type of path provided, filepath,moose obj, moose path are taken,
if source is empty then nothing to copy,
if destination was empty list is update with new object
-Oct 11: missing group are copied instead of creating one in new path which also copies Annotator info
+Oct 11 : missing group are copied instead of creating one in new path which also copies Annotator info
earlier if user asked to save the model, it was saving default to kkit format, now user need to run the command to save (if this is run in command)
To check: When Gui is allowed to merge 2 models, need to see what happens
@@ -113,34 +117,33 @@ def mergeChemModel(src,des):
grpNotcopiedyet = []
dictComptA = dict( [ (i.name,i) for i in moose.wildcardFind(modelA+'/##[ISA=ChemCompt]') ] )
dictComptB = dict( [ (i.name,i) for i in moose.wildcardFind(modelB+'/##[ISA=ChemCompt]') ] )
+
poolNotcopiedyet = []
if len(dictComptA):
for key in list(dictComptA.keys()):
if key not in dictComptB:
- # if compartmentname from modelB does not exist in modelA, then copy
+ # if compartment name from modelB does not exist in modelA, then copy
copy = moose.copy(dictComptA[key],moose.element(modelB))
+
dictComptB[key]=moose.element(copy)
else:
- #if compartmentname from modelB exist in modelA,
+ #if compartment name from modelB exist in modelA,
#volume is not same, then change volume of ModelB same as ModelA
if abs(dictComptB[key].volume - dictComptA[key].volume):
#hack for now
while (abs(dictComptB[key].volume - dictComptA[key].volume) != 0.0):
dictComptA[key].volume = float(dictComptB[key].volume)
dictComptB = dict( [ (i.name,i) for i in moose.wildcardFind(modelB+'/##[ISA=ChemCompt]') ] )
-
+
#Mergering pool
- poolMerge(dictComptB[key],dictComptA[key],poolNotcopiedyet)
-
- if grpNotcopiedyet:
- # objA = moose.element(comptA).parent.name
- # if not moose.exists(objA+'/'+comptB.name+'/'+bpath.name):
- # print bpath
- # moose.copy(bpath,moose.element(objA+'/'+comptB.name))
- pass
+ poolMerge(dictComptA[key],dictComptB[key],poolNotcopiedyet)
+
+
comptBdict = comptList(modelB)
+
poolListinb = {}
poolListinb = updatePoolList(comptBdict)
+
R_Duplicated, R_Notcopiedyet,R_Daggling = [], [], []
E_Duplicated, E_Notcopiedyet,E_Daggling = [], [], []
for key in list(dictComptA.keys()):
@@ -148,27 +151,27 @@ def mergeChemModel(src,des):
funcExist,funcNotallowed = functionMerge(dictComptB,dictComptA,key)
poolListinb = updatePoolList(dictComptB)
- R_Duplicated,R_Notcopiedyet,R_Daggling = reacMerge(dictComptB,dictComptA,key,poolListinb)
+ R_Duplicated,R_Notcopiedyet,R_Daggling = reacMerge(dictComptA,dictComptB,key,poolListinb)
poolListinb = updatePoolList(dictComptB)
E_Duplicated,E_Notcopiedyet,E_Daggling = enzymeMerge(dictComptB,dictComptA,key,poolListinb)
- '''
- if isinstance(src, str):
- if os.path.isfile(src) == True:
- spath, sfile = os.path.split(src)
- else:
- sfile = src
- else:
- sfile = src
- if isinstance(des, str):
- print " A str",des
- if os.path.isfile(des) == True:
- dpath, dfile = os.path.split(des)
- else:
- dfile = des
- else:
- dfile = des
- '''
+
+ # if isinstance(src, str):
+ # if os.path.isfile(src) == True:
+ # spath, sfile = os.path.split(src)
+ # else:
+ # sfile = src
+ # else:
+ # sfile = src
+ # if isinstance(des, str):
+ # print " A str",des
+ # if os.path.isfile(des) == True:
+ # dpath, dfile = os.path.split(des)
+ # else:
+ # dfile = des
+ # else:
+ # dfile = des
+
print("\nThe content of %s (src) model is merged to %s (des)." %(sfile, dfile))
# Here any error or warning during Merge is written it down
if funcExist:
@@ -216,6 +219,7 @@ def mergeChemModel(src,des):
print ("Enzyme:")
for ed in list(E_Daggling):
print ("%s " %str(ed.name))
+
## Model is saved
print ("\n ")
print ('\nMerged model is available under moose.element(\'%s\')' %(modelB))
@@ -249,6 +253,7 @@ def mergeChemModel(src,des):
# print (' If you are in python terminal you could save \n >moose.mooseWriteKkit(\'%s\',\'filename.g\')' %(modelB))
# print (' If you are in python terminal you could save \n >moose.mooseWriteSBML(\'%s\',\'filename.g\')' %(modelB))
#return modelB
+
else:
print ('\nSource file has no objects to copy(\'%s\')' %(modelA))
return moose.element('/')
@@ -375,62 +380,72 @@ def deleteSolver(modelRoot):
if moose.exists((st_ksolve).path):
moose.delete(st_ksolve)
-def poolMerge(comptA,comptB,poolNotcopiedyet):
-
- aCmptGrp = moose.wildcardFind(comptA.path+'/#[TYPE=Neutral]')
- aCmptGrp = aCmptGrp +(moose.element(comptA.path),)
-
- bCmptGrp = moose.wildcardFind(comptB.path+'/#[TYPE=Neutral]')
- bCmptGrp = bCmptGrp +(moose.element(comptB.path),)
-
- objA = moose.element(comptA.path).parent.name
- objB = moose.element(comptB.path).parent.name
- for bpath in bCmptGrp:
- grp_cmpt = ((bpath.path).replace(objB,objA)).replace('[0]','')
- if moose.exists(grp_cmpt) :
- if moose.element(grp_cmpt).className != bpath.className:
+def poolMerge(comptS,comptD,poolNotcopiedyet):
+ #Here from source file all the groups are check if exist, if doesn't exist then create those groups
+ #Then for that group pools are copied
+ SCmptGrp = moose.wildcardFind(comptS.path+'/#[TYPE=Neutral]')
+ SCmptGrp = SCmptGrp +(moose.element(comptS.path),)
+
+ DCmptGrp = moose.wildcardFind(comptD.path+'/#[TYPE=Neutral]')
+ DCmptGrp = DCmptGrp +(moose.element(comptD.path),)
+
+ objS = moose.element(comptS.path).parent.name
+ objD = moose.element(comptD.path).parent.name
+
+ for spath in SCmptGrp:
+ grp_cmpt = ((spath.path).replace(objS,objD)).replace('[0]','')
+ if moose.exists(grp_cmpt):
+ #If path exist, but its not the Neutral obj then creating a neutral obj with _grp
+ #It has happened that pool, reac, enz name might exist with the same name, which when tried to create a group
+ # it silently ignored the path and object copied under that pool instead of Group
+ if moose.element(grp_cmpt).className != spath.className:
grp_cmpt = grp_cmpt+'_grp'
- bpath.name = bpath.name+"_grp"
- l = moose.Neutral(grp_cmpt)
+ moose.Neutral(grp_cmpt)
+ # If group name is changed while creating in destination, then in source file the same is changed,
+ # so that later path issue doens't come
+ spath.name = spath.name+'_grp'
else:
- #moose.Neutral(grp_cmpt)
- src = bpath
- srcpath = (bpath.parent).path
- des = srcpath.replace(objB,objA)
- moose.copy(bpath,moose.element(des))
-
- apath = moose.element(bpath.path.replace(objB,objA))
-
- bpoollist = moose.wildcardFind(bpath.path+'/#[ISA=PoolBase]')
- apoollist = moose.wildcardFind(apath.path+'/#[ISA=PoolBase]')
- for bpool in bpoollist:
- if bpool.name not in [apool.name for apool in apoollist]:
- copied = copy_deleteUnlyingPoolObj(bpool,apath)
+ #Neutral obj from src if doesn't exist in destination,then create src's Neutral obj in des
+ src = spath
+ srcpath = (spath.parent).path
+ des = srcpath.replace(objS,objD)
+ moose.Neutral(moose.element(des).path+'/'+spath.name)
+
+ dpath = moose.element(spath.path.replace(objS,objD))
+ spoollist = moose.wildcardFind(spath.path+'/#[ISA=PoolBase]')
+ dpoollist = moose.wildcardFind(dpath.path+'/#[ISA=PoolBase]')
+ #check made, for a given Neutral or group if pool doesn't exist then copied
+ # but some pool if enzyme cplx then there are holded untill that enzyme is copied in
+ # `enzymeMerge` function
+ for spool in spoollist:
+ if spool.name not in [dpool.name for dpool in dpoollist]:
+ copied = copy_deleteUnlyingPoolObj(spool,dpath)
if copied == False:
#hold it for later, this pool may be under enzyme, as cplx
- poolNotcopiedyet.append(bpool)
-
+ poolNotcopiedyet.append(spool)
+
def copy_deleteUnlyingPoolObj(pool,path):
# check if this pool is under compartement or under enzyme?(which is enzyme_cplx)
# if enzyme_cplx then don't copy untill this perticular enzyme is copied
- # case: This enzyme_cplx might exist in modelA if enzyme exist
- # which will automatically copie's the pool
+ # case: This enzyme_cplx might exist in modelA if enzyme exist then this
+ # will automatically copie's the pool
copied = False
if pool.parent.className not in ["Enz","ZombieEnz","MMenz","ZombieMMenz"]:
- poolcopied = moose.copy(pool,path)
- copied = True
- # deleting function and enzyme which gets copied if exist under pool
- # This is done to ensure daggling function / enzyme not copied.
- funclist = []
- for types in ['setConc','setN','increment']:
- funclist.extend(moose.element(poolcopied).neighbors[types])
-
- for fl in funclist:
- moose.delete(fl)
- enzlist = moose.element(poolcopied).neighbors['reac']
- for el in list(set(enzlist)):
- moose.delete(el.path)
+ if path.className in ["Neutral","CubeMesh","CyclMesh"]:
+ poolcopied = moose.copy(pool,path)
+ copied = True
+ # deleting function and enzyme which gets copied if exist under pool
+ # This is done to ensure daggling function / enzyme not copied.
+ funclist = []
+ for types in ['setConc','setN','increment']:
+ funclist.extend(moose.element(poolcopied).neighbors[types])
+
+ for fl in funclist:
+ moose.delete(fl)
+ enzlist = moose.element(poolcopied).neighbors['reac']
+ for el in list(set(enzlist)):
+ moose.delete(el.path)
return copied
def updatePoolList(comptAdict):
@@ -439,88 +454,90 @@ def updatePoolList(comptAdict):
poolListina[key] = plist
return poolListina
-def enzymeMerge(comptA,comptB,key,poolListina):
+def enzymeMerge(comptD,comptS,key,poolListind):
war_msg = ""
RE_Duplicated, RE_Notcopiedyet, RE_Daggling = [], [], []
- comptApath = moose.element(comptA[key]).path
- comptBpath = moose.element(comptB[key]).path
- objA = moose.element(comptApath).parent.name
- objB = moose.element(comptBpath).parent.name
- enzyListina = moose.wildcardFind(comptApath+'/##[ISA=EnzBase]')
- enzyListinb = moose.wildcardFind(comptBpath+'/##[ISA=EnzBase]')
- for eb in enzyListinb:
- eBsubname, eBprdname = [],[]
- eBsubname = subprdList(eb,"sub")
- eBprdname = subprdList(eb,"prd")
+ comptDpath = moose.element(comptD[key]).path
+ comptSpath = moose.element(comptS[key]).path
+ objD = moose.element(comptDpath).parent.name
+ objS = moose.element(comptSpath).parent.name
+ #nzyListina => enzyListind
+
+ enzyListind = moose.wildcardFind(comptDpath+'/##[ISA=EnzBase]')
+ enzyListins = moose.wildcardFind(comptSpath+'/##[ISA=EnzBase]')
+ for es in enzyListins:
+ eSsubname, eSprdname = [],[]
+ eSsubname = subprdList(es,"sub")
+ eSprdname = subprdList(es,"prd")
allexists, allexistp = False, False
allclean = False
- poolinAlist = poolListina[findCompartment(eb).name]
- for pA in poolinAlist:
- if eb.parent.name == pA.name:
- eapath = eb.parent.path.replace(objB,objA)
+ poolinDlist = poolListind[findCompartment(es).name]
+ for pD in poolinDlist:
+ if es.parent.name == pD.name:
+ edpath = es.parent.path.replace(objS,objD)
- if not moose.exists(eapath+'/'+eb.name):
+ if not moose.exists(edpath+'/'+es.name):
#This will take care
# -- If same enzparent name but different enzyme name
# -- or different parent/enzyme name
- if eBsubname and eBprdname:
- allexists = checkexist(eBsubname,objB,objA)
- allexistp = checkexist(eBprdname,objB,objA)
+ if eSsubname and eSprdname:
+ allexists = checkexist(eSsubname,objS,objD)
+ allexistp = checkexist(eSprdname,objS,objD)
if allexists and allexistp:
- enzPool = moose.element(pA.path)
- eapath = eb.parent.path.replace(objB,objA)
- enz = moose.element(moose.copy(eb,moose.element(eapath)))
+ enzPool = moose.element(pD.path)
+ edpath = es.parent.path.replace(objS,objD)
+ enz = moose.element(moose.copy(es,moose.element(edpath)))
enzPool = enz.parent
- if eb.className in ["ZombieEnz","Enz"]:
+ if es.className in ["ZombieEnz","Enz"]:
moose.connect(moose.element(enz),"enz",enzPool,"reac")
- if eb.className in ["ZombieMMenz","MMenz"]:
+ if es.className in ["ZombieMMenz","MMenz"]:
moose.connect(enzPool,"nOut",enz,"enzDest")
- connectObj(enz,eBsubname,"sub",comptA,war_msg)
- connectObj(enz,eBprdname,"prd",comptA,war_msg)
+ connectObj(enz,eSsubname,"sub",comptD,war_msg)
+ connectObj(enz,eSprdname,"prd",comptD,war_msg)
allclean = True
else:
# didn't find sub or prd for this Enzyme
- RE_Notcopiedyet.append(eb)
+ RE_Notcopiedyet.append(es)
else:
# -- it is dagging reaction
- RE_Daggling.append(eb)
+ RE_Daggling.append(es)
else:
#Same Enzyme name
# -- Same substrate and product including same volume then don't copy
# -- different substrate/product or if sub/prd's volume is different then DUPLICATE the Enzyme
allclean = False
- ea = moose.element(eb.path.replace(objB,objA))
- eAsubname = subprdList(ea,"sub")
- eBsubname = subprdList(eb,"sub")
- hasSamenoofsublen,hasSameS,hasSamevols = same_len_name_vol(eAsubname,eBsubname)
-
- eAprdname = subprdList(ea,"prd")
- eBprdname = subprdList(eb,"prd")
- hasSamenoofprdlen,hasSameP,hasSamevolp = same_len_name_vol(eAprdname,eBprdname)
+ ed = moose.element(es.path.replace(objS,objD))
+ eDsubname = subprdList(ed,"sub")
+ eSsubname = subprdList(es,"sub")
+ hasSamenoofsublen,hasSameS,hasSamevols = same_len_name_vol(eDsubname,eSsubname)
+
+ eDprdname = subprdList(ed,"prd")
+ eSprdname = subprdList(es,"prd")
+ hasSamenoofprdlen,hasSameP,hasSamevolp = same_len_name_vol(eDprdname,eSprdname)
if not all((hasSamenoofsublen,hasSameS,hasSamevols,hasSamenoofprdlen,hasSameP,hasSamevolp)):
# May be different substrate or product or volume of Sub/prd may be different,
# Duplicating the enzyme
- if eBsubname and eBprdname:
+ if eSsubname and eSprdname:
allexists,allexistp = False,False
- allexists = checkexist(eBsubname,objB,objA)
- allexistp = checkexist(eBprdname,objB,objA)
+ allexists = checkexist(eSsubname,objS,objD)
+ allexistp = checkexist(eSprdname,objS,objD)
if allexists and allexistp:
- eb.name = eb.name+"_duplicated"
- if eb.className in ["ZombieEnz","Enz"]:
- eapath = eb.parent.path.replace(objB,objA)
- enz = moose.copy(eb,moose.element(eapath))
- moose.connect(enz, 'enz', eapath, 'reac' )
-
- if eb.className in ["ZombieMMenz","MMenz"]:
- eapath = eb.parent.path.replace(objB,objA)
- enz = moose.copy(eb,moose.element(eapath))
+ es.name = es.name+"_duplicated"
+ if es.className in ["ZombieEnz","Enz"]:
+ edpath = es.parent.path.replace(objS,objD)
+ enz = moose.copy(es,moose.element(edpath))
+ moose.connect(enz, 'enz', edpath, 'reac' )
+
+ if es.className in ["ZombieMMenz","MMenz"]:
+ edpath = es.parent.path.replace(objS,objD)
+ enz = moose.copy(es,moose.element(edpath))
enzinfo = moose.Annotator(enz.path+'/info')
moose.connect(moose.element(enz).parent,"nOut",moose.element(enz),"enzDest")
#moose.connect(moose.element(enz),"enz",moose.element(enz).parent,"reac")
- connectObj(enz,eBsubname,"sub",comptA,war_msg)
- connectObj(enz,eBprdname,"prd",comptA,war_msg)
+ connectObj(enz,eSsubname,"sub",comptD,war_msg)
+ connectObj(enz,eSprdname,"prd",comptD,war_msg)
RE_Duplicated.append(enz)
allclean = True
else:
@@ -531,51 +548,53 @@ def enzymeMerge(comptA,comptB,key,poolListina):
if not allclean:
# didn't find sub or prd for this enzyme
# -- it may be connected Enzyme cplx
- if eBsubname and eBprdname:
- RE_Notcopiedyet.append(eb)
+ if eSsubname and eSprdname:
+ RE_Notcopiedyet.append(es)
else:
- RE_Daggling.append(eb)
+ RE_Daggling.append(es)
return RE_Duplicated,RE_Notcopiedyet,RE_Daggling
-def reacMerge(comptA,comptB,key,poolListina):
+def reacMerge(comptS,comptD,key,poolListina):
RE_Duplicated, RE_Notcopiedyet, RE_Daggling = [], [], []
war_msg = ""
- comptApath = moose.element(comptA[key]).path
- comptBpath = moose.element(comptB[key]).path
- objA = moose.element(comptApath).parent.name
- objB = moose.element(comptBpath).parent.name
- reacListina = moose.wildcardFind(comptApath+'/##[ISA=ReacBase]')
- reacListinb = moose.wildcardFind(comptBpath+'/##[ISA=ReacBase]')
- for rb in reacListinb:
- rBsubname, rBprdname = [],[]
- rBsubname = subprdList(rb,"sub")
- rBprdname = subprdList(rb,"prd")
+ comptSpath = moose.element(comptS[key]).path
+ comptDpath = moose.element(comptD[key]).path
+ objS = moose.element(comptSpath).parent.name
+ objD = moose.element(comptDpath).parent.name
+
+ reacListins = moose.wildcardFind(comptSpath+'/##[ISA=ReacBase]')
+ reacListind = moose.wildcardFind(comptDpath+'/##[ISA=ReacBase]')
+
+ for rs in reacListins:
+ rSsubname, rSprdname = [],[]
+ rSsubname = subprdList(rs,"sub")
+ rSprdname = subprdList(rs,"prd")
allexists, allexistp = False, False
allclean = False
- if rb.name not in [ra.name for ra in reacListina]:
+ if rs.name not in [rd.name for rd in reacListind]:
# reaction name not found then copy
# And assuming that pools are copied earlier EXPECT POOL CPLX
- #To be assured the it takes correct compartment name incase reaction sub's
- #belongs to different compt
- key = findCompartment(rb).name
- if rBsubname and rBprdname:
- allexists = checkexist(rBsubname,objB,objA)
- allexistp = checkexist(rBprdname,objB,objA)
+ # To be assured the it takes correct compartment name incase reaction sub's
+ # belongs to different compt
+ key = findCompartment(rs).name
+ if rSsubname and rSprdname:
+ allexists = checkexist(rSsubname,objS,objD)
+ allexistp = checkexist(rSprdname,objS,objD)
if allexists and allexistp:
- rapath = rb.parent.path.replace(objB,objA)
- reac = moose.copy(rb,moose.element(rapath))
- connectObj(reac,rBsubname,"sub",comptA,war_msg)
- connectObj(reac,rBprdname,"prd",comptA,war_msg)
+ rdpath = rs.parent.path.replace(objS,objD)
+ reac = moose.copy(rs,moose.element(rdpath))
+ connectObj(reac,rSsubname,"sub",comptD,war_msg)
+ connectObj(reac,rSprdname,"prd",comptD,war_msg)
allclean = True
else:
# didn't find sub or prd for this reaction
# -- it may be connected Enzyme cplx
- RE_Notcopiedyet.append(rb)
+ RE_Notcopiedyet.append(rs)
else:
# -- it is dagging reaction
- RE_Daggling.append(rb)
+ RE_Daggling.append(rs)
#print ("This reaction \""+rb.path+"\" has no substrate/product daggling reaction are not copied")
#war_msg = war_msg+"\nThis reaction \""+rb.path+"\" has no substrate/product daggling reaction are not copied"
@@ -583,30 +602,31 @@ def reacMerge(comptA,comptB,key,poolListina):
#Same reaction name
# -- Same substrate and product including same volume then don't copy
# -- different substrate/product or if sub/prd's volume is different then DUPLICATE the reaction
+
allclean = False
- for ra in reacListina:
- if rb.name == ra.name:
- rAsubname = subprdList(ra,"sub")
- rBsubname = subprdList(rb,"sub")
- hasSamenoofsublen,hasSameS,hasSamevols = same_len_name_vol(rAsubname,rBsubname)
-
- rAprdname = subprdList(ra,"prd")
- rBprdname = subprdList(rb,"prd")
- hasSamenoofprdlen,hasSameP,hasSamevolp = same_len_name_vol(rAprdname,rBprdname)
+ for rd in reacListind:
+ if rs.name == rd.name:
+ rSsubname = subprdList(rs,"sub")
+ rDsubname = subprdList(rd,"sub")
+ hasSamenoofsublen,hasSameS,hasSamevols = same_len_name_vol(rSsubname,rDsubname)
+
+ rSprdname = subprdList(rs,"prd")
+ rDprdname = subprdList(rd,"prd")
+ hasSamenoofprdlen,hasSameP,hasSamevolp = same_len_name_vol(rSprdname,rDprdname)
if not all((hasSamenoofsublen,hasSameS,hasSamevols,hasSamenoofprdlen,hasSameP,hasSamevolp)):
# May be different substrate or product or volume of Sub/prd may be different,
# Duplicating the reaction
- if rBsubname and rBprdname:
+ if rSsubname and rSprdname:
allexists,allexistp = False,False
- allexists = checkexist(rBsubname,objB,objA)
- allexistp = checkexist(rBprdname,objB,objA)
+ allexists = checkexist(rSsubname,objS,objD)
+ allexistp = checkexist(rSprdname,objS,objD)
if allexists and allexistp:
- rb.name = rb.name+"_duplicated"
+ rs.name = rs.name+"_duplicated"
#reac = moose.Reac(comptA[key].path+'/'+rb.name+"_duplicated")
- rapath = rb.parent.path.replace(objB,objA)
- reac = moose.copy(rb,moose.element(rapath))
- connectObj(reac,rBsubname,"sub",comptA,war_msg)
- connectObj(reac,rBprdname,"prd",comptA,war_msg)
+ rdpath = rs.parent.path.replace(objS,objD)
+ reac = moose.copy(rs,moose.element(rdpath))
+ connectObj(reac,rSsubname,"sub",comptD,war_msg)
+ connectObj(reac,rSprdname,"prd",comptD,war_msg)
RE_Duplicated.append(reac)
allclean = True
else:
@@ -617,11 +637,11 @@ def reacMerge(comptA,comptB,key,poolListina):
if not allclean:
# didn't find sub or prd for this reaction
# -- it may be connected Enzyme cplx
- if rBsubname and rBprdname:
- RE_Notcopiedyet.append(rb)
+ if rSsubname and rSprdname:
+ RE_Notcopiedyet.append(rs)
else:
- RE_Daggling.append(rb)
-
+ RE_Daggling.append(rs)
+
return RE_Duplicated,RE_Notcopiedyet,RE_Daggling
def subprdList(reac,subprd):
@@ -704,6 +724,7 @@ def mooseIsInstance(element, classNames):
if __name__ == "__main__":
+
try:
sys.argv[1]
except IndexError:
@@ -712,7 +733,7 @@ if __name__ == "__main__":
else:
src = sys.argv[1]
if not os.path.exists(src):
- print("Filename or path does not exist",src)
+ print("Filename or path does not exist %s." %(src))
else:
try:
sys.argv[2]
@@ -722,10 +743,12 @@ if __name__ == "__main__":
else:
des = sys.argv[2]
if not os.path.exists(src):
- print("Filename or path does not exist",des)
+ print("Filename or path does not exist %s." %(des))
exit(0)
else:
+ print ("src and des %s, %s." %(src, des))
mergered = mergeChemModel(src,des)
+
'''
try:
sys.argv[3]
diff --git a/moose-core/python/moose/chemUtil/add_Delete_ChemicalSolver.py b/moose-core/python/moose/chemUtil/add_Delete_ChemicalSolver.py
index 05bd4a0a233f792f177bcc833beed7292efe3dc6..cc350ec6087ac1b731b0a274b75eeeebb8a3de77 100644
--- a/moose-core/python/moose/chemUtil/add_Delete_ChemicalSolver.py
+++ b/moose-core/python/moose/chemUtil/add_Delete_ChemicalSolver.py
@@ -1,11 +1,35 @@
# -*- coding: utf-8 -*-
import moose
+from fixXreacs import fixXreacs
+
+def positionCompt( compt ):
+ i = 0
+ while (i != len(compt)-1):
+ #print "PositionCompt ", compt[i+1],compt[i+1].volume, compt[i], compt[i].volume
+ compt[i+1].x1 += compt[i].x1
+ compt[i+1].x0 += compt[i].x1
+ i += 1
def moosedeleteChemSolver(modelRoot):
"""Delete solvers from Chemical Compartment
"""
compts = moose.wildcardFind(modelRoot + '/##[ISA=ChemCompt]')
+ for compt in compts:
+ if moose.exists(compt.path + '/stoich'):
+ st = moose.element(compt.path + '/stoich')
+ st_ksolve = st.ksolve
+ st_dsolve = st.dsolve
+
+ moose.delete(st)
+ if moose.exists((st_ksolve).path):
+ moose.delete(st_ksolve)
+ print("KSolver is deleted for modelpath %s " % modelRoot)
+ if moose.exists((st_dsolve).path):
+ moose.delete(st_dsolve)
+ print("DSolver is deleted for modelpath %s " % modelRoot)
+ '''
+ compts = moose.wildcardFind(modelRoot + '/##[ISA=ChemCompt]')
for compt in compts:
if moose.exists(compt.path + '/stoich'):
st = moose.element(compt.path + '/stoich')
@@ -15,6 +39,7 @@ def moosedeleteChemSolver(modelRoot):
moose.delete(st_ksolve)
print("Solver is deleted for modelpath %s " % modelRoot)
+ '''
def mooseaddChemSolver(modelRoot, solver):
"""
@@ -58,6 +83,56 @@ def mooseaddChemSolver(modelRoot, solver):
def setCompartmentSolver(modelRoot, solver):
+ comptlist = dict((c, c.volume) for c in moose.wildcardFind(modelRoot + '/##[ISA=ChemCompt]'))
+ comptVol = {}
+ compts = []
+ vol = [v for k,v in comptlist.items()]
+ volumeSort = sorted(vol)
+ for k,v in comptlist.items():
+ comptVol[k]= v
+ for volSor in volumeSort:
+ for a,b in comptVol.items():
+ if b == volSor:
+ compts.append(a)
+
+ #compts = [key for key, value in sorted(comptlist.items(), key=lambda (k,v): (v,k))]
+ if ( len(compts) == '0'):
+ print ("Atleast one compartment is required ")
+ return
+ else:
+ if ( len(compts) > 3 ):
+ print ("Warning: setSolverOnCompt Cannot handle " , len(compts) , " chemical compartments\n")
+ return;
+
+ elif (len(compts) >1 ):
+ positionCompt(compts)
+
+ fixXreacs( modelRoot )
+
+ for compt in compts:
+ if solver != 'ee':
+ if (solver == 'gsl') or (solver == 'Runge Kutta'):
+ ksolve = moose.Ksolve(compt.path + '/ksolve')
+ if (solver == 'gssa') or (solver == 'Gillespie'):
+ ksolve = moose.Gsolve(compt.path + '/gsolve')
+
+ dsolve = moose.Dsolve(compt.path+'/dsolve')
+ stoich = moose.Stoich(compt.path + '/stoich')
+ stoich.compartment = compt
+ stoich.ksolve = ksolve
+ stoich.dsolve = dsolve
+ stoich.path = compt.path + "/##"
+ ksolveList = moose.wildcardFind(modelRoot+'/##[ISA=Ksolve]')
+ dsolveList = moose.wildcardFind(modelRoot+'/##[ISA=Dsolve]')
+ stoichList = moose.wildcardFind(modelRoot+'/##[ISA=Stoich]')
+
+ i = 0
+ while(i < len(dsolveList)-1):
+ dsolveList[i+1].buildMeshJunctions(dsolveList[i])
+ i += 1
+
+ print( " Solver is added to model path %s" % modelRoot )
+ '''
compts = moose.wildcardFind(modelRoot + '/##[ISA=ChemCompt]')
for compt in compts:
if (solver == 'gsl') or (solver == 'Runge Kutta'):
@@ -80,3 +155,4 @@ def setCompartmentSolver(modelRoot, solver):
for i in stoichList:
i.filterXreacs()
print( " Solver is added to model path %s" % modelRoot )
+ '''
diff --git a/moose-core/python/moose/utils.py b/moose-core/python/moose/utils.py
index 37bf0deb9f82cc233d4f8f96de103b2c0c83cc15..9abf7ae5d2636cbccfc80cb0c19535233c4f6def 100644
--- a/moose-core/python/moose/utils.py
+++ b/moose-core/python/moose/utils.py
@@ -23,7 +23,6 @@ import re
from moose.moose_constants import *
from moose.print_utils import *
-
# Print and Plot utilities.
try:
from moose.plot_utils import *
@@ -341,7 +340,10 @@ def autoposition(root):
stack.extend([childcomp for childcomp in map(moose.element, comp.neighbors['raxial']) if childcomp.z == 0])
return ret
-
+def loadModel(filename, target,method='ee'):
+ moose.loadModel(filename,target)
+ moose.mooseaddChemSolver(target,method)
+
def readcell_scrambled(filename, target, method='ee'):
"""A special version for handling cases where a .p file has a line
with specified parent yet to be defined.
diff --git a/moose-core/python/rdesigneur/moogul.py b/moose-core/python/rdesigneur/moogul.py
new file mode 100644
index 0000000000000000000000000000000000000000..9f2ecba50b89d1808e711e471474307ca8e70197
--- /dev/null
+++ b/moose-core/python/rdesigneur/moogul.py
@@ -0,0 +1,314 @@
+# Moogul.py: MOOSE Graphics Using Lines
+# This is a fallback graphics interface for displaying neurons using
+# regular matplotlib routines.
+# Put in because the GL versions like moogli need all sorts of difficult
+# libraries and dependencies.
+# Copyright (C) Upinder S. Bhalla NCBS 2018
+# This program is licensed under the GNU Public License version 3.
+#
+
+import moose
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d.art3d import Line3DCollection
+
+class MoogulError( Exception ):
+ def __init__( self, value ):
+ self.value = value
+ def __str__( self ):
+ return repr( self.value )
+
+class MooView:
+ ''' The MooView class is a window in which to display one or more
+ moose cells, using the MooCell class.'''
+ def __init__( self, swx = 10, swy = 12, hideAxis = True
+ ):
+ plt.ion()
+ self.fig_ = plt.figure( figsize = (swx, swy) )
+ self.ax = self.fig_.add_subplot(111, projection='3d' )
+ self.drawables_ = []
+ self.fig_.canvas.mpl_connect("key_press_event", self.moveView )
+ plt.rcParams['keymap.xscale'] = ''
+ plt.rcParams['keymap.yscale'] = ''
+ plt.rcParams['keymap.zoom'] = ''
+ plt.rcParams['keymap.back'] = ''
+ plt.rcParams['keymap.home'] = ''
+ plt.rcParams['keymap.forward'] = ''
+ plt.rcParams['keymap.all_axes'] = ''
+ self.hideAxis = hideAxis
+ if self.hideAxis:
+ self.ax.set_axis_off()
+ #self.ax.margins( tight = True )
+ self.ax.margins()
+ self.sensitivity = 7.0 # degrees rotation
+ self.zoom = 1.05
+
+ def addDrawable( self, n ):
+ self.drawables_.append( n )
+
+ def firstDraw( self, rotation = 0.0, elev = 0.0, azim = 0.0 ):
+ self.coordMax = 0.0
+ self.coordMin = 0.0
+ if rotation == 0.0:
+ self.doRotation = False
+ self.rotation = 7 # default rotation per frame, in degrees.
+ else:
+ self.doRotation = True
+ self.rotation = rotation * 180/np.pi # arg units: radians/frame
+
+ self.azim = azim * 180/np.pi
+ self.elev = elev * 180/np.pi
+ for i in self.drawables_:
+ cmax, cmin = i.drawForTheFirstTime( self.ax )
+ self.coordMax = max( cmax, self.coordMax )
+ self.coordMin = min( cmin, self.coordMin )
+ if self.coordMin == self.coordMax:
+ self.coordMax = 1+self.coordMin
+
+
+ self.ax.set_xlim3d( self.coordMin, self.coordMax )
+ self.ax.set_ylim3d( self.coordMin, self.coordMax )
+ self.ax.set_zlim3d( self.coordMin, self.coordMax )
+ self.ax.view_init( elev = self.elev, azim = self.azim )
+ #self.ax.view_init( elev = -80.0, azim = 90.0 )
+ #self.colorbar = plt.colorbar( self.drawables_[0].segments )
+ self.colorbar = self.fig_.colorbar( self.drawables_[0].segments )
+ self.colorbar.set_label( self.drawables_[0].fieldInfo[3])
+ self.timeStr = self.ax.text2D( 0.05, 0.05,
+ "Time= 0.0", transform=self.ax.transAxes)
+ self.fig_.canvas.draw()
+ plt.show()
+
+ def updateValues( self ):
+ time = moose.element( '/clock' ).currentTime
+ self.timeStr.set_text( "Time= " + str(time) )
+ for i in self.drawables_:
+ i.updateValues()
+ if self.doRotation and abs( self.rotation ) < 120:
+ self.ax.azim += self.rotation
+ #self.fig_.canvas.draw()
+ plt.pause(0.001)
+
+ def moveView(self, event):
+ x0 = self.ax.get_xbound()[0]
+ x1 = self.ax.get_xbound()[1]
+ xk = (x0 - x1) / self.sensitivity
+ y0 = self.ax.get_ybound()[0]
+ y1 = self.ax.get_ybound()[1]
+ yk = (y0 - y1) / self.sensitivity
+ z0 = self.ax.get_zbound()[0]
+ z1 = self.ax.get_zbound()[1]
+ zk = (z0 - z1) / self.sensitivity
+
+ if event.key == "up" or event.key == "k":
+ self.ax.set_ybound( y0 + yk, y1 + yk )
+ if event.key == "down" or event.key == "j":
+ self.ax.set_ybound( y0 - yk, y1 - yk )
+ if event.key == "left" or event.key == "h":
+ self.ax.set_xbound( x0 + xk, x1 + xk )
+ if event.key == "right" or event.key == "l":
+ self.ax.set_xbound( x0 - xk, x1 - xk )
+ if event.key == "ctrl+up":
+ self.ax.set_zbound( z0 + zk, z1 + zk )
+ if event.key == "ctrl+down":
+ self.ax.set_zbound( z0 - zk, z1 - zk )
+ if event.key == "." or event.key == ">": # zoom in, bigger
+ self.ax.set_xbound( x0/self.zoom, x1/self.zoom )
+ self.ax.set_ybound( y0/self.zoom, y1/self.zoom )
+ self.ax.set_zbound( z0/self.zoom, z1/self.zoom )
+ if event.key == "," or event.key == "<": # zoom out, smaller
+ self.ax.set_xbound( x0*self.zoom, x1*self.zoom )
+ self.ax.set_ybound( y0*self.zoom, y1*self.zoom )
+ self.ax.set_zbound( z0*self.zoom, z1*self.zoom )
+ if event.key == "a": # autoscale to fill view.
+ self.ax.set_xlim3d( self.coordMin, self.coordMax )
+ self.ax.set_ylim3d( self.coordMin, self.coordMax )
+ self.ax.set_zlim3d( self.coordMin, self.coordMax )
+ if event.key == "p": # pitch
+ self.ax.elev += self.sensitivity
+ if event.key == "P":
+ self.ax.elev -= self.sensitivity
+ if event.key == "y": # yaw
+ self.ax.azim += self.sensitivity
+ if event.key == "Y":
+ self.ax.azim -= self.sensitivity
+ # Don't have anything for roll
+ if event.key == "g":
+ self.hideAxis = not self.hideAxis
+ if self.hideAxis:
+ self.ax.set_axis_off()
+ else:
+ self.ax.set_axis_on()
+ if event.key == "t": # Turn on/off twisting/autorotate
+ self.doRotation = not self.doRotation
+ if event.key == "?": # Print out help for these commands
+ self.printMoogulHelp()
+
+ self.fig_.canvas.draw()
+
+ def printMoogulHelp( self ):
+ print( '''
+ Key bindings for Moogul:
+ Up or k: pan object up
+ Down or j: pan object down
+ left or h: pan object left. Bug: direction depends on azimuth.
+ right or l: pan object right Bug: direction depends on azimuth
+ . or >: Zoom in: make object appear bigger
+ , or <: Zoom out: make object appear smaller
+ a: Autoscale to fill view
+ p: Pitch down
+ P: Pitch up
+ y: Yaw counterclockwise
+ Y: Yaw counterclockwise
+ g: Toggle visibility of grid
+ t: Toggle turn (rotation along long axis of cell)
+ ?: Print this help page.
+ ''')
+
+#####################################################################
+
+class MooDrawable:
+ ''' Base class for drawing things'''
+ def __init__( self,
+ fieldInfo, field, relativeObj, maxLineWidth,
+ colormap,
+ lenScale, diaScale, autoscale,
+ valMin, valMax
+ ):
+ self.field = field
+ self.relativeObj = relativeObj
+ self.maxLineWidth = maxLineWidth
+ self.lenScale = lenScale
+ self.diaScale = diaScale
+ self.colormap = colormap
+ self.autoscale = autoscale
+ self.valMin = valMin
+ self.valMax = valMax
+ self.fieldInfo = fieldInfo
+ self.fieldScale = fieldInfo[2]
+ #FieldInfo = [baseclass, fieldGetFunc, scale, axisText, min, max]
+
+ def updateValues( self ):
+ ''' Obtains values from the associated cell'''
+ self.val = np.array([moose.getField( i, self.field, 'double' ) for i in self.activeObjs]) * self.fieldScale
+ cmap = plt.get_cmap( self.colormap )
+ if self.autoscale:
+ valMin = min( self.val )
+ valMax = max( self.val )
+ else:
+ valMin = self.valMin
+ valMax = self.valMax
+ scaleVal = (self.val - valMin) / (valMax - valMin)
+ self.rgba = [ cmap(i) for i in scaleVal ]
+ self.segments.set_color( self.rgba )
+ return
+
+ def drawForTheFirstTime( self, ax ):
+ self.segments = Line3DCollection( self.activeCoords,
+ linewidths = self.linewidth, cmap = plt.get_cmap(self.colormap) )
+ self.cax = ax.add_collection3d( self.segments )
+ self.segments.set_array( self.valMin + np.array( range( len(self.activeCoords) ) ) * (self.valMax-self.valMin)/len(self.activeCoords) )
+ return self.coordMax, self.coordMin
+
+
+
+
+#####################################################################
+
+class MooNeuron( MooDrawable ):
+ ''' Draws collection of line segments of defined dia and color'''
+ def __init__( self,
+ neuronId,
+ fieldInfo,
+ field = 'Vm',
+ relativeObj = '.',
+ maxLineWidth = 20,
+ colormap = 'jet',
+ lenScale = 1e6, diaScale = 1e6, autoscale = False,
+ valMin = -0.1, valMax = 0.05,
+ ):
+ #self.isFieldOnCompt =
+ #field in ( 'Vm', 'Im', 'Rm', 'Cm', 'Ra', 'inject', 'diameter' )
+
+ MooDrawable.__init__( self, fieldInfo, field = field,
+ relativeObj = relativeObj, maxLineWidth = maxLineWidth,
+ colormap = colormap, lenScale = lenScale,
+ diaScale = diaScale, autoscale = autoscale,
+ valMin = valMin, valMax = valMax )
+ self.neuronId = neuronId
+ self.updateCoords()
+
+ def updateCoords( self ):
+ ''' Obtains coords from the associated cell'''
+ self.compts_ = moose.wildcardFind( self.neuronId.path + "/#[ISA=CompartmentBase]" )
+ # Matplotlib3d isn't able to do full rotations about an y axis,
+ # which is what the NeuroMorpho models use, so
+ # here we shuffle the axes around. Should be an option.
+ #coords = np.array([[[i.x0,i.y0,i.z0],[i.x,i.y,i.z]]
+ #for i in self.compts_])
+ coords = np.array([[[i.z0,i.x0,i.y0],[i.z,i.x,i.y]]
+ for i in self.compts_])
+ dia = np.array([i.diameter for i in self.compts_])
+ if self.relativeObj == '.':
+ self.activeCoords = coords
+ self.activeDia = dia
+ self.activeObjs = self.compts_
+ else:
+ self.activeObjs = []
+ self.activeCoords = []
+ self.activeDia = []
+ for i,j,k in zip( self.compts_, coords, dia ):
+ if moose.exists( i.path + '/' + self.relativeObj ):
+ elm = moose.element( i.path + '/' + self.relativeObj )
+ self.activeObjs.append( elm )
+ self.activeCoords.append( j )
+ self.activeDia.append( k )
+
+ self.activeCoords = np.array( self.activeCoords ) * self.lenScale
+ self.coordMax = np.amax( self.activeCoords )
+ self.coordMin = np.amin( self.activeCoords )
+ self.linewidth = np.array( [ min(self.maxLineWidth, 1 + int(i * self.diaScale )) for i in self.activeDia ] )
+
+ return
+
+#####################################################################
+class MooReacSystem( MooDrawable ):
+ ''' Draws collection of line segments of defined dia and color'''
+ def __init__( self,
+ mooObj, fieldInfo,
+ field = 'conc',
+ relativeObj = '.',
+ maxLineWidth = 100,
+ colormap = 'jet',
+ lenScale = 1e6, diaScale = 20e6, autoscale = False,
+ valMin = 0.0, valMax = 1.0
+ ):
+
+ MooDrawable.__init__( self, fieldInfo, field = field,
+ relativeObj = relativeObj, maxLineWidth = maxLineWidth,
+ colormap = colormap, lenScale = lenScale,
+ diaScale = diaScale, autoscale = autoscale,
+ valMin = valMin, valMax = valMax )
+ self.mooObj = mooObj
+ self.updateCoords()
+
+ def updateCoords( self ):
+ ''' For now a dummy cylinder '''
+ dx = 1e-6
+ dummyDia = 20e-6
+ numObj = len( self.mooObj )
+ x = np.arange( 0, (numObj+1) * dx, dx )
+ y = np.zeros( numObj + 1)
+ z = np.zeros( numObj + 1)
+
+ coords = np.array([[[i*dx,0,0],[(i+1)*dx,0,0]] for i in range( numObj )] )
+ dia = np.ones( numObj ) * dummyDia
+ self.activeCoords = coords
+ self.activeDia = dia
+ self.activeObjs = self.mooObj
+ self.activeCoords = np.array( self.activeCoords ) * self.lenScale
+ self.coordMax = np.amax( self.activeCoords )
+ self.coordMin = np.amin( self.activeCoords )
+ self.linewidth = np.array( [ min(self.maxLineWidth, 1 + int(i * self.diaScale )) for i in self.activeDia ] )
+
+ return
diff --git a/moose-core/python/rdesigneur/rdesigneur.py b/moose-core/python/rdesigneur/rdesigneur.py
index 0ef8d5c45c81c01f2a230e9f6cdcce9ffc50cd9d..66a2ea68996ee4dae577254b438c40a2b9c72dc4 100644
--- a/moose-core/python/rdesigneur/rdesigneur.py
+++ b/moose-core/python/rdesigneur/rdesigneur.py
@@ -28,8 +28,12 @@ import sys
import time
import rdesigneur.rmoogli as rmoogli
-import rdesigneur.fixXreacs as fixXreacs
from rdesigneur.rdesigneurProtos import *
+import fixXreacs
+#from . import fixXreacs
+#from rdesigneur.rmoogli import *
+#import rmoogli
+#from rdesigneurProtos import *
from moose.neuroml.NeuroML import NeuroML
from moose.neuroml.ChannelML import ChannelML
@@ -194,9 +198,9 @@ class rdesigneur:
self.buildChemDistrib()
self._configureSolvers()
self.buildAdaptors()
+ self._buildStims()
self._buildPlots()
self._buildMoogli()
- self._buildStims()
self._configureHSolve()
self._configureClocks()
if self.verbose:
@@ -332,6 +336,19 @@ class rdesigneur:
# Here are the functions to build the type-specific prototypes.
################################################################
def buildCellProto( self ):
+ # cellProtoList args:
+ # Option 1: zero args: make standard soma, len and dia 500 um.
+ # Option 2: [name, library_proto_name]: uses library proto
+ # Option 3: [fname.suffix, cellname ]: Loads cell from file
+ # Option 4: [moose<Classname>, cellname]: Makes proto of class
+ # Option 5: [funcname, cellname]: Calls named function with specified name of cell to be made.
+ # Option 6: [path, cellname]: Copies path to library as proto
+ # Option 7: [libraryName, cellname]: Renames library entry as proto
+ # Below two options only need the first two args, rest are optional
+ # Defailt values are given.
+ # Option 8: [somaProto,name, somaDia=5e-4, somaLen=5e-4]
+ # Option 9: [ballAndStick,name, somaDia=10e-6, somaLen=10e-6,
+ # dendDia=4e-6, dendLen=200e-6, numDendSeg=1]
if len( self.cellProtoList ) == 0:
''' Make HH squid model sized compartment:
len and dia 500 microns. CM = 0.01 F/m^2, RA =
@@ -339,6 +356,7 @@ class rdesigneur:
self.elecid = makePassiveHHsoma( name = 'cell' )
assert( moose.exists( '/library/cell/soma' ) )
self.soma = moose.element( '/library/cell/soma' )
+ return
'''
self.elecid = moose.Neuron( '/library/cell' )
@@ -350,7 +368,11 @@ class rdesigneur:
'''
for i in self.cellProtoList:
- if self.checkAndBuildProto( "cell", i, \
+ if i[0] == 'somaProto':
+ self._buildElecSoma( i )
+ elif i[0] == 'ballAndStick':
+ self._buildElecBallAndStick( i )
+ elif self.checkAndBuildProto( "cell", i, \
["Compartment", "SymCompartment"], ["swc", "p", "nml", "xml"] ):
self.elecid = moose.element( '/library/' + i[1] )
else:
@@ -394,11 +416,73 @@ class rdesigneur:
["Pool"], ["g", "sbml", "xml" ] ):
self._loadChem( i[0], i[1] )
self.chemid = moose.element( '/library/' + i[1] )
+ ################################################################
+ def _buildElecSoma( self, args ):
+ parms = [ 'somaProto', 'soma', 5e-4, 5e-4 ] # somaDia, somaLen
+ for i in range( len(args) ):
+ parms[i] = args[i]
+ cell = moose.Neuron( '/library/' + parms[1] )
+ buildCompt( cell, 'soma', dia = parms[2], dx = parms[3] )
+ self.elecid = cell
+ return cell
+
+ ################################################################
+ def _buildElecBallAndStick( self, args ):
+ parms = [ 'ballAndStick', 'cell', 10e-6, 10e-6, 4e-6, 200e-6, 1 ] # somaDia, somaLen, dendDia, dendLen, dendNumSeg
+ for i in range( len(args) ):
+ parms[i] = args[i]
+ if parms[6] <= 0:
+ return _self.buildElecSoma( parms[:4] )
+ cell = moose.Neuron( '/library/' + parms[1] )
+ prev = buildCompt( cell, 'soma', dia = args[2], dx = args[3] )
+ dx = parms[5]/parms[6]
+ x = prev.x
+ for i in range( parms[6] ):
+ compt = buildCompt( cell, 'dend' + str(i), x = x, dx = dx, dia = args[4] )
+ moose.connect( prev, 'axial', compt, 'raxial' )
+ prev = compt
+ x += dx
+ self.elecid = cell
+ return cell
+ ################################################################
+ def _buildVclampOnCompt( self, dendCompts, spineCompts, stimInfo ):
+ # stimInfo = [path, geomExpr, relPath, field, expr_string]
+ stimObj = []
+ for i in dendCompts + spineCompts:
+ vclamp = make_vclamp( name = 'vclamp', parent = i.path )
+ moose.connect( i, 'VmOut', vclamp, 'sensedIn' )
+ moose.connect( vclamp, 'currentOut', i, 'injectMsg' )
+ stimObj.append( vclamp )
+
+ return stimObj
+
+ def _buildSynInputOnCompt( self, dendCompts, spineCompts, stimInfo, doPeriodic = False ):
+ # stimInfo = [path, geomExpr, relPath, field, expr_string]
+ # Here we hack geomExpr to use it for the syn weight. We assume it
+ # is just a number. In due course
+ # it should be possible to actually evaluate it according to geom.
+ synWeight = float( stimInfo[1] )
+ stimObj = []
+ for i in dendCompts + spineCompts:
+ path = i.path + '/' + stimInfo[2] + '/sh/synapse[0]'
+ if moose.exists( path ):
+ synInput = make_synInput( name='synInput', parent=path )
+ synInput.doPeriodic = doPeriodic
+ moose.element(path).weight = synWeight
+ moose.connect( synInput, 'spikeOut', path, 'addSpike' )
+ stimObj.append( synInput )
+ return stimObj
+
################################################################
# Here we set up the distributions
################################################################
def buildPassiveDistrib( self ):
+ # [. path field expr [field expr]...]
+ # RM, RA, CM set specific values, per unit area etc.
+ # Ra, Ra, Cm set absolute values.
+ # Also does Em, Ek, initVm
+ # Expression can use p, g, L, len, dia, maxP, maxG, maxL.
temp = []
for i in self.passiveDistrib:
temp.extend( i )
@@ -444,6 +528,14 @@ class rdesigneur:
self.elecid.spineDistribution = temp
def buildChemDistrib( self ):
+ # Orig format [chem, elecPath, install, expr]
+ # where chem and install were not being used.
+ # Modified format [chemLibPath, elecPath, newChemName, expr]
+ # chemLibPath is name of chemCompt or even group on library
+ # If chemLibPath has multiple compts on it, then the smaller ones
+ # become endoMeshes, scaled as per original.
+ # As a backward compatibility hack, if the newChemName == 'install'
+ # we use the default naming.
for i in self.chemDistrib:
pair = i[1] + " " + i[3]
# Assign any other params. Possibly the first param should
@@ -515,11 +607,11 @@ class rdesigneur:
# Put in stuff to go through fields if the target is a chem object
field = plotSpec[3]
if not field in knownFields:
- print("Warning: Rdesigneur::_parseComptField: Unknown field '", field, "'")
+ print("Warning: Rdesigneur::_parseComptField: Unknown field '{}'".format( field ) )
return (), ""
kf = knownFields[field] # Find the field to decide type.
- if ( kf[0] == 'CaConcBase' or kf[0] == 'ChanBase' or kf[0] == 'NMDAChan' ):
+ if ( kf[0] == 'CaConcBase' or kf[0] == 'ChanBase' or kf[0] == 'NMDAChan' or kf[0] == 'VClamp' ):
objList = self._collapseElistToPathAndClass( comptList, plotSpec[2], kf[0] )
return objList, kf[1]
elif (field == 'n' or field == 'conc' or field == 'volume' ):
@@ -561,6 +653,9 @@ class rdesigneur:
'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ),
'spikeTime':('CompartmentBase', 'getVm', 1, 'Spike Times (s)'),
'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ),
+ 'Cm':('CompartmentBase', 'getCm', 1e12, 'Memb. capacitance (pF)' ),
+ 'Rm':('CompartmentBase', 'getRm', 1e-9, 'Memb. Res (GOhm)' ),
+ 'Ra':('CompartmentBase', 'getRa', 1e-6, 'Axial Res (MOhm)' ),
'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ),
'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ),
'modulation':('ChanBase', 'getModulation', 1, 'chan modulation (unitless)' ),
@@ -570,7 +665,8 @@ class rdesigneur:
'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ),
'n':('PoolBase', 'getN', 1, '# of molecules'),
'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ),
- 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' )
+ 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ),
+ 'current':('VClamp', 'getCurrent', 1e9, 'Holding Current (nA)')
}
graphs = moose.Neutral( self.modelPath + '/graphs' )
dummy = moose.element( '/' )
@@ -645,13 +741,25 @@ class rdesigneur:
################################################################
# Here we display the plots and moogli
################################################################
- def displayMoogli( self, moogliDt, runtime, rotation = math.pi/500.0, fullscreen = False):
- rmoogli.displayMoogli( self, moogliDt, runtime, rotation, fullscreen )
- def display( self ):
+ def displayMoogli( self, moogliDt, runtime, rotation = math.pi/500.0, fullscreen = False, block = True, azim = 0.0, elev = 0.0 ):
+ rmoogli.displayMoogli( self, moogliDt, runtime, rotation = rotation, fullscreen = fullscreen, azim = azim, elev = elev )
+ pr = moose.PyRun( '/model/updateMoogli' )
+
+ pr.runString = '''
+import rdesigneur.rmoogli
+rdesigneur.rmoogli.updateMoogliViewer()
+'''
+ moose.setClock( pr.tick, moogliDt )
+ moose.reinit()
+ moose.start( runtime )
+ if block:
+ self.display( len( self.moogNames ) + 1 )
+
+ def display( self, startIndex = 0 ):
import matplotlib.pyplot as plt
for i in self.plotNames:
- plt.figure( i[2] )
+ plt.figure( i[2] + startIndex )
plt.title( i[1] )
plt.xlabel( "Time (s)" )
plt.ylabel( i[4] )
@@ -670,7 +778,7 @@ class rdesigneur:
plt.plot( t, j.vector * i[3] )
if len( self.moogList ) > 0:
plt.ion()
- plt.show(block=True)
+ plt.show( block=True )
#This calls the _save function which saves only if the filenames have been specified
self._save()
@@ -689,6 +797,9 @@ class rdesigneur:
knownFields = {
'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ),
+ 'Cm':('CompartmentBase', 'getCm', 1e12, 'Memb. capacitance (pF)' ),
+ 'Rm':('CompartmentBase', 'getRm', 1e-9, 'Memb. Res (GOhm)' ),
+ 'Ra':('CompartmentBase', 'getRa', 1e-6, 'Axial Res (MOhm)' ),
'spikeTime':('CompartmentBase', 'getVm', 1, 'Spike Times (s)'),
'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ),
'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ),
@@ -700,7 +811,8 @@ class rdesigneur:
'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ),
'n':('PoolBase', 'getN', 1, '# of molecules'),
'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ),
- 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' )
+ 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ),
+ 'current':('VClamp', 'getCurrent', 1e9, 'Holding Current (nA)')
}
save_graphs = moose.Neutral( self.modelPath + '/save_graphs' )
@@ -898,24 +1010,41 @@ class rdesigneur:
'inject':('CompartmentBase', 'setInject'),
'Ca':('CaConcBase', 'getCa'),
'n':('PoolBase', 'setN'),
- 'conc':('PoolBase', 'setConc')
+ 'conc':('PoolBase', 'setConc'),
+ 'vclamp':('CompartmentBase', 'setInject'),
+ 'randsyn':('SynChan', 'addSpike'),
+ 'periodicsyn':('SynChan', 'addSpike')
}
stims = moose.Neutral( self.modelPath + '/stims' )
k = 0
+ # Stimlist = [path, geomExpr, relPath, field, expr_string]
for i in self.stimList:
pair = i[0] + " " + i[1]
dendCompts = self.elecid.compartmentsFromExpression[ pair ]
spineCompts = self.elecid.spinesFromExpression[ pair ]
- stimObj, stimField = self._parseComptField( dendCompts, i, knownFields )
- stimObj2, stimField2 = self._parseComptField( spineCompts, i, knownFields )
- assert( stimField == stimField2 )
- stimObj3 = stimObj + stimObj2
+ #print( "pair = {}, numcompts = {},{} ".format( pair, len( dendCompts), len( spineCompts ) ) )
+ if i[3] == 'vclamp':
+ stimObj3 = self._buildVclampOnCompt( dendCompts, spineCompts, i )
+ stimField = 'commandIn'
+ elif i[3] == 'randsyn':
+ stimObj3 = self._buildSynInputOnCompt( dendCompts, spineCompts, i )
+ stimField = 'setRate'
+ elif i[3] == 'periodicsyn':
+ stimObj3 = self._buildSynInputOnCompt( dendCompts, spineCompts, i, doPeriodic = True )
+ stimField = 'setRate'
+ else:
+ stimObj, stimField = self._parseComptField( dendCompts, i, knownFields )
+ stimObj2, stimField2 = self._parseComptField( spineCompts, i, knownFields )
+ assert( stimField == stimField2 )
+ stimObj3 = stimObj + stimObj2
numStim = len( stimObj3 )
if numStim > 0:
funcname = stims.path + '/stim' + str(k)
k += 1
func = moose.Function( funcname )
func.expr = i[4]
+ if i[3] == 'vclamp': # Hack to clean up initial condition
+ func.doEvalAtReinit = 1
for q in stimObj3:
moose.connect( func, 'valueOut', q, stimField )
@@ -1189,7 +1318,7 @@ class rdesigneur:
return
if not hasattr( self, 'dendCompt' ):
raise BuildError( "configureSolvers: no chem meshes defined." )
- fixXreacs.fixXreacs( self.modelPath )
+ fixXreacs.fixXreacs( self.chemid.path )
dmksolve = moose.Ksolve( self.dendCompt.path + '/ksolve' )
dmdsolve = moose.Dsolve( self.dendCompt.path + '/dsolve' )
dmstoich = moose.Stoich( self.dendCompt.path + '/stoich' )
@@ -1241,13 +1370,12 @@ class rdesigneur:
return
# Sort comptlist in decreasing order of volume
sortedComptlist = sorted( comptlist, key=lambda x: -x.volume )
- if ( len( sortedComptlist ) != 3 ):
- print("loadChem: Require 3 chem compartments, have: ",\
- len( sortedComptlist ))
- return False
- sortedComptlist[0].name = 'dend'
- sortedComptlist[1].name = 'spine'
- sortedComptlist[2].name = 'psd'
+ if ( len( sortedComptlist ) >= 1 ):
+ sortedComptlist[0].name = 'dend'
+ if ( len( sortedComptlist ) >= 2 ):
+ sortedComptlist[1].name = 'spine'
+ if ( len( sortedComptlist ) >= 3 ):
+ sortedComptlist[2].name = 'psd'
################################################################
diff --git a/moose-core/python/rdesigneur/rdesigneurProtos.py b/moose-core/python/rdesigneur/rdesigneurProtos.py
index 1c22ea0cbc68a695f7c653cf0f4a0bba3e023231..ec5d705936d16a3fe624afa12dead0f83dc99e7d 100644
--- a/moose-core/python/rdesigneur/rdesigneurProtos.py
+++ b/moose-core/python/rdesigneur/rdesigneurProtos.py
@@ -49,7 +49,10 @@ import moose
import math
from moose import utils
-EREST_ACT = -70e-3
+EREST_ACT = -0.060
+ECA = 0.080
+EK = -0.075
+SOMA_A = 3.32e-9
per_ms = 1e3
PI = 3.14159265359
FaradayConst = 96485.3365 # Coulomb/mol
@@ -106,6 +109,329 @@ def make_HH_K(name = 'HH_K', parent='/library', vmin=-120e-3, vmax=40e-3, vdivs=
k.tick = -1
return k
+#/========================================================================
+#/ Tabchannel Na Hippocampal cell channel
+#/========================================================================
+def make_Na( name ):
+ if moose.exists( '/library/' + name ):
+ return
+ Na = moose.HHChannel( '/library/' + name )
+ Na.Ek = 0.055 # V
+ Na.Gbar = 300 * SOMA_A # S
+ Na.Gk = 0 # S
+ Na.Xpower = 2
+ Na.Ypower = 1
+ Na.Zpower = 0
+
+ xgate = moose.element( Na.path + '/gateX' )
+ xA = np.array( [ 320e3 * (0.0131 + EREST_ACT),
+ -320e3, -1.0, -1.0 * (0.0131 + EREST_ACT), -0.004,
+ -280e3 * (0.0401 + EREST_ACT), 280e3, -1.0,
+ -1.0 * (0.0401 + EREST_ACT), 5.0e-3,
+ 3000, -0.1, 0.05 ] )
+ xgate.alphaParms = xA
+
+ ygate = moose.element( Na.path + '/gateY' )
+ yA = np.array( [ 128.0, 0.0, 0.0, -1.0 * (0.017 + EREST_ACT), 0.018,
+ 4.0e3, 0.0, 1.0, -1.0 * (0.040 + EREST_ACT), -5.0e-3,
+ 3000, -0.1, 0.05 ] )
+ ygate.alphaParms = yA
+ return Na
+
+#========================================================================
+# Tabchannel K(DR) Hippocampal cell channel
+#========================================================================
+def make_K_DR( name ):
+ if moose.exists( '/library/' + name ):
+ return
+ K_DR = moose.HHChannel( '/library/' + name )
+ K_DR.Ek = EK # V
+ K_DR.Gbar = 150 * SOMA_A # S
+ K_DR.Gk = 0 # S
+ K_DR.Xpower = 1
+ K_DR.Ypower = 0
+ K_DR.Zpower = 0
+
+ xgate = moose.element( K_DR.path + '/gateX' )
+ xA = np.array( [ 16e3 * (0.0351 + EREST_ACT),
+ -16e3, -1.0, -1.0 * (0.0351 + EREST_ACT), -0.005,
+ 250, 0.0, 0.0, -1.0 * (0.02 + EREST_ACT), 0.04,
+ 3000, -0.1, 0.05 ] )
+ xgate.alphaParms = xA
+ return K_DR
+
+#========================================================================
+# Tabchannel K(A) Hippocampal cell channel
+#========================================================================
+def make_K_A( name ):
+ if moose.exists( '/library/' + name ):
+ return
+ K_A = moose.HHChannel( '/library/' + name )
+ K_A.Ek = EK # V
+ K_A.Gbar = 50 * SOMA_A # S
+ K_A.Gk = 0 # S
+ K_A.Xpower = 1
+ K_A.Ypower = 1
+ K_A.Zpower = 0
+
+ xgate = moose.element( K_A.path + '/gateX' )
+ xA = np.array( [ 20e3 * (0.0131 + EREST_ACT),
+ -20e3, -1.0, -1.0 * (0.0131 + EREST_ACT), -0.01,
+ -17.5e3 * (0.0401 + EREST_ACT),
+ 17.5e3, -1.0, -1.0 * (0.0401 + EREST_ACT), 0.01,
+ 3000, -0.1, 0.05 ] )
+ xgate.alphaParms = xA
+
+ ygate = moose.element( K_A.path + '/gateY' )
+ yA = np.array( [ 1.6, 0.0, 0.0, 0.013 - EREST_ACT, 0.018,
+ 50.0, 0.0, 1.0, -1.0 * (0.0101 + EREST_ACT), -0.005,
+ 3000, -0.1, 0.05 ] )
+ ygate.alphaParms = yA
+
+#//========================================================================
+#// Tabulated Ca Channel
+#//========================================================================
+
+def make_Ca( name ):
+ if moose.exists( '/library/' + name):
+ return
+ Ca = moose.HHChannel( '/library/' + name )
+ Ca.Ek = ECA
+ Ca.Gbar = 40 * SOMA_A
+ Ca.Gk = 0
+ Ca.Xpower = 2
+ Ca.Ypower = 1
+ Ca.Zpower = 0
+
+ xgate = moose.element( Ca.path + '/gateX' )
+ xA = np.array( [ 1.6e3, 0, 1.0, -1.0 * (0.065 + EREST_ACT), -0.01389, -20e3 * (0.0511 + EREST_ACT), 20e3, -1.0, -1.0 * (0.0511 + EREST_ACT), 5.0e-3, 3000, -0.1, 0.05 ] )
+# xgate.min = -0.1
+# xgate.max = 0.05
+# xgate.divs = 3000
+#// Converting Traub's expressions for the gCa/s alpha and beta functions
+#// to SI units and entering the A, B, C, D and F parameters, we get:
+# xgate.alpha( 1.6e3, 0, 1.0, -1.0 * (0.065 + EREST_ACT), -0.01389 )
+# xgate.beta( -20e3 * (0.0511 + EREST_ACT), 20e3, -1.0, -1.0 * (0.0511 + EREST_ACT), 5.0e-3 )
+ #xgate.setupAlpha( xA )
+ xgate.alphaParms = xA
+
+
+# The Y gate (gCa/r) is not quite of this form. For V > EREST_ACT, alpha =
+# 5*{exp({-50*(V - EREST_ACT)})}. Otherwise, alpha = 5. Over the entire
+# range, alpha + beta = 5. To create the Y_A and Y_B tables, we use some
+# of the pieces of the setupalpha function.
+ ygate = moose.element( Ca.path + '/gateY' )
+ ygate.min = -0.1
+ ygate.max = 0.05
+ ygate.divs = 3000
+ yA = np.zeros( (ygate.divs + 1), dtype=float)
+ yB = np.zeros( (ygate.divs + 1), dtype=float)
+
+
+#Fill the Y_A table with alpha values and the Y_B table with (alpha+beta)
+ dx = (ygate.max - ygate.min)/ygate.divs
+ x = ygate.min
+ for i in range( ygate.divs + 1 ):
+ if ( x > EREST_ACT):
+ yA[i] = 5.0 * math.exp( -50 * (x - EREST_ACT) )
+ else:
+ yA[i] = 0.0
+ #yB[i] = 6.0 - yA[i]
+ yB[i] = 5.0
+ x += dx
+ ygate.tableA = yA
+ ygate.tableB = yB
+# Tell the cell reader that the current from this channel must be fed into
+# the Ca_conc pool of calcium.
+ addmsg1 = moose.Mstring( Ca.path + '/addmsg1' )
+ addmsg1.value = '. IkOut ../Ca_conc current'
+# in some compartments, whe have an NMDA_Ca_conc object to put the current
+# into.
+ addmsg2 = moose.Mstring( Ca.path + '/addmsg2' )
+ addmsg2.value = '. IkOut ../NMDA_Ca_conc current'
+# As we typically use the cell reader to create copies of these prototype
+#elements in one or more compartments, we need some way to be sure that the
+#needed messages are established. Although the cell reader has enough
+#information to create the messages which link compartments to their channels
+#and to other adjacent compartments, it most be provided with the information
+#needed to establish additional messages. This is done by placing the
+#message string in a user-defined field of one of the elements which is
+#involved in the message. The cell reader recognizes the added object names
+#"addmsg1", "addmsg2", etc. as indicating that they are to be
+#evaluated and used to set up messages. The paths are relative to the
+#element which contains the message string in its added field. Thus,
+#"../Ca_conc" refers to the sibling element Ca_conc and "."
+#refers to the Ca element itself.
+
+
+#/*************************************************************************
+#Next, we need an element to take the Calcium current calculated by the Ca
+#channel and convert it to the Ca concentration. The "Ca_concen" object
+#solves the equation dC/dt = B*I_Ca - C/tau, and sets Ca = Ca_base + C. As
+#it is easy to make mistakes in units when using this Calcium diffusion
+#equation, the units used here merit some discussion.
+
+#With Ca_base = 0, this corresponds to Traub's diffusion equation for
+#concentration, except that the sign of the current term here is positive, as
+#GENESIS uses the convention that I_Ca is the current flowing INTO the
+#compartment through the channel. In SI units, the concentration is usually
+#expressed in moles/m^3 (which equals millimoles/liter), and the units of B
+#are chosen so that B = 1/(ion_charge * Faraday * volume). Current is
+#expressed in amperes and one Faraday = 96487 coulombs. However, in this
+#case, Traub expresses the concentration in arbitrary units, current in
+#microamps and uses tau = 13.33 msec. If we use the same concentration units,
+#but express current in amperes and tau in seconds, our B constant is then
+#10^12 times the constant (called "phi") used in the paper. The actual value
+#used will be typically be determined by the cell reader from the cell
+#parameter file. However, for the prototype channel we wlll use Traub's
+#corrected value for the soma. (An error in the paper gives it as 17,402
+#rather than 17.402.) In our units, this will be 17.402e12.
+
+#*************************************************************************/
+
+
+#//========================================================================
+#// Ca conc
+#//========================================================================
+
+def make_Ca_conc( name ):
+ if moose.exists( '/library/' + name ):
+ return
+ conc = moose.CaConc( '/library/tempName' )
+ conc.name = name
+ conc.tau = 0.013333 # sec
+ conc.B = 17.402e12 # Curr to conc conversion for soma
+ conc.Ca_base = 0.00000
+
+#This Ca_concen element should receive a message from any calcium channels
+# with the current going through the channel. Here we have this specified
+# in the Ca channel, with the idea that more than one channel might
+# contribute Ca ions to this calcium pool. In the original GENESIS file
+# this was specified here in make_Ca_conc.
+
+#========================================================================
+# Tabulated Ca-dependent K AHP Channel
+#========================================================================
+
+# This is a tabchannel which gets the calcium concentration from Ca_conc
+# in order to calculate the activation of its Z gate. It is set up much
+# like the Ca channel, except that the A and B tables have values which are
+# functions of concentration, instead of voltage.
+
+def make_K_AHP( name ):
+ if moose.exists( '/library/' + name ):
+ return
+ K_AHP = moose.HHChannel( '/library/' + name )
+ K_AHP.Ek = EK # V
+ K_AHP.Gbar = 8 * SOMA_A # S
+ K_AHP.Gk = 0 # S
+ K_AHP.Xpower = 0
+ K_AHP.Ypower = 0
+ K_AHP.Zpower = 1
+
+ zgate = moose.element( K_AHP.path + '/gateZ' )
+ xmax = 500.0
+ zgate.min = 0
+ zgate.max = xmax
+ zgate.divs = 3000
+ zA = np.zeros( (zgate.divs + 1), dtype=float)
+ zB = np.zeros( (zgate.divs + 1), dtype=float)
+ dx = (zgate.max - zgate.min)/zgate.divs
+ x = zgate.min
+ for i in range( zgate.divs + 1 ):
+ zA[i] = min( 0.02 * CA_SCALE * x, 10 )
+ zB[i] = 1.0
+ x = x + dx
+
+ zgate.tableA = zA
+ zgate.tableB = zB
+ addmsg1 = moose.Mstring( K_AHP.path + '/addmsg1' )
+ addmsg1.value = '../Ca_conc concOut . concen'
+# Use an added field to tell the cell reader to set up a message from the
+# Ca_Conc with concentration info, to the current K_AHP object.
+
+
+#//========================================================================
+#// Ca-dependent K Channel - K(C) - (vdep_channel with table and tabgate)
+#//========================================================================
+
+#The expression for the conductance of the potassium C-current channel has a
+#typical voltage and time dependent activation gate, where the time dependence
+#arises from the solution of a differential equation containing the rate
+#parameters alpha and beta. It is multiplied by a function of calcium
+#concentration that is given explicitly rather than being obtained from a
+#differential equation. Therefore, we need a way to multiply the activation
+#by a concentration dependent value which is determined from a lookup table.
+#This is accomplished by using the Z gate with the new tabchannel "instant"
+#field, introduced in GENESIS 2.2, to implement an "instantaneous" gate for
+#the multiplicative Ca-dependent factor in the conductance.
+
+def make_K_C( name ):
+ if moose.exists( '/library/' + name ):
+ return
+ K_C = moose.HHChannel( '/library/' + name )
+ K_C.Ek = EK # V
+ K_C.Gbar = 100.0 * SOMA_A # S
+ K_C.Gk = 0 # S
+ K_C.Xpower = 1
+ K_C.Zpower = 1
+ K_C.instant = 4 # Flag: 0x100 means Z gate is instant.
+ K_C.useConcentration = 1
+
+ # Now make a X-table for the voltage-dependent activation parameter.
+ xgate = moose.element( K_C.path + '/gateX' )
+ xgate.min = -0.1
+ xgate.max = 0.05
+ xgate.divs = 3000
+ xA = np.zeros( (xgate.divs + 1), dtype=float)
+ xB = np.zeros( (xgate.divs + 1), dtype=float)
+ dx = (xgate.max - xgate.min)/xgate.divs
+ x = xgate.min
+ for i in range( xgate.divs + 1 ):
+ alpha = 0.0
+ beta = 0.0
+ if (x < EREST_ACT + 0.05):
+ alpha = math.exp( 53.872 * (x - EREST_ACT) - 0.66835 ) / 0.018975
+ beta = 2000* (math.exp ( (EREST_ACT + 0.0065 - x)/0.027)) - alpha
+ else:
+ alpha = 2000 * math.exp( ( EREST_ACT + 0.0065 - x)/0.027 )
+ beta = 0.0
+ xA[i] = alpha
+ xB[i] = alpha + beta
+ x = x + dx
+ xgate.tableA = xA
+ xgate.tableB = xB
+
+# Create a table for the function of concentration, allowing a
+# concentration range of 0 to 200, with 3000 divisions. This is done
+# using the Z gate, which can receive a CONCEN message. By using
+# the "instant" flag, the A and B tables are evaluated as lookup tables,
+# rather than being used in a differential equation.
+ zgate = moose.element( K_C.path + '/gateZ' )
+ zgate.min = 0.0
+ xmax = 150.0
+ zgate.max = xmax
+ zgate.divs = 3000
+ zA = np.zeros( (zgate.divs + 1), dtype=float)
+ zB = np.zeros( (zgate.divs + 1), dtype=float)
+ dx = ( zgate.max - zgate.min)/ zgate.divs
+ x = zgate.min
+ #CaScale = 100000.0 / 250.0e-3
+ for i in range( zgate.divs + 1 ):
+ zA[i] = min( 1000.0, x * CA_SCALE / (250 * xmax ) )
+ zB[i] = 1000.0
+ x += dx
+ zgate.tableA = zA
+ zgate.tableB = zB
+
+# Now we need to provide for messages that link to external elements.
+# The message that sends the Ca concentration to the Z gate tables is stored
+# in an added field of the channel, so that it may be found by the cell
+# reader.
+ addmsg1 = moose.Mstring( K_C.path + '/addmsg1' )
+ addmsg1.value = '../Ca_conc concOut . concen'
+
+
#========================================================================
# SynChan: Glu receptor
#========================================================================
@@ -131,14 +457,15 @@ def make_GABA( name ):
if moose.exists( '/library/' + name ):
return
GABA = moose.SynChan( '/library/' + name )
- GABA.Ek = EK + 10.0e-3
- GABA.tau1 = 4.0e-3
- GABA.tau2 = 9.0e-3
+ GABA.Ek = -0.065 # V
+ GABA.tau1 = 4.0e-3 # s
+ GABA.tau2 = 9.0e-3 # s
sh = moose.SimpleSynHandler( GABA.path + '/sh' )
moose.connect( sh, 'activationOut', GABA, 'activation' )
sh.numSynapses = 1
sh.synapse[0].weight = 1
+#========================================================================
def makeChemOscillator( name = 'osc', parent = '/library' ):
model = moose.Neutral( parent + '/' + name )
@@ -248,7 +575,7 @@ def transformNMDAR( path ):
# compartment, it builds it on 'pa'. It places the compartment
# on the end of 'prev', and at 0,0,0 otherwise.
-def buildCompt( pa, name, RM = 1.0, RA = 1.0, CM = 0.01, dia = 1.0e-6, x = 0.0, y = 0.0, z = 0.0, dx = 10e-6, dy = 0.0, dz = 0.0 ):
+def buildCompt( pa, name, RM = 1.0, RA = 1.0, CM = 0.01, dia = 1.0e-6, x = 0.0, y = 0.0, z = 0.0, dx = 10e-6, dy = 0.0, dz = 0.0, Em = -0.065, initVm = -0.065 ):
length = np.sqrt( dx * dx + dy * dy + dz * dz )
compt = moose.Compartment( pa.path + '/' + name )
compt.x0 = x
@@ -264,6 +591,8 @@ def buildCompt( pa, name, RM = 1.0, RA = 1.0, CM = 0.01, dia = 1.0e-6, x = 0.0,
compt.Ra = length * RA / xa
compt.Rm = RM / sa
compt.Cm = CM * sa
+ compt.Em = Em
+ compt.initVm = initVm
return compt
def buildComptWrapper( pa, name, length, dia, xoffset, RM, RA, CM ):
@@ -289,8 +618,6 @@ def buildSyn( name, compt, Ek, tau1, tau2, Gbar, CM ):
# Utility function, borrowed from proto18.py, for making an LCa channel.
# Based on Traub's 91 model, I believe.
def make_LCa( name = 'LCa', parent = '/library' ):
- EREST_ACT = -0.060 #/* hippocampal cell resting potl */
- ECA = 0.140 + EREST_ACT #// 0.080
if moose.exists( parent + '/' + name ):
return
Ca = moose.HHChannel( parent + '/' + name )
@@ -327,80 +654,39 @@ def make_LCa( name = 'LCa', parent = '/library' ):
return Ca
######################################################################
-# Derived from : squid/electronics.py
-# Description:
-# Author: Subhasis Ray
-# Maintainer:
-# Created: Wed Feb 22 00:53:38 2012 (+0530)
-# Version:
-# Last-Updated: Fri May 04 16:35:40 2018 (+0530)
-# By: Upi
-# Update #: 221
+def make_vclamp( name = 'Vclamp', parent = '/library' ):
+ if moose.exists( '/library/' + name ):
+ return
+ vclamp = moose.VClamp( parent + '/' + name )
+ vclamp.mode = 0 # Default. could try 1, 2 as well
+ vclamp.tau = 0.2e-3 # lowpass filter for command voltage input
+ vclamp.ti = 20e-6 # Integral time
+ vclamp.td = 5e-6 # Differential time. Should it be >= dt?
+ vclamp.gain = 0.00005 # Gain of vclamp ckt.
-# Change log:
-#
-# 2012-02-22 23:22:30 (+0530) Subha - the circuitry put in a class.
-# 2018-05-04 23:22:30 (+0530) Upi - Adapted for Rdesigneur
-#
+ # Connect voltage clamp circuitry
+ addmsg1 = moose.Mstring( vclamp.path + '/addmsg1' )
+ addmsg1.value = '. currentOut .. injectMsg'
+ addmsg2 = moose.Mstring( vclamp.path + '/addmsg2' )
+ addmsg2.value = '.. VmOut . sensedIn'
-# Code:
+ return vclamp
-class ClampCircuit(moose.Neutral):
- """Container for a Voltage-Clamp/Current clamp circuit."""
- defaults = {
- 'level1': 25.0e-3,
- 'width1': 50.0e-3,
- 'delay1': 2.0e-3,
- 'delay2': 1e3,
- 'trigMode': 0,
- 'delay3': 1e6
- }
- def __init__(self, path ):
- moose.Neutral.__init__(self, path)
- '''
- self.pulsegen = moose.PulseGen(path+"/pulse") # holding voltage/current generator
- self.pulsegen.count = 2
- self.pulsegen.baseLevel = -65.0e-3
- self.pulsegen.firstLevel = -40.0e-3
- self.pulsegen.firstWidth = 50.0e-3
- self.pulsegen.firstDelay = 2.0e-3
- self.pulsegen.secondDelay = 0.0
- self.pulsegen.trigMode = 2
- self.gate = moose.PulseGen(path+"/gate") # holding voltage/current generator
- self.gate.level[0] = 1.0
- self.gate.delay[0] = 0.0
- self.gate.width[0] = 1e3
- moose.connect(self.gate, 'output', self.pulsegen, 'input')
- '''
- self.lowpass = moose.RC(path+"/lowpass") # lowpass filter
- self.lowpass.R = 1.0
- self.lowpass.C = 0.03
- self.vclamp = moose.DiffAmp(path+"/vclamp")
- self.vclamp.gain = 1.0
- self.vclamp.saturation = 1e10
- self.iclamp = moose.DiffAmp(path+"/iclamp")
- self.iclamp.gain = 0.0
- self.iclamp.saturation = 1e10
- self.pid = moose.PIDController(path+"/pid")
- self.pid.gain = 0.5
- self.pid.tauI = 0.02e-3
- self.pid.tauD = 0.005e-3
- self.pid.saturation = 1e7
- # Connect voltage clamp circuitry
- #moose.connect(self.pulsegen, "output", self.lowpass, "injectIn")
- moose.connect(self.lowpass, "output", self.vclamp, "plusIn")
- moose.connect(self.vclamp, "output", self.pid, "commandIn")
- #moose.connect(compartment, "VmOut", self.pid, "sensedIn")
- #moose.connect(self.pid, "output", compartment, "injectMsg")
- addmsg1 = moose.Mstring( path + '/addmsg1' )
- addmsg1.value = './pid output .. injectMsg'
- addmsg2 = moose.Mstring( path + '/addmsg2' )
- addmsg2.value = '.. VmOut ./pid sensedIn'
+######################################################################
-def make_vclamp( name = 'Vclamp', parent = '/library' ):
+def make_synInput( name = 'RandSpike', parent = '/library' ):
if moose.exists( '/library/' + name ):
return
- vclamp = ClampCircuit( parent + '/' + name )
+ rs = moose.RandSpike( parent + '/' + name + '_rs' )
+ rs.rate = 0 # mean firing rate
+ rs.refractT = 5e-3 # 5 ms.
+
+
+ # Connect rand spike to channel that it is sitting on.
+ addmsg1 = moose.Mstring( rs.path + '/addmsg1' )
+ addmsg1.value = '. spikeOut ../sh/synapse[0] addSpike'
+
+ return rs
######################################################################
diff --git a/moose-core/python/rdesigneur/rmoogli.py b/moose-core/python/rdesigneur/rmoogli.py
index fcb1a1209f03b4f542f63b0fee38f169588dabb0..f98549c0b2f9b4109774da64d57a68ebbb9396a5 100644
--- a/moose-core/python/rdesigneur/rmoogli.py
+++ b/moose-core/python/rdesigneur/rmoogli.py
@@ -1,178 +1,53 @@
-# -*- coding: utf-8 -*-
-#########################################################################
-## rdesigneur0_4.py ---
-## This program is part of 'MOOSE', the
-## Messaging Object Oriented Simulation Environment.
-## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
-## It is made available under the terms of the
-## GNU General Public License version 2 or later.
-## See the file COPYING.LIB for the full notice.
-#########################################################################
-
-import math
-import sys
-import moose
-import os
-
-# Check if DISPLAY environment variable is properly set. If not, warn the user
-# and continue.
-hasDisplay = True
-display = os.environ.get('DISPLAY', '' )
-if not display:
- hasDisplay = False
- print( "Warning: Environment variable DISPLAY is not set."
- " Did you forget to pass -X or -Y switch to ssh command?\n"
- "Anyway, MOOSE will continue without graphics.\n"
- )
-
-hasMoogli = True
-
-if hasDisplay:
- try:
- import matplotlib
- from PyQt4 import QtGui
- import moogli
- import moogli.extensions.moose
- app = QtGui.QApplication(sys.argv)
- except Exception as e:
- print( 'Warning: Moogli not found. All moogli calls will use dummy functions' )
- hasMoogli = False
-
-
-runtime = 0.0
-moogliDt = 1.0
-rotation = math.pi / 500.0
-
-def getComptParent( obj ):
- k = moose.element(obj)
- while not k.isA[ "CompartmentBase" ]:
- if k == moose.element( '/' ):
- return k.path
- k = moose.element( k.parent )
- return k.path
-
-#######################################################################
-## Here we set up the callback functions for the viewer
-def prelude( view ):
- view.home()
- view.pitch( math.pi / 2.0 )
- view.zoom( 0.05 )
- #network.groups["soma"].set( "color", moogli.colors.RED )
-
-# This func is used for the first viewer, it has to handle advancing time.
-def interlude( view ):
- moose.start( moogliDt )
- val = [ moose.getField( i, view.mooField, "double" ) * view.mooScale for i in view.mooObj ]
- view.mooGroup.set("color", val, view.mapper)
- view.yaw( rotation )
- #print moogliDt, len( val ), runtime
- if moose.element("/clock").currentTime >= runtime:
- view.stop()
-
-# This func is used for later viewers, that don't handle advancing time.
-def interlude2( view ):
- val = [ moose.getField( i, view.mooField, "double" ) * view.mooScale for i in view.mooObj ]
-
- view.mooGroup.set("color", val, view.mapper)
- view.yaw( rotation )
- if moose.element("/clock").currentTime >= runtime:
- view.stop()
-
-def postlude( view ):
- view.rd.display()
-
-def makeMoogli( rd, mooObj, moogliEntry, fieldInfo ):
- if not hasMoogli:
- return None
- mooField = moogliEntry[3]
+from __future__ import absolute_import, print_function
+# rmoogli.py: rdesigneur Moogli interface
+# This is a fallback version designed to work with moogul but using almost
+# the same interface as far as rdesigneur is concerned.
+# Put in because the GL versions like moogli need all sorts of difficult
+# libraries and dependencies.
+# Copyright (C) Upinder S. Bhalla NCBS 2018
+# This program is licensed under the GNU Public License version 3.
+
+import rdesigneur.moogul as moogul
+mooViews = []
+
+def makeMoogli( rd, mooObj, args, fieldInfo ):
+ #mooObj is currently poorly handled. Ideally it should simply be a
+ # vector of objects to plot, each with coords. This could be readily
+ # used to build up the display in a neutral manner.
+ # Another version is to have a vector of objects to plot, only as a
+ # way to get at their fields. We would separately need mapping to
+ # neuron compartments and index along length.
+ # Cleaner still would be to have the C code give a vector of values
+ # For now it means something different for chem and elec displays.
+ #moogliEntry = [elecPath,bool,whichObjToDisplay,FieldToDisplay,titleForDisplay,rangeMin,rangeMax]
+ mooField = args[3]
+ relObjPath = args[2]
numMoogli = len( mooObj )
- network = moogli.extensions.moose.read( path = rd.elecid.path, vertices=15)
- #print len( network.groups["spine"].shapes )
- #print len( network.groups["dendrite"].shapes )
- #print len( network.groups["soma"].shapes )
- #soma = network.groups["soma"].shapes[ '/model/elec/soma']
- #print network.groups["soma"].shapes
-
- soma = network.groups["soma"].shapes[ rd.elecid.path + '/soma[0]']
- if ( mooField == 'n' or mooField == 'conc' ):
- updateGroup = soma.subdivide( numMoogli )
- displayObj = mooObj
- else:
- shell = moose.element( '/' )
- displayObj = [i for i in mooObj if i != shell ]
- cpa = [getComptParent( i ) for i in displayObj ]
- updateGroup = moogli.Group( "update" )
- updateShapes = [network.shapes[i] for i in cpa]
- #print "########### Len( cpa, mooObj ) = ", len( cpa ), len( mooObj ), len( updateShapes )
- updateGroup.attach_shapes( updateShapes )
-
- normalizer = moogli.utilities.normalizer(
- moogliEntry[5], moogliEntry[6],
- clipleft =True,
- clipright = True )
- colormap = moogli.colors.MatplotlibColorMap(matplotlib.cm.rainbow)
- mapper = moogli.utilities.mapper(colormap, normalizer)
-
- viewer = moogli.Viewer("Viewer")
- viewer.setWindowTitle( moogliEntry[4] )
- if ( mooField == 'n' or mooField == 'conc' ):
- viewer.attach_shapes( updateGroup.shapes.values())
- viewer.detach_shape(soma)
- else:
- viewer.attach_shapes(network.shapes.values())
-
- if len( rd.moogNames ) == 0:
- view = moogli.View("main-view",
- prelude=prelude,
- interlude=interlude,
- postlude = postlude)
+ viewer = moogul.MooView()
+ if mooField == 'n' or mooField == 'conc':
+ #moogul.updateDiffCoords( mooObj )
+ reacSystem = moogul.MooReacSystem( mooObj, fieldInfo,
+ field = mooField, relativeObj = relObjPath,
+ valMin = args[5], valMax = args[6] )
+ viewer.addDrawable( reacSystem )
else:
- view = moogli.View("main-view",
- prelude=prelude,
- interlude=interlude2)
+ neuron = moogul.MooNeuron( rd.elecid, fieldInfo,
+ field = mooField, relativeObj = relObjPath,
+ valMin = args[5], valMax = args[6] )
+ viewer.addDrawable( neuron )
- cb = moogli.widgets.ColorBar(id="cb",
- title=fieldInfo[3],
- text_color=moogli.colors.BLACK,
- position=moogli.geometry.Vec3f(0.975, 0.5, 0.0),
- size=moogli.geometry.Vec3f(0.30, 0.05, 0.0),
- text_font="/usr/share/fonts/truetype/ubuntu-font-family/Ubuntu-R.ttf",
- orientation=math.pi / 2.0,
- text_character_size=16,
- label_formatting_precision=0,
- colormap=moogli.colors.MatplotlibColorMap(matplotlib.cm.rainbow),
- color_resolution=100,
- scalar_range=moogli.geometry.Vec2f(
- moogliEntry[5],
- moogliEntry[6]))
- cb.set_num_labels(3)
- view.attach_color_bar(cb)
- view.rd = rd
- view.mooObj = displayObj
- view.mooGroup = updateGroup
- view.mooField = mooField
- view.mooScale = fieldInfo[2]
- view.mapper = mapper
- viewer.attach_view(view)
return viewer
-def displayMoogli( rd, _dt, _runtime, _rotation, fullscreen = False ):
- if not hasMoogli:
- return None
- global runtime
- global moogliDt
- global rotation
- runtime = _runtime
- moogliDt = _dt
- rotation = _rotation
+def updateMoogliViewer():
+ for i in mooViews:
+ i.updateValues()
+
+
+def displayMoogli( rd, _dt, _runtime, rotation = 0.0, fullscreen = False, azim = 0.0, elev = 0.0 ):
+ global mooViews
+ mooViews = rd.moogNames
for i in rd.moogNames:
- if fullscreen:
- i.showMaximized()
- else:
- i.show()
- i.start()
- #viewer.showMaximized()
- #viewer.show()
- #viewer.start()
- return app.exec_()
+ i.firstDraw( rotation = rotation, azim = azim, elev = elev )
+ # rotation in radians/frame, azim, elev in radians.
+
diff --git a/moose-core/tests/python/test_kkit.py b/moose-core/tests/python/test_kkit.py
index 8ba4fd2b3f57eb643e4bf537695c5f4012b4ea5a..ba194d63c44f58abb5f0c1d8a3235bed5f42f92e 100644
--- a/moose-core/tests/python/test_kkit.py
+++ b/moose-core/tests/python/test_kkit.py
@@ -23,7 +23,9 @@ def main():
runtime = float( sys.argv[2] )
if ( len( sys.argv ) == 4 ):
solver = sys.argv[3]
- modelId = moose.loadModel( mfile, 'model', solver )
+ modelId = moose.loadModel( mfile, 'model')
+ moose.mooseaddChemSolver('model',solver)
+
# Increase volume so that the stochastic solver gssa
# gives an interesting output
#compt = moose.element( '/model/kinetics' )
diff --git a/moose-core/tests/python/test_negative_value_flag.py b/moose-core/tests/python/test_negative_value_flag.py
index 35c37ad7c83c01dc35faeaa28e6abb568cda2fc6..7b073975de42881f5abab9f8cb4f40c6f029e12f 100644
--- a/moose-core/tests/python/test_negative_value_flag.py
+++ b/moose-core/tests/python/test_negative_value_flag.py
@@ -35,7 +35,8 @@ def main():
if ( len( sys.argv ) == 4 ):
solver = sys.argv[3]
- modelId = moose.loadModel( mfile, 'model', solver )
+ modelId = moose.loadModel( mfile, 'model')
+ moose.mooseaddChemSolver('model',solver)
moose.element( '/model/kinetics/neuroNOS/nNOS.arg' ).concInit = 0.1
moose.reinit()
moose.start( runtime )