From db571897a7ad9506157681a4d151aaa289d7d3bf Mon Sep 17 00:00:00 2001
From: HarshaRani <ranishashi@gmail.com>
Date: Wed, 1 Aug 2018 11:58:41 +0530
Subject: [PATCH] Squashed 'moose-core/' changes from cddc3bd..1271dcd
1271dcd Merge pull request #278 from upibhalla/master
c400af3 Minor extension to rdesigneur to handle display of the passive electrical parameters of compartments: Rm, Ra, Cm. These will be scaled if spine geometry changes.
a77c611 Fixes to rdesigneur for NeuroMesh models. Key bugfix to Dsolve for NeuroMesh models. Some code cleanup to eliminate compiler warnings about uninitialized variables.
840021f Added lots of channel prototypes to rdesigneurProtos. Small fixes to rdesigneur for multiscale models.
c4f9fa5 Updates to Rdesigneur to handle synaptic input stimulus. Modified RandSpike with a doPeriodic flag, so it can be used for periodic stimuli as well.
cb2dd4b Merge branch 'master' of https://github.com/BhallaLab/moose-core
fc6ca10 Updates to rdesigneur to handle voltage clamp specification as a stimulus. Added doEvalAtReinit flag to Function so that it can initialize at a precomputed value rather than zero.
0445da5 Merge pull request #276 from upibhalla/master
48d1d4c Updates to rdesigneur to handle soma parameters and ball-and-stick model in cellProto argument
694fe5e Merge branch 'master' of https://github.com/BhallaLab/moose-core
180109a Changed Enz number-rate scaling to use the volume of the enzyme molecule, rather than the volume of the substrate. This usage gives results consistent with Genesis/kkit
f9de66e Merge pull request #273 from upibhalla/master
c2d1a6f Fixed grotesque bug with Dsolve.cpp which manifested only if pool names were of length 13 and they were diffusing.
b9448c0 Merge branch 'master' of https://github.com/BhallaLab/moose-core
bcd9ddb Fixes to Enz so that kcat of zero doesn't set up NaNs in Km and K1. Display updates to rdesigneur, rmoogli, moogul.
eeea183 Merge pull request #270 from upibhalla/master
5e840d2 Merge branch 'master' into master
9ac8d79 Merge pull request #272 from BhallaLab/dilawar-patch-2
8f276ef Update travis_prepare_osx.sh
2d00294 Update rmoogli.py
5e793c3 Update moogul.py
5108d7e Update rmoogli.py
42a4f3b Merge branch 'master' of https://github.com/BhallaLab/moose-core
8ea8b2e Fixes to moogli when run from outside the rdesigneur directory
c411242 minor cleanup on rdesigneur
ea99dcc Updated version of moogul, integrated now with rdesigneur
d81e796 Merge pull request #264 from hrani/master
bf0119b Update test_negative_value_flag.py
6e71d8e Update test_kkit.py
fdc7d66 cleaned up print stm
51dfc2c cleaned up print stm
b055a5c cleaned up wrt to python3 compilation problem
09da84b indentation corrected
00bc557 indentation corrected
b55acd9 Cleaned wrt to diffusion,clean up to remove in ReadKKit removed the setting up of solver, in mooseaddsolver() cleaned up in setting and deleting the solver, added utils.loadmodel function to call loadmodel function with setting for solver
9441bf9 Merge branch 'master' into master
08f0804 Merge pull request #266 from upibhalla/master
0ad961c Update rdesigneur.py
3ec65d7 Update rdesigneur.py
b6b58a4 Moved fixXreacs.py from python/rdesigneur directory to python directory
736f7f5 Merge branch 'master' of https://github.com/BhallaLab/moose-core
38a1d03 Bugfix for fixXreacs so that it works properly if there was already some text in the notes.
b8fcd11 Merge branch 'master' of http://github.com/hrani/moose-core
ad6925f Merge branch 'master' of https://github.com/BhallaLab/moose-core
9d5be9d Merge branch 'master' into master
cc13c50 travis failing with libgfortran not found, hopefully this should fix it
b246358 rename modelA and modelB to modelS and modelD, cleanup while copying the object at Neutral level
ad7236b enzcomplex pool cleaned-up
8fce41c indentation corrected while writting annotation for enzymecomplex
git-subtree-dir: moose-core
git-subtree-split: 1271dcd74baf95a973debcdef0fc8b348010bb78
---
.travis.yml | 4 +
.travis/travis_prepare_osx.sh | 2 +-
biophysics/MarkovSolver.cpp | 6 +-
biophysics/MatrixOps.cpp | 12 +-
biophysics/RandSpike.cpp | 51 +-
biophysics/RandSpike.h | 4 +
builtins/Function.cpp | 33 +-
builtins/Function.h | 5 +
builtins/HDF5WriterBase.cpp | 3 +-
diffusion/Dsolve.cpp | 17 +-
diffusion/Dsolve.h | 5 +-
kinetics/Enz.cpp | 18 +-
kinetics/MMenz.cpp | 2 +
kinetics/ReadKkit.cpp | 25 +-
kinetics/WriteKkit.cpp | 6 +-
ksolve/ZombieEnz.cpp | 7 +-
ksolve/ZombieMMenz.cpp | 4 +
python/{rdesigneur => }/fixXreacs.py | 2 +-
python/moose/SBML/readSBML.py | 33 +-
python/moose/SBML/writeSBML.py | 54 ++-
python/moose/chemMerge/merge.py | 383 +++++++--------
.../chemUtil/add_Delete_ChemicalSolver.py | 76 +++
python/moose/utils.py | 6 +-
python/rdesigneur/moogul.py | 314 +++++++++++++
python/rdesigneur/rdesigneur.py | 178 ++++++-
python/rdesigneur/rdesigneurProtos.py | 436 +++++++++++++++---
python/rdesigneur/rmoogli.py | 217 ++-------
tests/python/test_kkit.py | 4 +-
tests/python/test_negative_value_flag.py | 3 +-
29 files changed, 1371 insertions(+), 539 deletions(-)
rename python/{rdesigneur => }/fixXreacs.py (99%)
create mode 100644 python/rdesigneur/moogul.py
diff --git a/.travis.yml b/.travis.yml
index b5f3f688..a89ee791 100644
--- a/.travis.yml
+++ b/.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/.travis/travis_prepare_osx.sh b/.travis/travis_prepare_osx.sh
index 65bff743..d8aa29a1 100755
--- a/.travis/travis_prepare_osx.sh
+++ b/.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/biophysics/MarkovSolver.cpp b/biophysics/MarkovSolver.cpp
index 1ab4adc1..b5e3ddbb 100644
--- a/biophysics/MarkovSolver.cpp
+++ b/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/biophysics/MatrixOps.cpp b/biophysics/MatrixOps.cpp
index a983ef74..f0f478c3 100644
--- a/biophysics/MatrixOps.cpp
+++ b/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/biophysics/RandSpike.cpp b/biophysics/RandSpike.cpp
index 0c2e09ff..4678970f 100644
--- a/biophysics/RandSpike.cpp
+++ b/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/biophysics/RandSpike.h b/biophysics/RandSpike.h
index ea0cf595..2a7da3de 100644
--- a/biophysics/RandSpike.h
+++ b/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/builtins/Function.cpp b/builtins/Function.cpp
index 79413c1b..ea0c8ae8 100644
--- a/builtins/Function.cpp
+++ b/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/builtins/Function.h b/builtins/Function.h
index ed6c4517..83f2b903 100644
--- a/builtins/Function.h
+++ b/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/builtins/HDF5WriterBase.cpp b/builtins/HDF5WriterBase.cpp
index 80c47cbe..747a2d8d 100644
--- a/builtins/HDF5WriterBase.cpp
+++ b/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/diffusion/Dsolve.cpp b/diffusion/Dsolve.cpp
index 61ffb1bc..d2f27122 100644
--- a/diffusion/Dsolve.cpp
+++ b/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/diffusion/Dsolve.h b/diffusion/Dsolve.h
index d6028e62..b8197ef4 100644
--- a/diffusion/Dsolve.h
+++ b/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/kinetics/Enz.cpp b/kinetics/Enz.cpp
index 2362c0b3..f41d109a 100644
--- a/kinetics/Enz.cpp
+++ b/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/kinetics/MMenz.cpp b/kinetics/MMenz.cpp
index 8b396ceb..cd03b00f 100644
--- a/kinetics/MMenz.cpp
+++ b/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/kinetics/ReadKkit.cpp b/kinetics/ReadKkit.cpp
index d3b5bd41..8172ae90 100644
--- a/kinetics/ReadKkit.cpp
+++ b/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/kinetics/WriteKkit.cpp b/kinetics/WriteKkit.cpp
index 0f8533c0..c32b1fa0 100644
--- a/kinetics/WriteKkit.cpp
+++ b/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/ksolve/ZombieEnz.cpp b/ksolve/ZombieEnz.cpp
index 361ac111..e73c27ac 100644
--- a/ksolve/ZombieEnz.cpp
+++ b/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/ksolve/ZombieMMenz.cpp b/ksolve/ZombieMMenz.cpp
index 95f163a0..dc400897 100644
--- a/ksolve/ZombieMMenz.cpp
+++ b/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/python/rdesigneur/fixXreacs.py b/python/fixXreacs.py
similarity index 99%
rename from python/rdesigneur/fixXreacs.py
rename to python/fixXreacs.py
index 19379f5f..de7bec63 100644
--- a/python/rdesigneur/fixXreacs.py
+++ b/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/python/moose/SBML/readSBML.py b/python/moose/SBML/readSBML.py
index 06de82bc..d9b0e90f 100644
--- a/python/moose/SBML/readSBML.py
+++ b/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/python/moose/SBML/writeSBML.py b/python/moose/SBML/writeSBML.py
index 8491b46a..8124611b 100644
--- a/python/moose/SBML/writeSBML.py
+++ b/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/python/moose/chemMerge/merge.py b/python/moose/chemMerge/merge.py
index 49d7d716..d3d80058 100644
--- a/python/moose/chemMerge/merge.py
+++ b/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/python/moose/chemUtil/add_Delete_ChemicalSolver.py b/python/moose/chemUtil/add_Delete_ChemicalSolver.py
index 05bd4a0a..cc350ec6 100644
--- a/python/moose/chemUtil/add_Delete_ChemicalSolver.py
+++ b/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/python/moose/utils.py b/python/moose/utils.py
index 37bf0deb..9abf7ae5 100644
--- a/python/moose/utils.py
+++ b/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/python/rdesigneur/moogul.py b/python/rdesigneur/moogul.py
new file mode 100644
index 00000000..9f2ecba5
--- /dev/null
+++ b/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/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py
index 0ef8d5c4..66a2ea68 100644
--- a/python/rdesigneur/rdesigneur.py
+++ b/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/python/rdesigneur/rdesigneurProtos.py b/python/rdesigneur/rdesigneurProtos.py
index 1c22ea0c..ec5d7059 100644
--- a/python/rdesigneur/rdesigneurProtos.py
+++ b/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/python/rdesigneur/rmoogli.py b/python/rdesigneur/rmoogli.py
index fcb1a120..f98549c0 100644
--- a/python/rdesigneur/rmoogli.py
+++ b/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/tests/python/test_kkit.py b/tests/python/test_kkit.py
index 8ba4fd2b..ba194d63 100644
--- a/tests/python/test_kkit.py
+++ b/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/tests/python/test_negative_value_flag.py b/tests/python/test_negative_value_flag.py
index 35c37ad7..7b073975 100644
--- a/tests/python/test_negative_value_flag.py
+++ b/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 )
--
GitLab