diff --git a/moose-core/CMakeLists.txt b/moose-core/CMakeLists.txt
index 9bf8255c1dc2d3f3cafe068ef33daa5b84392b5f..3be4235b3d563ff57adfe0d952ebecc84c44d993 100644
--- a/moose-core/CMakeLists.txt
+++ b/moose-core/CMakeLists.txt
@@ -130,6 +130,7 @@ include_directories( ${CMAKE_CURRENT_SOURCE_DIR} )
if(WITH_BOOST)
set(WITH_BOOST_ODE ON)
+ set(WITH_GSL OFF)
endif(WITH_BOOST)
# If using BOOST ODE2 library to solve ODE system, then don't use GSL.
diff --git a/moose-core/biophysics/Neuron.cpp b/moose-core/biophysics/Neuron.cpp
index aa86759175fc376489128abe3f4a14d1caae32b5..3fbab6d0a237cd878e4874e1accd3046cf403c61 100644
--- a/moose-core/biophysics/Neuron.cpp
+++ b/moose-core/biophysics/Neuron.cpp
@@ -315,11 +315,17 @@ const Cinfo* Neuron::initCinfo()
"The expression for the *spacing* field must evaluate to > 0 for "
"the spine to be installed. For example, if the expresssion is\n"
" H(1 - L) \n"
- "then the systemwill only put spines closer than "
+ "then the system will only put spines closer than "
"one length constant from the soma, and zero elsewhere. \n"
"Available spine parameters are: \n"
"spacing, minSpacing, size, sizeDistrib "
- "angle, angleDistrib \n",
+ "angle, angleDistrib \n"
+ "minSpacing sets the granularity of sampling (typically about 0.1*"
+ "spacing) for the usual case where spines are spaced randomly. "
+ "If minSpacing < 0 then the spines are spaced equally at "
+ "'spacing', unless the dendritic segment length is smaller than "
+ "'spacing'. In that case it falls back to the regular random "
+ "placement method.",
&Neuron::setSpineDistribution,
&Neuron::getSpineDistribution
);
@@ -1788,6 +1794,19 @@ static void addPos( unsigned int segIndex, unsigned int eIndex,
vector< unsigned int >& elistIndex,
vector< double >& pos )
{
+ if ( minSpacing < 0.0 ) {
+ // Use uniform spacing
+ for ( double position = spacing * 0.5;
+ position < dendLength; position += spacing ) {
+ seglistIndex.push_back( segIndex );
+ elistIndex.push_back( eIndex );
+ pos.push_back( position );
+ }
+ if ( dendLength > spacing * 0.5 )
+ return;
+ // If the dend length is too small for regular placement,
+ // fall back to using probability to decide if segment gets spine
+ }
if ( minSpacing < spacing * 0.1 && minSpacing < 1e-7 )
minSpacing = spacing * 0.1;
if ( minSpacing > spacing * 0.5 )
@@ -1860,13 +1879,6 @@ void Neuron::makeSpacingDistrib( const vector< ObjId >& elist,
{
double spacing = val[ j + nuParser::EXPR ];
double spacingDistrib = parser.eval( val.begin() + j );
- if ( spacingDistrib > spacing || spacingDistrib < 0 )
- {
- cout << "Warning: Neuron::makeSpacingDistrib: " <<
- "0 < " << spacingDistrib << " < " << spacing <<
- " fails on " << elist[i].path() << ". Using 0.\n";
- spacingDistrib = 0.0;
- }
map< Id, unsigned int>::const_iterator
lookupDend = segIndex_.find( elist[i] );
if ( lookupDend != segIndex_.end() )
diff --git a/moose-core/ksolve/Ksolve.cpp b/moose-core/ksolve/Ksolve.cpp
index 2277e5ffb37d5041eb56487cfa09766568b9e980..81bc118fc29118f5a746740710879a286468ffea 100644
--- a/moose-core/ksolve/Ksolve.cpp
+++ b/moose-core/ksolve/Ksolve.cpp
@@ -235,7 +235,7 @@ static const Cinfo* ksolveCinfo = Ksolve::initCinfo();
Ksolve::Ksolve()
:
-#if USE_GSL
+#ifdef USE_GSL
method_( "rk5" ),
#elif USE_BOOST_ODE
method_( "rk5a" ),
@@ -248,19 +248,12 @@ Ksolve::Ksolve()
dsolve_(),
dsolvePtr_( 0 )
{
+ ;
}
Ksolve::~Ksolve()
{
-#if 0
- char* p = getenv( "MOOSE_SHOW_SOLVER_PERF" );
- if( p != NULL )
- {
- cout << "Info: Ksolve (+Dsolve) took " << totalTime_ << " seconds and took " << numSteps_
- << " steps." << endl;
-
- }
-#endif
+ ;
}
//////////////////////////////////////////////////////////////
@@ -378,6 +371,7 @@ void Ksolve::setStoich( Id stoich )
assert( stoich.element()->cinfo()->isA( "Stoich" ) );
stoich_ = stoich;
stoichPtr_ = reinterpret_cast< Stoich* >( stoich.eref().data() );
+
if ( !isBuilt_ )
{
OdeSystem ode;
@@ -386,6 +380,7 @@ void Ksolve::setStoich( Id stoich )
// ode.initStepSize = getEstimatedDt();
ode.initStepSize = 0.01; // This will be overridden at reinit.
ode.method = method_;
+
#ifdef USE_GSL
ode.gslSys.dimension = stoichPtr_->getNumAllPools();
if ( ode.gslSys.dimension == 0 )
@@ -416,6 +411,7 @@ void Ksolve::setStoich( Id stoich )
#endif
isBuilt_ = true;
}
+
}
Id Ksolve::getDsolve() const
diff --git a/moose-core/ksolve/OdeSystem.h b/moose-core/ksolve/OdeSystem.h
index 6a3ede6dcf26d589ada9b4ffd0dfa15a16cfb557..511513b57cccfe8cdb7e2911e58e769bdb9b4eab 100644
--- a/moose-core/ksolve/OdeSystem.h
+++ b/moose-core/ksolve/OdeSystem.h
@@ -26,15 +26,16 @@ class OdeSystem {
{;}
std::string method;
- // GSL stuff
+
+ double initStepSize;
+ double epsAbs; // Absolute error
+ double epsRel; // Relative error
#ifdef USE_GSL
+ // GSL stuff
gsl_odeiv2_system gslSys;
const gsl_odeiv2_step_type* gslStep;
#endif
- double initStepSize;
- double epsAbs; // Absolute error
- double epsRel; // Relative error
#if USE_BOOST_ODE
size_t dimension;
diff --git a/moose-core/ksolve/SteadyStateBoost.cpp b/moose-core/ksolve/SteadyStateBoost.cpp
index a96439baec7d42630d9fe2441641a2c5c5521bd0..5d31165d66d52edc2d493aab558b429f151aff1b 100644
--- a/moose-core/ksolve/SteadyStateBoost.cpp
+++ b/moose-core/ksolve/SteadyStateBoost.cpp
@@ -85,27 +85,6 @@ struct reac_info
const Cinfo* SteadyState::initCinfo()
{
- /**
- * This picks up the entire Stoich data structure
- static Finfo* gslShared[] =
- {
- new SrcFinfo( "reinitSrc", Ftype0() ),
- new DestFinfo( "assignStoich",
- Ftype1< void* >(),
- RFCAST( &SteadyState::assignStoichFunc )
- ),
- new DestFinfo( "setMolN",
- Ftype2< double, unsigned int >(),
- RFCAST( &SteadyState::setMolN )
- ),
- new SrcFinfo( "requestYsrc", Ftype0() ),
- new DestFinfo( "assignY",
- Ftype1< double* >(),
- RFCAST( &SteadyState::assignY )
- ),
- };
- */
-
/**
* These are the fields of the SteadyState class
*/
@@ -207,59 +186,59 @@ const Cinfo* SteadyState::initCinfo()
// MsgDest definitions
///////////////////////////////////////////////////////
static DestFinfo setupMatrix( "setupMatrix",
- "This function initializes and rebuilds the matrices used "
- "in the calculation.",
- new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
- );
+ "This function initializes and rebuilds the matrices used "
+ "in the calculation.",
+ new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
+ );
static DestFinfo settle( "settle",
- "Finds the nearest steady state to the current initial "
- "conditions. This function rebuilds the entire calculation "
- "only if the object has not yet been initialized.",
- new OpFunc0< SteadyState >( &SteadyState::settleFunc )
- );
+ "Finds the nearest steady state to the current initial "
+ "conditions. This function rebuilds the entire calculation "
+ "only if the object has not yet been initialized.",
+ new OpFunc0< SteadyState >( &SteadyState::settleFunc )
+ );
static DestFinfo resettle( "resettle",
- "Finds the nearest steady state to the current initial "
- "conditions. This function rebuilds the entire calculation ",
- new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
- );
+ "Finds the nearest steady state to the current initial "
+ "conditions. This function rebuilds the entire calculation ",
+ new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
+ );
static DestFinfo showMatrices( "showMatrices",
- "Utility function to show the matrices derived for the calculations on the reaction system. Shows the Nr, gamma, and total matrices",
- new OpFunc0< SteadyState >( &SteadyState::showMatrices )
- );
+ "Utility function to show the matrices derived for the calculations on the reaction system. Shows the Nr, gamma, and total matrices",
+ new OpFunc0< SteadyState >( &SteadyState::showMatrices )
+ );
static DestFinfo randomInit( "randomInit",
- "Generate random initial conditions consistent with the mass"
- "conservation rules. Typically invoked in order to scan"
- "states",
- new EpFunc0< SteadyState >(
- &SteadyState::randomizeInitialCondition )
- );
+ "Generate random initial conditions consistent with the mass"
+ "conservation rules. Typically invoked in order to scan"
+ "states",
+ new EpFunc0< SteadyState >(
+ &SteadyState::randomizeInitialCondition )
+ );
///////////////////////////////////////////////////////
// Shared definitions
///////////////////////////////////////////////////////
static Finfo * steadyStateFinfos[] =
{
- &stoich, // Value
- &badStoichiometry, // ReadOnlyValue
- &isInitialized, // ReadOnlyValue
- &nIter, // ReadOnlyValue
- &status, // ReadOnlyValue
- &maxIter, // Value
- &convergenceCriterion, // ReadOnlyValue
- &numVarPools, // ReadOnlyValue
- &rank, // ReadOnlyValue
- &stateType, // ReadOnlyValue
- &nNegEigenvalues, // ReadOnlyValue
- &nPosEigenvalues, // ReadOnlyValue
- &solutionStatus, // ReadOnlyValue
- &total, // LookupValue
- &eigenvalues, // ReadOnlyLookupValue
- &setupMatrix, // DestFinfo
- &settle, // DestFinfo
- &resettle, // DestFinfo
- &showMatrices, // DestFinfo
- &randomInit, // DestFinfo
+ &stoich, // Value
+ &badStoichiometry, // ReadOnlyValue
+ &isInitialized, // ReadOnlyValue
+ &nIter, // ReadOnlyValue
+ &status, // ReadOnlyValue
+ &maxIter, // Value
+ &convergenceCriterion, // ReadOnlyValue
+ &numVarPools, // ReadOnlyValue
+ &rank, // ReadOnlyValue
+ &stateType, // ReadOnlyValue
+ &nNegEigenvalues, // ReadOnlyValue
+ &nPosEigenvalues, // ReadOnlyValue
+ &solutionStatus, // ReadOnlyValue
+ &total, // LookupValue
+ &eigenvalues, // ReadOnlyLookupValue
+ &setupMatrix, // DestFinfo
+ &settle, // DestFinfo
+ &resettle, // DestFinfo
+ &showMatrices, // DestFinfo
+ &randomInit, // DestFinfo
};
@@ -278,7 +257,7 @@ const Cinfo* SteadyState::initCinfo()
"Note that the method finds unstable as well as stable fixed "
"points.\n "
"The SteadyState class also provides a utility function "
- "*randomInit()* to "
+ "*randomInit()* to "
"randomly initialize the concentrations, within the constraints "
"of stoichiometry. This is useful if you are trying to find "
"the major fixed points of the system. Note that this is "
@@ -368,12 +347,10 @@ void SteadyState::setStoich( Id value )
numVarPools_ = Field< unsigned int >::get( stoich_, "numVarPools" );
nReacs_ = Field< unsigned int >::get( stoich_, "numRates" );
setupSSmatrix();
- double vol = LookupField< unsigned int, double >::get(
- stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
+ double vol = LookupField< unsigned int, double >::get(stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
pool_.setVolume( vol );
pool_.setStoich( stoichPtr, 0 );
- pool_.updateAllRateTerms( stoichPtr->getRateTerms(),
- stoichPtr->getNumCoreRates() );
+ pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
isInitialized_ = 1;
}
@@ -513,9 +490,9 @@ void SteadyState::showMatrices()
return;
}
int numConsv = numVarPools_ - rank_;
- cout << "Totals: ";
+ cout << "Totals: ";
for ( int i = 0; i < numConsv; ++i )
- cout << total_[i] << " ";
+ cout << total_[i] << " ";
cout << endl;
cout << "gamma " << gamma_ << endl;
cout << "Nr " << Nr_ << endl;
@@ -549,9 +526,7 @@ void SteadyState::setupSSmatrix()
{
double x = 0;
if ( j == colIndex[k] && k < rowStart[i+1] )
- {
x = entry[k++];
- }
N(i,j) = x;
LU_(i,j) = x;
}
@@ -714,7 +689,7 @@ static bool isSolutionValid( const vector< double >& x )
for( size_t i = 0; i < x.size(); i++ )
{
double v = x[i];
- if ( std::isnan( v ) or std::isinf( v ) )
+ if ( std::isnan( v ) || std::isinf( v ) )
{
cout << "Warning: SteadyState iteration gave nan/inf concs\n";
return false;
diff --git a/moose-core/ksolve/SteadyStateGsl.cpp b/moose-core/ksolve/SteadyStateGsl.cpp
index 295b5256d5edcd4e574ef7b9e99222db4aef9e23..5aa0941dacf666a66978c5b1ac66872dfe8731f9 100644
--- a/moose-core/ksolve/SteadyStateGsl.cpp
+++ b/moose-core/ksolve/SteadyStateGsl.cpp
@@ -87,20 +87,20 @@ const Cinfo* SteadyState::initCinfo()
* This picks up the entire Stoich data structure
static Finfo* gslShared[] =
{
- new SrcFinfo( "reinitSrc", Ftype0::global() ),
- new DestFinfo( "assignStoich",
- Ftype1< void* >::global(),
- RFCAST( &SteadyState::assignStoichFunc )
- ),
- new DestFinfo( "setMolN",
- Ftype2< double, unsigned int >::global(),
- RFCAST( &SteadyState::setMolN )
- ),
- new SrcFinfo( "requestYsrc", Ftype0::global() ),
- new DestFinfo( "assignY",
- Ftype1< double* >::global(),
- RFCAST( &SteadyState::assignY )
- ),
+ new SrcFinfo( "reinitSrc", Ftype0::global() ),
+ new DestFinfo( "assignStoich",
+ Ftype1< void* >::global(),
+ RFCAST( &SteadyState::assignStoichFunc )
+ ),
+ new DestFinfo( "setMolN",
+ Ftype2< double, unsigned int >::global(),
+ RFCAST( &SteadyState::setMolN )
+ ),
+ new SrcFinfo( "requestYsrc", Ftype0::global() ),
+ new DestFinfo( "assignY",
+ Ftype1< double* >::global(),
+ RFCAST( &SteadyState::assignY )
+ ),
};
*/
@@ -205,61 +205,61 @@ const Cinfo* SteadyState::initCinfo()
// MsgDest definitions
///////////////////////////////////////////////////////
static DestFinfo setupMatrix( "setupMatrix",
- "This function initializes and rebuilds the matrices used "
- "in the calculation.",
- new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
- );
+ "This function initializes and rebuilds the matrices used "
+ "in the calculation.",
+ new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
+ );
static DestFinfo settle( "settle",
- "Finds the nearest steady state to the current initial "
- "conditions. This function rebuilds the entire calculation "
- "only if the object has not yet been initialized.",
- new OpFunc0< SteadyState >( &SteadyState::settleFunc )
- );
+ "Finds the nearest steady state to the current initial "
+ "conditions. This function rebuilds the entire calculation "
+ "only if the object has not yet been initialized.",
+ new OpFunc0< SteadyState >( &SteadyState::settleFunc )
+ );
static DestFinfo resettle( "resettle",
- "Finds the nearest steady state to the current initial "
- "conditions. This function rebuilds the entire calculation ",
- new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
- );
+ "Finds the nearest steady state to the current initial "
+ "conditions. This function rebuilds the entire calculation ",
+ new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
+ );
static DestFinfo showMatrices( "showMatrices",
- "Utility function to show the matrices derived for the calculations on the reaction system. Shows the Nr, gamma, and total matrices",
- new OpFunc0< SteadyState >( &SteadyState::showMatrices )
- );
+ "Utility function to show the matrices derived for the "
+ "calculations on the reaction system. Shows the Nr, gamma, and total matrices",
+ new OpFunc0< SteadyState >( &SteadyState::showMatrices )
+ );
static DestFinfo randomInit( "randomInit",
- "Generate random initial conditions consistent with the mass"
- "conservation rules. Typically invoked in order to scan"
- "states",
- new EpFunc0< SteadyState >(
- &SteadyState::randomizeInitialCondition )
- );
+ "Generate random initial conditions consistent with the mass"
+ "conservation rules. Typically invoked in order to scan"
+ "states",
+ new EpFunc0< SteadyState >(
+ &SteadyState::randomizeInitialCondition )
+ );
+
///////////////////////////////////////////////////////
// Shared definitions
///////////////////////////////////////////////////////
static Finfo * steadyStateFinfos[] =
{
- &stoich, // Value
- &badStoichiometry, // ReadOnlyValue
- &isInitialized, // ReadOnlyValue
- &nIter, // ReadOnlyValue
- &status, // ReadOnlyValue
- &maxIter, // Value
- &convergenceCriterion, // ReadOnlyValue
- &numVarPools, // ReadOnlyValue
- &rank, // ReadOnlyValue
- &stateType, // ReadOnlyValue
- &nNegEigenvalues, // ReadOnlyValue
- &nPosEigenvalues, // ReadOnlyValue
- &solutionStatus, // ReadOnlyValue
- &total, // LookupValue
- &eigenvalues, // ReadOnlyLookupValue
- &setupMatrix, // DestFinfo
- &settle, // DestFinfo
- &resettle, // DestFinfo
- &showMatrices, // DestFinfo
- &randomInit, // DestFinfo
-
-
+ &stoich, // Value
+ &badStoichiometry, // ReadOnlyValue
+ &isInitialized, // ReadOnlyValue
+ &nIter, // ReadOnlyValue
+ &status, // ReadOnlyValue
+ &maxIter, // Value
+ &convergenceCriterion, // ReadOnlyValue
+ &numVarPools, // ReadOnlyValue
+ &rank, // ReadOnlyValue
+ &stateType, // ReadOnlyValue
+ &nNegEigenvalues, // ReadOnlyValue
+ &nPosEigenvalues, // ReadOnlyValue
+ &solutionStatus, // ReadOnlyValue
+ &total, // LookupValue
+ &eigenvalues, // ReadOnlyLookupValue
+ &setupMatrix, // DestFinfo
+ &settle, // DestFinfo
+ &resettle, // DestFinfo
+ &showMatrices, // DestFinfo
+ &randomInit, // DestFinfo
};
static string doc[] =
@@ -276,7 +276,7 @@ const Cinfo* SteadyState::initCinfo()
"Note that the method finds unstable as well as stable fixed "
"points.\n "
"The SteadyState class also provides a utility function "
- "*randomInit()* to "
+ "*randomInit()* to "
"randomly initialize the concentrations, within the constraints "
"of stoichiometry. This is useful if you are trying to find "
"the major fixed points of the system. Note that this is "
@@ -316,7 +316,6 @@ static const Cinfo* steadyStateCinfo = SteadyState::initCinfo();
///////////////////////////////////////////////////
// Class function definitions
///////////////////////////////////////////////////
-
SteadyState::SteadyState()
:
nIter_( 0 ),
@@ -382,9 +381,10 @@ void SteadyState::setStoich( Id value )
double vol = LookupField< unsigned int, double >::get(
stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
pool_.setVolume( vol );
- pool_.setStoich( stoichPtr, 0 );
- pool_.updateAllRateTerms( stoichPtr->getRateTerms(),
- stoichPtr->getNumCoreRates() );
+
+ pool_.setStoich( stoichPtr, nullptr );
+
+ pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
isInitialized_ = 1;
}
@@ -548,9 +548,9 @@ void SteadyState::showMatrices()
return;
}
int numConsv = numVarPools_ - rank_;
- cout << "Totals: ";
+ cout << "Totals: ";
for ( int i = 0; i < numConsv; ++i )
- cout << total_[i] << " ";
+ cout << total_[i] << " ";
cout << endl;
#ifdef USE_GSL
print_gsl_mat( gamma_, "gamma" );
@@ -584,7 +584,7 @@ void SteadyState::setupSSmatrix()
{
gsl_matrix_set (LU_, i, i + nReacs_, 1 );
unsigned int k = rowStart[i];
- // cout << endl << i << ": ";
+ // cout << endl << i << ": ";
for ( unsigned int j = 0; j < nReacs_; ++j )
{
double x = 0;
@@ -592,7 +592,7 @@ void SteadyState::setupSSmatrix()
{
x = entry[k++];
}
- // cout << " " << x;
+ // cout << " " << x;
gsl_matrix_set (N, i, j, x);
gsl_matrix_set (LU_, i, j, x );
}
@@ -637,10 +637,10 @@ void SteadyState::setupSSmatrix()
/*
cout << "S = (";
for ( unsigned int j = 0; j < numVarPools_; ++j )
- cout << s_->S()[ j ] << ", ";
+ cout << s_->S()[ j ] << ", ";
cout << "), Sinit = ( ";
for ( unsigned int j = 0; j < numVarPools_; ++j )
- cout << s_->Sinit()[ j ] << ", ";
+ cout << s_->Sinit()[ j ] << ", ";
cout << ")\n";
*/
Id ksolve = Field< Id >::get( stoich_, "ksolve" );
diff --git a/moose-core/ksolve/VoxelPools.cpp b/moose-core/ksolve/VoxelPools.cpp
index ea575304362e063afda23ec222e09c227a93019f..25e8ac7e875aed68b074c59d4b8a64dcb88a0593 100644
--- a/moose-core/ksolve/VoxelPools.cpp
+++ b/moose-core/ksolve/VoxelPools.cpp
@@ -6,7 +6,8 @@
** GNU Lesser General Public License version 2.1
** See the file COPYING.LIB for the full notice.
**********************************************************************/
-#include "header.h"
+
+#include "../basecode/header.h"
#ifdef USE_GSL
#include <gsl/gsl_errno.h>
@@ -66,9 +67,12 @@ void VoxelPools::reinit( double dt )
void VoxelPools::setStoich( Stoich* s, const OdeSystem* ode )
{
stoichPtr_ = s;
- absTol_ = ode->epsAbs;
- relTol_ = ode->epsRel;
- method_ = ode->method;
+ if( ode )
+ {
+ epsAbs_ = ode->epsAbs;
+ epsRel_ = ode->epsRel;
+ method_ = ode->method;
+ }
#ifdef USE_GSL
if ( ode )
@@ -77,10 +81,9 @@ void VoxelPools::setStoich( Stoich* s, const OdeSystem* ode )
if ( driver_ )
gsl_odeiv2_driver_free( driver_ );
- driver_ = gsl_odeiv2_driver_alloc_y_new(
- &sys_, ode->gslStep, ode->initStepSize,
- ode->epsAbs, ode->epsRel
- );
+ driver_ = gsl_odeiv2_driver_alloc_y_new( &sys_, ode->gslStep
+ , ode->initStepSize, ode->epsAbs, ode->epsRel
+ );
}
#endif
VoxelPoolsBase::reinit();
@@ -89,6 +92,7 @@ void VoxelPools::setStoich( Stoich* s, const OdeSystem* ode )
void VoxelPools::advance( const ProcInfo* p )
{
double t = p->currTime - p->dt;
+
#ifdef USE_GSL
int status = gsl_odeiv2_driver_apply( driver_, &t, p->currTime, varS());
if ( status != GSL_SUCCESS )
@@ -135,114 +139,41 @@ void VoxelPools::advance( const ProcInfo* p )
* user should provide the stepping size when using fixed dt. This feature
* can be incredibly useful on large system.
*/
- const double fixedDt = 0.1;
- if( method_ == "rk2" )
+ // Variout stepper times are listed here:
+ // https://www.boost.org/doc/libs/1_68_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.explicit_steppers
+
+ auto sys = [this](const vector_type_& dy, vector_type_& dydt, const double t) {
+ VoxelPools::evalRates(this, dy, dydt); };
+
+ // This is usually the default method for boost: Runge Kutta Fehlberg
+ if( method_ == "rk5")
+ odeint::integrate_const( rk_karp_stepper_type_()
+ , sys , Svec() , p->currTime - p->dt, p->currTime, p->dt);
+ else if( method_ == "rk2" )
odeint::integrate_const( rk_midpoint_stepper_type_()
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
- );
+ , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
else if( method_ == "rk4" )
odeint::integrate_const( rk4_stepper_type_()
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
- );
- else if( method_ == "rk5")
- odeint::integrate_const( rk_karp_stepper_type_()
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
- );
+ , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt );
else if( method_ == "rk5a")
- odeint::integrate_adaptive(
- odeint::make_controlled<rk_karp_stepper_type_>( absTol_, relTol_ )
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt
- , p->currTime
- , p->dt
- );
+ odeint::integrate_adaptive( odeint::make_controlled<rk_dopri_stepper_type_>( epsAbs_, epsRel_ )
+ , sys , Svec() , p->currTime - p->dt , p->currTime, p->dt );
else if ("rk54" == method_ )
odeint::integrate_const( rk_karp_stepper_type_()
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
- );
+ , sys , Svec() , p->currTime - p->dt, p->currTime, p->dt);
else if ("rk54a" == method_ )
- odeint::integrate_adaptive(
- odeint::make_controlled<rk_karp_stepper_type_>( absTol_, relTol_ )
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt
- , p->currTime
- , p->dt
- );
- else if ("rk5" == method_ )
- odeint::integrate_const( rk_dopri_stepper_type_()
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt
- , p->currTime
- , std::min( p->dt, fixedDt )
- );
- else if ("rk5a" == method_ )
- odeint::integrate_adaptive(
- odeint::make_controlled<rk_dopri_stepper_type_>( absTol_, relTol_ )
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt
- , p->currTime
- , p->dt
- );
+ odeint::integrate_adaptive( odeint::make_controlled<rk_karp_stepper_type_>( epsAbs_, epsRel_ )
+ , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
else if( method_ == "rk8" )
odeint::integrate_const( rk_felhberg_stepper_type_()
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt )
- );
+ , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
else if( method_ == "rk8a" )
- odeint::integrate_adaptive(
- odeint::make_controlled<rk_felhberg_stepper_type_>( absTol_, relTol_ )
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt
- , p->currTime
- , p->dt
- );
-
+ odeint::integrate_adaptive( odeint::make_controlled<rk_felhberg_stepper_type_>( epsAbs_, epsRel_ )
+ , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
else
- odeint::integrate_adaptive(
- odeint::make_controlled<rk_karp_stepper_type_>( absTol_, relTol_ )
- , [this](const vector_type_& dy, vector_type_& dydt, const double t) {
- VoxelPools::evalRates(this, dy, dydt );
- }
- , Svec()
- , p->currTime - p->dt
- , p->currTime
- , p->dt
- );
+ odeint::integrate_adaptive( odeint::make_controlled<rk_karp_stepper_type_>( epsAbs_, epsRel_ )
+ , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
#endif
if ( !stoichPtr_->getAllowNegative() ) // clean out negatives
@@ -353,8 +284,7 @@ void VoxelPools::updateRates( const double* s, double* yprime ) const
stoichPtr_->getNumProxyPools();
// totVar should include proxyPools if this voxel does not use them
unsigned int totInvar = stoichPtr_->getNumBufPools();
- assert( N.nColumns() == 0 ||
- N.nRows() == stoichPtr_->getNumAllPools() );
+ assert( N.nColumns() == 0 || N.nRows() == stoichPtr_->getNumAllPools() );
assert( N.nColumns() == rates_.size() );
for ( vector< RateTerm* >::const_iterator
diff --git a/moose-core/ksolve/VoxelPools.h b/moose-core/ksolve/VoxelPools.h
index a39b151c4ed8185c29829c980a38bc5782169da9..04c5bb8c7ed03f96dc1743e8ac5f7bd299448f89 100644
--- a/moose-core/ksolve/VoxelPools.h
+++ b/moose-core/ksolve/VoxelPools.h
@@ -99,8 +99,8 @@ private:
gsl_odeiv2_system sys_;
#endif
- double absTol_;
- double relTol_;
+ double epsAbs_;
+ double epsRel_;
string method_;
};
diff --git a/moose-core/python/rdesigneur/rdesigneur.py b/moose-core/python/rdesigneur/rdesigneur.py
index 568bac944455073b43e71ad2cf2dbaa3faf69e9c..447faa110d3b6f9469e9f1719f33bd785ab8d33a 100644
--- a/moose-core/python/rdesigneur/rdesigneur.py
+++ b/moose-core/python/rdesigneur/rdesigneur.py
@@ -73,6 +73,7 @@ class rdesigneur:
combineSegments = True,
stealCellFromLibrary = False,
verbose = True,
+ benchmark = False,
addSomaChemCompt = False, # Put a soma chemCompt on neuroMesh
addEndoChemCompt = False, # Put an endo compartment, typically for ER, on each of the NeuroMesh compartments.
diffusionLength= 2e-6,
@@ -96,8 +97,8 @@ class rdesigneur:
chemDistrib = [],
adaptorList= [],
stimList = [],
- plotList = [],
- moogList = [],
+ plotList = [], # elecpath, geom_expr, object, field, title ['wave' [min max]]
+ moogList = [],
params = None
):
""" Constructor of the rdesigner. This just sets up internal fields
@@ -110,6 +111,7 @@ class rdesigneur:
self.combineSegments = combineSegments
self.stealCellFromLibrary = stealCellFromLibrary
self.verbose = verbose
+ self.benchmark = benchmark
self.addSomaChemCompt = addSomaChemCompt
self.addEndoChemCompt = addEndoChemCompt
self.diffusionLength= diffusionLength
@@ -138,11 +140,16 @@ class rdesigneur:
self.params = params
self.adaptorList = adaptorList
- self.stimList = stimList
- self.plotList = plotList
- self.saveList = plotList #ADDED BY Sarthak
+ try:
+ self.stimList = [ rstim.convertArg(i) for i in stimList ]
+ self.plotList = [ rplot.convertArg(i) for i in plotList ]
+ self.moogList = [ rmoog.convertArg(i) for i in moogList ]
+ except BuildError as msg:
+ print("Error: rdesigneur: " + msg)
+ quit()
+
+ #self.saveList = plotList #ADDED BY Sarthak
self.saveAs = []
- self.moogList = moogList
self.plotNames = []
self.wavePlotNames = []
self.saveNames = []
@@ -166,7 +173,9 @@ class rdesigneur:
################################################################
def _printModelStats( self ):
- print("\n\tRdesigneur: Elec model has",
+ if not self.verbose:
+ return
+ print("\nRdesigneur: Elec model has",
self.elecid.numCompartments, "compartments and",
self.elecid.numSpines, "spines on",
len( self.cellPortionElist ), "compartments.")
@@ -184,17 +193,15 @@ class rdesigneur:
return
self.model = moose.Neutral( modelPath )
self.modelPath = modelPath
- print( "[INFO ] rdesigneur: Building model. " )
funcs = [ self.installCellFromProtos, self.buildPassiveDistrib
, self.buildChanDistrib, self.buildSpineDistrib, self.buildChemDistrib
, self._configureSolvers, self.buildAdaptors, self._buildStims
, self._buildPlots, self._buildMoogli, self._configureHSolve
- , self._configureClocks, self._printModelStats, self._savePlots
+ , self._configureClocks, self._printModelStats
]
for i, _func in enumerate( funcs ):
- if self.verbose:
+ if self.benchmark:
print( " + (%d/%d) executing %25s"%(i, len(funcs), _func.__name__), end=' ' )
- sys.stdout.flush()
t0 = time.time()
try:
_func( )
@@ -203,11 +210,12 @@ class rdesigneur:
moose.delete( self.model )
return
t = time.time() - t0
- if self.verbose:
+ if self.benchmark:
msg = r' ... DONE'
if t > 1:
msg += ' %.3f sec' % t
print( msg )
+ sys.stdout.flush()
def installCellFromProtos( self ):
if self.stealCellFromLibrary:
@@ -463,10 +471,10 @@ class rdesigneur:
# 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] )
+ synWeight = float( stimInfo.geom_expr )
stimObj = []
for i in dendCompts + spineCompts:
- path = i.path + '/' + stimInfo[2] + '/sh/synapse[0]'
+ path = i.path + '/' + stimInfo.relpath + '/sh/synapse[0]'
if moose.exists( path ):
synInput = make_synInput( name='synInput', parent=path )
synInput.doPeriodic = doPeriodic
@@ -486,6 +494,10 @@ class rdesigneur:
# Expression can use p, g, L, len, dia, maxP, maxG, maxL.
temp = []
for i in self.passiveDistrib:
+ if (len( i ) < 3) or (len(i) %2 != 1):
+ raise BuildError( "buildPassiveDistrib: Need 3 + N*2 arguments, have {}".format( len(i) ) )
+
+ temp.append( '.' )
temp.extend( i )
temp.extend( [""] )
self.elecid.passiveDistribution = temp
@@ -608,18 +620,17 @@ class rdesigneur:
# [ region_wildcard, region_expr, path, field, title]
def _parseComptField( self, comptList, plotSpec, knownFields ):
# Put in stuff to go through fields if the target is a chem object
- field = plotSpec[3]
+ field = plotSpec.field
if not field in knownFields:
print("Warning: Rdesigneur::_parseComptField: Unknown field '{}'".format( field ) )
return (), ""
kf = knownFields[field] # Find the field to decide type.
if kf[0] in ['CaConcBase', 'ChanBase', 'NMDAChan', 'VClamp']:
- objList = self._collapseElistToPathAndClass( comptList, plotSpec[2], kf[0] )
+ objList = self._collapseElistToPathAndClass( comptList, plotSpec.relpath, kf[0] )
return objList, kf[1]
-
elif field in [ 'n', 'conc', 'volume']:
- path = plotSpec[2]
+ path = plotSpec.relpath
pos = path.find( '/' )
if pos == -1: # Assume it is in the dend compartment.
path = 'dend/' + path
@@ -640,13 +651,13 @@ class rdesigneur:
voxelVec = [i for i in range(len( em ) ) if em[i] in comptSet ]
# Here we collapse the voxelVec into objects to plot.
- allObj = moose.vec( self.modelPath + '/chem/' + plotSpec[2] )
+ allObj = moose.vec( self.modelPath + '/chem/' + plotSpec.relpath )
#print "####### allObj=", self.modelPath + '/chem/' + plotSpec[2]
if len( allObj ) >= len( voxelVec ):
objList = [ allObj[int(j)] for j in voxelVec]
else:
objList = []
- print( "Warn: Rdesigneur::_parseComptField: unknown Object: '%s'" % plotSpec[2] )
+ print( "Warning: Rdesigneur::_parseComptField: unknown Object: '", plotSpec.relpath, "'" )
#print "############", chemCompt, len(objList), kf[1]
return objList, kf[1]
@@ -678,7 +689,7 @@ class rdesigneur:
dummy = moose.element( '/' )
k = 0
for i in self.plotList:
- pair = i[0] + " " + i[1]
+ pair = i.elecpath + ' ' + i.geom_expr
dendCompts = self.elecid.compartmentsFromExpression[ pair ]
spineCompts = self.elecid.spinesFromExpression[ pair ]
plotObj, plotField = self._parseComptField( dendCompts, i, knownFields )
@@ -686,26 +697,26 @@ class rdesigneur:
assert( plotField == plotField2 )
plotObj3 = plotObj + plotObj2
numPlots = sum( q != dummy for q in plotObj3 )
- if numPlots == 0:
- return
-
- tabname = graphs.path + '/plot' + str(k)
- scale = knownFields[i[3]][2]
- units = knownFields[i[3]][3]
- ymin = i[6] if len(i) > 7 else 0
- ymax = i[7] if len(i) > 7 else 0
- if len( i ) > 5 and i[5] == 'wave':
- self.wavePlotNames.append( [ tabname, i[4], k, scale, units, i[3], ymin, ymax ] )
- else:
- self.plotNames.append( [ tabname, i[4], k, scale, units, i[3], ymin, ymax ] )
- k += 1
- if i[3] in [ 'n', 'conc', 'volume', 'Gbar' ]:
- tabs = moose.Table2( tabname, numPlots )
- else:
- tabs = moose.Table( tabname, numPlots )
- if i[3] == 'spikeTime':
- tabs.vec.threshold = -0.02 # Threshold for classifying Vm as a spike.
- tabs.vec.useSpikeMode = True # spike detect mode on
+ #print( "PlotList: {0}: numobj={1}, field ={2}, nd={3}, ns={4}".format( pair, numPlots, plotField, len( dendCompts ), len( spineCompts ) ) )
+ if numPlots > 0:
+ tabname = graphs.path + '/plot' + str(k)
+ scale = knownFields[i.field][2]
+ units = knownFields[i.field][3]
+ if i.mode == 'wave':
+ self.wavePlotNames.append( [ tabname, i.title, k, scale, units, i ] )
+ else:
+ self.plotNames.append( [ tabname, i.title, k, scale, units, i.field, i.ymin, i.ymax ] )
+ if len( i.saveFile ) > 4 and i.saveFile[-4] == '.xml' or i.saveFile:
+ self.saveNames.append( [ tabname, len(self.saveNames), scale, units, i ] )
+
+ k += 1
+ if i.field == 'n' or i.field == 'conc' or i.field == 'volume' or i.field == 'Gbar':
+ tabs = moose.Table2( tabname, numPlots )
+ else:
+ tabs = moose.Table( tabname, numPlots )
+ if i.field == 'spikeTime':
+ tabs.vec.threshold = -0.02 # Threshold for classifying Vm as a spike.
+ tabs.vec.useSpikeMode = True # spike detect mode on
vtabs = moose.vec( tabs )
q = 0
@@ -731,8 +742,8 @@ class rdesigneur:
moogliBase = moose.Neutral( self.modelPath + '/moogli' )
k = 0
for i in self.moogList:
- kf = knownFields[i[3]]
- pair = i[0] + " " + i[1]
+ kf = knownFields[i.field]
+ pair = i.elecpath + " " + i.geom_expr
dendCompts = self.elecid.compartmentsFromExpression[ pair ]
spineCompts = self.elecid.spinesFromExpression[ pair ]
dendObj, mooField = self._parseComptField( dendCompts, i, knownFields )
@@ -740,13 +751,6 @@ class rdesigneur:
assert( mooField == mooField2 )
mooObj3 = dendObj + spineObj
numMoogli = len( mooObj3 )
- #dendComptMap = self.dendCompt.elecComptMap
- #self.moogliViewer = rmoogli.makeMoogli( self, mooObj3, mooField )
- if len( i ) == 5:
- i.extend( kf[4:6] )
- elif len( i ) == 6:
- i.extend( [kf[5]] )
- #self.moogliViewer = rmoogli.makeMoogli( self, mooObj3, i, kf )
self.moogNames.append( rmoogli.makeMoogli( self, mooObj3, i, kf ) )
@@ -815,10 +819,15 @@ rdesigneur.rmoogli.updateMoogliViewer()
plt.title( i[1] )
plt.xlabel( "position (voxels)" )
plt.ylabel( i[4] )
- mn = np.min(vpts)
- mx = np.max(vpts)
- if mn/mx < 0.3:
- mn = 0
+ plotinfo = i[5]
+ if plotinfo.ymin == plotinfo.ymax:
+ mn = np.min(vpts)
+ mx = np.max(vpts)
+ if mn/mx < 0.3:
+ mn = 0
+ else:
+ mn = plotinfo.ymin
+ mx = plotinfo.ymax
ax.set_ylim( mn, mx )
line, = plt.plot( range( len( vtab ) ), vpts[0] )
timeLabel = plt.text( len(vtab ) * 0.05, mn + 0.9*(mx-mn), 'time = 0' )
@@ -847,134 +856,13 @@ rdesigneur.rmoogli.updateMoogliViewer()
Email address: sarthaks442@gmail.com
Heavily modified by U.S. Bhalla
'''
-
- def _savePlots( self ):
- if self.verbose:
- print( 'rdesigneur: Saving plots ...', end = ' ' )
- sys.stdout.flush()
-
- 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)' ),
- 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ),
- 'modulation':('ChanBase', 'getModulation', 1, 'chan modulation (unitless)' ),
- 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)' ),
- 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)' ),
- 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)' ),
- '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)' ),
- 'current':('VClamp', 'getCurrent', 1e9, 'Holding Current (nA)')
- }
-
- save_graphs = moose.Neutral( self.modelPath + '/save_graphs' )
- dummy = moose.element( '/' )
- k = 0
-
- for i in self.saveList:
- pair = i[0] + " " + i[1]
- dendCompts = self.elecid.compartmentsFromExpression[ pair ]
- spineCompts = self.elecid.spinesFromExpression[ pair ]
- plotObj, plotField = self._parseComptField( dendCompts, i, knownFields )
- plotObj2, plotField2 = self._parseComptField( spineCompts, i, knownFields )
- assert( plotField == plotField2 )
- plotObj3 = plotObj + plotObj2
- numPlots = sum( i != dummy for i in plotObj3 )
- if numPlots > 0:
- save_tabname = save_graphs.path + '/save_plot' + str(k)
- scale = knownFields[i[3]][2]
- units = knownFields[i[3]][3]
- self.saveNames.append( ( save_tabname, i[4], k, scale, units ) )
- k += 1
- if i[3] in [ 'n', 'conc', 'volume', 'Gbar' ]:
- save_tabs = moose.Table2( save_tabname, numPlots )
- save_vtabs = moose.vec( save_tabs )
- else:
- save_tabs = moose.Table( save_tabname, numPlots )
- save_vtabs = moose.vec( save_tabs )
- if i[3] == 'spikeTime':
- save_vtabs.threshold = -0.02 # Threshold for classifying Vm as a spike.
- save_vtabs.useSpikeMode = True # spike detect mode on
- q = 0
- for p in [ x for x in plotObj3 if x != dummy ]:
- moose.connect( save_vtabs[q], 'requestOut', p, plotField )
- q += 1
-
- if self.verbose:
- print( ' ... DONE.' )
-
- def _getTimeSeriesTable( self ):
-
- '''
- This function gets the list with all the details of the simulation
- required for plotting.
- This function adds flexibility in terms of the details
- we wish to store.
- '''
-
- knownFields = {
- 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ),
- 'spikeTime':('CompartmentBase', 'getVm', 1, 'Spike Times (s)'),
- 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ),
- 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ),
- 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ),
- 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)' ),
- 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)' ),
- 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)' ),
- '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)' )
- }
-
- '''
- This takes data from plotList
- saveList is exactly like plotList but with a few additional arguments:
- ->It will have a resolution option, i.e., the number of decimal figures to which the value should be rounded
- ->There is a list of "saveAs" formats
- With saveList, the user will able to set what all details he wishes to be saved.
- '''
-
- for i,ind in enumerate(self.saveNames):
- pair = self.saveList[i][0] + " " + self.saveList[i][1]
- dendCompts = self.elecid.compartmentsFromExpression[ pair ]
- spineCompts = self.elecid.spinesFromExpression[ pair ]
- # Here we get the object details from plotList
- savePlotObj, plotField = self._parseComptField( dendCompts, self.saveList[i], knownFields )
- savePlotObj2, plotField2 = self._parseComptField( spineCompts, self.saveList[i], knownFields )
- savePlotObj3 = savePlotObj + savePlotObj2
-
- rowList = list(ind)
- save_vtab = moose.vec( ind[0] )
- t = np.arange( 0, save_vtab[0].vector.size, 1 ) * save_vtab[0].dt
-
- rowList.append(save_vtab[0].dt)
- rowList.append(t)
- rowList.append([jvec.vector * ind[3] for jvec in save_vtab]) #get values
- rowList.append(self.saveList[i][3])
- rowList.append(filter(lambda obj: obj.path != '/', savePlotObj3)) #this filters out dummy elements
-
- if (type(self.saveList[i][-1])==int):
- rowList.append(self.saveList[i][-1])
- else:
- rowList.append(12)
-
- self.tabForXML.append(rowList)
- rowList = []
-
- timeSeriesTable = self.tabForXML # the list with all the details of plot
- return timeSeriesTable
-
- def _writeXML( self, filename, timeSeriesData ): #to write to XML file
-
- plotData = timeSeriesData
- print("[CAUTION] The '%s' file might be very large if all the compartments are to be saved." % filename)
+ def _writeXML( self, plotData, time, vtab ):
+ tabname = plotData[0]
+ idx = plotData[1]
+ scale = plotData[2]
+ units = plotData[3]
+ rp = plotData[4]
+ filename = rp.saveFile[:-4] + str(idx) + '.xml'
root = etree.Element("TimeSeriesPlot")
parameters = etree.SubElement( root, "parameters" )
if self.params == None:
@@ -987,36 +875,33 @@ rdesigneur.rmoogli.updateMoogliViewer()
#plotData contains all the details of a single plot
title = etree.SubElement( root, "timeSeries" )
- title.set( 'title', str(plotData[1]))
- title.set( 'field', str(plotData[8]))
- title.set( 'scale', str(plotData[3]))
- title.set( 'units', str(plotData[4]))
- title.set( 'dt', str(plotData[5]))
+ title.set( 'title', rp.title)
+ title.set( 'field', rp.field)
+ title.set( 'scale', str(scale) )
+ title.set( 'units', units)
+ title.set( 'dt', str(vtab[0].dt) )
+ res = rp.saveResolution
p = []
- assert(len(plotData[7]) == len(plotData[9]))
-
- res = plotData[10]
- for ind, jvec in enumerate(plotData[7]):
+ for t, v in zip( time, vtab ):
p.append( etree.SubElement( title, "data"))
- p[-1].set( 'path', str(plotData[9][ind].path))
- p[-1].text = ''.join( str(round(value,res)) + ' ' for value in jvec )
+ p[-1].set( 'path', v.path )
+ p[-1].text = ''.join( str(round(y,res)) + ' ' for y in v.vector )
tree = etree.ElementTree(root)
tree.write(filename)
- def _writeCSV(self, filename, timeSeriesData):
-
- plotData = timeSeriesData
- dataList = []
- header = []
- time = plotData[6]
- res = plotData[10]
-
- for ind, jvec in enumerate(plotData[7]):
- header.append(plotData[9][ind].path)
- dataList.append([round(value,res) for value in jvec.tolist()])
- dl = [tuple(lst) for lst in dataList]
- rows = zip(tuple(time), *dl)
- header.insert(0, "time")
+ def _writeCSV( self, plotData, time, vtab ):
+ tabname = plotData[0]
+ idx = plotData[1]
+ scale = plotData[2]
+ units = plotData[3]
+ rp = plotData[4]
+ filename = rp.saveFile[:-4] + str(idx) + '.csv'
+
+ header = ["time",]
+ valMatrix = [time,]
+ header.extend( [ v.path for v in vtab ] )
+ valMatrix.extend( [ v.vector for v in vtab ] )
+ nv = np.array( valMatrix ).T
with open(filename, 'wb') as f:
writer = csv.writer(f, quoting=csv.QUOTE_MINIMAL)
writer.writerow(header)
@@ -1024,42 +909,25 @@ rdesigneur.rmoogli.updateMoogliViewer()
writer.writerow(row)
##########****SAVING*****###############
- def _saveFormats(self, timeSeriesData, k, *filenames):
- "This takes in the filenames and writes to corresponding format."
- if filenames:
- for filename in filenames:
- for name in filename:
- print (name)
- if name[-4:] == '.xml':
- self._writeXML(name, timeSeriesData)
- print(name, " written")
- elif name[-4:] == '.csv':
- self._writeCSV(name, timeSeriesData)
- print(name, " written")
- else:
- print("Save format not known")
- pass
- else:
- pass
def _save( self ):
- timeSeriesTable = self._getTimeSeriesTable()
- for i,sList in enumerate(self.saveList):
-
- if (len(sList) >= 6) and (type(sList[5]) != int):
- self.saveAs.extend(filter(lambda fmt: type(fmt)!=int, sList[5:]))
- try:
- timeSeriesData = timeSeriesTable[i]
- except IndexError:
- print("The object to be plotted has all dummy elements.")
- pass
- self._saveFormats(timeSeriesData, i, self.saveAs)
- self.saveAs=[]
+ for i in self.saveNames:
+ tabname = i[0]
+ idx = i[1]
+ scale = i[2]
+ units = i[3]
+ rp = i[4] # The rplot data structure, it has the setup info.
+
+ vtab = moose.vec( tabname )
+ t = np.arange( 0, vtab[0].vector.size, 1 ) * vtab[0].dt
+ ftype = rp.filename[-4:]
+ if ftype == '.xml':
+ self._writeXML( i, t, vtab )
+ elif ftype == '.csv':
+ self._writeCSV( i, t, vtab )
else:
- pass
- else:
- pass
+ print("Save format '{}' not known, please use .csv or .xml".format( ftype ) )
################################################################
# Here we set up the stims
@@ -1080,17 +948,17 @@ rdesigneur.rmoogli.updateMoogliViewer()
k = 0
# Stimlist = [path, geomExpr, relPath, field, expr_string]
for i in self.stimList:
- pair = i[0] + " " + i[1]
+ pair = i.elecpath + " " + i.geom_expr
dendCompts = self.elecid.compartmentsFromExpression[ pair ]
spineCompts = self.elecid.spinesFromExpression[ pair ]
#print( "pair = {}, numcompts = {},{} ".format( pair, len( dendCompts), len( spineCompts ) ) )
- if i[3] == 'vclamp':
+ if i.field == 'vclamp':
stimObj3 = self._buildVclampOnCompt( dendCompts, spineCompts, i )
stimField = 'commandIn'
- elif i[3] == 'randsyn':
+ elif i.field == 'randsyn':
stimObj3 = self._buildSynInputOnCompt( dendCompts, spineCompts, i )
stimField = 'setRate'
- elif i[3] == 'periodicsyn':
+ elif i.field == 'periodicsyn':
stimObj3 = self._buildSynInputOnCompt( dendCompts, spineCompts, i, doPeriodic = True )
stimField = 'setRate'
else:
@@ -1103,7 +971,7 @@ rdesigneur.rmoogli.updateMoogliViewer()
funcname = stims.path + '/stim' + str(k)
k += 1
func = moose.Function( funcname )
- func.expr = i[4]
+ func.expr = i.expr
#if i[3] == 'vclamp': # Hack to clean up initial condition
func.doEvalAtReinit = 1
for q in stimObj3:
@@ -1619,3 +1487,95 @@ rdesigneur.rmoogli.updateMoogliViewer()
for j in range( i[1], i[2] ):
moose.connect( i[3], 'requestOut', chemVec[j], chemFieldSrc)
msg = moose.connect( i[3], 'output', elObj, elecFieldDest )
+
+
+#######################################################################
+# Some helper classes, used to define argument lists.
+#######################################################################
+
+class baseplot:
+ def __init__( self,
+ elecpath='soma', geom_expr='1', relpath='.', field='Vm' ):
+ self.elecpath = elecpath
+ self.geom_expr = geom_expr
+ self.relpath = relpath
+ self.field = field
+
+class rplot( baseplot ):
+ def __init__( self,
+ elecpath = 'soma', geom_expr = '1', relpath = '.', field = 'Vm',
+ title = 'Membrane potential',
+ mode = 'time',
+ ymin = 0.0, ymax = 0.0,
+ saveFile = "", saveResolution = 3, show = True ):
+ baseplot.__init__( self, elecpath, geom_expr, relpath, field )
+ self.title = title
+ self.mode = mode # Options: time, wave, wave_still, raster
+ self.ymin = ymin # If ymin == ymax, it autoscales.
+ self.ymax = ymax
+ if len( saveFile ) < 5:
+ self.saveFile = ""
+ else:
+ f = saveFile.split('.')
+ if len(f) < 2 or ( f[-1] != 'xml' and f[-1] != 'csv' ):
+ raise BuildError( "rplot: Filetype is '{}', must be of type .xml or .csv.".format( f[-1] ) )
+ self.saveFile = saveFile
+ self.show = show
+
+ def printme( self ):
+ print( "{}, {}, {}, {}, {}, {}, {}, {}, {}, {}".format(
+ self.elecpath,
+ self.geom_expr, self.relpath, self.field, self.title,
+ self.mode, self.ymin, self.ymax, self.saveFile, self.show ) )
+
+ @staticmethod
+ def convertArg( arg ):
+ if isinstance( arg, rplot ):
+ return arg
+ elif isinstance( arg, list ):
+ return rplot( *arg )
+ else:
+ raise BuildError( "rplot initialization failed" )
+
+class rmoog( baseplot ):
+ def __init__( self,
+ elecpath = 'soma', geom_expr = '1', relpath = '.', field = 'Vm',
+ title = 'Membrane potential',
+ ymin = 0.0, ymax = 0.0,
+ show = True ): # Could put in other display options.
+ baseplot.__init__( self, elecpath, geom_expr, relpath, field )
+ self.title = title
+ self.ymin = ymin # If ymin == ymax, it autoscales.
+ self.ymax = ymax
+ self.show = show
+
+ @staticmethod
+ def convertArg( arg ):
+ if isinstance( arg, rmoog ):
+ return arg
+ elif isinstance( arg, list ):
+ return rmoog( *arg )
+ else:
+ raise BuildError( "rmoog initialization failed" )
+
+ # Stimlist = [path, geomExpr, relPath, field, expr_string]
+class rstim( baseplot ):
+ def __init__( self,
+ elecpath = 'soma', geom_expr = '1', relpath = '.', field = 'inject', expr = '0'):
+ baseplot.__init__( self, elecpath, geom_expr, relpath, field )
+ self.expr = expr
+
+ def printme( self ):
+ print( "{}, {}, {}, {}, {}, {}, {}, {}, {}, {}".format(
+ self.elecpath,
+ self.geom_expr, self.relpath, self.field, self.expr ) )
+
+ @staticmethod
+ def convertArg( arg ):
+ if isinstance( arg, rstim ):
+ return arg
+ elif isinstance( arg, list ):
+ return rstim( *arg )
+ else:
+ raise BuildError( "rstim initialization failed" )
+
diff --git a/moose-core/python/rdesigneur/rmoogli.py b/moose-core/python/rdesigneur/rmoogli.py
index 16317b5686aca2abb90e58faad963dcd17b21061..f673dad7ca0e975286e7e0c6dab39da9cc15967b 100644
--- a/moose-core/python/rdesigneur/rmoogli.py
+++ b/moose-core/python/rdesigneur/rmoogli.py
@@ -29,7 +29,7 @@ def makeMoogli( rd, mooObj, args, fieldInfo ):
else:
ymin = fieldInfo[4]
ymax = fieldInfo[5]
- print( "fieldinfo = {}, ymin = {}, ymax = {}".format( fieldInfo, ymin, ymax ))
+ #print( "fieldinfo = {}, ymin = {}, ymax = {}".format( fieldInfo, ymin, ymax ))
viewer = moogul.MooView()
if mooField == 'n' or mooField == 'conc':
diff --git a/moose-core/tests/python/chem_models/19085.cspace b/moose-core/tests/python/chem_models/19085.cspace
new file mode 100644
index 0000000000000000000000000000000000000000..1d2cf355c02f13483bdf5ed379ba58e4b823d30e
--- /dev/null
+++ b/moose-core/tests/python/chem_models/19085.cspace
@@ -0,0 +1 @@
+M101: |DabX|Jbca| 5.59269 0.0157641 0.172865 0.361005 4.72728 1.08558 0.0982933
\ No newline at end of file
diff --git a/moose-core/tests/python/testdisabled_dose_response.py b/moose-core/tests/python/testdisabled_dose_response.py
new file mode 100644
index 0000000000000000000000000000000000000000..a840565aaccac9cf4bd64110da69dca07ff9c9a0
--- /dev/null
+++ b/moose-core/tests/python/testdisabled_dose_response.py
@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+## Makes and plots the dose response curve for bistable models
+## Author: Sahil Moza
+## June 26, 2014
+## Update:
+## Friday 14 September 2018 05:48:42 PM IST
+## Tunrned into a light-weight test by Dilawar Singh
+
+import os
+import sys
+import moose
+import numpy as np
+
+sdir_ = os.path.dirname( os.path.realpath( __file__ ) )
+
+def setupSteadyState(simdt,plotDt):
+
+ ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
+ stoich = moose.Stoich( '/model/kinetics/stoich' )
+ stoich.compartment = moose.element('/model/kinetics')
+ stoich.ksolve = ksolve
+ stoich.path = "/model/kinetics/##"
+ state = moose.SteadyState( '/model/kinetics/state' )
+ moose.reinit()
+ state.stoich = stoich
+ state.showMatrices()
+ state.convergenceCriterion = 1e-8
+ return ksolve, state
+
+def parseModelName(fileName):
+ pos1=fileName.rfind('/')
+ pos2=fileName.rfind('.')
+ directory=fileName[:pos1]
+ prefix=fileName[pos1+1:pos2]
+ suffix=fileName[pos2+1:len(fileName)]
+ return directory, prefix, suffix
+
+# Solve for the steady state
+def getState( ksolve, state, vol):
+ scale = 1.0 / ( vol * 6.022e23 )
+ moose.reinit()
+ state.randomInit() # Removing random initial condition to systematically make Dose reponse curves.
+ moose.start( 2.0 ) # Run the model for 2 seconds.
+ state.settle()
+
+ vector = []
+ a = moose.element( '/model/kinetics/a' ).conc
+ for x in ksolve.nVec[0]:
+ vector.append( x * scale)
+ failedSteadyState = any([np.isnan(x) for x in vector])
+ if not (failedSteadyState):
+ return state.stateType, state.solutionStatus, a, vector
+
+
+def main():
+ # Setup parameters for simulation and plotting
+ simdt= 1e-2
+ plotDt= 1
+
+ # Factors to change in the dose concentration in log scale
+ factorExponent = 10 ## Base: ten raised to some power.
+ factorBegin = -10
+ factorEnd = 11
+ factorStepsize = 2
+ factorScale = 10.0 ## To scale up or down the factors
+
+ # Load Model and set up the steady state solver.
+ # model = sys.argv[1] # To load model from a file.
+ model = os.path.join( sdir_, 'chem_models/19085.cspace' )
+ modelPath, modelName, modelType = parseModelName(model)
+ outputDir = modelPath
+
+ modelId = moose.loadModel(model, 'model', 'ee')
+ dosePath = '/model/kinetics/b/DabX' # The dose entity
+
+ ksolve, state = setupSteadyState( simdt, plotDt)
+ vol = moose.element( '/model/kinetics' ).volume
+ iterInit = 100
+ solutionVector = []
+ factorArr = []
+
+ enz = moose.element(dosePath)
+ init = float(enz.kcat) # Dose parameter
+
+ # Change Dose here to .
+ for factor in range(factorBegin, factorEnd, factorStepsize ):
+ scale = factorExponent ** (factor/factorScale)
+ enz.kcat = init * scale
+ print( "scale={:.3f}\tkcat={:.3f}".format( scale, enz.kcat) )
+ for num in range(iterInit):
+ stateType, solStatus, a, vector = getState( ksolve, state, vol)
+ if solStatus == 0:
+ #solutionVector.append(vector[0]/sum(vector))
+ solutionVector.append(a)
+ factorArr.append(scale)
+
+ joint = np.array([factorArr, solutionVector])
+ joint = joint[:,joint[1,:].argsort()]
+ got = np.mean( joint ), np.std( joint )
+ expected = (1.2247, 2.46)
+ # Close upto 2 decimal place is good enough.
+ assert np.isclose(got, expected, atol=1e-2).all(), "Got %s, expected %s" % (got, expected)
+ print( joint )
+
+if __name__ == '__main__':
+ main()
diff --git a/moose-examples/.travis.yml b/moose-examples/.travis.yml
index 0f209164e1c6fd9c79f41d8d6eef77a22550fbc6..46cb8a38d0d45115070f835661545650f12f7ae4 100644
--- a/moose-examples/.travis.yml
+++ b/moose-examples/.travis.yml
@@ -1,4 +1,9 @@
sudo : required
+language: python
+python:
+ - "2.7"
+ - "3.6"
+ - "3.7"
notifications:
email:
@@ -7,16 +12,8 @@ notifications:
- bhalla@ncbs.res.in
- hrani@ncbs.res.in
-install:
- - sudo ./.travis_prepare.sh
-
script:
+ - pip install pymoose --user --pre
- ./.travis_run.sh
exclude: [vendor]
-
-deploy:
- provider: pages
- skip-cleanup: true
- github-token: $GHI_TOKEN #
- keep-history: true
diff --git a/moose-examples/.travis_prepare.sh b/moose-examples/.travis_prepare.sh
deleted file mode 100755
index a2ad83cc3a0099473bde794b6a81647d87016bf2..0000000000000000000000000000000000000000
--- a/moose-examples/.travis_prepare.sh
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/usr/bin/env bash
-# ROOT should run this script.
-VERSION=$(lsb_release -r | cut -f2)
-apt-get -y update
-apt-get install cmake coreutils --force-yes
-apt-get -y --force-yes install python-qt4 python-pip graphviz
-apt-get -y --force-yes install python-h5py python-scipy python-pygraphviz
-wget -nv https://download.opensuse.org/repositories/home:moose/xUbuntu_$VERSION/Release.key -O Release.key
-apt-key add - < Release.key
-cat <<EOF > /etc/apt/sources.list.d/home:moose.list
-deb http://download.opensuse.org/repositories/home:/moose/xUbuntu_${VERSION}/ /
-EOF
-apt-get update
-apt-get install python-numpy python-matplotlib python-networkx
-apt-get install pymoose
-python -m pip install python-libsbml --user
diff --git a/moose-examples/neuroml2/NML2_SingleCompHHCell.nml b/moose-examples/neuroml2/NML2_SingleCompHHCell.nml
new file mode 100644
index 0000000000000000000000000000000000000000..1b1d43b92606d2c1e0b78e88315c5d3b11b11271
--- /dev/null
+++ b/moose-examples/neuroml2/NML2_SingleCompHHCell.nml
@@ -0,0 +1,90 @@
+<?xml version="1.0" encoding="UTF-8"?>
+
+<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"
+ xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+ xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 ../Schemas/NeuroML2/NeuroML_v2beta4.xsd"
+ id="NML2_SingleCompHHCell">
+
+ <!-- Single compartment cell with HH channels -->
+
+ <!-- This is a "pure" NeuroML 2 file. It can be included in a LEMS file for use in a simulaton
+ by the LEMS interpreter, see LEMS_NML2_Ex5_DetCell.xml -->
+
+ <ionChannelHH id="passiveChan" conductance="10pS">
+ <notes>Leak conductance</notes>
+ </ionChannelHH>
+
+
+ <ionChannelHH id="naChan" conductance="10pS" species="na">
+ <notes>Na channel</notes>
+
+ <gateHHrates id="m" instances="3">
+ <forwardRate type="HHExpLinearRate" rate="1per_ms" midpoint="-40mV" scale="10mV"/>
+ <reverseRate type="HHExpRate" rate="4per_ms" midpoint="-65mV" scale="-18mV"/>
+ </gateHHrates>
+
+ <gateHHrates id="h" instances="1">
+ <forwardRate type="HHExpRate" rate="0.07per_ms" midpoint="-65mV" scale="-20mV"/>
+ <reverseRate type="HHSigmoidRate" rate="1per_ms" midpoint="-35mV" scale="10mV"/>
+ </gateHHrates>
+
+ </ionChannelHH>
+
+
+ <ionChannelHH id="kChan" conductance="10pS" species="k">
+
+ <gateHHrates id="n" instances="4">
+ <forwardRate type="HHExpLinearRate" rate="0.1per_ms" midpoint="-55mV" scale="10mV"/>
+ <reverseRate type="HHExpRate" rate="0.125per_ms" midpoint="-65mV" scale="-80mV"/>
+ </gateHHrates>
+
+ </ionChannelHH>
+
+
+
+ <cell id="hhcell">
+
+ <morphology id="morph1">
+ <segment id="0" name="soma">
+ <proximal x="0" y="0" z="0" diameter="17.841242"/> <!--Gives a convenient surface area of 1000.0 um^2-->
+ <distal x="0" y="0" z="0" diameter="17.841242"/>
+ </segment>
+
+ <segmentGroup id="soma_group">
+ <member segment="0"/>
+ </segmentGroup>
+
+ </morphology>
+
+ <biophysicalProperties id="bioPhys1">
+
+ <membraneProperties>
+
+ <channelDensity id="leak" ionChannel="passiveChan" condDensity="3.0 S_per_m2" erev="-54.3mV" ion="non_specific"/>
+ <channelDensity id="naChans" ionChannel="naChan" condDensity="120.0 mS_per_cm2" erev="50.0 mV" ion="na"/>
+ <channelDensity id="kChans" ionChannel="kChan" condDensity="360 S_per_m2" erev="-77mV" ion="k"/>
+
+ <spikeThresh value="-20mV"/>
+ <specificCapacitance value="1.0 uF_per_cm2"/>
+ <initMembPotential value="-65mV"/>
+
+ </membraneProperties>
+
+ <intracellularProperties>
+ <resistivity value="0.03 kohm_cm"/> <!-- Note: not used in single compartment simulations -->
+ </intracellularProperties>
+
+ </biophysicalProperties>
+
+ </cell>
+
+ <pulseGenerator id="pulseGen1" delay="100ms" duration="100ms" amplitude="0.08nA"/>
+
+
+ <network id="net1">
+ <population id="hhpop" component="hhcell" size="1"/>
+ <explicitInput target="hhpop[0]" input="pulseGen1"/>
+ </network>
+
+</neuroml>
+
diff --git a/moose-examples/neuroml2/converter.py b/moose-examples/neuroml2/converter.py
new file mode 100644
index 0000000000000000000000000000000000000000..0aec54f3ca940d66eabd83fd2d986208a29a2b03
--- /dev/null
+++ b/moose-examples/neuroml2/converter.py
@@ -0,0 +1,209 @@
+# -*- coding: utf-8 -*-
+# converter.py ---
+#
+# Filename: mtoneuroml.py
+# Description:
+# Author:
+# Maintainer:
+# Created: Mon Apr 22 12:15:23 2013 (+0530)
+# Version:
+# Last-Updated: Wed Jul 10 16:36:14 2013 (+0530)
+# By: subha
+# Update #: 819
+# URL:
+# Keywords:
+# Compatibility:
+#
+#
+
+# Commentary:
+#
+# Utility for converting a MOOSE model into NeuroML2. This uses Python
+# libNeuroML.
+#
+#
+
+# Change log:
+#
+# Tue May 21 16:58:03 IST 2013 - Subha moved the code for function
+# fitting to hhfit.py.
+
+#
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; see the file COPYING. If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
+# Floor, Boston, MA 02110-1301, USA.
+#
+#
+
+# Code:
+
+#!!!!! TODO: unit conversion !!!!
+
+try:
+ from future_builtins import zip
+except ImportError:
+ pass
+import traceback
+import warnings
+from collections import deque
+import numpy as np
+from scipy.optimize import curve_fit
+from matplotlib import pyplot as plt
+
+import moose
+from moose.utils import autoposition
+import neuroml
+import hhfit
+
+
+def convert_morphology(root, positions='auto'):
+ """Convert moose neuron morphology contained under `root` into a
+ NeuroML object. The id of the return object is
+ {root.name}_morphology. Each segment object gets the numeric value
+ of the moose id of the object. The name of the segments are same
+ as the corresponding moose compartment.
+
+ Parameters
+ ----------
+ root : a moose element containing a single cell model.
+
+ positions : string
+ flag to indicate if the positions of the end points of the
+ compartments are explicitly available in the compartments or
+ should be automatically generated. Possible values:
+
+ `auto` - automatically generate z coordinates using length of the
+ compartments.
+
+ `explicit` - model has explicit coordinates for all compartments.
+
+ Return
+ ------
+ a neuroml.Morphology instance.
+
+ """
+ if positions == 'auto':
+ queue = deque([autoposition(root)])
+ elif positions == 'explicit':
+ compartments = moose.wildcardFind('%s/##[TYPE=Compartment]' % (root.path))
+ queue = deque([compartment for compartment in map(moose.element, compartments)
+ if len(compartment.neighbours['axial']) == 0])
+ if len(queue) != 1:
+ raise Exception('There must be one and only one top level compartment. Found %d' % (len(topcomp_list)))
+ else:
+ raise Exception('allowed values for keyword argument positions=`auto` or `explicit`')
+ comp_seg = {}
+ parent = None
+ while len(queue) > 0:
+ compartment = queue.popleft()
+ proximal = neuroml.Point3DWithDiam(x=compartment.x0,
+ y=compartment.y0,
+ z=compartment.z0,
+ diameter=compartment.diameter)
+ distal = neuroml.Point3DWithDiam(x=compartment.x,
+ y=compartment.y,
+ z=compartment.z,
+ diameter=compartment.diameter)
+ plist = list(map(moose.element, compartment.neighbours['axial']))
+ try:
+ parent = neuroml.SegmentParent(segments=comp_seg[moose.element(plist[0])].id)
+ except (KeyError, IndexError) as e:
+ parent = None
+ segment = neuroml.Segment(id=compartment.id_.value,
+ proximal=proximal,
+ distal=distal,
+ parent=parent)
+ # TODO: For the time being using numerical value of the moose
+ # id for neuroml id.This needs to be updated for handling
+ # array elements
+ segment.name = compartment.name
+ comp_seg[compartment] = segment
+ queue.extend([comp for comp in map(moose.element, compartment.neighbours['raxial'])])
+ morph = neuroml.Morphology(id='%s_morphology' % (root.name))
+ morph.segments.extend(comp_seg.values())
+ return morph
+
+
+def define_vdep_rate(fn, name):
+ """Define new component type with generic expressions for voltage
+ dependent rate.
+
+ """
+ ctype = neuroml.ComponentType(name)
+ # This is going to be ugly ...
+
+
+
+def convert_hhgate(gate):
+ """Convert a MOOSE gate into GateHHRates in NeuroML"""
+ hh_rates = neuroml.GateHHRates(id=gate.id_.value, name=gate.name)
+ alpha = gate.tableA.copy()
+ beta = gate.tableB - alpha
+ vrange = np.linspace(gate.min, gate.max, len(alpha))
+ afn, ap = hhfit.find_ratefn(vrange, alpha)
+ bfn, bp = hhfit.find_ratefn(vrange, beta)
+ if afn is None:
+ raise Exception('could not find a fitting function for `alpha`')
+ if bfn is None:
+ raise Exception('could not find a fitting function for `alpha`')
+ afn_type = fn_rate_map[afn]
+ afn_component_type = None
+ if afn_type is None:
+ afn_type, afn_component_type = define_component_type(afn)
+ hh_rates.forward_rate = neuroml.HHRate(type=afn_type,
+ midpoint='%gmV' % (ap[2]),
+ scale='%gmV' % (ap[1]),
+ rate='%gper_ms' % (ap[0]))
+ bfn_type = fn_rate_map[bfn]
+ bfn_component_type = None
+ if bfn_type is None:
+ bfn_type, bfn_component_type = define_component_type(bfn)
+ hh_rates.reverse_rate = neuroml.HHRate(type=bfn_type,
+ midpoint='%gmV' % (bp[2]),
+ scale='%gmV' % (bp[1]),
+ rate='%gper_ms' % (bp[0]))
+ return hh_rates, afn_component_type, bfn_component_type
+
+
+def convert_hhchannel(channel):
+ """Convert a moose HHChannel object into a neuroml element.
+
+ TODO: need to check useConcentration option for Ca2+ and V
+ dependent gates. How to handle generic expressions???
+
+ """
+ nml_channel = neuroml.IonChannel(id=channel.id_.value
+ , name=channel.name, type='ionChannelHH'
+ , conductance=channel.Gbar
+ )
+ if channel.Xpower > 0:
+ hh_rate_x = convert_hhgate(channel.gateX[0])
+ hh_rate_x.instances = channel.Xpower
+ nml_channel.gate.append(hh_rate_x)
+
+ if channel.Ypower > 0:
+ hh_rate_y = convert_hhgate(channel.gateY[0])
+ hh_rate_y.instances = channel.Ypower
+ nml_channel.gate.append(hh_rate_y)
+
+ if channel.Zpower > 0:
+ hh_rate_z = convert_hhgate(channel.gateZ[0])
+ hh_rate_y.instances = channel.Zpower
+ nml_channel.gate.append(hh_rate_z)
+ return nml_channel
+
+
+#
+# converter.py ends here
diff --git a/moose-examples/neuroml2/passiveCell.nml b/moose-examples/neuroml2/passiveCell.nml
new file mode 100644
index 0000000000000000000000000000000000000000..9b91f2530dc08dd5a2d50b9e01f133b94552cf7a
--- /dev/null
+++ b/moose-examples/neuroml2/passiveCell.nml
@@ -0,0 +1,52 @@
+<?xml version="1.0" encoding="UTF-8"?>
+
+<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"
+ xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+ xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 ../Schemas/NeuroML2/NeuroML_v2beta4.xsd"
+ id="NML2_SingleCompHHCell">
+
+
+ <ionChannel type="ionChannelPassive" id="passiveChan" conductance="10pS">
+ <notes>Leak conductance</notes>
+ </ionChannel>
+
+ <cell id="passiveCell">
+
+ <morphology id="morph1">
+ <segment id="0" name="soma">
+ <proximal x="0" y="0" z="0" diameter="17.841242"/> <!--Gives a convenient surface area of 1000.0 ?m�-->
+ <distal x="0" y="0" z="0" diameter="17.841242"/>
+ </segment>
+
+ <segmentGroup id="soma_group">
+ <member segment="0"/>
+ </segmentGroup>
+
+ </morphology>
+
+ <biophysicalProperties id="bioPhys1">
+
+ <membraneProperties>
+ <channelDensity id="leak" ionChannel="passiveChan" condDensity="3.0S_per_m2" erev="-54.3mV" ion="non_specific"/>
+ <spikeThresh value="-20mV"/>
+ <specificCapacitance value="1.0uF_per_cm2"/>
+ <initMembPotential value="-66.6mV"/>
+ </membraneProperties>
+
+ <intracellularProperties>
+ <resistivity value="0.03kohm_cm"/> <!-- Note: not used in single compartment simulations -->
+ </intracellularProperties>
+
+ </biophysicalProperties>
+
+ </cell>
+
+ <pulseGenerator id="pulseGen1" delay="50ms" duration="50ms" amplitude="0.08nA"/>
+
+
+ <network id="net1">
+ <population id="pop0" component="passiveCell" size="1"/>
+ <explicitInput target="pop0[0]" input="pulseGen1"/>
+ </network>
+
+</neuroml>
diff --git a/moose-examples/neuroml2/run_cell.py b/moose-examples/neuroml2/run_cell.py
new file mode 100644
index 0000000000000000000000000000000000000000..174358d0b426fc0b2e3fa77232de8a4eb1c036b7
--- /dev/null
+++ b/moose-examples/neuroml2/run_cell.py
@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+# run_cell.py ---
+#
+# Filename: run_cell.py
+# Description:
+# Author:
+# Maintainer: P Gleeson
+# Version:
+# URL:
+# Keywords:
+# Compatibility:
+#
+#
+
+# Commentary:
+#
+#
+#
+#
+
+# Change log:
+# Sunday 16 September 2018 10:04:24 AM IST
+# - Tweaked file to to make it compatible with moose.
+#
+#
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; see the file COPYING. If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
+# Floor, Boston, MA 02110-1301, USA.
+#
+#
+
+# Code:
+
+import moose
+import sys
+import numpy as np
+
+
+def run(nogui):
+
+ filename = 'passiveCell.nml'
+ print('Loading: %s'%filename)
+ reader = moose.mooseReadNML2( filename )
+ assert reader
+ reader.read(filename)
+
+ msoma = reader.getComp(reader.doc.networks[0].populations[0].id,0,0)
+ print(msoma)
+
+
+ data = moose.Neutral('/data')
+
+ pg = reader.getInput('pulseGen1')
+
+ inj = moose.Table('%s/pulse' % (data.path))
+ moose.connect(inj, 'requestOut', pg, 'getOutputValue')
+
+
+ vm = moose.Table('%s/Vm' % (data.path))
+ moose.connect(vm, 'requestOut', msoma, 'getVm')
+
+ simdt = 1e-6
+ plotdt = 1e-4
+ simtime = 150e-3
+
+ if (1):
+ #moose.showmsg( '/clock' )
+ for i in range(8):
+ moose.setClock( i, simdt )
+ moose.setClock( 8, plotdt )
+ moose.reinit()
+ else:
+ utils.resetSim([model.path, data.path], simdt, plotdt, simmethod='ee')
+ moose.showmsg( '/clock' )
+
+ moose.start(simtime)
+
+ print("Finished simulation!")
+
+ t = np.linspace(0, simtime, len(vm.vector))
+
+ if not nogui:
+ import matplotlib.pyplot as plt
+
+ plt.subplot(211)
+ plt.plot(t, vm.vector * 1e3, label='Vm (mV)')
+ plt.legend()
+ plt.title('Vm')
+ plt.subplot(212)
+ plt.title('Input')
+ plt.plot(t, inj.vector * 1e9, label='injected (nA)')
+ #plt.plot(t, gK.vector * 1e6, label='K')
+ #plt.plot(t, gNa.vector * 1e6, label='Na')
+ plt.legend()
+ plt.show()
+ plt.close()
+
+if __name__ == '__main__':
+ nogui = '-nogui' in sys.argv
+ run(nogui)
diff --git a/moose-examples/neuroml2/run_hhcell.py b/moose-examples/neuroml2/run_hhcell.py
new file mode 100644
index 0000000000000000000000000000000000000000..2b38cd67d9f9ff8db70444582f4476d3f265e189
--- /dev/null
+++ b/moose-examples/neuroml2/run_hhcell.py
@@ -0,0 +1,151 @@
+# -*- coding: utf-8 -*-
+# run_hhcell.py ---
+#
+# Filename: run_hhcell.py
+# Description:
+# Author:
+# Maintainer: P Gleeson
+# Version:
+# URL:
+# Keywords:
+# Compatibility:
+#
+#
+
+# Commentary:
+#
+#
+#
+#
+
+# Change log:
+#
+#
+#
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 3, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; see the file COPYING. If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
+# Floor, Boston, MA 02110-1301, USA.
+#
+#
+
+# Code:
+
+import moose
+import sys
+import numpy as np
+
+def test_channel_gates():
+ """Creates prototype channels under `/library` and plots the time
+ constants (tau) and activation (minf, hinf, ninf) parameters for the
+ channel gates.
+
+ """
+ import matplotlib.pyplot as plt
+ lib = moose.Neutral('/library')
+ m = moose.element('/library[0]/naChan[0]/gateX')
+ h = moose.element('/library[0]/naChan[0]/gateY')
+ n = moose.element('/library[0]/kChan[0]/gateX')
+ v = np.linspace(n.min,n.max, n.divs+1)
+
+ plt.subplot(221)
+ plt.plot(v, 1/m.tableB, label='tau_m')
+ plt.plot(v, 1/h.tableB, label='tau_h')
+ plt.plot(v, 1/n.tableB, label='tau_n')
+ plt.legend()
+
+ plt.subplot(222)
+ plt.plot(v, m.tableA/m.tableB, label='m_inf')
+ plt.plot(v, h.tableA/h.tableB, label='h_inf')
+ plt.plot(v, n.tableA/n.tableB, label='n_inf')
+ plt.legend()
+
+ plt.subplot(223)
+ plt.plot(v, m.tableA, label='mA(alpha)')
+ plt.plot(v, h.tableA, label='hA(alpha)')
+ plt.plot(v, n.tableA, label='nA(alpha)')
+ plt.legend()
+ plt.subplot(224)
+
+ plt.plot(v, m.tableB, label='mB')
+ plt.plot(v, m.tableB-m.tableA, label='mB-A(beta)')
+
+ plt.plot(v, h.tableB, label='hB')
+ plt.plot(v, h.tableB-h.tableA, label='hB-A(beta)')
+
+ plt.plot(v, n.tableB, label='nB')
+ plt.plot(v, n.tableB-n.tableA, label='nB-nA(beta)')
+ plt.legend()
+
+ plt.show()
+
+
+def run(nogui):
+ filename = 'NML2_SingleCompHHCell.nml'
+ print('Loading: %s'%filename)
+ reader = moose.mooseReadNML2(filename)
+ msoma = reader.getComp(reader.doc.networks[0].populations[0].id,0,0)
+ print(msoma)
+ data = moose.Neutral('/data')
+ pg = reader.getInput('pulseGen1')
+ inj = moose.Table('%s/pulse' % (data.path))
+ moose.connect(inj, 'requestOut', pg, 'getOutputValue')
+ vm = moose.Table('%s/Vm' % (data.path))
+ moose.connect(vm, 'requestOut', msoma, 'getVm')
+
+ simdt = 1e-6
+ plotdt = 1e-4
+ simtime = 300e-3
+ if (1):
+ #moose.showmsg( '/clock' )
+ for i in range(8):
+ moose.setClock( i, simdt )
+ moose.setClock( 8, plotdt )
+ moose.reinit()
+ else:
+ utils.resetSim([model.path, data.path], simdt, plotdt, simmethod='ee')
+ moose.showmsg( '/clock' )
+ moose.start(simtime)
+
+ print("Finished simulation!")
+
+ t = np.linspace(0, simtime, len(vm.vector))
+
+ if not nogui:
+ import matplotlib.pyplot as plt
+
+ vfile = open('moose_v_hh.dat','w')
+
+ for i in range(len(t)):
+ vfile.write('%s\t%s\n'%(t[i],vm.vector[i]))
+ vfile.close()
+ plt.subplot(211)
+ plt.plot(t, vm.vector * 1e3, label='Vm (mV)')
+ plt.legend()
+ plt.title('Vm')
+ plt.subplot(212)
+ plt.title('Input')
+ plt.plot(t, inj.vector * 1e9, label='injected (nA)')
+ #plt.plot(t, gK.vector * 1e6, label='K')
+ #plt.plot(t, gNa.vector * 1e6, label='Na')
+ plt.legend()
+ plt.figure()
+ test_channel_gates()
+ plt.show()
+ plt.close()
+
+
+if __name__ == '__main__':
+ nogui = '-nogui' in sys.argv
+ run(nogui)
diff --git a/moose-examples/squid/squid.py b/moose-examples/squid/squid.py
index dae8ff9dc5e83e37bd6249a074965f140900cdc7..de3efc46ea16f41d1de9f58c26dda95e1f6c4058 100644
--- a/moose-examples/squid/squid.py
+++ b/moose-examples/squid/squid.py
@@ -28,8 +28,6 @@
# Code:
import sys
-sys.path.append('../../python')
-
import numpy
import moose
@@ -99,22 +97,24 @@ class IonChannel(moose.HHChannel):
def alpha_m(self):
if self.Xpower == 0:
return numpy.array([])
- return numpy.array(moose.HHGate('%s/gateX' % (self.path)).tableA)
+ return numpy.array(moose.element('%s/gateX' % (self.path)).tableA)
@property
def beta_m(self):
if self.Xpower == 0:
return numpy.array([])
- return numpy.array(moose.HHGate('%s/gateX' % (self.path)).tableB) - numpy.array(moose.HHGate('%s/gateX' % (self.path)).tableA)
+ return numpy.array(moose.element('%s/gateX' % (self.path)).tableB) - \
+ numpy.array(moose.element('%s/gateX' % (self.path)).tableA)
@property
def alpha_h(self):
if self.Ypower == 0:
return numpy.array([])
- return numpy.array(moose.HHGate('%s/gateY' % (self.path)).tableA)
+ return numpy.array(moose.element('%s/gateY' % (self.path)).tableA)
@property
def beta_h(self):
if self.Ypower == 0:
return numpy.array([])
- return numpy.array(moose.HHGate('%s/gateY' % (self.path)).tableB) - numpy.array(moose.HHGate('%s/gateY' % (self.path)).tableA)
+ return numpy.array(moose.element('%s/gateY' % (self.path)).tableB) \
+ - numpy.array(moose.element('%s/gateY' % (self.path)).tableA)
class SquidAxon(moose.Compartment):
EREST_ACT = 0.0 # can be -70 mV if not following original HH convention
diff --git a/moose-examples/squid/squid_demo.py b/moose-examples/squid/squid_demo.py
index d243b0f5b03c88b0b84d20029cae8999008702c1..d87d615c922129db69130d691bdc73f35472af6c 100644
--- a/moose-examples/squid/squid_demo.py
+++ b/moose-examples/squid/squid_demo.py
@@ -1,5 +1,4 @@
-
- # -*- coding: utf-8 -*-
+# -*- coding: utf-8 -*-
# squidgui.py ---
#
# Filename: squidgui.py
@@ -48,18 +47,21 @@
# Code:
import sys
-sys.path.append('../../python')
import os
-os.environ['NUMPTHREADS'] = '1'
-
from collections import defaultdict
import time
-from PyQt4 import QtGui
-from PyQt4 import QtCore
+pyqt_ver_ = 4
+try:
+ from PyQt5 import QtGui, QtCore
+ pyqt_ver_ = 5
+except ImportError as e:
+ from PyQt4 import QtGui
+ from PyQt4 import QtCore
import numpy
from matplotlib.figure import Figure
-from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas, NavigationToolbar2QT as NavigationToolbar
+from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
+from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
import moose
@@ -163,7 +165,6 @@ def set_default_line_edit_size(widget):
widget.setMinimumSize(default_line_edit_size)
widget.setMaximumSize(default_line_edit_size)
-
class SquidGui(QtGui.QMainWindow):
defaults = {}
defaults.update(SquidAxon.defaults)
diff --git a/moose-examples/squid/squid_demo_qt5.py b/moose-examples/squid/squid_demo_qt5.py
new file mode 100644
index 0000000000000000000000000000000000000000..59979883d8ff111d951ee70eb95c25f5d2173e47
--- /dev/null
+++ b/moose-examples/squid/squid_demo_qt5.py
@@ -0,0 +1,861 @@
+# -*- coding: utf-8 -*-
+# Description: Squid Model
+# Author: Subha
+# Maintainer: Dilawar Singh <dilawars@ncbs.res.in>
+# Created: Mon Jul 9 18:23:55 2012 (+0530)
+# Version:
+# Last-Updated: Wednesday 12 September 2018 04:23:52 PM IST
+# PyQt5 version
+
+import sys
+import os
+from collections import defaultdict
+import time
+
+from PyQt5 import QtGui, QtCore
+from PyQt5.QtWidgets import QMainWindow, QApplication, QGroupBox, QSizePolicy
+from PyQt5.QtWidgets import QLabel, QLineEdit, QGridLayout, QDockWidget
+from PyQt5.QtWidgets import QCheckBox, QTabWidget, QComboBox, QWidget
+from PyQt5.QtWidgets import QVBoxLayout, QFrame, QHBoxLayout, QAction
+from PyQt5.QtWidgets import QToolButton, QScrollArea, QTextBrowser
+
+import numpy
+from matplotlib.figure import Figure
+from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
+from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
+
+import moose
+
+from squid import *
+from squid_setup import SquidSetup
+from electronics import ClampCircuit
+
+
+tooltip_Nernst = """<h3>Ionic equilibrium potential</h3>
+<p/>
+The equilibrium potential for ion C is given by Nernst equation:
+<p>
+E<sub>C</sub> = (RT/zF) * ln([C]<sub>out</sub> / [C]<sub>in</sub>)
+</p>
+where R is the ideal gas constant (8.3145 J/mol K),<br>
+ T is absolute temperature,<br>
+ z is the valence of the ion,<br>
+ F is Faraday's constant 96480 C/mol,<br>
+ [C]<sub>out</sub> is concentration of C outside the membrane,<br>
+ [C]<sub>in</sub> is concentration of C inside the membrane."""
+
+tooltip_Erest = """<h3>Resting membrane potential</h3>
+<p/>
+The resting membrane potential is determined by the ionic
+concentrations inside and outside the cell membrane and is given by
+the Goldman-Hodgkin-Katz equation:
+<p>
+
+V = (RT/F) * ln((P<sub>K</sub>[K<sup>+</sup>]<sub>out</sub> + P<sub>Na</sub>[Na<sup>+</sup>]<sub>out</sub> + P<sub>Cl</sub>[Cl<sup>-</sup>]<sub>in</sub>) / (P<sub>K</sub>[K<sup>+</sup>]in + P<sub>Na</sub>[Na<sup>+</sup>]<sub>in</sub> + P<sub>Cl</sub>[Cl<sup>-</sup>]<sub>out</sub>))
+
+</p>
+where P<sub>C</sub> is the permeability of the membrane to ion C.
+
+"""
+
+tooltip_NaChan = """<h3>Na+ channel conductance</h3>
+<p/>
+The Na<sup>+</sup> channel conductance in squid giant axon is given by:
+
+<p> G<sub>Na</sub> = GÌ„<sub>Na</sub> * m<sup>3</sup> * h </p>
+
+and the current through this channel is:
+<p>
+I<sub>Na</sub> = G<sub>Na</sub> * (V - E<sub>Na</sub>) = GÌ„<sub>Na</sub> * m<sup>3</sup> * h * (V - E<sub>Na</sub>)
+</p>
+
+where GÌ„<sub>Na</sub> is the peak conductance of Na<sup>+</sup> channel, m is
+the fraction of activation gates open and h is the fraction of
+deactivation gates open. The transition from open to closed state has
+first order kinetics:
+<p> dm/dt = α<sub>m</sub> * ( 1 - m) - β<sub>m</sub> * m </p>
+and similarly for h.
+
+The steady state values are:
+<p> m<sub>∞</sub> = α<sub>m</sub>/(α<sub>m</sub> + β<sub>m</sub>) </p>
+and time constant for steady state is:
+<p>τ<sub>m</sub> = 1/ (α<sub>m</sub> + β<sub>m</sub>) </p>
+and similarly for h.
+"""
+
+tooltip_KChan = """<h3>K+ channel conductance</h3>
+<p/>The K+ channel conductance in squid giant axon is given by:
+
+<p> G<sub>K</sub> = GÌ„<sub>K</sub> * n<sup>4</sup></p>
+
+and the current through this channel is:
+<p>
+I<sub>K</sub> = G<sub>K</sub> * (V - E<sub>K</sub>) = GÌ„<sub>K</sub> * n<sup>4</sup> * (V - E<sub>K</sub>)
+</p>
+where GÌ„<sub>K</sub> is the peak conductance of K<sup>+</sup> channel,
+n is the fraction of activation gates open. The transition from open
+to closed state has first order kinetics: <p> dn/dt = α<sub>n</sub> *
+( 1 - n) - β<sub>n</sub> * n </p>.
+
+The steady state values are:
+<p>
+n<sub>∞</sub> = α<sub>n</sub>/(α<sub>n</sub> + β<sub>n</sub>)
+</p>
+and time constant for steady state is:
+<p>
+τ<sub>n</sub> = 1/ (α<sub>n</sub> + β<sub>n</sub>)
+
+</p>
+and similarly for h.
+"""
+
+tooltip_Im = """<h3>Membrane current</h3>
+<p/>
+The current through the membrane is given by:
+<p>
+I<sub>m</sub> = C<sub>m</sub> dV/dt + I<sub>K</sub> + I<sub>Na</sub> + I<sub>L</sub>
+</p><p>
+ = C<sub>m</sub> dV/dt + G<sub>K</sub>(V, t) * (V - E<sub>K</sub>) + G<sub>Na</sub> * (V - E<sub>Na</sub>) + G<sub>L</sub> * (V - E<sub>L</sub>)
+</p>
+where G<sub>L</sub> is the leak current and E<sub>L</sub> is the leak reversal potential.
+
+"""
+
+default_line_edit_size = QtCore.QSize(80, 25)
+def set_default_line_edit_size(widget):
+ widget.setMinimumSize(default_line_edit_size)
+ widget.setMaximumSize(default_line_edit_size)
+
+class SquidGui( QMainWindow ):
+ defaults = {}
+ defaults.update(SquidAxon.defaults)
+ defaults.update(ClampCircuit.defaults)
+ defaults.update({'runtime': 50.0,
+ 'simdt': 0.01,
+ 'plotdt': 0.1,
+ 'vclamp.holdingV': 0.0,
+ 'vclamp.holdingT': 10.0,
+ 'vclamp.prepulseV': 0.0,
+ 'vclamp.prepulseT': 0.0,
+ 'vclamp.clampV': 50.0,
+ 'vclamp.clampT': 20.0,
+ 'iclamp.baseI': 0.0,
+ 'iclamp.firstI': 0.1,
+ 'iclamp.firstT': 40.0,
+ 'iclamp.firstD': 5.0,
+ 'iclamp.secondI': 0.0,
+ 'iclamp.secondT': 0.0,
+ 'iclamp.secondD': 0.0
+ })
+ def __init__(self, *args):
+ QMainWindow.__init__(self, *args)
+ self.squid_setup = SquidSetup()
+ self._plotdt = SquidGui.defaults['plotdt']
+ self._plot_dict = defaultdict(list)
+ self.setWindowTitle('Squid Axon simulation')
+ self.setDockNestingEnabled(True)
+ self._createRunControl()
+ self.addDockWidget(QtCore.Qt.LeftDockWidgetArea, self._runControlDock)
+ self._runControlDock.setFeatures(QDockWidget.AllDockWidgetFeatures)
+ self._createChannelControl()
+ self._channelCtrlBox.setWindowTitle('Channel properties')
+ self._channelControlDock.setFeatures(QDockWidget.AllDockWidgetFeatures)
+ self.addDockWidget(QtCore.Qt.LeftDockWidgetArea, self._channelControlDock)
+ self._createElectronicsControl()
+ self._electronicsDock.setFeatures(QDockWidget.AllDockWidgetFeatures)
+ self._electronicsDock.setWindowTitle('Electronics')
+ self.addDockWidget(QtCore.Qt.LeftDockWidgetArea, self._electronicsDock)
+ self._createPlotWidget()
+ self.setCentralWidget(self._plotWidget)
+ self._createStatePlotWidget()
+ self._createHelpMessage()
+ self._helpWindow.setVisible(False)
+ self._statePlotWidget.setWindowFlags(QtCore.Qt.Window)
+ self._statePlotWidget.setWindowTitle('State plot')
+ self._initActions()
+ self._createRunToolBar()
+ self._createPlotToolBar()
+
+ def getFloatInput(self, widget, name):
+ try:
+ return float(str(widget.text()))
+ except ValueError:
+ QMessageBox.critical(self, 'Invalid input', 'Please enter a valid number for {}'.format(name))
+ raise
+
+
+ def _createPlotWidget(self):
+ self._plotWidget = QWidget()
+ self._plotFigure = Figure()
+ self._plotCanvas = FigureCanvas(self._plotFigure)
+ self._plotCanvas.setSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding)
+ self._plotCanvas.updateGeometry()
+ self._plotCanvas.setParent(self._plotWidget)
+ self._plotCanvas.mpl_connect('scroll_event', self._onScroll)
+ self._plotFigure.set_canvas(self._plotCanvas)
+ # Vm and command voltage go in the same subplot
+ self._vm_axes = self._plotFigure.add_subplot(2,2,1, title='Membrane potential')
+ self._vm_axes.set_ylim(-20.0, 120.0)
+ # Channel conductances go to the same subplot
+ self._g_axes = self._plotFigure.add_subplot(2,2,2, title='Channel conductance')
+ self._g_axes.set_ylim(0.0, 0.5)
+ # Injection current for Vclamp/Iclamp go to the same subplot
+ self._im_axes = self._plotFigure.add_subplot(2,2,3, title='Injection current')
+ self._im_axes.set_ylim(-0.5, 0.5)
+ # Channel currents go to the same subplot
+ self._i_axes = self._plotFigure.add_subplot(2,2,4, title='Channel current')
+ self._i_axes.set_ylim(-10, 10)
+ for axis in self._plotFigure.axes:
+ axis.set_autoscale_on(False)
+ layout = QVBoxLayout()
+ layout.addWidget(self._plotCanvas)
+ self._plotNavigator = NavigationToolbar(self._plotCanvas, self._plotWidget)
+ layout.addWidget(self._plotNavigator)
+ self._plotWidget.setLayout(layout)
+
+ def _createStatePlotWidget(self):
+ self._statePlotWidget = QWidget()
+ self._statePlotFigure = Figure()
+ self._statePlotCanvas = FigureCanvas(self._statePlotFigure)
+ self._statePlotCanvas.setSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding)
+ self._statePlotCanvas.updateGeometry()
+ self._statePlotCanvas.setParent(self._statePlotWidget)
+ self._statePlotFigure.set_canvas(self._statePlotCanvas)
+ self._statePlotFigure.subplots_adjust(hspace=0.5)
+ self._statePlotAxes = self._statePlotFigure.add_subplot(2,1,1, title='State plot')
+ self._state_plot, = self._statePlotAxes.plot([], [], label='state')
+ self._activationParamAxes = self._statePlotFigure.add_subplot(2,1,2, title='H-H activation parameters vs time')
+ self._activationParamAxes.set_xlabel('Time (ms)')
+ #for axis in self._plotFigure.axes:
+ # axis.autoscale(False)
+ self._stateplot_xvar_label = QLabel('Variable on X-axis')
+ self._stateplot_xvar_combo = QComboBox()
+ self._stateplot_xvar_combo.addItems(['V', 'm', 'n', 'h'])
+ self._stateplot_xvar_combo.setCurrentIndex(0)
+ self._stateplot_xvar_combo.setEditable(False)
+ # self.connect(self._stateplot_xvar_combo,
+ # QtCore.SIGNAL('currentIndexChanged(const QString&)'),
+ # self._statePlotXSlot)
+ self._stateplot_xvar_combo.currentIndexChanged.connect( self._statePlotXSlot )
+ self._stateplot_yvar_label = QLabel('Variable on Y-axis')
+ self._stateplot_yvar_combo = QComboBox()
+ self._stateplot_yvar_combo.addItems(['V', 'm', 'n', 'h'])
+ self._stateplot_yvar_combo.setCurrentIndex(2)
+ self._stateplot_yvar_combo.setEditable(False)
+ self._stateplot_yvar_combo.currentIndexChanged.connect(self._statePlotYSlot)
+ self._statePlotNavigator = NavigationToolbar(self._statePlotCanvas, self._statePlotWidget)
+ frame = QFrame()
+ frame.setFrameStyle(QFrame.StyledPanel + QFrame.Raised)
+ layout = QHBoxLayout()
+ layout.addWidget(self._stateplot_xvar_label)
+ layout.addWidget(self._stateplot_xvar_combo)
+ layout.addWidget(self._stateplot_yvar_label)
+ layout.addWidget(self._stateplot_yvar_combo)
+ frame.setLayout(layout)
+ self._closeStatePlotAction = QAction('Close', self)
+ self._closeStatePlotAction.triggered.connect(self._statePlotWidget.close)
+ self._closeStatePlotButton = QToolButton()
+ self._closeStatePlotButton.setDefaultAction(self._closeStatePlotAction)
+ layout = QVBoxLayout()
+ layout.addWidget(frame)
+ layout.addWidget(self._statePlotCanvas)
+ layout.addWidget(self._statePlotNavigator)
+ layout.addWidget(self._closeStatePlotButton)
+ self._statePlotWidget.setLayout(layout)
+ # Setting the close event so that when the help window is
+ # closed the ``State plot`` button becomes unchecked
+ self._statePlotWidget.closeEvent = lambda event: self._showStatePlotAction.setChecked(False)
+
+ def _createRunControl(self):
+ self._runControlBox = QGroupBox(self)
+ self._runControlBox.setSizePolicy(QSizePolicy.Preferred, QSizePolicy.Preferred)
+ self._runTimeLabel = QLabel("Run time (ms)", self._runControlBox)
+ self._simTimeStepLabel = QLabel("Simulation time step (ms)", self._runControlBox)
+ self._runTimeEdit = QLineEdit('%g' % (SquidGui.defaults['runtime']), self._runControlBox)
+ set_default_line_edit_size(self._runTimeEdit)
+ self._simTimeStepEdit = QLineEdit('%g' % (SquidGui.defaults['simdt']), self._runControlBox)
+ set_default_line_edit_size(self._simTimeStepEdit)
+ layout = QGridLayout()
+ layout.addWidget(self._runTimeLabel, 0, 0)
+ layout.addWidget(self._runTimeEdit, 0, 1)
+ layout.addWidget(self._simTimeStepLabel, 1, 0)
+ layout.addWidget(self._simTimeStepEdit, 1, 1)
+ layout.setColumnStretch(2, 1.0)
+ layout.setRowStretch(2, 1.0)
+ self._runControlBox.setLayout(layout)
+ self._runControlDock = QDockWidget('Simulation', self)
+ self._runControlDock.setWidget(self._runControlBox)
+
+ def _createChannelControl(self):
+ self._channelControlDock = QDockWidget('Channels', self)
+ self._channelCtrlBox = QGroupBox(self)
+ self._naConductanceToggle = QCheckBox('Block Na+ channel', self._channelCtrlBox)
+ self._naConductanceToggle.setToolTip('<html>%s</html>' % (tooltip_NaChan))
+ self._kConductanceToggle = QCheckBox('Block K+ channel', self._channelCtrlBox)
+ self._kConductanceToggle.setToolTip('<html>%s</html>' % (tooltip_KChan))
+ self._kOutLabel = QLabel('[K+]out (mM)', self._channelCtrlBox)
+ self._kOutEdit = QLineEdit('%g' % (self.squid_setup.squid_axon.K_out),
+ self._channelCtrlBox)
+ self._kOutLabel.setToolTip('<html>%s</html>' % (tooltip_Nernst))
+ self._kOutEdit.setToolTip('<html>%s</html>' % (tooltip_Nernst))
+ set_default_line_edit_size(self._kOutEdit)
+ self._naOutLabel = QLabel('[Na+]out (mM)', self._channelCtrlBox)
+ self._naOutEdit = QLineEdit('%g' % (self.squid_setup.squid_axon.Na_out),
+ self._channelCtrlBox)
+ self._naOutLabel.setToolTip('<html>%s</html>' % (tooltip_Nernst))
+ self._naOutEdit.setToolTip('<html>%s</html>' % (tooltip_Nernst))
+ set_default_line_edit_size(self._naOutEdit)
+ self._kInLabel = QLabel('[K+]in (mM)', self._channelCtrlBox)
+ self._kInEdit = QLineEdit('%g' % (self.squid_setup.squid_axon.K_in),
+ self._channelCtrlBox)
+ self._kInEdit.setToolTip(tooltip_Nernst)
+ self._naInLabel = QLabel('[Na+]in (mM)', self._channelCtrlBox)
+ self._naInEdit = QLineEdit('%g' % (self.squid_setup.squid_axon.Na_in),
+ self._channelCtrlBox)
+ self._naInEdit.setToolTip('<html>%s</html>' % (tooltip_Nernst))
+ self._temperatureLabel = QLabel('Temperature (C)', self._channelCtrlBox)
+ self._temperatureEdit = QLineEdit('%g' % (self.defaults['temperature'] - CELSIUS_TO_KELVIN),
+ self._channelCtrlBox)
+ self._temperatureEdit.setToolTip('<html>%s</html>' % (tooltip_Nernst))
+ set_default_line_edit_size(self._temperatureEdit)
+ for child in self._channelCtrlBox.children():
+ if isinstance(child, QLineEdit):
+ set_default_line_edit_size(child)
+ layout = QGridLayout(self._channelCtrlBox)
+ layout.addWidget(self._naConductanceToggle, 0, 0)
+ layout.addWidget(self._kConductanceToggle, 1, 0)
+ layout.addWidget(self._naOutLabel, 2, 0)
+ layout.addWidget(self._naOutEdit, 2, 1)
+ layout.addWidget(self._naInLabel, 3, 0)
+ layout.addWidget(self._naInEdit, 3, 1)
+ layout.addWidget(self._kOutLabel, 4, 0)
+ layout.addWidget(self._kOutEdit, 4, 1)
+ layout.addWidget(self._kInLabel, 5, 0)
+ layout.addWidget(self._kInEdit, 5, 1)
+ layout.addWidget(self._temperatureLabel, 6, 0)
+ layout.addWidget(self._temperatureEdit, 6, 1)
+ layout.setRowStretch(7, 1.0)
+ self._channelCtrlBox.setLayout(layout)
+ self._channelControlDock.setWidget(self._channelCtrlBox)
+ return self._channelCtrlBox
+
+ def __get_stateplot_data(self, name):
+ data = []
+ if name == 'V':
+ data = self.squid_setup.vm_table.vector
+ elif name == 'm':
+ data = self.squid_setup.m_table.vector
+ elif name == 'h':
+ data = self.squid_setup.h_table.vector
+ elif name == 'n':
+ data = self.squid_setup.n_table.vector
+ else:
+ raise ValueError('Unrecognized selection: %s' % (name))
+ return numpy.asarray(data)
+
+ def _statePlotYSlot(self, selectedItem):
+ ydata = self.__get_stateplot_data(str(selectedItem))
+ self._state_plot.set_ydata(ydata)
+ self._statePlotAxes.set_ylabel(selectedItem)
+ if str(selectedItem) == 'V':
+ self._statePlotAxes.set_ylim(-20, 120)
+ else:
+ self._statePlotAxes.set_ylim(0, 1)
+ self._statePlotCanvas.draw()
+
+ def _statePlotXSlot(self, selectedItem):
+ xdata = self.__get_stateplot_data(str(selectedItem))
+ self._state_plot.set_xdata(xdata)
+ self._statePlotAxes.set_xlabel(selectedItem)
+ if str(selectedItem) == 'V':
+ self._statePlotAxes.set_xlim(-20, 120)
+ else:
+ self._statePlotAxes.set_xlim(0, 1)
+ self._statePlotCanvas.draw()
+
+ def _createElectronicsControl(self):
+ """Creates a tabbed widget of voltage clamp and current clamp controls"""
+ self._electronicsTab = QTabWidget(self)
+ self._electronicsTab.addTab(self._getIClampCtrlBox(), 'Current clamp')
+ self._electronicsTab.addTab(self._getVClampCtrlBox(), 'Voltage clamp')
+ self._electronicsDock = QDockWidget(self)
+ self._electronicsDock.setWidget(self._electronicsTab)
+
+ def _getVClampCtrlBox(self):
+ vClampPanel = QGroupBox(self)
+ self._vClampCtrlBox = vClampPanel
+ self._holdingVLabel = QLabel("Holding Voltage (mV)", vClampPanel)
+ self._holdingVEdit = QLineEdit('%g' % (SquidGui.defaults['vclamp.holdingV']), vClampPanel)
+ self._holdingTimeLabel = QLabel("Holding Time (ms)", vClampPanel)
+ self._holdingTimeEdit = QLineEdit('%g' % (SquidGui.defaults['vclamp.holdingT']), vClampPanel)
+ self._prePulseVLabel = QLabel("Pre-pulse Voltage (mV)", vClampPanel)
+ self._prePulseVEdit = QLineEdit('%g' % (SquidGui.defaults['vclamp.prepulseV']), vClampPanel)
+ self._prePulseTimeLabel = QLabel("Pre-pulse Time (ms)", vClampPanel)
+ self._prePulseTimeEdit = QLineEdit('%g' % (SquidGui.defaults['vclamp.prepulseT']), vClampPanel)
+ self._clampVLabel = QLabel("Clamp Voltage (mV)", vClampPanel)
+ self._clampVEdit = QLineEdit('%g' % (SquidGui.defaults['vclamp.clampV']), vClampPanel)
+ self._clampTimeLabel = QLabel("Clamp Time (ms)", vClampPanel)
+ self._clampTimeEdit = QLineEdit('%g' % (SquidGui.defaults['vclamp.clampT']), vClampPanel)
+ for child in vClampPanel.children():
+ if isinstance(child, QLineEdit):
+ set_default_line_edit_size(child)
+ layout = QGridLayout(vClampPanel)
+ layout.addWidget(self._holdingVLabel, 0, 0)
+ layout.addWidget(self._holdingVEdit, 0, 1)
+ layout.addWidget(self._holdingTimeLabel, 1, 0)
+ layout.addWidget(self._holdingTimeEdit, 1, 1)
+ layout.addWidget(self._prePulseVLabel, 2, 0)
+ layout.addWidget(self._prePulseVEdit, 2, 1)
+ layout.addWidget(self._prePulseTimeLabel,3,0)
+ layout.addWidget(self._prePulseTimeEdit, 3, 1)
+ layout.addWidget(self._clampVLabel, 4, 0)
+ layout.addWidget(self._clampVEdit, 4, 1)
+ layout.addWidget(self._clampTimeLabel, 5, 0)
+ layout.addWidget(self._clampTimeEdit, 5, 1)
+ layout.setRowStretch(6, 1.0)
+ vClampPanel.setLayout(layout)
+ return self._vClampCtrlBox
+
+ def _getIClampCtrlBox(self):
+ iClampPanel = QGroupBox(self)
+ self._iClampCtrlBox = iClampPanel
+ self._baseCurrentLabel = QLabel("Base Current Level (uA)",iClampPanel)
+ self._baseCurrentEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.baseI']),iClampPanel)
+ self._firstPulseLabel = QLabel("First Pulse Current (uA)", iClampPanel)
+ self._firstPulseEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.firstI']), iClampPanel)
+ self._firstDelayLabel = QLabel("First Onset Delay (ms)", iClampPanel)
+ self._firstDelayEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.firstD']),iClampPanel)
+ self._firstPulseWidthLabel = QLabel("First Pulse Width (ms)", iClampPanel)
+ self._firstPulseWidthEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.firstT']), iClampPanel)
+ self._secondPulseLabel = QLabel("Second Pulse Current (uA)", iClampPanel)
+ self._secondPulseEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.secondI']), iClampPanel)
+ self._secondDelayLabel = QLabel("Second Onset Delay (ms)", iClampPanel)
+ self._secondDelayEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.secondD']),iClampPanel)
+ self._secondPulseWidthLabel = QLabel("Second Pulse Width (ms)", iClampPanel)
+ self._secondPulseWidthEdit = QLineEdit('%g' % (SquidGui.defaults['iclamp.secondT']), iClampPanel)
+ self._pulseMode = QComboBox(iClampPanel)
+ self._pulseMode.addItem("Single Pulse")
+ self._pulseMode.addItem("Pulse Train")
+ for child in iClampPanel.children():
+ if isinstance(child, QLineEdit):
+ set_default_line_edit_size(child)
+ layout = QGridLayout(iClampPanel)
+ layout.addWidget(self._baseCurrentLabel, 0, 0)
+ layout.addWidget(self._baseCurrentEdit, 0, 1)
+ layout.addWidget(self._firstPulseLabel, 1, 0)
+ layout.addWidget(self._firstPulseEdit, 1, 1)
+ layout.addWidget(self._firstDelayLabel, 2, 0)
+ layout.addWidget(self._firstDelayEdit, 2, 1)
+ layout.addWidget(self._firstPulseWidthLabel, 3, 0)
+ layout.addWidget(self._firstPulseWidthEdit, 3, 1)
+ layout.addWidget(self._secondPulseLabel, 4, 0)
+ layout.addWidget(self._secondPulseEdit, 4, 1)
+ layout.addWidget(self._secondDelayLabel, 5, 0)
+ layout.addWidget(self._secondDelayEdit, 5, 1)
+ layout.addWidget(self._secondPulseWidthLabel, 6, 0)
+ layout.addWidget(self._secondPulseWidthEdit, 6, 1)
+ layout.addWidget(self._pulseMode, 7, 0, 1, 2)
+ layout.setRowStretch(8, 1.0)
+ # layout.setSizeConstraint(QLayout.SetFixedSize)
+ iClampPanel.setLayout(layout)
+ return self._iClampCtrlBox
+
+ def _overlayPlots(self, overlay):
+ if not overlay:
+ for axis in (self._plotFigure.axes + self._statePlotFigure.axes):
+ title = axis.get_title()
+ axis.clear()
+ axis.set_title(title)
+ suffix = ''
+ else:
+ suffix = '_%d' % (len(self._plot_dict['vm']))
+ self._vm_axes.set_xlim(0.0, self._runtime)
+ self._g_axes.set_xlim(0.0, self._runtime)
+ self._im_axes.set_xlim(0.0, self._runtime)
+ self._i_axes.set_xlim(0.0, self._runtime)
+ self._vm_plot, = self._vm_axes.plot([], [], label='Vm%s'%(suffix))
+ self._plot_dict['vm'].append(self._vm_plot)
+ self._command_plot, = self._vm_axes.plot([], [], label='command%s'%(suffix))
+ self._plot_dict['command'].append(self._command_plot)
+ # Channel conductances go to the same subplot
+ self._gna_plot, = self._g_axes.plot([], [], label='Na%s'%(suffix))
+ self._plot_dict['gna'].append(self._gna_plot)
+ self._gk_plot, = self._g_axes.plot([], [], label='K%s'%(suffix))
+ self._plot_dict['gk'].append(self._gk_plot)
+ # Injection current for Vclamp/Iclamp go to the same subplot
+ self._iclamp_plot, = self._im_axes.plot([], [], label='Iclamp%s'%(suffix))
+ self._vclamp_plot, = self._im_axes.plot([], [], label='Vclamp%s'%(suffix))
+ self._plot_dict['iclamp'].append(self._iclamp_plot)
+ self._plot_dict['vclamp'].append(self._vclamp_plot)
+ # Channel currents go to the same subplot
+ self._ina_plot, = self._i_axes.plot([], [], label='Na%s'%(suffix))
+ self._plot_dict['ina'].append(self._ina_plot)
+ self._ik_plot, = self._i_axes.plot([], [], label='K%s'%(suffix))
+ self._plot_dict['ik'].append(self._ik_plot)
+ # self._i_axes.legend()
+ # State plots
+ self._state_plot, = self._statePlotAxes.plot([], [], label='state%s'%(suffix))
+ self._plot_dict['state'].append(self._state_plot)
+ self._m_plot, = self._activationParamAxes.plot([],[], label='m%s'%(suffix))
+ self._h_plot, = self._activationParamAxes.plot([], [], label='h%s'%(suffix))
+ self._n_plot, = self._activationParamAxes.plot([], [], label='n%s'%(suffix))
+ self._plot_dict['m'].append(self._m_plot)
+ self._plot_dict['h'].append(self._h_plot)
+ self._plot_dict['n'].append(self._n_plot)
+ if self._showLegendAction.isChecked():
+ for axis in (self._plotFigure.axes + self._statePlotFigure.axes):
+ axis.legend()
+
+ def _updateAllPlots(self):
+ self._updatePlots()
+ self._updateStatePlot()
+
+ def _updatePlots(self):
+ if len(self.squid_setup.vm_table.vector) <= 0:
+ return
+ vm = numpy.asarray(self.squid_setup.vm_table.vector)
+ cmd = numpy.asarray(self.squid_setup.cmd_table.vector)
+ ik = numpy.asarray(self.squid_setup.ik_table.vector)
+ ina = numpy.asarray(self.squid_setup.ina_table.vector)
+ iclamp = numpy.asarray(self.squid_setup.iclamp_table.vector)
+ vclamp = numpy.asarray(self.squid_setup.vclamp_table.vector)
+ gk = numpy.asarray(self.squid_setup.gk_table.vector)
+ gna = numpy.asarray(self.squid_setup.gna_table.vector)
+ time_series = numpy.linspace(0, self._plotdt * len(vm), len(vm))
+ self._vm_plot.set_data(time_series, vm)
+ time_series = numpy.linspace(0, self._plotdt * len(cmd), len(cmd))
+ self._command_plot.set_data(time_series, cmd)
+ time_series = numpy.linspace(0, self._plotdt * len(ik), len(ik))
+ self._ik_plot.set_data(time_series, ik)
+ time_series = numpy.linspace(0, self._plotdt * len(ina), len(ina))
+ self._ina_plot.set_data(time_series, ina)
+ time_series = numpy.linspace(0, self._plotdt * len(iclamp), len(iclamp))
+ self._iclamp_plot.set_data(time_series, iclamp)
+ time_series = numpy.linspace(0, self._plotdt * len(vclamp), len(vclamp))
+ self._vclamp_plot.set_data(time_series, vclamp)
+ time_series = numpy.linspace(0, self._plotdt * len(gk), len(gk))
+ self._gk_plot.set_data(time_series, gk)
+ time_series = numpy.linspace(0, self._plotdt * len(gna), len(gna))
+ self._gna_plot.set_data(time_series, gna)
+ # self._vm_axes.margins(y=0.1)
+ # self._g_axes.margin(y=0.1)
+ # self._im_axes.margins(y=0.1)
+ # self._i_axes.margins(y=0.1)
+ if self._autoscaleAction.isChecked():
+ for axis in self._plotFigure.axes:
+ axis.relim()
+ axis.margins(0.1, 0.1)
+ axis.autoscale_view(tight=True)
+ else:
+ self._vm_axes.set_ylim(-20.0, 120.0)
+ self._g_axes.set_ylim(0.0, 0.5)
+ self._im_axes.set_ylim(-0.5, 0.5)
+ self._i_axes.set_ylim(-10, 10)
+ self._vm_axes.set_xlim(0.0, time_series[-1])
+ self._g_axes.set_xlim(0.0, time_series[-1])
+ self._im_axes.set_xlim(0.0, time_series[-1])
+ self._i_axes.set_xlim(0.0, time_series[-1])
+ self._plotCanvas.draw()
+
+ def _updateStatePlot(self):
+ if len(self.squid_setup.vm_table.vector) <= 0:
+ return
+ sx = str(self._stateplot_xvar_combo.currentText())
+ sy = str(self._stateplot_yvar_combo.currentText())
+ xdata = self.__get_stateplot_data(sx)
+ ydata = self.__get_stateplot_data(sy)
+ minlen = min(len(xdata), len(ydata))
+ self._state_plot.set_data(xdata[:minlen], ydata[:minlen])
+ self._statePlotAxes.set_xlabel(sx)
+ self._statePlotAxes.set_ylabel(sy)
+ if sx == 'V':
+ self._statePlotAxes.set_xlim(-20, 120)
+ else:
+ self._statePlotAxes.set_xlim(0, 1)
+ if sy == 'V':
+ self._statePlotAxes.set_ylim(-20, 120)
+ else:
+ self._statePlotAxes.set_ylim(0, 1)
+ self._activationParamAxes.set_xlim(0, self._runtime)
+ m = self.__get_stateplot_data('m')
+ n = self.__get_stateplot_data('n')
+ h = self.__get_stateplot_data('h')
+ time_series = numpy.linspace(0, self._plotdt*len(m), len(m))
+ self._m_plot.set_data(time_series, m)
+ time_series = numpy.linspace(0, self._plotdt*len(h), len(h))
+ self._h_plot.set_data(time_series, h)
+ time_series = numpy.linspace(0, self._plotdt*len(n), len(n))
+ self._n_plot.set_data(time_series, n)
+ if self._autoscaleAction.isChecked():
+ for axis in self._statePlotFigure.axes:
+ axis.relim()
+ axis.set_autoscale_on(True)
+ axis.autoscale_view(True)
+ self._statePlotCanvas.draw()
+
+ def _runSlot(self):
+ if moose.isRunning():
+ print('Stopping simulation in progress ...')
+ moose.stop()
+ self._runtime = self.getFloatInput(self._runTimeEdit, self._runTimeLabel.text())
+ self._overlayPlots(self._overlayAction.isChecked())
+ self._simdt = self.getFloatInput(self._simTimeStepEdit, self._simTimeStepLabel.text())
+ clampMode = None
+ singlePulse = True
+ if self._electronicsTab.currentWidget() == self._vClampCtrlBox:
+ clampMode = 'vclamp'
+ baseLevel = self.getFloatInput(self._holdingVEdit, self._holdingVLabel.text())
+ firstDelay = self.getFloatInput(self._holdingTimeEdit, self._holdingTimeLabel.text())
+ firstWidth = self.getFloatInput(self._prePulseTimeEdit, self._prePulseTimeLabel.text())
+ firstLevel = self.getFloatInput(self._prePulseVEdit, self._prePulseVLabel.text())
+ secondDelay = firstWidth
+ secondWidth = self.getFloatInput(self._clampTimeEdit, self._clampTimeLabel.text())
+ secondLevel = self.getFloatInput(self._clampVEdit, self._clampVLabel.text())
+ if not self._autoscaleAction.isChecked():
+ self._im_axes.set_ylim(-10.0, 10.0)
+ else:
+ clampMode = 'iclamp'
+ baseLevel = self.getFloatInput(self._baseCurrentEdit, self._baseCurrentLabel.text())
+ firstDelay = self.getFloatInput(self._firstDelayEdit, self._firstDelayLabel.text())
+ firstWidth = self.getFloatInput(self._firstPulseWidthEdit, self._firstPulseWidthLabel.text())
+ firstLevel = self.getFloatInput(self._firstPulseEdit, self._firstPulseLabel.text())
+ secondDelay = self.getFloatInput(self._secondDelayEdit, self._secondDelayLabel.text())
+ secondLevel = self.getFloatInput(self._secondPulseEdit, self._secondPulseLabel.text())
+ secondWidth = self.getFloatInput(self._secondPulseWidthEdit, self._secondPulseWidthLabel.text())
+ singlePulse = (self._pulseMode.currentIndex() == 0)
+ if not self._autoscaleAction.isChecked():
+ self._im_axes.set_ylim(-0.4, 0.4)
+ self.squid_setup.clamp_ckt.configure_pulses(baseLevel=baseLevel,
+ firstDelay=firstDelay,
+ firstWidth=firstWidth,
+ firstLevel=firstLevel,
+ secondDelay=secondDelay,
+ secondWidth=secondWidth,
+ secondLevel=secondLevel,
+ singlePulse=singlePulse)
+ if self._kConductanceToggle.isChecked():
+ self.squid_setup.squid_axon.specific_gK = 0.0
+ else:
+ self.squid_setup.squid_axon.specific_gK = SquidAxon.defaults['specific_gK']
+ if self._naConductanceToggle.isChecked():
+ self.squid_setup.squid_axon.specific_gNa = 0.0
+ else:
+ self.squid_setup.squid_axon.specific_gNa = SquidAxon.defaults['specific_gNa']
+ self.squid_setup.squid_axon.celsius = self.getFloatInput(self._temperatureEdit, self._temperatureLabel.text())
+ self.squid_setup.squid_axon.K_out = self.getFloatInput(self._kOutEdit, self._kOutLabel.text())
+ self.squid_setup.squid_axon.Na_out = self.getFloatInput(self._naOutEdit, self._naOutLabel.text())
+ self.squid_setup.squid_axon.K_in = self.getFloatInput(self._kInEdit, self._kInLabel.text())
+ self.squid_setup.squid_axon.Na_in = self.getFloatInput(self._naInEdit, self._naInLabel.text())
+ self.squid_setup.squid_axon.updateEk()
+ self.squid_setup.schedule(self._simdt, self._plotdt, clampMode)
+ # The following line is for use with Qthread
+ self.squid_setup.run(self._runtime)
+ self._updateAllPlots()
+
+ def _toggleDocking(self, on):
+ self._channelControlDock.setFloating(on)
+ self._electronicsDock.setFloating(on)
+ self._runControlDock.setFloating(on)
+
+ def _restoreDocks(self):
+ self._channelControlDock.setVisible(True)
+ self._electronicsDock.setVisible(True)
+ self._runControlDock.setVisible(True)
+
+ def _initActions(self):
+ self._runAction = QAction(self.tr('Run'), self)
+ self._runAction.setShortcut(self.tr('F5'))
+ self._runAction.setToolTip('Run simulation (F5)')
+ self._runAction.triggered.connect( self._runSlot)
+ self._resetToDefaultsAction = QAction(self.tr('Restore defaults'), self)
+ self._resetToDefaultsAction.setToolTip('Reset all settings to their default values')
+ self._resetToDefaultsAction.triggered.connect( self._useDefaults)
+ self._showLegendAction = QAction(self.tr('Display legend'), self)
+ self._showLegendAction.setCheckable(True)
+ self._showLegendAction.toggled.connect(self._showLegend)
+ self._showStatePlotAction = QAction(self.tr('State plot'), self)
+ self._showStatePlotAction.setCheckable(True)
+ self._showStatePlotAction.setChecked(False)
+ self._showStatePlotAction.toggled.connect(self._statePlotWidget.setVisible)
+ self._autoscaleAction = QAction(self.tr('Auto-scale plots'), self)
+ self._autoscaleAction.setCheckable(True)
+ self._autoscaleAction.setChecked(False)
+ self._autoscaleAction.toggled.connect(self._autoscale)
+ self._overlayAction = QAction('Overlay plots', self)
+ self._overlayAction.setCheckable(True)
+ self._overlayAction.setChecked(False)
+ self._dockAction = QAction('Undock all', self)
+ self._dockAction.setCheckable(True)
+ self._dockAction.setChecked(False)
+ # self._dockAction.toggle.connect( self._toggleDocking)
+ self._restoreDocksAction = QAction('Show all', self)
+ self._restoreDocksAction.triggered.connect( self._restoreDocks)
+ self._quitAction = QAction(self.tr('&Quit'), self)
+ self._quitAction.setShortcut(self.tr('Ctrl+Q'))
+ self._quitAction.triggered.connect(qApp.closeAllWindows)
+
+
+
+ def _createRunToolBar(self):
+ self._simToolBar = self.addToolBar(self.tr('Simulation control'))
+ self._simToolBar.addAction(self._quitAction)
+ self._simToolBar.addAction(self._runAction)
+ self._simToolBar.addAction(self._resetToDefaultsAction)
+ self._simToolBar.addAction(self._dockAction)
+ self._simToolBar.addAction(self._restoreDocksAction)
+
+ def _createPlotToolBar(self):
+ self._plotToolBar = self.addToolBar(self.tr('Plotting control'))
+ self._plotToolBar.addAction(self._showLegendAction)
+ self._plotToolBar.addAction(self._autoscaleAction)
+ self._plotToolBar.addAction(self._overlayAction)
+ self._plotToolBar.addAction(self._showStatePlotAction)
+ self._plotToolBar.addAction(self._helpAction)
+ self._plotToolBar.addAction(self._helpBiophysicsAction)
+
+ def _showLegend(self, on):
+ if on:
+ for axis in (self._plotFigure.axes + self._statePlotFigure.axes):
+ axis.legend().set_visible(True)
+ else:
+ for axis in (self._plotFigure.axes + self._statePlotFigure.axes):
+ axis.legend().set_visible(False)
+ self._plotCanvas.draw()
+ self._statePlotCanvas.draw()
+
+ def _autoscale(self, on):
+ if on:
+ for axis in (self._plotFigure.axes + self._statePlotFigure.axes):
+ axis.relim()
+ axis.set_autoscale_on(True)
+ axis.autoscale_view(True)
+ else:
+ for axis in self._plotFigure.axes:
+ axis.set_autoscale_on(False)
+ self._vm_axes.set_ylim(-20.0, 120.0)
+ self._g_axes.set_ylim(0.0, 0.5)
+ self._im_axes.set_ylim(-0.5, 0.5)
+ self._i_axes.set_ylim(-10, 10)
+ self._plotCanvas.draw()
+ self._statePlotCanvas.draw()
+
+ def _useDefaults(self):
+ self._runTimeEdit.setText('%g' % (self.defaults['runtime']))
+ self._simTimeStepEdit.setText('%g' % (self.defaults['simdt']))
+ self._overlayAction.setChecked(False)
+ self._naConductanceToggle.setChecked(False)
+ self._kConductanceToggle.setChecked(False)
+ self._kOutEdit.setText('%g' % (SquidGui.defaults['K_out']))
+ self._naOutEdit.setText('%g' % (SquidGui.defaults['Na_out']))
+ self._kInEdit.setText('%g' % (SquidGui.defaults['K_in']))
+ self._naInEdit.setText('%g' % (SquidGui.defaults['Na_in']))
+ self._temperatureEdit.setText('%g' % (SquidGui.defaults['temperature'] - CELSIUS_TO_KELVIN))
+ self._holdingVEdit.setText('%g' % (SquidGui.defaults['vclamp.holdingV']))
+ self._holdingTimeEdit.setText('%g' % (SquidGui.defaults['vclamp.holdingT']))
+ self._prePulseVEdit.setText('%g' % (SquidGui.defaults['vclamp.prepulseV']))
+ self._prePulseTimeEdit.setText('%g' % (SquidGui.defaults['vclamp.prepulseT']))
+ self._clampVEdit.setText('%g' % (SquidGui.defaults['vclamp.clampV']))
+ self._clampTimeEdit.setText('%g' % (SquidGui.defaults['vclamp.clampT']))
+ self._baseCurrentEdit.setText('%g' % (SquidGui.defaults['iclamp.baseI']))
+ self._firstPulseEdit.setText('%g' % (SquidGui.defaults['iclamp.firstI']))
+ self._firstDelayEdit.setText('%g' % (SquidGui.defaults['iclamp.firstD']))
+ self._firstPulseWidthEdit.setText('%g' % (SquidGui.defaults['iclamp.firstT']))
+ self._secondPulseEdit.setText('%g' % (SquidGui.defaults['iclamp.secondI']))
+ self._secondDelayEdit.setText('%g' % (SquidGui.defaults['iclamp.secondD']))
+ self._secondPulseWidthEdit.setText('%g' % (SquidGui.defaults['iclamp.secondT']))
+ self._pulseMode.setCurrentIndex(0)
+
+ def _onScroll(self, event):
+ if event.inaxes is None:
+ return
+ axes = event.inaxes
+ zoom = 0.0
+ if event.button == 'up':
+ zoom = -1.0
+ elif event.button == 'down':
+ zoom = 1.0
+ if zoom != 0.0:
+ self._plotNavigator.push_current()
+ axes.get_xaxis().zoom(zoom)
+ axes.get_yaxis().zoom(zoom)
+ self._plotCanvas.draw()
+
+ def closeEvent(self, event):
+ qApp.closeAllWindows()
+
+ def _showBioPhysicsHelp(self):
+ self._createHelpMessage()
+ self._helpMessageText.setText('<html><p>%s</p><p>%s</p><p>%s</p><p>%s</p><p>%s</p></html>' %
+ (tooltip_Nernst,
+ tooltip_Erest,
+ tooltip_KChan,
+ tooltip_NaChan,
+ tooltip_Im))
+ self._helpWindow.setVisible(True)
+
+ def _showRunningHelp(self):
+ self._createHelpMessage()
+ self._helpMessageText.setSource(QtCore.QUrl(self._helpBaseURL))
+ self._helpWindow.setVisible(True)
+
+ def _createHelpMessage(self):
+ if hasattr(self, '_helpWindow'):
+ return
+ self._helpWindow = QWidget()
+ self._helpWindow.setWindowFlags(QtCore.Qt.Window)
+ layout = QVBoxLayout()
+ self._helpWindow.setLayout(layout)
+ self._helpMessageArea = QScrollArea()
+ self._helpMessageText = QTextBrowser()
+ self._helpMessageText.setOpenExternalLinks(True)
+ self._helpMessageArea.setWidget(self._helpMessageText)
+ layout.addWidget(self._helpMessageText)
+ self._squidGuiPath = os.path.dirname(os.path.abspath(__file__))
+ self._helpBaseURL = os.path.join(self._squidGuiPath,'help.html')
+ self._helpMessageText.setSource(QtCore.QUrl(self._helpBaseURL))
+ self._helpMessageText.setSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding)
+ self._helpMessageArea.setSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding)
+ self._helpMessageText.setMinimumSize(800, 600)
+ self._closeHelpAction = QAction('Close', self)
+ self._closeHelpAction.triggered.connect(self._helpWindow.close)
+ # Setting the close event so that the ``Help`` button is
+ # unchecked when the help window is closed
+ self._helpWindow.closeEvent = lambda event: self._helpAction.setChecked(False)
+ self._helpTOCAction = QAction('Help running demo', self)
+ self._helpTOCAction.triggered.connect( self._jumpToHelpTOC)
+ # This panel is for putting two buttons using horizontal
+ # layout
+ panel = QFrame()
+ panel.setFrameStyle(QFrame.StyledPanel + QFrame.Raised)
+ layout.addWidget(panel)
+ layout = QHBoxLayout()
+ panel.setLayout(layout)
+ self._helpAction = QAction('Help running', self)
+ self._helpAction.triggered.connect(self._showRunningHelp)
+ self._helpBiophysicsAction = QAction('Help biophysics', self)
+ self._helpBiophysicsAction.triggered.connect(self._showBioPhysicsHelp)
+ self._helpTOCButton = QToolButton()
+ self._helpTOCButton.setDefaultAction(self._helpTOCAction)
+ self._helpBiophysicsButton = QToolButton()
+ self._helpBiophysicsButton.setDefaultAction(self._helpBiophysicsAction)
+ layout.addWidget(self._helpTOCButton)
+ layout.addWidget(self._helpBiophysicsButton)
+ self._closeHelpButton = QToolButton()
+ self._closeHelpButton.setDefaultAction(self._closeHelpAction)
+ layout.addWidget(self._closeHelpButton)
+
+ def _jumpToHelpTOC(self):
+ self._helpMessageText.setSource(QtCore.QUrl(self._helpBaseURL))
+
+if __name__ == '__main__':
+ app = QApplication(sys.argv)
+ # app.connect(app, QtCore.SIGNAL('lastWindowClosed()'), app, QtCore.SLOT('quit()'))
+ app.quitOnLastWindowClosed( )
+ qApp = app
+ squid_gui = SquidGui()
+ squid_gui.show()
+ print(squid_gui.size())
+ sys.exit(app.exec_())
+
+#
+# squidgui.py ends here
diff --git a/moose-examples/squid/test_squid.py b/moose-examples/squid/test_squid.py
index 005cb1626077b6164a00f7eb17ec0c788fa3e670..5bbef8bcec843677903be28e162c44e514633238 100644
--- a/moose-examples/squid/test_squid.py
+++ b/moose-examples/squid/test_squid.py
@@ -30,6 +30,7 @@
import unittest
import pylab
+import numpy
from squid import SquidModel
class SquidAxonTest(unittest.TestCase):
@@ -112,6 +113,5 @@ class SquidAxonTest(unittest.TestCase):
self.assertLessEqual(difference, numpy.mean(beta_m)*1e-6)
-
#
# test_squid.py ends here
diff --git a/moose-examples/tutorials/Electrophys/CableInjectEquivCkt.png b/moose-examples/tutorials/Electrophys/CableInjectEquivCkt.png
new file mode 100644
index 0000000000000000000000000000000000000000..d1938db8fdf25b373a13fce83b704d5edd031040
Binary files /dev/null and b/moose-examples/tutorials/Electrophys/CableInjectEquivCkt.png differ
diff --git a/moose-examples/tutorials/Electrophys/README.txt b/moose-examples/tutorials/Electrophys/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..585fd34eb33aea4c68345cc080ba48190fdd7090
--- /dev/null
+++ b/moose-examples/tutorials/Electrophys/README.txt
@@ -0,0 +1,80 @@
+This is a list of files in the Electrophys tutorials directory
+
+ephys1_cable.py:
+
+This is a model of a simple uniform passive cable. It has two plots: one
+of membrane potential Vm as a function of time, sampled at four locations
+along the cable, and the other of Vm as a function of position, sampled at
+5 times (intervals of 5 ms) during the simulation.
+The user can use the sliders to scale the following parameters up and down:
+RM
+RA
+CM
+length
+diameter
+
+In addition, the user can toggle between a long (steady-state) current
+injection stimulus, and a brief pulse of 1 ms.
+
+Finally, when the user hits "Quit" the time-series of the soma compartment
+is printed out.
+
+Things to do:
+1. Vary diameter, length, RM, CM, RA. Check that the decay behaves as per
+ cable theory. Get an intuition for these.
+2. Check that lambda (length constant) is according to equations. Check
+ how much decay there is for different L (electrotonic length).
+2a. Note that if the cable is shorter than L, the signal doesn't decay so
+ well. Why?
+2b. Convince yourself that you'll never see a passive signal along a long
+ axon.
+3. Examine propagation of the depolarization for a brief current pulse.
+ Look for the three signatures of propagation along a dendrite.
+4. Run a simulation of the long stimulus. When you hit Quit, the program will
+ dump the soma potential charging time-course into a file. Check that
+ the tau (time-constant) of the soma is according to equations when the
+ L is very small, but deviates when it is longer. Use Rall's expression
+ for L as a function of the time-courses of soma depolarization.
+
+
+
+ephys2_Rall_law.py:
+
+This explores the implication of Rall's Law and cable branching.
+This is a model of a branched cell, compared with the model of a uniform
+cylindrical cell. The sliders vary the parameters of the branches. Two
+time-points are displayed; 10 ms and 50 ms.
+
+Things to do:
+1. Vary all the passive parameters of branches, see how they affect the
+ propagation past the branch point.
+2. Match up the branched cell (blue plot/dots) to the cylinder (red line).
+ See if the resultant parameters obey Rall's Law.
+
+Todo: - Stimuli in any end or both dend ends.
+
+Neuronal summation
+ - Synaptic input at Y tips and branch point, both I and E.
+ - Vary weights
+ - Vary time since start for each.
+ - Vary taus of I.
+ - Have spiking soma, to set up thresholding as an option.
+
+Channel mixer:
+ - Modulate the conductance of battery of channels
+ - Modulate Ca_conc tau
+ - Plot Ca and Vm
+ - Set different stim pulse amplitude and duration.
+
+NMDA receptor and associativity:
+ - Plot Vm and Ca
+ - Give glu and NMDA input different time and different ampl
+
+AP propagation:
+ - Vary dia
+ - Vary RM, CM, RA
+ - Vary channel densities
+
+Squid demo:
+ - Already have it.
+
diff --git a/moose-examples/tutorials/Electrophys/RallsLaw.png b/moose-examples/tutorials/Electrophys/RallsLaw.png
new file mode 100644
index 0000000000000000000000000000000000000000..21fae9010f07863d14dd3a3cf25720f8dfe6846d
Binary files /dev/null and b/moose-examples/tutorials/Electrophys/RallsLaw.png differ
diff --git a/moose-examples/tutorials/Electrophys/ephys1_cable.py b/moose-examples/tutorials/Electrophys/ephys1_cable.py
new file mode 100644
index 0000000000000000000000000000000000000000..515ea4fda7ca2b9fbbbd33e4023d264eb902a8b2
--- /dev/null
+++ b/moose-examples/tutorials/Electrophys/ephys1_cable.py
@@ -0,0 +1,212 @@
+########################################################################
+# This example demonstrates a cable
+# Copyright (C) Upinder S. Bhalla NCBS 2018
+# Released under the terms of the GNU Public License V3.
+########################################################################
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+from matplotlib.widgets import Slider, Button
+import numpy as np
+import warnings
+import moose
+import rdesigneur as rd
+
+numDendSeg = 50
+interval = 0.005
+lines = []
+tplot = []
+RM = 1.0
+RA = 1.0
+CM = 0.01
+length = 0.002
+dia = 1e-6
+stimStr = "2e-10"
+runtime = 50e-3
+elecPlotDt = 0.001
+
+def makeModel():
+ segLen = length / numDendSeg
+ rdes = rd.rdesigneur(
+ stealCellFromLibrary = True,
+ elecPlotDt = elecPlotDt,
+ verbose = False,
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]
+ # The numerical arguments are all optional
+ cellProto =
+ [['ballAndStick', 'soma', dia, segLen, dia, length, numDendSeg]],
+ passiveDistrib = [[ '#', 'RM', str(RM), 'CM', str(CM), 'RA', str(RA) ]],
+ stimList = [['soma', '1', '.', 'inject', stimStr ]],
+ plotList = [['dend3,dend18,dend33,dend49', '1', '.', 'Vm' ]],
+ )
+ rdes.buildModel()
+
+def main():
+ global vec
+ warnings.filterwarnings("ignore", category=UserWarning, module="matplotlib")
+ makeDisplay()
+ quit()
+
+def updateRM(RM_):
+ global RM
+ RM = RM_
+ updateDisplay()
+
+def updateCM(CM_):
+ global CM
+ CM = CM_
+ updateDisplay()
+
+def updateRA(RA_):
+ global RA
+ RA = RA_
+ updateDisplay()
+
+def updateLen(val): # This does the length of the entire cell.
+ global length
+ length = val * 0.001
+ updateDisplay()
+
+def updateDia(val):
+ global dia
+ dia = val * 1e-6
+ updateDisplay()
+
+class stimToggle():
+ def __init__( self, toggle, ax ):
+ self.duration = 1
+ self.toggle = toggle
+ self.ax = ax
+
+ def click( self, event ):
+ global stimStr
+ if self.duration < 0.5:
+ self.duration = 1.0
+ self.toggle.label.set_text( "Long Stim" )
+ self.toggle.color = "yellow"
+ self.toggle.hovercolor = "yellow"
+ stimStr = "2e-10"
+ else:
+ self.duration = 0.001
+ self.toggle.label.set_text( "Short Stim" )
+ self.toggle.color = "orange"
+ self.toggle.hovercolor = "orange"
+ stimStr = "1e-9*(t<0.001)"
+ #stimStr = "10e-9*(t<0.001)-9e-9*(t>0.001&&t<0.002)"
+ updateDisplay()
+
+def updateDisplay():
+ makeModel()
+ vec = moose.wildcardFind( '/model/elec/#[ISA=CompartmentBase]' )
+ tabvec = moose.wildcardFind( '/model/graphs/plot#' )
+ #vec[0].inject = 1e-10
+ moose.reinit()
+ initdt = 0.0001
+ dt = initdt
+ currtime = 0.0
+ for i in lines:
+ moose.start( dt )
+ currtime += dt
+ i.set_ydata( [v.Vm * 1000 for v in vec] )
+ dt = interval
+ #print( "inRunSim v0={}, v10={}".format( vec[0].Vm, vec[10].Vm ) )
+ moose.start( runtime - currtime )
+ for i,tab in zip( tplot, tabvec ):
+ i.set_ydata( tab.vector * 1000 )
+
+ moose.delete( '/model' )
+ moose.delete( '/library' )
+ # Put in something here for the time-series on soma
+
+def doQuit( event ):
+ tab = []
+ makeModel()
+ soma = moose.element( '/model/elec/soma' )
+ moose.reinit()
+ with open('output.txt', 'w') as file:
+ file.write( "0.000 {:.2f}\n".format( soma.Vm * 1000 ) )
+ for t in np.arange( 0, 0.1, 0.001 ):
+ moose.start(0.001)
+ file.write( "{:.3f} {:.2f}\n".format( t+0.001, soma.Vm*1000 ) )
+ quit()
+
+def makeDisplay():
+ global tplot
+ global lines
+ #img = mpimg.imread( 'CableEquivCkt.png' )
+ img = mpimg.imread( 'CableInjectEquivCkt.png' )
+ #plt.ion()
+ fig = plt.figure( figsize=(10,12) )
+ png = fig.add_subplot(321)
+ imgplot = plt.imshow( img )
+ plt.axis('off')
+ ax1 = fig.add_subplot(322)
+ plt.ylabel( 'Vm (mV)' )
+ plt.ylim( -80, 40 )
+ plt.xlabel( 'time (ms)' )
+ plt.title( "Membrane potential vs time at 4 positions." )
+ t = np.arange( 0.0, runtime + elecPlotDt / 2.0, elecPlotDt ) * 1000 #ms
+ for i,col,name in zip( range( 4 ), ['b-', 'g-', 'r-', 'm-' ], ['a', 'b', 'c', 'd'] ):
+ ln, = ax1.plot( t, np.zeros(len(t)), col, label='pos= ' + name )
+ tplot.append(ln)
+ plt.legend()
+
+ ax2 = fig.add_subplot(312)
+ #ax2.margins( 0.05 )
+ #ax.set_ylim( 0, 0.1 )
+ plt.ylabel( 'Vm (mV)' )
+ plt.ylim( -80, 50 )
+ plt.xlabel( 'Position (microns)' )
+ #ax2.autoscale( enable = True, axis = 'y' )
+ plt.title( "Membrane potential vs position, at 5 times." )
+ t = np.arange( 0, numDendSeg+1, 1 ) #sec
+ for i,col in zip( range( 5 ), ['k-', 'b-', 'g-', 'y-', 'm-' ] ):
+ ln, = ax2.plot( t, np.zeros(numDendSeg+1), col, label="t={}".format(i*interval) )
+ lines.append(ln)
+ plt.legend()
+ ax = fig.add_subplot(313)
+ #ax.margins( 0.05 )
+ plt.axis('off')
+ axcolor = 'palegreen'
+ axStim = plt.axes( [0.02,0.05, 0.20,0.03], facecolor='green' )
+ axReset = plt.axes( [0.25,0.05, 0.30,0.03], facecolor='blue' )
+ axQuit = plt.axes( [0.60,0.05, 0.30,0.03], facecolor='blue' )
+ axRM = plt.axes( [0.25,0.1, 0.65,0.03], facecolor=axcolor )
+ axCM = plt.axes( [0.25,0.15, 0.65,0.03], facecolor=axcolor )
+ axRA = plt.axes( [0.25,0.20, 0.65,0.03], facecolor=axcolor )
+ axLen = plt.axes( [0.25,0.25, 0.65,0.03], facecolor=axcolor )
+ axDia = plt.axes( [0.25,0.30, 0.65,0.03], facecolor=axcolor )
+ #aInit = Slider( axAinit, 'A init conc', 0, 10, valinit=1.0, valstep=0.2)
+ stim = Button( axStim, 'Long Stim', color = 'yellow' )
+ stimObj = stimToggle( stim, axStim )
+
+ reset = Button( axReset, 'Reset', color = 'cyan' )
+ q = Button( axQuit, 'Quit', color = 'pink' )
+ RM = Slider( axRM, 'RM ( ohm.m^2 )', 0.1, 10, valinit=1.0 )
+ CM = Slider( axCM, 'CM ( Farad/m^2)', 0.001, 0.1, valinit=0.01 )
+ RA = Slider( axRA, 'RA ( ohm.m', 0.1, 10, valinit=1.0 )
+ length = Slider( axLen, 'Total length of cell (mm)', 0.1, 10, valinit=2.0 )
+ dia = Slider( axDia, 'Diameter of cell (um)', 0.1, 10, valinit=1.0 )
+ def resetParms( event ):
+ RM.reset()
+ CM.reset()
+ RA.reset()
+ length.reset()
+ dia.reset()
+
+
+ stim.on_clicked( stimObj.click )
+ reset.on_clicked( resetParms )
+ q.on_clicked( doQuit )
+ RM.on_changed( updateRM )
+ CM.on_changed( updateCM )
+ RA.on_changed( updateRA )
+ length.on_changed( updateLen )
+ dia.on_changed( updateDia )
+ plt.tight_layout()
+
+ updateDisplay()
+ plt.show()
+
+# Run the 'main' if this script is executed standalone.
+if __name__ == '__main__':
+ main()
diff --git a/moose-examples/tutorials/Electrophys/ephys2_Rall_law.py b/moose-examples/tutorials/Electrophys/ephys2_Rall_law.py
new file mode 100644
index 0000000000000000000000000000000000000000..a5f6e90d25f5716473cc474319097a7fc80150d9
--- /dev/null
+++ b/moose-examples/tutorials/Electrophys/ephys2_Rall_law.py
@@ -0,0 +1,268 @@
+########################################################################
+# This example demonstrates a cable
+# Copyright (C) Upinder S. Bhalla NCBS 2018
+# Released under the terms of the GNU Public License V3.
+########################################################################
+import moose
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+from matplotlib.widgets import Slider, Button
+import numpy as np
+import rdesigneur as rd
+
+numDendSeg = 10 # Applies both to dend and to branches.
+interval1 = 0.010
+interval2 = 0.05 - interval1
+lines = []
+# These 5 params vary only for the branches.
+RM = 1.0
+RA = 1.0
+CM = 0.01
+length = 0.001
+dia = 2e-6
+# The params below are fixed, apply to the soma and dend.
+dendRM = 2.0
+dendRA = 1.5
+dendCM = 0.02
+dendLength = 0.001
+dendDia = 2e-6
+
+# Stimulus in Amps applied to soma.
+stimStr = "2e-10"
+
+class lineWrapper():
+ def __init__( self ):
+ self.YdendLines = 0
+ self.Ybranch1Lines = 0
+ self.Ybranch2Lines = 0
+ self.CylLines = 0
+
+def makeYmodel():
+ segLen = dendLength / numDendSeg
+ rdes = rd.rdesigneur(
+ stealCellFromLibrary = True,
+ verbose = False,
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]
+ # The numerical arguments are all optional
+ cellProto =
+ [['ballAndStick', 'cellBase', dendDia, segLen, dendDia, dendLength, numDendSeg]],
+ passiveDistrib = [[ '#', 'RM', str(dendRM), 'CM', str(dendCM), 'RA', str(dendRA) ]],
+ stimList = [['soma', '1', '.', 'inject', stimStr ]],
+ )
+ # Build the arms of the Y for a branching cell.
+ pa = moose.element( '/library/cellBase' )
+ x = length
+ y = 0.0
+ dx = length / ( numDendSeg * np.sqrt(2.0) )
+ dy = 0.0
+ prevc1 = moose.element( '/library/cellBase/dend{}'.format( numDendSeg-1 ) )
+ prevc2 = prevc1
+ for i in range( numDendSeg ):
+ c1 = rd.buildCompt( pa, 'branch1_{}'.format(i), RM = RM, CM = CM, RA = RA, dia = dia, x=x, y=y, dx = dx, dy = dy )
+ c2 = rd.buildCompt( pa, 'branch2_{}'.format(i), RM = RM, CM = CM, RA = RA, dia = dia, x=x, y=-y, dx = dx, dy = -dy )
+ moose.connect( prevc1, 'axial', c1, 'raxial' )
+ moose.connect( prevc2, 'axial', c2, 'raxial' )
+ prevc1 = c1
+ prevc2 = c2
+ x += dx
+ y += dy
+ rdes.buildModel()
+
+def makeCylModel():
+ segLen = dendLength / numDendSeg
+ rdes = rd.rdesigneur(
+ stealCellFromLibrary = True,
+ verbose = False,
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]
+ # The numerical arguments are all optional
+ cellProto =
+ [['ballAndStick', 'soma', dendDia, segLen, dendDia, 2*dendLength, 2*numDendSeg]],
+ passiveDistrib = [[ '#', 'RM', str(dendRM), 'CM', str(dendCM), 'RA', str(dendRA) ]],
+ stimList = [['soma', '1', '.', 'inject', stimStr ]],
+ )
+ rdes.buildModel()
+
+def main():
+ global vec
+ makeDisplay()
+ quit()
+
+def updateRM(RM_):
+ global RM
+ RM = RM_
+ updateDisplay()
+
+def updateCM(CM_):
+ global CM
+ CM = CM_
+ updateDisplay()
+
+def updateRA(RA_):
+ global RA
+ RA = RA_
+ updateDisplay()
+
+def updateLen(val): # This does the length of the entire cell.
+ global length
+ length = val * 0.001
+ updateDisplay()
+
+def updateDia(val):
+ global dia
+ dia = val * 1e-6
+ updateDisplay()
+
+class stimToggle():
+ def __init__( self, toggle, ax ):
+ self.duration = 1
+ self.toggle = toggle
+ self.ax = ax
+
+ def click( self, event ):
+ global stimStr
+ if self.duration < 0.5:
+ self.duration = 1.0
+ self.toggle.label.set_text( "Long Stim" )
+ self.toggle.color = "yellow"
+ self.toggle.hovercolor = "yellow"
+ stimStr = "2e-10"
+ else:
+ self.duration = 0.001
+ self.toggle.label.set_text( "Short Stim" )
+ self.toggle.color = "orange"
+ self.toggle.hovercolor = "orange"
+ stimStr = "40e-9*(t<0.001)-36e-9*(t>0.001&&t<0.002)"
+ updateDisplay()
+ #self.ax.update()
+ #self.ax.redraw_in_frame()
+ #self.ax.draw()
+
+def printSomaVm():
+ print("This is somaVm" )
+
+def updateDisplay():
+ makeCylModel()
+ moose.element( '/model/elec/' ).name = 'Cyl'
+ moose.element( '/model' ).name = 'model1'
+ makeYmodel()
+ moose.element( '/model/elec/' ).name = 'Y'
+ moose.move( '/model1/Cyl', '/model' )
+ #moose.le( '/model/Y' )
+ #print "################################################"
+ #moose.le( '/model/Cyl' )
+ vecYdend = moose.wildcardFind( '/model/Y/soma,/model/Y/dend#' )
+ vecYbranch1 = moose.wildcardFind( '/model/Y/branch1#' )
+ vecYbranch2 = moose.wildcardFind( '/model/Y/branch2#' )
+ vecCyl = moose.wildcardFind( '/model/Cyl/#[ISA=CompartmentBase]' )
+ #vec[0].inject = 1e-10
+ moose.reinit()
+ dt = interval1
+ for i in lines:
+ moose.start( dt )
+ #print( len(vecCyl), len(vecYdend), len(vecYbranch1), len(vecYbranch2) )
+ i.CylLines.set_ydata( [v.Vm*1000 for v in vecCyl] )
+ i.YdendLines.set_ydata( [v.Vm*1000 for v in vecYdend] )
+ i.Ybranch1Lines.set_ydata( [v.Vm*1000 for v in vecYbranch1] )
+ i.Ybranch2Lines.set_ydata( [v.Vm*1000 for v in vecYbranch2] )
+ dt = interval2
+
+ moose.delete( '/model' )
+ moose.delete( '/model1' )
+ moose.delete( '/library' )
+ # Put in something here for the time-series on soma
+
+def doQuit( event ):
+ cylLength = 2*dendLength*float(2*numDendSeg+1)/(2*numDendSeg)
+ print( "The cylinder parameters were:\n"
+ "Length = {} microns\n"
+ "Diameter = {} microns\n"
+ "RM = {} Ohms.m^2\n"
+ "RA = {} Ohms.m\n"
+ "CM = {} Farads/m^2\n"
+ "Your branch parameters were:\n"
+ "Length = {:.2f} microns\n"
+ "Diameter = {:.2f} microns\n"
+ "RM = {:.2f} Ohms.m^2\n"
+ "RA = {:.2f} Ohms.m\n"
+ "CM = {:.3f} Farads/m^2\n".format(
+ cylLength*1e6, dendDia*1e6, dendRM, dendRA, dendCM,
+ length*1e6, dia*1e6, RM, RA, CM ) )
+ print( "The branch point was at 1100 microns from the left" )
+
+ quit()
+
+def makeDisplay():
+ global lines
+ #img = mpimg.imread( 'CableEquivCkt.png' )
+ img = mpimg.imread( 'RallsLaw.png' )
+ #plt.ion()
+ fig = plt.figure( figsize=(10,12) )
+ png = fig.add_subplot(311)
+ imgplot = plt.imshow( img )
+ plt.axis('off')
+ ax2 = fig.add_subplot(312)
+ #ax.set_ylim( 0, 0.1 )
+ plt.ylabel( 'Vm (mV)' )
+ plt.ylim( -70, 0.0 )
+ plt.xlabel( 'Position (microns)' )
+ #ax2.autoscale( enable = True, axis = 'y' )
+ plt.title( "Membrane potential as a function of position along cell." )
+ #for i,col in zip( range( 5 ), ['k', 'b', 'g', 'y', 'm' ] ):
+ for i,col in zip( range( 2 ), ['b', 'k' ] ):
+ lw = lineWrapper()
+ lw.YdendLines, = ax2.plot( np.arange(0, numDendSeg+1, 1 ),
+ np.zeros(numDendSeg+1), col + '.' )
+ lw.Ybranch1Lines, = ax2.plot( np.arange(0, numDendSeg, 1) + numDendSeg + 1,
+ np.zeros(numDendSeg), col + ':' )
+ lw.Ybranch2Lines, = ax2.plot( np.arange(0, numDendSeg, 1) + numDendSeg + 1,
+ np.zeros(numDendSeg) + numDendSeg + 1, col + '.' )
+ lw.CylLines, = ax2.plot( np.arange(0, numDendSeg*2+1, 1),
+ np.zeros(numDendSeg*2+1), 'r-' )
+ lines.append( lw )
+
+ ax = fig.add_subplot(313)
+ plt.axis('off')
+ axcolor = 'palegreen'
+ axStim = plt.axes( [0.02,0.05, 0.20,0.03], facecolor='green' )
+ axReset = plt.axes( [0.25,0.05, 0.30,0.03], facecolor='blue' )
+ axQuit = plt.axes( [0.60,0.05, 0.30,0.03], facecolor='blue' )
+ axRM = plt.axes( [0.25,0.1, 0.65,0.03], facecolor=axcolor )
+ axCM = plt.axes( [0.25,0.15, 0.65,0.03], facecolor=axcolor )
+ axRA = plt.axes( [0.25,0.20, 0.65,0.03], facecolor=axcolor )
+ axLen = plt.axes( [0.25,0.25, 0.65,0.03], facecolor=axcolor )
+ axDia = plt.axes( [0.25,0.30, 0.65,0.03], facecolor=axcolor )
+ #aInit = Slider( axAinit, 'A init conc', 0, 10, valinit=1.0, valstep=0.2)
+ stim = Button( axStim, 'Long Stim', color = 'yellow' )
+ stimObj = stimToggle( stim, axStim )
+
+ reset = Button( axReset, 'Reset', color = 'cyan' )
+ q = Button( axQuit, 'Quit', color = 'pink' )
+ RM = Slider( axRM, 'RM ( ohm.m^2 )', 0.1, 10, valinit=1.0 )
+ CM = Slider( axCM, 'CM ( Farad/m^2)', 0.001, 0.05, valinit=0.01, valfmt = '%0.3f' )
+ RA = Slider( axRA, 'RA ( ohm.m', 0.1, 10, valinit=1.0 )
+ length = Slider( axLen, 'Length of branches (mm)', 0.1, 10, valinit=2.0 )
+ dia = Slider( axDia, 'Diameter of branches (um)', 0.1, 10, valinit=1.0 )
+ def resetParms( event ):
+ RM.reset()
+ CM.reset()
+ RA.reset()
+ length.reset()
+ dia.reset()
+
+
+ stim.on_clicked( stimObj.click )
+ reset.on_clicked( resetParms )
+ q.on_clicked( doQuit )
+ RM.on_changed( updateRM )
+ CM.on_changed( updateCM )
+ RA.on_changed( updateRA )
+ length.on_changed( updateLen )
+ dia.on_changed( updateDia )
+
+ updateDisplay()
+
+ plt.show()
+
+# Run the 'main' if this script is executed standalone.
+if __name__ == '__main__':
+ main()
diff --git a/moose-examples/tutorials/Rdesigneur/README.txt b/moose-examples/tutorials/Rdesigneur/README.txt
index b64d89a987600062e60963a0484823aa31a8f3a6..e76845689429a3080d8b04131c3ad7798ddf29b3 100644
--- a/moose-examples/tutorials/Rdesigneur/README.txt
+++ b/moose-examples/tutorials/Rdesigneur/README.txt
@@ -291,6 +291,16 @@ ex9.2_spines_in_neuronal_morpho.py: Add spines to a neuron built from a
- See if you can deliver the current injection to the spine. Hint: the
name of the spine compartments is 'head#' where # is the index of the
spine.
+
+
+ex9.3_spiral_spines.py: Just for fun. Illustrates how to place spines in a
+spiral around the dendrite. For good measure the spines get bigger the further
+they are from the soma. Note that the uniform spacing of spines is signified
+by the negative minSpacing term, the fourth argument to spineDistrib.
+ Suggestions:
+ - Play with expressions for spine size and angular placement.
+ - See what happens if the segment size gets smaller than the
+ spine spacing.
To come:
rdes_ex10.py: Build a spiny neuron, and insert the oscillatory chemical model
diff --git a/moose-examples/tutorials/Rdesigneur/ex3.2_squid_axon_propgn.py b/moose-examples/tutorials/Rdesigneur/ex3.2_squid_axon_propgn.py
index 9729a0d653942428734e7966205fe06f789ff876..cb3e0a5587b67c1c9abd78eda16ab478ee41e1d6 100644
--- a/moose-examples/tutorials/Rdesigneur/ex3.2_squid_axon_propgn.py
+++ b/moose-examples/tutorials/Rdesigneur/ex3.2_squid_axon_propgn.py
@@ -43,7 +43,7 @@ rdes = rd.rdesigneur(
chanDistrib = [
['Na', '#', 'Gbar', '1200' ],
['K', '#', 'Gbar', '360' ]],
- stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.2) * 2e-11' ]],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.005 && t<0.2) * 2e-11' ]],
plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']],
moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
)
@@ -51,4 +51,4 @@ rdes = rd.rdesigneur(
rdes.buildModel()
moose.reinit()
-rdes.displayMoogli( 0.00005, 0.05, 0.0 )
+rdes.displayMoogli( 0.00005, 0.04, 0.0 )
diff --git a/moose-examples/tutorials/Rdesigneur/ex3.3_AP_collision.py b/moose-examples/tutorials/Rdesigneur/ex3.3_AP_collision.py
index b6a8ce2a80a0ad0435703eb4159f646f43e36ae8..dfe4e8cda578cf12863d7d04a2419c06c7e04b8b 100644
--- a/moose-examples/tutorials/Rdesigneur/ex3.3_AP_collision.py
+++ b/moose-examples/tutorials/Rdesigneur/ex3.3_AP_collision.py
@@ -52,4 +52,4 @@ rdes = rd.rdesigneur(
rdes.buildModel()
moose.reinit()
-rdes.displayMoogli( 0.00005, 0.05, 0.0 )
+rdes.displayMoogli( 0.00005, 0.03, 0.0 )
diff --git a/moose-examples/tutorials/Rdesigneur/ex3.4_myelinated_axon.py b/moose-examples/tutorials/Rdesigneur/ex3.4_myelinated_axon.py
index e91f7898249c78688016dcb98c0e547b7101b430..c609529d0eed4ced43c5ca76345860818e1df40a 100644
--- a/moose-examples/tutorials/Rdesigneur/ex3.4_myelinated_axon.py
+++ b/moose-examples/tutorials/Rdesigneur/ex3.4_myelinated_axon.py
@@ -58,12 +58,6 @@ rdes = rd.rdesigneur(
moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
)
-
rdes.buildModel()
-
-for i in moose.wildcardFind( "/model/elec/#/Na" ):
- print(i.parent.name, i.Gbar)
-
moose.reinit()
-
rdes.displayMoogli( 0.00005, 0.05, 0.0 )
diff --git a/moose-examples/tutorials/Rdesigneur/ex7.2_CICR.py b/moose-examples/tutorials/Rdesigneur/ex7.2_CICR.py
index a44e4a1a7833327b788fd7e98dc730eb89357fdb..38447a8d0a2c040874cc6c8f1b91e88481a1c1fb 100644
--- a/moose-examples/tutorials/Rdesigneur/ex7.2_CICR.py
+++ b/moose-examples/tutorials/Rdesigneur/ex7.2_CICR.py
@@ -15,6 +15,7 @@ rdes = rd.rdesigneur(
turnOffElec = True,
chemDt = 0.005,
chemPlotDt = 0.02,
+ numWaveFrames = 200,
diffusionLength = 1e-6,
useGssa = False,
addSomaChemCompt = False,
diff --git a/moose-examples/tutorials/Rdesigneur/ex9.3_spiral_spines.py b/moose-examples/tutorials/Rdesigneur/ex9.3_spiral_spines.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe5c39c544f28a24c70025c2a8ad63bb2959edce
--- /dev/null
+++ b/moose-examples/tutorials/Rdesigneur/ex9.3_spiral_spines.py
@@ -0,0 +1,19 @@
+##########################################################################
+# This illustrates some of the capabilities for spine placement.
+# It has spines whose size increase with distance from the soma.
+# Further, the angular direction of the spines spirals around the dendrite.
+##########################################################################
+import moose
+import rdesigneur as rd
+rdes = rd.rdesigneur(
+ cellProto = [['ballAndStick', 'elec', 10e-6, 10e-6, 2e-6, 300e-6, 50]],
+ spineProto = [['makePassiveSpine()', 'spine']],
+ spineDistrib = [['spine', '#dend#', '3e-6', '-1e-6', '1+p*2e4', '0', 'p*6.28e7', '0']],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.02) * 1e-9' ]],
+ moogList = [['#', '1', '.', 'Vm', 'Soma potential']]
+)
+
+rdes.buildModel()
+
+moose.reinit()
+rdes.displayMoogli( 0.0002, 0.025, 0.02 )
diff --git a/moose-gui/objectedit.py b/moose-gui/objectedit.py
index cbc8c22bffb77b358b864c094e490aea7372c9d6..d5d66b60c96fa506af50e620260a109f5deb8e29 100644
--- a/moose-gui/objectedit.py
+++ b/moose-gui/objectedit.py
@@ -175,10 +175,14 @@ class ObjectEditModel(QtCore.QAbstractTableModel):
self.fields.append(fieldName)
#harsha: For signalling models will be pulling out notes field from Annotator
# can updates if exist for other types also
- if ( isinstance(self.mooseObject, moose.PoolBase)
- or isinstance(self.mooseObject,moose.EnzBase)
- or isinstance(self.mooseObject,moose.Neutral)) :
- self.fields.append("Color")
+ if (isinstance (self.mooseObject,moose.ChemCompt) or \
+ isinstance(self.mooseObject,moose.ReacBase) or \
+ isinstance(moose.element(moose.element(self.mooseObject).parent),moose.EnzBase) \
+ ):
+ pass
+ else:
+ self.fields.append("Color")
+
flag = QtCore.Qt.ItemIsEnabled | QtCore.Qt.ItemIsSelectable | QtCore.Qt.ItemIsEditable
self.fieldFlags[fieldName] = flag
diff --git a/moose-gui/plugins/kkit.py b/moose-gui/plugins/kkit.py
index 4870a71d167ec86b8370bf493a81a717a17493d2..3cd531bb4498af9fc86fc75bc11a9a4376f602d0 100644
--- a/moose-gui/plugins/kkit.py
+++ b/moose-gui/plugins/kkit.py
@@ -6,13 +6,14 @@ __version__ = "1.0.0"
__maintainer__ = "HarshaRani"
__email__ = "hrani@ncbs.res.in"
__status__ = "Development"
-__updated__ = "Sep 7 2018"
+__updated__ = "Sep 11 2018"
#Change log:
# 2018
-#Jun 18: update the color of the group from objecteditor
+#sep 11: comparment size is calculated based on group sceneBoundingRect size
#Sep 07: in positionChange all the group's boundingRect is calculated
# and when group is moved the children's position are stored
+#Jun 18: update the color of the group from objecteditor
#code
import math
@@ -908,7 +909,16 @@ class KineticsWidget(EditorWidgetBase):
# rectcompt = calculateChildBoundingRect(grpcompt)
rectgrp = calculateChildBoundingRect(v)
v.setRect(rectgrp.x()-10,rectgrp.y()-10,(rectgrp.width()+20),(rectgrp.height()+20))
-
+ for k, v in self.qGraCompt.items():
+ #rectcompt = v.childrenBoundingRect()
+ rectcompt = calculateChildBoundingRect(v)
+ comptBoundingRect = v.boundingRect()
+ if not comptBoundingRect.contains(rectcompt):
+ self.updateCompartmentSize(v)
+
+ else:
+ rectcompt = calculateChildBoundingRect(v)
+ v.setRect(rectcompt.x()-10,rectcompt.y()-10,(rectcompt.width()+20),(rectcompt.height()+20))
else:
mobj = self.mooseId_GObj[element(mooseObject)]
self.updateArrow(mobj)
@@ -938,11 +948,11 @@ class KineticsWidget(EditorWidgetBase):
comptBoundingRect = v.boundingRect()
if not comptBoundingRect.contains(rectcompt):
self.updateCompartmentSize(v)
+
else:
rectcompt = calculateChildBoundingRect(v)
v.setRect(rectcompt.x()-10,rectcompt.y()-10,(rectcompt.width()+20),(rectcompt.height()+20))
- pass
-
+
def updateGrpSize(self,grp):
compartmentBoundary = grp.rect()
diff --git a/moose-gui/plugins/kkitUtil.py b/moose-gui/plugins/kkitUtil.py
index 7185988005fc7c1c15e2795e878870bd1fd9b73f..b7dde1b7cc061f53c9075db806ac4fce9ba32117 100644
--- a/moose-gui/plugins/kkitUtil.py
+++ b/moose-gui/plugins/kkitUtil.py
@@ -5,12 +5,16 @@ __version__ = "1.0.0"
__maintainer__ = "HarshaRani"
__email__ = "hrani@ncbs.res.in"
__status__ = "Development"
-__updated__ = "Oct 18 2017"
+__updated__ = "Sep 17 2018"
'''
+2018
+Sep 17: when vertical or horizontal layout is applied for group, compartment size is recalculated
+Sep 11: group size is calculated based on sceneBoundingRect for compartment size
+2017
Oct 18 some of the function moved to this file from kkitOrdinateUtils
'''
-from moose import Annotator,element
+from moose import Annotator,element,ChemCompt
from kkitQGraphics import PoolItem, ReacItem,EnzItem,CplxItem,GRPItem,ComptItem
from PyQt4 import QtCore,QtGui,QtSvg
from PyQt4.QtGui import QColor
@@ -127,9 +131,10 @@ def moveX(reference, collider, layoutPt, margin):
layoutPt.drawLine_arrow(itemignoreZooming=False)
def handleCollisions(compartments, moveCallback, layoutPt,margin = 5.0):
-
+ print " handelCollision"
if len(compartments) is 0 : return
compartments = sorted(compartments, key = lambda c: c.sceneBoundingRect().center().x())
+ print " compartment ",compartments
reference = compartments.pop(0);
print (reference.name)
referenceRect = reference.sceneBoundingRect()
@@ -138,6 +143,18 @@ def handleCollisions(compartments, moveCallback, layoutPt,margin = 5.0):
)
for collider in colliders:
moveCallback(reference, collider, layoutPt,margin)
+ #print (reference.mobj).parent
+ if isinstance(element(((reference.mobj).parent)),ChemCompt):
+ v = layoutPt.qGraCompt[element(((reference.mobj).parent))]
+ #layoutPt.updateCompartmentSize(x)
+ rectcompt = calculateChildBoundingRect(v)
+ comptBoundingRect = v.boundingRect()
+ if not comptBoundingRect.contains(rectcompt):
+ layoutPt.updateCompartmentSize(v)
+
+ else:
+ rectcompt = calculateChildBoundingRect(v)
+ v.setRect(rectcompt.x()-10,rectcompt.y()-10,(rectcompt.width()+20),(rectcompt.height()+20))
return handleCollisions(compartments, moveCallback, layoutPt,margin)
def calculateChildBoundingRect(compt):
@@ -149,7 +166,6 @@ def calculateChildBoundingRect(compt):
ypos = []
xpos = []
for l in compt.childItems():
-
''' All the children including pool,reac,enz,polygon(arrow),table '''
if not isinstance(l,QtSvg.QGraphicsSvgItem):
if (not isinstance(l,QtGui.QGraphicsPolygonItem)):
@@ -158,11 +174,18 @@ def calculateChildBoundingRect(compt):
xpos.append(l.pos().x())
ypos.append(l.pos().y()+l.boundingRect().bottomRight().y())
ypos.append(l.pos().y())
+
else:
+ xpos.append(l.sceneBoundingRect().x())
+ xpos.append(l.sceneBoundingRect().bottomRight().x())
+ ypos.append(l.sceneBoundingRect().y())
+ ypos.append(l.sceneBoundingRect().bottomRight().y())
+ '''
xpos.append(l.rect().x())
xpos.append(l.boundingRect().bottomRight().x())
ypos.append(l.rect().y())
ypos.append(l.boundingRect().bottomRight().y())
+ '''
if (isinstance(l,PoolItem) or isinstance(l,EnzItem)):
''' For Enz cplx height and for pool function height needs to be taken'''
for ll in l.childItems():
diff --git a/moose-gui/plugins/kkitViewcontrol.py b/moose-gui/plugins/kkitViewcontrol.py
index a813bffe41e544d21228295f1b2328f260d18af0..06fa8b7f9cca58f97379d0af03360736abe6e311 100644
--- a/moose-gui/plugins/kkitViewcontrol.py
+++ b/moose-gui/plugins/kkitViewcontrol.py
@@ -144,13 +144,14 @@ class GraphicalView(QtGui.QGraphicsView):
elif gsolution[1] == GROUP_INTERIOR:
groupInteriorfound = True
groupList.append(gsolution)
-
if item.name == COMPARTMENT:
csolution = (item, self.resolveCompartmentInteriorAndBoundary(item, position))
if csolution[1] == COMPARTMENT_BOUNDARY:
- comptInteriorfound = True
- comptBoundary.append(csolution)
-
+ return csolution
+ # elif csolution[1] == COMPARTMENT_INTERIOR:
+ # comptInteriorfound = True
+ # comptBoundary.append(csolution)
+
if groupInteriorfound:
if comptInteriorfound:
return comptBoundary[0]
@@ -181,7 +182,6 @@ class GraphicalView(QtGui.QGraphicsView):
##This is kept for reference, so that if object (P,R,E,Tab,Fun) is moved outside the compartment,
#then it need to be pull back to original position
self.state["press"]["scenepos"] = item.parent().scenePos()
-
if itemType == COMPARTMENT_INTERIOR or itemType == GROUP_BOUNDARY or itemType == GROUP_INTERIOR:
self.removeConnector()
@@ -194,6 +194,9 @@ class GraphicalView(QtGui.QGraphicsView):
if itemType == GROUP_BOUNDARY:
popupmenu = QtGui.QMenu('PopupMenu', self)
popupmenu.addAction("DeleteGroup", lambda : self.deleteGroup(item,self.layoutPt))
+ popupmenu.addAction("LinearLayout", lambda : handleCollisions(list(self.layoutPt.qGraGrp.values()), moveX, self.layoutPt))
+ popupmenu.addAction("VerticalLayout" ,lambda : handleCollisions(list(self.layoutPt.qGraGrp.values()), moveMin, self.layoutPt ))
+
#popupmenu.addAction("CloneGroup" ,lambda : handleCollisions(comptList, moveMin, self.layoutPt ))
popupmenu.exec_(self.mapToGlobal(event.pos()))
@@ -361,6 +364,12 @@ class GraphicalView(QtGui.QGraphicsView):
grpCmpt = self.findGraphic_groupcompt(item)
if movedGraphObj.parentItem() != grpCmpt:
'''Not same compartment/group to which it belonged to '''
+ if isinstance(movedGraphObj,FuncItem):
+ funcPool = moose.element((movedGraphObj.mobj.neighbors['valueOut'])[0])
+ parentGrapItem = self.layoutPt.mooseId_GObj[moose.element(funcPool)]
+ if parentGrapItem.parentItem() != grpCmpt:
+ self.objectpullback("Functionparent",grpCmpt,movedGraphObj,xx,yy)
+
if isinstance(movedGraphObj,(EnzItem,MMEnzItem)):
parentPool = moose.element((movedGraphObj.mobj.neighbors['enzDest'])[0])
if isinstance(parentPool,PoolBase):
@@ -558,16 +567,21 @@ class GraphicalView(QtGui.QGraphicsView):
if messgtype.lower() == "all":
messgstr = "The object name \'{0}\' exist in \'{1}\' {2}".format(movedGraphObj.mobj.name,item.mobj.name,desObj)
elif messgtype.lower() =="enzymeparent":
- messgstr = "The Enzyme parent \'{0}\' doesn't exist in \'{2}\' {1} \n If you need to move the enzyme to {1} first parent pool needs to be moved".format(movedGraphObj.mobj.parent.name,desObj,item.mobj.name)
+ messgstr = "The Enzyme parent \'{0}\' doesn't exist in \'{2}\' {1} \n If you need to move the enzyme to {1} first move the parent pool".format(movedGraphObj.mobj.parent.name,desObj,item.mobj.name)
elif messgtype.lower() == "enzyme":
messgstr = "The Enzyme \'{0}\' already exist in \'{2}\' {1}".format(movedGraphObj.mobj.name,desObj,item.mobj.name)
elif messgtype.lower() == "empty":
messgstr = "The object can't be moved to empty space"
+ elif messgtype.lower() =="functionparent":
+ messgstr = "The Function parent \'{0}\' doesn't exist in \'{2}\' {1} \n If you need to move the function to {1} first move the parent pool".format(movedGraphObj.mobj.parent.name,desObj,item.mobj.name)
+
QtGui.QMessageBox.warning(None,'Could not move the object', messgstr )
QtGui.QApplication.setOverrideCursor(QtGui.QCursor(Qt.Qt.ArrowCursor))
def moveObjSceneParent(self,item,movedGraphObj,itempos,eventpos):
''' Scene parent object needs to be updated '''
+ if isinstance(movedGraphObj,FuncItem):
+ return
prevPar = movedGraphObj.parentItem()
movedGraphObj.setParentItem(item)