From 6733daac7a05453ee3b1036c06741d6aeb4e45d7 Mon Sep 17 00:00:00 2001
From: Dilawar Singh <dilawar@users.noreply.github.com>
Date: Thu, 20 Sep 2018 16:03:05 +0530
Subject: [PATCH] 2018 sep20 (#243)
* Pull subtree to master branch.
* Squashed 'moose-core/' changes from d229eba6bb..8463cc73e5
8463cc73e5 Update setup.cmake.py (#307)
8c38fc6d60 Removed deprecated warnings from neuroml2 reader (#305)
6df4332d2d Fixes to BOOST based steadystate solver (#306)
99dd7d2503 Merge pull request #302 from BhallaLab/chhennapoda
86c4244522 Merge branch 'master' into chhennapoda
ec06b242ae HotFix: Fix regression in StreadyState solver caused by #293. (#304)
6b79f6701c prefer .so suffix over CPython version specific suffix . this way a wheel compiled by any version of python3 will work with any python3.
92da6a7fb8 removed python-libsbml from install dependencies; added matplotlib. sbml and neuroml packages should be installed by user and not automatically. They may not be available for all distribution of python on all OSes.
e2bc19c0e4 Merge branch 'master' into chhennapoda
385a5cf0a1 Merge pull request #303 from upibhalla/master
080767c5aa Fixes to rdesigneur due to bad merge. Cleanup on rmoogli
62f8bc89e0 Tweak to Neuron.cpp so it can handle insertion of spines with uniform spacing.
b6088109bf Fixes for linux.
406fee2a06 Temp commit to test on OSX.
cc3c604b9f Merge branch 'master' into chhennapoda
342092829c Update moose_test.py (#301)
61f4f66843 show timeout in moose.test() .
abf822c4a6 with python2, just test with GSL. Build python3 first.
088c61673c Made changes to build system so that pymoose directory can be built in isolation. Build system is bit more module. Ideally I should make sure that each subdirectory can built in isolation.
3c0ff27903 Few changes to build with python3.6.
5d688f63d8 Use option --relative in setup.py and tell cmake not to look into usr prefix anymore.
1f3fe9fc35 Merge branch 'chhennapoda' of github.com:BhallaLab/moose-core into chhennapoda
3e19983c83 removed unneeded library which is never used.
c2fbab8726 cmake related changes to build on centos5 (manylinux docker image).
f1230f1518 Merge branch 'master' of https://github.com/upibhalla/moose-core
23a2892cc8 Merge branch 'master' of https://github.com/BhallaLab/moose-core
8a835c1332 minor bugfix to rdesigneur
6ef1d567ec Merge branch 'master' into master
a9a2e758ff Merge branch 'master' of https://github.com/BhallaLab/moose-core
5e1294d991 Further updates to rdesigneur for argument handling
a702bded54 Put in kwargs based specification of plotList, stimList and moogList. Folded savePlots into plotList but not yet tested.
git-subtree-dir: moose-core
git-subtree-split: 8463cc73e5d837ade8e191404359bc278f377b9d
* Fixes to build system.
* Some test only runs on travis only.
* Don't run moose.test()
---
.travis_build_linux.sh | 8 +
.travis_build_osx.sh | 1 -
CMakeLists.txt | 5 +-
moose-core/.travis.yml | 4 +-
moose-core/.travis/travis_build_linux.sh | 18 +-
moose-core/.travis/travis_prepare_linux.sh | 19 +-
moose-core/CMakeLists.txt | 129 +++++------
moose-core/biophysics/ChanCommon.cpp | 66 +++---
moose-core/biophysics/ChanCommon.h | 143 ++++++-------
moose-core/biophysics/Leakage.cpp | 19 +-
moose-core/biophysics/Leakage.h | 2 +-
moose-core/builtins/CMakeLists.txt | 17 +-
moose-core/builtins/testBuiltins.cpp | 4 +
moose-core/ksolve/CMakeLists.txt | 54 ++---
moose-core/ksolve/Ksolve.cpp | 4 -
moose-core/ksolve/NonlinearSystem.h | 200 ++++++++----------
moose-core/ksolve/SteadyStateBoost.cpp | 121 ++++++-----
moose-core/ksolve/SteadyStateGsl.cpp | 5 -
moose-core/ksolve/VoxelPools.cpp | 42 ++--
moose-core/pymoose/CMakeLists.txt | 12 +-
moose-core/python/moose/moose.py | 28 ++-
moose-core/python/moose/neuroml/ChannelML.py | 11 +-
moose-core/python/moose/neuroml/MorphML.py | 13 +-
moose-core/python/moose/neuroml/NetworkML.py | 2 +-
moose-core/python/moose/neuroml/NeuroML.py | 45 ++--
moose-core/python/moose/neuroml2/reader.py | 184 +++++++---------
moose-core/python/setup.cmake.py | 9 +-
moose-core/tests/python/test_Xreacs2.py | 6 +-
moose-core/tests/python/test_dose_response.py | 116 ++++++++++
moose-core/tests/python/test_neuroml2.py | 12 ++
.../tests/python/test_steady_state_solver.py | 152 +++++++++++++
scripts/pull_subtree.sh | 6 +-
32 files changed, 865 insertions(+), 592 deletions(-)
create mode 100644 moose-core/tests/python/test_dose_response.py
create mode 100644 moose-core/tests/python/test_steady_state_solver.py
diff --git a/.travis_build_linux.sh b/.travis_build_linux.sh
index a9397370..4906319c 100755
--- a/.travis_build_linux.sh
+++ b/.travis_build_linux.sh
@@ -4,3 +4,11 @@ set -e -x
PATH=/usr/bin:/usr/local/bin:$PATH
gbp buildpackage --git-ignore-branch --git-ignore-new -uc -us -d | tee _gbp.log
pwd
+ls ../*.deb
+
+if [ -n "$TRAVIS" ]; then
+ echo "We are on travis"
+ sudo dpkg -i ../*.deb
+ sudo apt-get install -f
+ python -c 'import moose; print(moose.__file__); print(moose.__version__)'
+fi
diff --git a/.travis_build_osx.sh b/.travis_build_osx.sh
index 3f08edd1..10d39acf 100755
--- a/.travis_build_osx.sh
+++ b/.travis_build_osx.sh
@@ -3,4 +3,3 @@ brew tap BhallaLab/moose
brew install --HEAD moose-nightly
python -m pip install networkx python-libsbml pyNeuroML matplotlib
python -c 'import moose; print( moose.__version__ )'
-python -c 'import moose; print( moose.test() )'
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 756c329d..cf83ee35 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -108,10 +108,7 @@ if(WITH_GUI)
add_dependencies(examples gui)
endif()
-
-# NOTE: leading / is important other whole path will be installed. Do not
-# install usr since it comes with CMAKE_INSTALL_PREFIX.
-install(DIRECTORY ${PYMOOSE_INSTALL_DIR}/usr/
+install(DIRECTORY ${PYMOOSE_INSTALL_DIR}/
DESTINATION ${CMAKE_INSTALL_PREFIX}
CONFIGURATIONS Release
PATTERN ".git" EXCLUDE
diff --git a/moose-core/.travis.yml b/moose-core/.travis.yml
index a89ee791..01188d6d 100644
--- a/moose-core/.travis.yml
+++ b/moose-core/.travis.yml
@@ -1,7 +1,7 @@
language: cpp
sudo: required
-group: edge
-
+group: travis_lts
+dist: xenial
os:
- linux
diff --git a/moose-core/.travis/travis_build_linux.sh b/moose-core/.travis/travis_build_linux.sh
index f715e6c4..6885c7c6 100755
--- a/moose-core/.travis/travis_build_linux.sh
+++ b/moose-core/.travis/travis_build_linux.sh
@@ -40,24 +40,34 @@ unset PYTHONPATH
$PYTHON2 -m compileall -q .
if type $PYTHON3 > /dev/null; then $PYTHON3 -m compileall -q . ; fi
-# This is only applicable on linux build.
+# Python2 and GSL.
+echo "Currently in `pwd`"
+(
+ mkdir -p _GSL_BUILD && cd _GSL_BUILD
+ cmake -DDEBUG=ON -DPYTHON_EXECUTABLE="$PYTHON2" ..
+ $MAKE && ctest --output-on-failure
+ sudo make install && cd /tmp
+ $PYTHON2 -c 'import moose;print(moose.__file__);print(moose.version())'
+)
+
+# Python3 with GSL and BOOST. This should be the only build after we drop
+# python2 support.
echo "Python3: Removed python2-networkx and install python3"
if type $PYTHON3 > /dev/null; then
sudo apt-get remove -qq python-networkx || echo "Error with apt"
sudo apt-get install -qq python3-networkx || echo "Error with apt"
-
# GSL.
(
mkdir -p _GSL_BUILD2 && cd _GSL_BUILD2 && \
cmake -DPYTHON_EXECUTABLE="$PYTHON3" -DDEBUG=ON ..
- $MAKE && ctest --output-on-failure
+ $MAKE && ctest --output-on-failure
)
# BOOST
(
mkdir -p _BOOST_BUILD2 && cd _BOOST_BUILD2 && \
cmake -DWITH_BOOST_ODE=ON -DPYTHON_EXECUTABLE="$PYTHON3" ..
- $MAKE && ctest --output-on-failure
+ $MAKE && ctest --output-on-failure
)
echo "All done"
else
diff --git a/moose-core/.travis/travis_prepare_linux.sh b/moose-core/.travis/travis_prepare_linux.sh
index 484f1a4e..0c0e604f 100755
--- a/moose-core/.travis/travis_prepare_linux.sh
+++ b/moose-core/.travis/travis_prepare_linux.sh
@@ -10,7 +10,7 @@
# OPTIONS: ---
# REQUIREMENTS: ---
# BUGS: ---
-# NOTES: ---
+# NOTES: Always run with sudo permission.
# AUTHOR: Dilawar Singh (), dilawars@ncbs.res.in
# ORGANIZATION: NCBS Bangalore
# CREATED: 01/02/2017 10:10:02 AM
@@ -21,15 +21,24 @@ set -o nounset # Treat unset variables as an error
set +e # Let installation fail in some command
apt-get install -qq libxml2-dev libbz2-dev
-apt-get install -qq libhdf5-serial-dev
apt-get install -qq make cmake
apt-get install -qq python-numpy python-matplotlib python-networkx python-pip
+apt-get install -qq python3-lxml python-lxml
apt-get install -qq python3-numpy python3-matplotlib python3-dev
-apt-get install -qq libboost-all-dev
-apt-get install -qq libgsl0-dev
apt-get install -qq python-pip python3-pip
apt-get install -qq libgraphviz-dev
+# Gsl
+apt-get install -qq libgsl0-dev || apt-get install -qq libgsl-dev
+
+# Boost related.
+apt-get install -qq liblapack-dev
+apt-get install -qq libboost-all-dev
+
# Dependencies for NML2
apt-get install -qq python-scipy python3-scipy
-#pip install pyNeuroML libNeuroML
+apt-get install -qq python-lxml python3-lxml
+apt-get install -qq python-setuptools python3-setuptools
+apt-get install -qq python-tornado python3-tornado
+python2 -m pip install pyNeuroML libNeuroML
+python3 -m pip install pyNeuroML libNeuroML
diff --git a/moose-core/CMakeLists.txt b/moose-core/CMakeLists.txt
index 3be4235b..a71896f6 100644
--- a/moose-core/CMakeLists.txt
+++ b/moose-core/CMakeLists.txt
@@ -76,6 +76,8 @@ else()
endif()
################################ CMAKE OPTIONS ##################################
+option(WITH_NSDF "Enable NSDF support. Requires hdf5" OFF )
+
option(DEBUG "Build with debug support" OFF)
option(GPROF "Build for profiling using gprof" OFF)
option(ENABLE_UNIT_TESTS "Enable unit tests (DEBUG should also be ON)" OFF)
@@ -187,68 +189,71 @@ if(WITH_BOOST_ODE)
add_definitions(-DUSE_BOOST_ODE -UUSE_GSL)
endif()
-
-find_package(HDF5 COMPONENTS CXX HL)
-if(NOT HDF5_FOUND)
- message(
- "==================================================================\n"
- " HDF5 not found. Disabling support. Required for nsdf. \n\n"
- " If you need hdf5 support, please install hdf5-dev or hdf5-devel\n"
- " package or equivalent.\n\n"
- " $ sudo apt-get install libhdf5-dev \n"
- " $ sudo yum install libhdf5-devel \n"
- " $ brew install hdf5 \n\n"
- " Otherwise, continue with 'make' and 'make install' \n"
- " If you install hdf5 to non-standard path, export environment \n"
- " variable HDF5_ROOT to the location. Rerun cmake \n"
- "================================================================ \n"
- )
-endif(NOT HDF5_FOUND)
-
-if(HDF5_FOUND)
- include_directories( ${HDF5_INCLUDE_DIRS} )
- add_definitions( -DUSE_HDF5 )
- if(HDF5_USE_STATIC_LIBRARIES)
- message(STATUS "Finding static HDF5 libraries in $ENV{HDF5_ROOT}")
- find_library(HDF5_CXX_LIBRARIES NAMES libhdf5.a
- PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
+# NSDF5 support. Disabled by default.
+if(WITH_NSDF)
+ find_package(HDF5 COMPONENTS CXX HL)
+ if(NOT HDF5_FOUND)
+ message(
+ "==================================================================\n"
+ " HDF5 not found. Disabling NSDF support.\n\n"
+ " If you need NSDF support, please install hdf5-dev or hdf5-devel\n"
+ " package or equivalent.\n\n"
+ " $ sudo apt-get install libhdf5-dev \n"
+ " $ sudo yum install libhdf5-devel \n"
+ " $ brew install hdf5 \n\n"
+ " Otherwise, continue with 'make' and 'make install' \n"
+ " If you install hdf5 to non-standard path, export environment \n"
+ " variable HDF5_ROOT to the location. Rerun cmake \n"
+ "================================================================ \n"
)
- find_library(HDF5_HL_LIBRARIES NAMES libhdf5_hl.a
- PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
- )
- set(HDF5_LIBRARIES ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES})
+ endif(NOT HDF5_FOUND)
+
+ if(HDF5_FOUND)
+ include_directories( ${HDF5_INCLUDE_DIRS} )
+ add_definitions( -DUSE_HDF5 -DENABLE_NSDF )
+ if(HDF5_USE_STATIC_LIBRARIES)
+ message(STATUS "Finding static HDF5 libraries in $ENV{HDF5_ROOT}")
+ find_library(HDF5_CXX_LIBRARIES NAMES libhdf5.a
+ PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
+ )
+ find_library(HDF5_HL_LIBRARIES NAMES libhdf5_hl.a
+ PATHS $ENV{HDF5_ROOT}/lib $ENV{HDF5_ROOT}/lib64
+ )
+ set(HDF5_LIBRARIES ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES})
+ endif()
+
+
+ # Make sure, HDF5_HL_LIBRARIES are set. The COMPONENTS in find_package may
+ # or may not work. See BhallaLab/moose-core#163.
+ if(NOT HDF5_HL_LIBRARIES)
+ set(HDF5_HL_LIBRARIES ${HDF5_HL_LIBRARIES})
+ endif(NOT HDF5_HL_LIBRARIES)
+ list(APPEND HDF5_LIBRARIES ${HDF5_HL_LIBRARIES})
+
+ # message(STATUS "MOOSE will use following HDF5 ${HDF5_LIBRARIES}" )
+ foreach(HDF5_LIB ${HDF5_LIBRARIES})
+ if(HDF5_LIB)
+ get_filename_component( HDF5_LIB_EXT ${HDF5_LIB} EXT )
+ if(HDF5_LIB_EXT)
+ if(${HDF5_LIB_EXT} STREQUAL ".a")
+ list(APPEND STATIC_LIBRARIES ${HDF5_LIB} )
+ else( )
+ list(APPEND SYSTEM_SHARED_LIBS ${HDF5_LIB} )
+ endif( )
+ endif()
+ endif( )
+ endforeach( )
+ else( HDF5_FOUND )
+ message(STATUS "HDF5 is not found" )
+ endif( HDF5_FOUND )
+ # This is a fix for new HDF5 package on Debian/Ubuntu which installs hdf5
+ # headers in non-standard path. issue #80.
+ if(HDF5_LIBRARY_DIRS)
+ set_target_properties( libmoose PROPERTIES LINK_FLAGS "-L${HDF5_LIBRARY_DIRS}" )
endif()
-
-
- # Make sure, HDF5_HL_LIBRARIES are set. The COMPONENTS in find_package may
- # or may not work. See BhallaLab/moose-core#163.
- if(NOT HDF5_HL_LIBRARIES)
- set(HDF5_HL_LIBRARIES ${HDF5_HL_LIBRARIES})
- endif(NOT HDF5_HL_LIBRARIES)
- list(APPEND HDF5_LIBRARIES ${HDF5_HL_LIBRARIES})
-
- # message(STATUS "MOOSE will use following HDF5 ${HDF5_LIBRARIES}" )
- foreach(HDF5_LIB ${HDF5_LIBRARIES})
- if(HDF5_LIB)
- get_filename_component( HDF5_LIB_EXT ${HDF5_LIB} EXT )
- if(HDF5_LIB_EXT)
- if(${HDF5_LIB_EXT} STREQUAL ".a")
- list(APPEND STATIC_LIBRARIES ${HDF5_LIB} )
- else( )
- list(APPEND SYSTEM_SHARED_LIBS ${HDF5_LIB} )
- endif( )
- endif()
- endif( )
- endforeach( )
-else( HDF5_FOUND )
- message(STATUS "HDF5 is not found" )
-endif( HDF5_FOUND )
-
-# This is a fix for new HDF5 package on Debian/Ubuntu which installs hdf5
-# headers in non-standard path. issue #80.
-if(HDF5_LIBRARY_DIRS)
- set_target_properties( libmoose PROPERTIES LINK_FLAGS "-L${HDF5_LIBRARY_DIRS}" )
-endif()
+else(WITH_NSDF)
+ message(STATUS "NSDF support is disabled" )
+endif(WITH_NSDF)
# Openmpi
if(WITH_MPI)
@@ -432,7 +437,7 @@ message(STATUS "binary distribution file ${PYMOOSE_BDIST_FILE}")
add_custom_target(bdist ALL DEPENDS ${PYMOOSE_BDIST_FILE} )
# Any command using setup.cmake.py must run in the same directory. Use option
-# `--relative` to prefix is aways fixed.
+# `--relative` to prefix is aways fixed.
add_custom_command( OUTPUT ${PYMOOSE_BDIST_FILE}
COMMAND ${PYTHON_EXECUTABLE} setup.cmake.py bdist_dumb -p ${_platform}
-d ${PYMOOSE_BDIST_DIR} --relative
@@ -464,7 +469,7 @@ add_dependencies( copy_pymoose _moose )
install(TARGETS moose.bin DESTINATION bin CONFIGURATIONS Debug)
install(TARGETS libmoose DESTINATION lib CONFIGURATIONS Debug)
-# install pymoose bdist.
+# install pymoose bdist.
install(DIRECTORY ${PYMOOSE_BDIST_INSTALL_DIR}/
DESTINATION ${CMAKE_INSTALL_PREFIX}
CONFIGURATIONS Release Debug
diff --git a/moose-core/biophysics/ChanCommon.cpp b/moose-core/biophysics/ChanCommon.cpp
index 9a015e4f..124fe67f 100644
--- a/moose-core/biophysics/ChanCommon.cpp
+++ b/moose-core/biophysics/ChanCommon.cpp
@@ -15,17 +15,19 @@
// Constructor
///////////////////////////////////////////////////
ChanCommon::ChanCommon()
- :
- Vm_( 0.0 ),
- Gbar_( 0.0 ), modulation_( 1.0 ),
- Ek_( 0.0 ),
- Gk_( 0.0 ), Ik_( 0.0 )
+ :
+ Vm_( 0.0 ),
+ Gbar_( 0.0 ), modulation_( 1.0 ),
+ Ek_( 0.0 ),
+ Gk_( 0.0 ), Ik_( 0.0 )
{
- ;
+ ;
}
ChanCommon::~ChanCommon()
-{;}
+{
+ ;
+}
///////////////////////////////////////////////////
// Field function definitions
@@ -33,54 +35,54 @@ ChanCommon::~ChanCommon()
void ChanCommon::vSetGbar( const Eref& e, double Gbar )
{
- Gbar_ = Gbar;
+ Gbar_ = Gbar;
}
double ChanCommon::vGetGbar( const Eref& e ) const
{
- return Gbar_;
+ return Gbar_;
}
void ChanCommon::vSetModulation( const Eref& e, double modulation )
{
- modulation_ = modulation;
+ modulation_ = modulation;
}
double ChanCommon::vGetModulation( const Eref& e ) const
{
- return modulation_;
+ return modulation_;
}
double ChanCommon::getModulation() const
{
- return modulation_;
+ return modulation_;
}
void ChanCommon::vSetEk( const Eref& e, double Ek )
{
- Ek_ = Ek;
+ Ek_ = Ek;
}
double ChanCommon::vGetEk( const Eref& e ) const
{
- return Ek_;
+ return Ek_;
}
void ChanCommon::vSetGk( const Eref& e, double Gk )
{
- Gk_ = Gk;
+ Gk_ = Gk;
}
double ChanCommon::vGetGk( const Eref& e ) const
{
- return Gk_;
+ return Gk_;
}
void ChanCommon::vSetIk( const Eref& e, double Ik )
{
- Ik_ = Ik;
+ Ik_ = Ik;
}
double ChanCommon::vGetIk( const Eref& e ) const
{
- return Ik_;
+ return Ik_;
}
///////////////////////////////////////////////////
@@ -89,7 +91,7 @@ double ChanCommon::vGetIk( const Eref& e ) const
void ChanCommon::vHandleVm( double Vm )
{
- Vm_ = Vm;
+ Vm_ = Vm;
}
///////////////////////////////////////////////////
@@ -99,34 +101,34 @@ void ChanCommon::vHandleVm( double Vm )
void ChanCommon::sendProcessMsgs( const Eref& e, const ProcPtr info )
{
- ChanBase::channelOut()->send( e, Gk_, Ek_ );
- // This is used if the channel connects up to a conc pool and
- // handles influx of ions giving rise to a concentration change.
- ChanBase::IkOut()->send( e, Ik_ );
- // Needed by GHK-type objects
- ChanBase::permeability()->send( e, Gk_ );
+ ChanBase::channelOut()->send( e, Gk_, Ek_ );
+ // This is used if the channel connects up to a conc pool and
+ // handles influx of ions giving rise to a concentration change.
+ ChanBase::IkOut()->send( e, Ik_ );
+ // Needed by GHK-type objects
+ ChanBase::permeability()->send( e, Gk_ );
}
void ChanCommon::sendReinitMsgs( const Eref& e, const ProcPtr info )
{
- ChanBase::channelOut()->send( e, Gk_, Ek_ );
- ChanBase::IkOut()->send( e, Ik_ );
- // Needed by GHK-type objects
- ChanBase::permeability()->send( e, Gk_ );
+ ChanBase::channelOut()->send( e, Gk_, Ek_ );
+ ChanBase::IkOut()->send( e, Ik_ );
+ // Needed by GHK-type objects
+ ChanBase::permeability()->send( e, Gk_ );
}
void ChanCommon::updateIk()
{
- Ik_ = ( Ek_ - Vm_ ) * Gk_;
+ Ik_ = ( Ek_ - Vm_ ) * Gk_;
}
double ChanCommon::getVm() const
{
- return Vm_;
+ return Vm_;
}
double ChanCommon::getGbar() const
{
- return Gbar_;
+ return Gbar_;
}
diff --git a/moose-core/biophysics/ChanCommon.h b/moose-core/biophysics/ChanCommon.h
index 007c2f37..50725bb6 100644
--- a/moose-core/biophysics/ChanCommon.h
+++ b/moose-core/biophysics/ChanCommon.h
@@ -20,78 +20,75 @@
class ChanCommon: public virtual ChanBase
{
- public:
- ChanCommon();
- ~ChanCommon();
-
- /////////////////////////////////////////////////////////////
- // Value field access function definitions
- /////////////////////////////////////////////////////////////
-
- void vSetGbar( const Eref& e, double Gbar );
- double vGetGbar( const Eref& e ) const;
- void vSetModulation( const Eref& e, double modulation );
- double vGetModulation( const Eref& e ) const;
- double getModulation() const;
- void vSetEk( const Eref& e, double Ek );
- double vGetEk( const Eref& e ) const;
- void vSetGk( const Eref& e, double Gk );
- double vGetGk( const Eref& e ) const;
- /// Ik is read-only for MOOSE, but we provide the set
- /// func for derived classes to update it.
- void vSetIk( const Eref& e, double Ic );
- double vGetIk( const Eref& e ) const;
-
- /////////////////////////////////////////////////////////////
- // Dest function definitions
- /////////////////////////////////////////////////////////////
-
- /**
- * Assign the local Vm_ to the incoming Vm from the compartment
- */
- void vHandleVm( double Vm );
-
- /////////////////////////////////////////////////////////////
- /**
- * This function sends out the messages expected of a channel,
- * after process. Used as a utility by various derived classes.
- */
- void sendProcessMsgs( const Eref& e, const ProcPtr info );
- void sendReinitMsgs( const Eref& e, const ProcPtr info );
-
- /////////////////////////////////////////////////////////////
-
- /**
- * Utility function for a common computation using local variables
- */
- void updateIk();
-
-
- /// Utility function to access Vm
- double getVm() const;
-
- /// Utility function to acces Gbar
- double getGbar() const;
-
- /// Specify the Class Info static variable for initialization.
- static const Cinfo* initCinfo();
- protected:
- /// Vm_ is input variable from compartment, used for most rates
- double Vm_;
-
- private:
- /// Channel maximal conductance
- double Gbar_;
- /// Channel modulation. Scales conductance.
- double modulation_;
- /// Reversal potential of channel
- double Ek_;
-
- /// Channel actual conductance depending on opening of gates.
- double Gk_;
- /// Channel current
- double Ik_;
+public:
+ ChanCommon();
+ ~ChanCommon();
+
+ /////////////////////////////////////////////////////////////
+ // Value field access function definitions
+ /////////////////////////////////////////////////////////////
+
+ void vSetGbar( const Eref& e, double Gbar );
+ double vGetGbar( const Eref& e ) const;
+ void vSetModulation( const Eref& e, double modulation );
+ double vGetModulation( const Eref& e ) const;
+ double getModulation() const;
+ void vSetEk( const Eref& e, double Ek );
+ double vGetEk( const Eref& e ) const;
+ void vSetGk( const Eref& e, double Gk );
+ double vGetGk( const Eref& e ) const;
+ /// Ik is read-only for MOOSE, but we provide the set
+ /// func for derived classes to update it.
+ void vSetIk( const Eref& e, double Ic );
+ double vGetIk( const Eref& e ) const;
+
+ /////////////////////////////////////////////////////////////
+ // Dest function definitions
+ /////////////////////////////////////////////////////////////
+
+ /**
+ * Assign the local Vm_ to the incoming Vm from the compartment
+ */
+ void vHandleVm( double Vm );
+
+ /////////////////////////////////////////////////////////////
+ /**
+ * This function sends out the messages expected of a channel,
+ * after process. Used as a utility by various derived classes.
+ */
+ void sendProcessMsgs( const Eref& e, const ProcPtr info );
+ void sendReinitMsgs( const Eref& e, const ProcPtr info );
+
+ /////////////////////////////////////////////////////////////
+
+ /**
+ * Utility function for a common computation using local variables
+ */
+ void updateIk();
+
+ /// Utility function to access Vm
+ double getVm() const;
+
+ /// Utility function to acces Gbar
+ double getGbar() const;
+
+ /// Specify the Class Info static variable for initialization.
+ static const Cinfo* initCinfo();
+protected:
+ /// Vm_ is input variable from compartment, used for most rates
+ double Vm_;
+
+private:
+ /// Channel maximal conductance
+ double Gbar_;
+ /// Channel modulation. Scales conductance.
+ double modulation_;
+ /// Reversal potential of channel
+ double Ek_;
+
+ /// Channel actual conductance depending on opening of gates.
+ double Gk_;
+ /// Channel current
+ double Ik_;
};
-
-
#endif // _ChanCommon_h
diff --git a/moose-core/biophysics/Leakage.cpp b/moose-core/biophysics/Leakage.cpp
index ae78b6ee..b3bf9b71 100644
--- a/moose-core/biophysics/Leakage.cpp
+++ b/moose-core/biophysics/Leakage.cpp
@@ -52,7 +52,8 @@
const Cinfo* Leakage::initCinfo()
{
- static string doc[] = {
+ static string doc[] =
+ {
"Name", "Leakage",
"Author", "Subhasis Ray, 2009, Upi Bhalla 2014 NCBS",
"Description", "Leakage: Passive leakage channel."
@@ -63,8 +64,8 @@ const Cinfo* Leakage::initCinfo()
static Cinfo LeakageCinfo(
"Leakage",
ChanBase::initCinfo(),
- 0,
- 0,
+ 0,
+ 0,
&dinfo,
doc,
sizeof( doc ) / sizeof( string ));
@@ -88,22 +89,22 @@ Leakage::~Leakage()
void Leakage::vProcess( const Eref & e, ProcPtr p )
{
- ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
- updateIk();
+ ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
+ updateIk();
sendProcessMsgs(e, p);
}
void Leakage::vReinit( const Eref & e, ProcPtr p )
{
- ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
- updateIk();
+ ChanCommon::vSetGk( e, this->vGetGbar( e ) * this->vGetModulation( e ));
+ updateIk();
sendReinitMsgs(e, p);
}
void Leakage::vSetGbar( const Eref& e, double gbar )
{
- ChanCommon::vSetGk( e, gbar * this->vGetModulation( e ) );
- ChanCommon::vSetGbar( e, gbar );
+ ChanCommon::vSetGk( e, gbar * this->vGetModulation( e ) );
+ ChanCommon::vSetGbar( e, gbar );
}
//
diff --git a/moose-core/biophysics/Leakage.h b/moose-core/biophysics/Leakage.h
index 959b5327..3b20eaf7 100644
--- a/moose-core/biophysics/Leakage.h
+++ b/moose-core/biophysics/Leakage.h
@@ -51,7 +51,7 @@
class Leakage: public ChanCommon
{
- public:
+public:
Leakage();
~Leakage();
void vProcess( const Eref & e, ProcPtr p );
diff --git a/moose-core/builtins/CMakeLists.txt b/moose-core/builtins/CMakeLists.txt
index bbc4ab70..64caa561 100644
--- a/moose-core/builtins/CMakeLists.txt
+++ b/moose-core/builtins/CMakeLists.txt
@@ -3,7 +3,7 @@ include_directories( ${CMAKE_SOURCE_DIR}/basecode )
include_directories( ${CMAKE_SOURCE_DIR}/external/muparser/include )
include_directories( ${CMAKE_SOURCE_DIR}/scheduling )
-add_library(moose_builtins
+set(SRCS
Arith.cpp
Group.cpp
Mstring.cpp
@@ -20,10 +20,17 @@ add_library(moose_builtins
Streamer.cpp
Stats.cpp
Interpol2D.cpp
- HDF5WriterBase.cpp
- NSDFWriter.cpp
- HDF5DataWriter.cpp
SpikeStats.cpp
testBuiltins.cpp
- testNSDF.cpp
)
+
+if(WITH_NSDF AND HDF5_FOUND)
+ list(APPEND SRCS
+ HDF5WriterBase.cpp
+ NSDFWriter.cpp
+ HDF5DataWriter.cpp
+ testNSDF.cpp
+ )
+endif()
+
+add_library(moose_builtins ${SRCS} )
diff --git a/moose-core/builtins/testBuiltins.cpp b/moose-core/builtins/testBuiltins.cpp
index c301a6a2..4c05470a 100644
--- a/moose-core/builtins/testBuiltins.cpp
+++ b/moose-core/builtins/testBuiltins.cpp
@@ -18,7 +18,9 @@
#include "../shell/Shell.h"
+#ifdef ENABLE_NSDF
extern void testNSDF();
+#endif
void testArith()
{
@@ -425,7 +427,9 @@ void testBuiltins()
{
testArith();
testTable();
+#if ENABLE_NSDF
testNSDF();
+#endif
}
void testBuiltinsProcess()
diff --git a/moose-core/ksolve/CMakeLists.txt b/moose-core/ksolve/CMakeLists.txt
index f86c69c0..c43875a0 100644
--- a/moose-core/ksolve/CMakeLists.txt
+++ b/moose-core/ksolve/CMakeLists.txt
@@ -2,9 +2,21 @@ cmake_minimum_required(VERSION 2.8)
include(CheckIncludeFileCXX)
include_directories(../builtins ../basecode ../utility ../kinetics)
-include_directories(../ )
include_directories(../mesh)
include_directories(../external/muparser/include )
+
+if(WITH_BOOST)
+ find_package(Boost REQUIRED COMPONENTS thread)
+ add_definitions( -DUSE_BOOST_ASYNC )
+ if("${CMAKE_BUILD_TYPE}" STREQUAL "Release")
+ add_definitions( -UBOOST_UBLAS_TYPE_CHECK )
+ else()
+ add_definitions( -DBOOST_UBLAS_TYPE_CHECK )
+ endif()
+ set(WITH_BOOST_ODE ON)
+ include_directories( ${Boost_INCLUDE_DIRS} )
+endif()
+
IF(WITH_BOOST_ODE)
check_include_file_cxx( ${Boost_INCLUDE_DIRS}/boost/numeric/odeint.hpp
ODEINT_EXISTS)
@@ -18,30 +30,24 @@ elseif(WITH_GSL)
include_directories( ${GSL_INCLUDE_DIRS} )
endif(WITH_BOOST_ODE)
-if(WITH_BOOST)
- find_package(Boost REQUIRED COMPONENTS thread)
- add_definitions( -DUSE_BOOST_ASYNC )
- include_directories( ${Boost_INCLUDE_DIRS} )
-endif()
-
set(KSOLVE_SRCS
- KinSparseMatrix.cpp
- ZombiePool.cpp
- ZombieFunction.cpp
- ZombieBufPool.cpp
- ZombieReac.cpp
- ZombieEnz.cpp
- ZombieMMenz.cpp
- VoxelPoolsBase.cpp
- VoxelPools.cpp
- GssaVoxelPools.cpp
- RateTerm.cpp
- FuncTerm.cpp
- Stoich.cpp
- Ksolve.cpp
- Gsolve.cpp
- ZombiePoolInterface.cpp
- testKsolve.cpp
+ KinSparseMatrix.cpp
+ ZombiePool.cpp
+ ZombieFunction.cpp
+ ZombieBufPool.cpp
+ ZombieReac.cpp
+ ZombieEnz.cpp
+ ZombieMMenz.cpp
+ VoxelPoolsBase.cpp
+ VoxelPools.cpp
+ GssaVoxelPools.cpp
+ RateTerm.cpp
+ FuncTerm.cpp
+ Stoich.cpp
+ Ksolve.cpp
+ Gsolve.cpp
+ ZombiePoolInterface.cpp
+ testKsolve.cpp
)
if(WITH_GSL)
diff --git a/moose-core/ksolve/Ksolve.cpp b/moose-core/ksolve/Ksolve.cpp
index 81bc118f..a072b4c0 100644
--- a/moose-core/ksolve/Ksolve.cpp
+++ b/moose-core/ksolve/Ksolve.cpp
@@ -235,11 +235,7 @@ static const Cinfo* ksolveCinfo = Ksolve::initCinfo();
Ksolve::Ksolve()
:
-#ifdef USE_GSL
method_( "rk5" ),
-#elif USE_BOOST_ODE
- method_( "rk5a" ),
-#endif
epsAbs_( 1e-7 ),
epsRel_( 1e-7 ),
numThreads_( 1 ),
diff --git a/moose-core/ksolve/NonlinearSystem.h b/moose-core/ksolve/NonlinearSystem.h
index 5792f287..e2ea270f 100644
--- a/moose-core/ksolve/NonlinearSystem.h
+++ b/moose-core/ksolve/NonlinearSystem.h
@@ -31,8 +31,6 @@
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
-
-
#include "VoxelPools.h"
using namespace std;
@@ -46,10 +44,10 @@ typedef ublas::matrix<value_type> matrix_type;
class ReacInfo
{
public:
- int rank;
- int num_reacs;
+ size_t rank;
+ size_t num_reacs;
size_t num_mols;
- int nIter;
+ size_t nIter;
double convergenceCriterion;
double* T;
VoxelPools* pool;
@@ -60,8 +58,10 @@ class ReacInfo
/* Matrix inversion routine.
- Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
- template<class T>
+ * Uses lu_factorize and lu_substitute in uBLAS to invert a matrix
+ */
+
+template<class T>
bool inverse(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
{
using namespace boost::numeric::ublas;
@@ -73,6 +73,7 @@ bool inverse(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
// perform LU-factorization
int res = lu_factorize(A,pm);
+
if( res != 0 ) return false;
// create identity matrix of "inverse"
@@ -103,14 +104,19 @@ public:
x1.resize( size_, 0);
ri.nVec.resize( size_ );
- dx_ = sqrt( numeric_limits<double>::epsilon() );
+
+ // Find machine epsilon to decide on dx. Machine eps is typically 2.22045e-16 on 64 bit machines/intel.
+ // Square root of this value is a very good value. The values matches with GSL solver nicely.
+ // Making this value too small or too big gonna cause error in blas methods especially when
+ // computing inverse.
+ // Gnu-GSL uses GSL_SQRT_DBL_EPSILON which is 1.4901161193847656e-08
+ // dx_ = 1.4901161193847656e-08;
+ dx_ = std::sqrt( numeric_limits<double>::epsilon() );
}
- vector_type compute_at(const vector_type& x)
+ bool compute_at(const vector_type& x, vector_type& result)
{
- vector_type result( size_ );
- system(x, result);
- return result;
+ return system(x, result);
}
int apply( )
@@ -118,22 +124,51 @@ public:
return system(x_, f_);
}
- int compute_jacobians( const vector_type& x, bool compute_inverse = true )
+ int compute_jacobians( const vector_type& x, bool compute_inverse = false )
{
for( size_t i = 0; i < size_; i++)
+ {
+ // This trick I leart by looking at GSL implmentation.
+ double dx = dx_ * std::fabs(x[i]);
+ if( dx == 0 )
+ dx = dx_;
+
for( size_t j = 0; j < size_; j++)
{
vector_type temp = x;
- temp[j] += dx_;
- J_(i, j) = (compute_at(temp)[i] - compute_at(x)[i]) / dx_;
+ temp[j] += dx;
+ vector_type res1, res2;
+ auto r1 = compute_at(temp, res1);
+ auto r2 = compute_at(x, res2);
+
+ if( 0 == r1 && 0 == r2 )
+ {
+ J_(i, j) = (res1[i] - res2[i]) / dx;
+ if( std::isnan(J_(i,j)) || std::isinf(J_(i,j)) )
+ {
+ /* Try increasing dx */
+ J_.clear();
+ return -1;
+ }
+ }
+ else
+ {
+ J_.clear();
+ return -1;
+ }
}
+ }
- // is_jacobian_valid_ = true;
- // Keep the inverted J_ ready
- //if(is_jacobian_valid_ and compute_inverse )
if( compute_inverse )
- inverse( J_, invJ_ );
-
+ {
+ try {
+ inverse( J_, invJ_ );
+ } catch ( exception & e ) {
+ J_.clear();
+ invJ_.clear();
+ return -1;
+ }
+ }
return 0;
}
@@ -147,9 +182,15 @@ public:
init[i] = x[i];
x_ = init;
- apply();
-
- compute_jacobians( init );
+ if( 0 == apply() )
+ {
+ if( 0 != compute_jacobians( init ) )
+ {
+ return;
+ }
+ }
+ else
+ return;
}
string to_string( )
@@ -167,22 +208,16 @@ public:
int system( const vector_type& x, vector_type& f )
{
- int num_consv = ri.num_mols - ri.rank;
+ size_t num_consv = ri.num_mols - ri.rank;
for ( size_t i = 0; i < ri.num_mols; ++i )
{
double temp = x[i] * x[i] ;
-#if 0
// if overflow
- if ( std::isnan( temp ) or std::isinf( temp ) )
- {
- cerr << "Failed: ";
- for( auto v : ri.nVec ) cerr << v << ", ";
- cerr << endl;
+ if ( std::isnan( temp ) || std::isinf( temp ) )
return -1;
- }
-#endif
- ri.nVec[i] = temp;
+ else
+ ri.nVec[i] = temp;
}
vector< double > vels;
@@ -192,21 +227,21 @@ public:
// y = Nr . v
// Note that Nr is row-echelon: diagonal and above.
- for ( int i = 0; i < ri.rank; ++i )
+ f.resize( ri.rank + num_consv);
+ for ( size_t i = 0; i < ri.rank; ++i )
{
double temp = 0;
- for ( int j = i; j < ri.num_reacs; ++j )
+ for ( size_t j = i; j < ri.num_reacs; ++j )
temp += ri.Nr(i, j ) * vels[j];
f[i] = temp ;
}
// dT = gamma.S - T
- for ( int i = 0; i < num_consv; ++i )
+ for ( size_t i = 0; i < num_consv; ++i )
{
double dT = - ri.T[i];
for ( size_t j = 0; j < ri.num_mols; ++j )
dT += ri.gamma(i, j) * x[j] * x[j];
-
f[ i + ri.rank] = dT ;
}
return 0;
@@ -222,99 +257,42 @@ public:
* @return If successful, return true. Check the variable `x_` at
* which the system f_ is close to zero (within the tolerance).
*/
- bool find_roots_gnewton( double tolerance = 1e-7 , size_t max_iter = 50)
+ bool find_roots_gnewton( double tolerance, size_t max_iter)
{
- //tolerance = sqrt( numeric_limits<double>::epsilon() );
double norm2OfDiff = 1.0;
size_t iter = 0;
- int status = apply();
+ if(0 != apply() )
+ return false;
- while( ublas::norm_2(f_) >= tolerance )
+ while( ublas::norm_2(f_) > tolerance )
{
iter += 1;
- compute_jacobians( x_, true );
+ ri.nIter = iter;
+
+ if( 0 != compute_jacobians( x_, true ) )
+ {
+ J_.clear();
+ invJ_.clear();
+ return false;
+ }
+
vector_type correction = ublas::prod( invJ_, f_ );
x_ -= correction;
// If could not compute the value of system successfully.
- status = apply();
- if( 0 != status )
+ if( 0 != apply() )
+ {
+ x_.clear();
return false;
+ }
if( iter >= max_iter )
break;
}
-
- ri.nIter = iter;
-
- if( iter >= max_iter )
- return false;
-
return true;
}
- /**
- * @brief Compute the slope of function in given dimension.
- *
- * @param which_dimen The index of dimension.
- *
- * @return Slope.
- */
- value_type slope( unsigned int which_dimen )
- {
- vector_type x = x_;
- x[which_dimen] += dx_;
- // x1 and x2 holds the f_ of system at x_ and x (which is x +
- // some step)
- system( x_, x1 );
- system( x, x2 );
- return ublas::norm_2( (x2 - x1)/dx_ );
- }
-
-
- /**
- * @brief Makes the first guess. After this call the Newton method.
- */
- void correction_step( )
- {
- // Get the jacobian at current point. Notice that in this method, we
- // don't have to compute inverse of jacobian
-
- vector_type direction( size_ );
-
- // Now take the largest step possible such that the value of system at
- // (x_ - step ) is lower than the value of system as x_.
- vector_type nextState( size_ );
-
- apply();
-
- unsigned int i = 0;
-
- double factor = 1e-2;
- while( true )
- {
- i += 1;
- compute_jacobians( x_, false );
- // Make a move in either side of direction. In whichever direction
- // the function decreases.
- direction = ublas::prod( J_, f_ );
- nextState = x_ - factor * direction;
- if( ublas::norm_2( compute_at( nextState ) ) >= ublas::norm_2(compute_at(x_)))
- factor = factor / 2.0;
- else
- {
- cerr << "Correction term applied ";
- x_ = nextState;
- apply();
- break;
- }
-
- if ( i > 20 )
- break;
- }
- }
-
public:
const size_t size_;
diff --git a/moose-core/ksolve/SteadyStateBoost.cpp b/moose-core/ksolve/SteadyStateBoost.cpp
index 5d31165d..98b67154 100644
--- a/moose-core/ksolve/SteadyStateBoost.cpp
+++ b/moose-core/ksolve/SteadyStateBoost.cpp
@@ -349,7 +349,7 @@ void SteadyState::setStoich( Id value )
setupSSmatrix();
double vol = LookupField< unsigned int, double >::get(stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
pool_.setVolume( vol );
- pool_.setStoich( stoichPtr, 0 );
+ pool_.setStoich( stoichPtr, nullptr );
pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
isInitialized_ = 1;
}
@@ -561,6 +561,7 @@ void SteadyState::setupSSmatrix()
total_.assign( nConsv, 0.0 );
Id ksolve = Field< Id >::get( stoich_, "ksolve" );
+
vector< double > nVec =
LookupField< unsigned int, vector< double > >::get(
ksolve,"nVec", 0 );
@@ -589,7 +590,7 @@ void SteadyState::classifyState( const double* T )
double tot = 0.0;
Stoich* s = reinterpret_cast< Stoich* >( stoich_.eref().data() );
vector< double > nVec = LookupField< unsigned int, vector< double > >::get(
- s->getKsolve(), "nVec", 0 );
+ s->getKsolve(), "nVec", 0 );
for ( unsigned int i = 0; i < numVarPools_; ++i )
tot += nVec[i];
@@ -600,14 +601,14 @@ void SteadyState::classifyState( const double* T )
for ( unsigned int i = 0; i < numVarPools_; ++i )
{
double orig = nVec[i];
- if ( std::isnan( orig ) or std::isinf( orig ) )
+ if ( std::isnan( orig ) || std::isinf( orig ) )
{
cout << "Warning: SteadyState::classifyState: orig=nan\n";
solutionStatus_ = 2; // Steady state OK, eig failed
J.clear();
return;
}
- if ( std::isnan( tot ) or std::isinf( tot ))
+ if ( std::isnan( tot ) || std::isinf( tot ))
{
cout << "Warning: SteadyState::classifyState: tot=nan\n";
solutionStatus_ = 2; // Steady state OK, eig failed
@@ -622,73 +623,72 @@ void SteadyState::classifyState( const double* T )
// Assign the rates for each mol.
for ( unsigned int j = 0; j < numVarPools_; ++j )
{
- if( std::isnan( yprime[j] ) or std::isinf( yprime[j] ) )
+ if( std::isnan( yprime[j] ) || std::isinf( yprime[j] ) )
{
- cout << "Warning: Overflow/underflow. Can't continue " << endl;
solutionStatus_ = 2;
+ J.clear();
return;
}
J(i, j) = yprime[j];
}
- }
- // Jacobian is now ready. Find eigenvalues.
- ublas::vector< std::complex< double > > eigenVec ( J.size1() );
+ // Jacobian is now ready. Find eigenvalues.
+ ublas::vector< std::complex< double > > eigenVec ( J.size1() );
- ublas::matrix< std::complex<double>, ublas::column_major >* vl, *vr;
- vl = NULL; vr = NULL;
+ ublas::matrix< std::complex<double>, ublas::column_major >* vl, *vr;
+ vl = NULL; vr = NULL;
- /*-----------------------------------------------------------------------------
- * INFO: Calling lapack routine geev to compute eigen vector of matrix J.
- *
- * Argument 3 and 4 are left- and right-eigenvectors. Since we do not need
- * them, they are set to NULL. Argument 2 holds eigen-vector and result is
- * copied to it (output ).
- *-----------------------------------------------------------------------------*/
- int status = lapack::geev( J, eigenVec, vl, vr, lapack::optimal_workspace() );
+ /*-----------------------------------------------------------------------------
+ * INFO: Calling lapack routine geev to compute eigen vector of matrix J.
+ *
+ * Argument 3 and 4 are left- and right-eigenvectors. Since we do not need
+ * them, they are set to NULL. Argument 2 holds eigen-vector and result is
+ * copied to it (output ).
+ *-----------------------------------------------------------------------------*/
+ int status = lapack::geev( J, eigenVec, vl, vr, lapack::optimal_workspace() );
- eigenvalues_.clear();
- eigenvalues_.resize( numVarPools_, 0.0 );
- if ( status != 0 )
- {
- cout << "Warning: SteadyState::classifyState failed to find eigenvalues. Status = " <<
- status << endl;
- solutionStatus_ = 2; // Steady state OK, eig classification failed
- }
- else // Eigenvalues are ready. Classify state.
- {
- nNegEigenvalues_ = 0;
- nPosEigenvalues_ = 0;
- for ( unsigned int i = 0; i < numVarPools_; ++i )
+ eigenvalues_.clear();
+ eigenvalues_.resize( numVarPools_, 0.0 );
+ if ( status != 0 )
{
- std::complex<value_type> z = eigenVec[ i ];
- double r = z.real();
- nNegEigenvalues_ += ( r < -EPSILON );
- nPosEigenvalues_ += ( r > EPSILON );
- eigenvalues_[i] = r;
- // We have a problem here because numVarPools_ usually > rank
- // This means we have several zero eigenvalues.
+ cout << "Warning: SteadyState::classifyState failed to find eigenvalues. Status = " <<
+ status << endl;
+ solutionStatus_ = 2; // Steady state OK, eig classification failed
+ }
+ else // Eigenvalues are ready. Classify state.
+ {
+ nNegEigenvalues_ = 0;
+ nPosEigenvalues_ = 0;
+ for ( unsigned int i = 0; i < numVarPools_; ++i )
+ {
+ std::complex<value_type> z = eigenVec[ i ];
+ double r = z.real();
+ nNegEigenvalues_ += ( r < -EPSILON );
+ nPosEigenvalues_ += ( r > EPSILON );
+ eigenvalues_[i] = r;
+ // We have a problem here because numVarPools_ usually > rank
+ // This means we have several zero eigenvalues.
+ }
+ if ( nNegEigenvalues_ == rank_ )
+ stateType_ = 0; // Stable
+ else if ( nPosEigenvalues_ == rank_ ) // Never see it.
+ stateType_ = 1; // Unstable
+ else if (nPosEigenvalues_ == 1)
+ stateType_ = 2; // Saddle
+ else if ( nPosEigenvalues_ >= 2 )
+ stateType_ = 3; // putative oscillatory
+ else if ( nNegEigenvalues_ == ( rank_ - 1) && nPosEigenvalues_ == 0 )
+ stateType_ = 4; // one zero or unclassified eigenvalue. Messy.
+ else
+ stateType_ = 5; // Other
}
- if ( nNegEigenvalues_ == rank_ )
- stateType_ = 0; // Stable
- else if ( nPosEigenvalues_ == rank_ ) // Never see it.
- stateType_ = 1; // Unstable
- else if (nPosEigenvalues_ == 1)
- stateType_ = 2; // Saddle
- else if ( nPosEigenvalues_ >= 2 )
- stateType_ = 3; // putative oscillatory
- else if ( nNegEigenvalues_ == ( rank_ - 1) && nPosEigenvalues_ == 0 )
- stateType_ = 4; // one zero or unclassified eigenvalue. Messy.
- else
- stateType_ = 5; // Other
}
}
static bool isSolutionValid( const vector< double >& x )
{
- for( size_t i = 0; i < x.size(); i++ )
+ for( auto &v : x )
{
- double v = x[i];
if ( std::isnan( v ) || std::isinf( v ) )
{
cout << "Warning: SteadyState iteration gave nan/inf concs\n";
@@ -703,6 +703,19 @@ static bool isSolutionValid( const vector< double >& x )
return true;
}
+static bool isSolutionPositive( const vector< double >& x )
+{
+ for( auto &v : x )
+ {
+ if( v < 0.0 )
+ {
+ cout << "Warning: SteadyState iteration gave negative concs" << endl;
+ return false;
+ }
+ }
+ return true;
+}
+
/**
* The settle function computes the steady state nearest the initial
* conditions.
@@ -773,7 +786,7 @@ void SteadyState::settle( bool forceSetup )
int status = 1;
// Find roots . If successful, set status to 0.
- if( ss->find_roots_gnewton( ) )
+ if( ss->find_roots_gnewton( convergenceCriterion_, maxIter_ ) )
status = 0;
if ( status == 0 && isSolutionValid( ss->ri.nVec ) )
diff --git a/moose-core/ksolve/SteadyStateGsl.cpp b/moose-core/ksolve/SteadyStateGsl.cpp
index 5aa0941d..eaec6d07 100644
--- a/moose-core/ksolve/SteadyStateGsl.cpp
+++ b/moose-core/ksolve/SteadyStateGsl.cpp
@@ -383,7 +383,6 @@ void SteadyState::setStoich( Id value )
pool_.setVolume( vol );
pool_.setStoich( stoichPtr, nullptr );
-
pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
isInitialized_ = 1;
}
@@ -936,13 +935,9 @@ int ss_func( const gsl_vector* x, void* params, gsl_vector* f )
{
double temp = op( gsl_vector_get( x, i ) );
if ( isNaN( temp ) || isInfinity( temp ) )
- {
return GSL_ERANGE;
- }
else
- {
ri->nVec[i] = temp;
- }
}
vector< double > vels;
ri->pool->updateReacVelocities( &ri->nVec[0], vels );
diff --git a/moose-core/ksolve/VoxelPools.cpp b/moose-core/ksolve/VoxelPools.cpp
index 25e8ac7e..53f5a173 100644
--- a/moose-core/ksolve/VoxelPools.cpp
+++ b/moose-core/ksolve/VoxelPools.cpp
@@ -1,11 +1,11 @@
-/**********************************************************************
-** This program is part of 'MOOSE', the
-** Messaging Object Oriented Simulation Environment.
-** Copyright (C) 2003-2010 Upinder S. Bhalla. and NCBS
-** It is made available under the terms of the
-** GNU Lesser General Public License version 2.1
-** See the file COPYING.LIB for the full notice.
-**********************************************************************/
+/*
+* This program is part of 'MOOSE', the
+* Messaging Object Oriented Simulation Environment.
+* Copyright (C) 2003-2010 Upinder S. Bhalla. and NCBS
+* It is made available under the terms of the
+* GNU Lesser General Public License version 2.1
+* See the file COPYING.LIB for the full notice.
+*/
#include "../basecode/header.h"
@@ -143,13 +143,20 @@ void VoxelPools::advance( const ProcInfo* p )
// 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
+ // Describe system to be used in boost solver calls.
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);
+ // This is usually the default method. It works well in practice. Tested
+ // with steady-state solver. Closest to GSL rk5 .
+ if( method_ == "rk5" || method_ == "gsl" || method_ == "boost" )
+ odeint::integrate_const(
+ make_dense_output( epsAbs_, epsRel_, odeint::runge_kutta_dopri5<vector_type_>() )
+ , sys , Svec() , p->currTime - p->dt , p->currTime , p->dt
+ );
+ else if( method_ == "rk5a" || method_ == "adaptive" )
+ odeint::integrate_adaptive( odeint::make_controlled<rk_dopri_stepper_type_>( epsAbs_, epsRel_ )
+ , sys , Svec() , p->currTime - p->dt , p->currTime, p->dt );
else if( method_ == "rk2" )
odeint::integrate_const( rk_midpoint_stepper_type_()
, sys, Svec(), p->currTime - p->dt, p->currTime, p->dt);
@@ -203,17 +210,6 @@ int VoxelPools::gslFunc( double t, const double* y, double *dydt,
VoxelPools* vp = reinterpret_cast< VoxelPools* >( params );
// Stoich* s = reinterpret_cast< Stoich* >( params );
double* q = const_cast< double* >( y ); // Assign the func portion.
-
- // Assign the buffered pools
- // Not possible because this is a static function
- // Not needed because dydt = 0;
- /*
- double* b = q + s->getNumVarPools();
- vector< double >::const_iterator sinit = Sinit_.begin() + s->getNumVarPools();
- for ( unsigned int i = 0; i < s->getNumBufPools(); ++i )
- *b++ = *sinit++;
- */
-
vp->stoichPtr_->updateFuncs( q, t );
vp->updateRates( y, dydt );
#ifdef USE_GSL
diff --git a/moose-core/pymoose/CMakeLists.txt b/moose-core/pymoose/CMakeLists.txt
index d25a0875..0a2a6388 100644
--- a/moose-core/pymoose/CMakeLists.txt
+++ b/moose-core/pymoose/CMakeLists.txt
@@ -4,18 +4,18 @@ set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../cmake_modules/")
find_package(PythonInterp REQUIRED)
-# Find Numpy
+# Find Numpy
find_package(NumPy REQUIRED)
include_directories(${NUMPY_INCLUDE_DIRS})
add_definitions( -std=c++11 )
add_definitions(-DUSE_NUMPY)
add_definitions(-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION)
-# set module extensiton. default is .so
+# set module extensiton. default is .so. Also check ../python/setup.cmake.py
execute_process( COMMAND
${PYTHON_EXECUTABLE} -c
"import importlib.machinery
-print(importlib.machinery.EXTENSION_SUFFIXES[0])"
+print(importlib.machinery.EXTENSION_SUFFIXES[-1])"
OUTPUT_VARIABLE PYTHON_SO_EXTENSION
OUTPUT_STRIP_TRAILING_WHITESPACE
)
@@ -25,7 +25,7 @@ if(NOT PYTHON_SO_EXTENSION)
endif()
message(STATUS "Python so extension ${PYTHON_SO_EXTENSION}" )
-# TARGET
+# TARGET
set(PYMOOSE_SRCS
moosemodule.cpp
vec.cpp
@@ -63,7 +63,7 @@ else()
COMPILE_DEFINITIONS "PYMOOSE"
COMPILE_FLAGS "${PYTHON_INCLUDE_FLAGS}"
)
-
+
endif()
# Remove prefix lib from python module.
@@ -82,7 +82,7 @@ if(MACOSX AND ("${PYTHON_VERSION_MAJOR}" STREQUAL "2"))
endif()
# see issue #80
-if(HDF5_LIBRARY_DIRS)
+if(HDF5_FOUND AND WITH_NSDF)
set_target_properties( _moose PROPERTIES LINK_FLAGS "-L${HDF5_LIBRARY_DIRS}" )
endif()
diff --git a/moose-core/python/moose/moose.py b/moose-core/python/moose/moose.py
index b9acb706..e6cf9bf7 100644
--- a/moose-core/python/moose/moose.py
+++ b/moose-core/python/moose/moose.py
@@ -29,10 +29,10 @@ try:
import moose.neuroml2 as _neuroml2
except Exception as e:
nml2Import_ = False
- nml2ImportError_ = '\n'.join( [
+ nml2ImportError_ = ' '.join( [
"NML2 support is disabled because `libneuroml` and "
- , "`pyneuroml` modules are not found."
- , " pip install pyneuroml libneuroml "
+ , "`pyneuroml` modules are not found.\n"
+ , " $ pip install pyneuroml libneuroml \n"
, " should fix it."
, " Actual error: %s " % e ]
)
@@ -229,27 +229,23 @@ def mergeChemModel(src, des):
# NML2 reader and writer function.
def mooseReadNML2( modelpath ):
- """Read NeuroML model (version 2).
-
+ """Read NeuroML model (version 2) and return reader object.
"""
global nml2Import_
- if nml2Import_:
- reader = _neuroml2.NML2Reader( )
- reader.read( modelpath )
- return reader
- else:
- mu.info( nml2ImportError_ )
- mu.warn( "Could not load NML2 support. Doing nothing" )
- return False
+ if not nml2Import_:
+ mu.warn( nml2ImportError_ )
+ raise RuntimeError( "Could not load NML2 support." )
+
+ reader = _neuroml2.NML2Reader( )
+ reader.read( modelpath )
+ return reader
def mooseWriteNML2( outfile ):
- mu.warn( "Writing to NML2 is not supported yet" )
+ raise NotImplementedError( "Writing to NML2 is not supported yet" )
################################################################
# Wrappers for global functions
################################################################
-
-
def pwe():
"""Print present working element. Convenience function for GENESIS
users. If you want to retrieve the element in stead of printing
diff --git a/moose-core/python/moose/neuroml/ChannelML.py b/moose-core/python/moose/neuroml/ChannelML.py
index 60a566fd..cc2c4443 100644
--- a/moose-core/python/moose/neuroml/ChannelML.py
+++ b/moose-core/python/moose/neuroml/ChannelML.py
@@ -15,14 +15,13 @@ readChannelML(...) / readSynapseML to load from an xml.etree xml element (could
from __future__ import print_function
from xml.etree import ElementTree as ET
-import string
-import os, sys
+import os
+import sys
import math
-
import moose
from moose.neuroml import utils
-from moose import utils as moose_utils
-from moose import print_utils as pu
+import moose.utils as mu
+import moose.print_utils as pu
class ChannelML():
@@ -497,7 +496,7 @@ def make_new_synapse(syn_name, postcomp, syn_name_full, nml_params):
## connect the SimpleSynHandler or the STDPSynHandler to the SynChan (double exp)
moose.connect( synhandler, 'activationOut', syn, 'activation' )
# mgblock connections if required
- childmgblock = moose_utils.get_child_Mstring(syn,'mgblockStr')
+ childmgblock = mu.get_child_Mstring(syn,'mgblockStr')
#### connect the post compartment to the synapse
if childmgblock.value=='True': # If NMDA synapse based on mgblock, connect to mgblock
mgblock = moose.Mg_block(syn.path+'/mgblock')
diff --git a/moose-core/python/moose/neuroml/MorphML.py b/moose-core/python/moose/neuroml/MorphML.py
index 0a5d0c15..d34a0c12 100644
--- a/moose-core/python/moose/neuroml/MorphML.py
+++ b/moose-core/python/moose/neuroml/MorphML.py
@@ -259,15 +259,18 @@ class MorphML():
if cableid in self.intFireCableIds:
mechanismname = self.intFireCableIds[cableid]
if mechanismname is not None: # this cableid is an intfire
- ## create LIF (subclass of Compartment) and set to default values
+ # create LIF (subclass of Compartment) and set to default values
moosecomp = moose.LIF(moosecomppath)
- moosechannel = moose.Neutral('/library/'+mechanismname)
+ mname = '/library/' + mechanismname
+ moosechannel = moose.element(mname) if moose.exists(mname) else moose.Neutral(mname)
+ # Mstring values are 'string'; make sure to convert them to
+ # float else it will seg-fault with python3+
moosechannelval = moose.Mstring(moosechannel.path+'/vReset')
- moosecomp.vReset = moosechannelval.value
+ moosecomp.vReset = float(moosechannelval.value)
moosechannelval = moose.Mstring(moosechannel.path+'/thresh')
- moosecomp.thresh = moosechannelval.value
+ moosecomp.thresh = float( moosechannelval.value )
moosechannelval = moose.Mstring(moosechannel.path+'/refracT')
- moosecomp.refractoryPeriod = moosechannelval.value
+ moosecomp.refractoryPeriod = eval(moosechannelval.value)
## refracG is currently not supported by moose.LIF
## when you implement it, check if refracG or g_refrac
## is a conductance density or a conductance, I think the former
diff --git a/moose-core/python/moose/neuroml/NetworkML.py b/moose-core/python/moose/neuroml/NetworkML.py
index d7f30692..14ae5f53 100644
--- a/moose-core/python/moose/neuroml/NetworkML.py
+++ b/moose-core/python/moose/neuroml/NetworkML.py
@@ -179,7 +179,7 @@ class NetworkML():
model_filenames = (cellname+'.xml', cellname+'.morph.xml')
success = False
for model_filename in model_filenames:
- model_path = find_first_file(model_filename,self.model_dir)
+ model_path = find_first_file(model_filename, self.model_dir)
if model_path is not None:
cellDict = mmlR.readMorphMLFromFile(model_path, self.params)
success = True
diff --git a/moose-core/python/moose/neuroml/NeuroML.py b/moose-core/python/moose/neuroml/NeuroML.py
index ceddf21d..5d9d8edf 100644
--- a/moose-core/python/moose/neuroml/NeuroML.py
+++ b/moose-core/python/moose/neuroml/NeuroML.py
@@ -1,8 +1,12 @@
# -*- coding: utf-8 -*-
-## Description: class NeuroML for loading NeuroML from single file into MOOSE
-## Version 1.0 by Aditya Gilra, NCBS, Bangalore, India, 2011 for serial MOOSE
-## Version 1.5 by Niraj Dudani, NCBS, Bangalore, India, 2012, ported to parallel MOOSE
-## Version 1.6 by Aditya Gilra, NCBS, Bangalore, India, 2012, further changes for parallel MOOSE
+from __future__ import print_function, division
+
+# Description: class NeuroML for loading NeuroML from single file into MOOSE
+# Version 1.0 by Aditya Gilra, NCBS, Bangalore, India, 2011 for serial MOOSE
+# Version 1.5 by Niraj Dudani, NCBS, Bangalore, India, 2012, ported to parallel MOOSE
+# Version 1.6 by Aditya Gilra, NCBS, Bangalore, India, 2012, further changes for parallel MOOSE
+# Maintainer: Dilawar Singh, dilawars@ncbs.res.in
+# - python3 related fixes.
"""
NeuroML.py is the preferred interface to read NeuroML files.
@@ -15,7 +19,8 @@ readNeuroMLFromFile(...) to load NeuroML from a file:
OR
-(b) the file could have only L3 (network) with L2 (channels/synapses) and L1 (cells) spread over multiple files;
+(b) the file could have only L3 (network) with L2 (channels/synapses) and L1
+(cells) spread over multiple files;
these multiple files should be in the same or lower-level directory
named as <chan/syn_name>.xml or <cell_name>.xml or <cell_name>.morph.xml
(essentially as generated by neuroConstruct's export).
@@ -44,7 +49,6 @@ In [2]: import moose.neuroml
In [3]: moose.neuroml.loadNeuroML_L123('Generated.net.xml')
"""
-from __future__ import print_function
import moose
from moose.utils import *
from xml.etree import ElementTree as ET
@@ -57,32 +61,27 @@ from moose.neuroml.utils import *
import sys
from os import path
-import logging
-console = logging.StreamHandler()
-console.setLevel(logging.INFO)
-formatter = logging.Formatter('%(name)-12s: %(levelname)-8s %(message)s')
-console.setFormatter(formatter)
-_logger = logging.getLogger('moose.nml.neuroml')
-_logger.addHandler(console)
-
-
class NeuroML():
def __init__(self):
pass
- def readNeuroMLFromFile(self,filename,params={},cellsDict={}):
+ def readNeuroMLFromFile(self, filename, params={}, cellsDict={}):
"""
For the format of params required to tweak what cells are loaded,
refer to the doc string of NetworkML.readNetworkMLFromFile().
Returns (populationDict,projectionDict),
see doc string of NetworkML.readNetworkML() for details.
"""
- _logger.info("Loading neuroml file %s " % filename)
+ mu.info("Loading neuroml file %s " % filename)
moose.Neutral('/library') # creates /library in MOOSE tree; elif present, wraps
tree = ET.parse(filename)
root_element = tree.getroot()
- self.model_dir = path.dirname( path.abspath( filename ) )
+
+ # if model_path is given in params, use it else use the directory of NML
+ # as model_dir.
+ self.model_dir = params.get('model_dir', path.dirname(path.abspath(filename)))
+
if 'lengthUnits' in list(root_element.attrib.keys()):
self.lengthUnits = root_element.attrib['lengthUnits']
else:
@@ -110,14 +109,14 @@ class NeuroML():
self.temperature = float(tag.find('.//{'+meta_ns+'}value').text)
self.temperature_default = False
if self.temperature_default:
- _logger.info("Using default temperature of %s degree Celsius" % self.temperature)
+ mu.info("Using default temperature of %s degree Celsius" % self.temperature)
self.nml_params = {
'temperature':self.temperature,
'model_dir':self.model_dir,
}
- _logger.debug("Loading channels and synapses into MOOSE /library ...")
+ mu.debug("Loading channels and synapses into MOOSE /library ...")
cmlR = ChannelML(self.nml_params)
for channels in root_element.findall('.//{'+neuroml_ns+'}channels'):
self.channelUnits = channels.attrib['units']
@@ -131,7 +130,7 @@ class NeuroML():
for ionConc in channels.findall('.//{'+cml_ns+'}ion_concentration'):
cmlR.readIonConcML(ionConc,units=self.channelUnits)
- _logger.debug("Loading cell definitions into MOOSE /library ...")
+ mu.debug("Loading cell definitions into MOOSE /library ...")
mmlR = MorphML(self.nml_params)
self.cellsDict = cellsDict
for cells in root_element.findall('.//{'+neuroml_ns+'}cells'):
@@ -145,7 +144,7 @@ class NeuroML():
and root_element.find('.//{'+nml_ns+'}populations') is None:
return (self.cellsDict,'no populations (L3 NetworkML) found.')
else:
- _logger.debug("Loading individual cells into MOOSE root ... ")
+ mu.debug("Loading individual cells into MOOSE root ... ")
nmlR = NetworkML(self.nml_params)
return nmlR.readNetworkML(root_element,self.cellsDict,\
params=params,lengthUnits=self.lengthUnits)
@@ -163,6 +162,6 @@ def loadNeuroML_L123(filename):
if __name__ == "__main__":
if len(sys.argv)<2:
- _logger.error("You need to specify the neuroml filename.")
+ mu.error("You need to specify the neuroml filename.")
sys.exit(1)
print(loadNeuroML_L123(sys.argv[1]))
diff --git a/moose-core/python/moose/neuroml2/reader.py b/moose-core/python/moose/neuroml2/reader.py
index 1fdf4fc1..75c9a7f2 100644
--- a/moose-core/python/moose/neuroml2/reader.py
+++ b/moose-core/python/moose/neuroml2/reader.py
@@ -9,48 +9,11 @@
# Version:
# Last-Updated: 15 Jan 2018, pgleeson
# 16 Jan 2018, dilawar, python3 compatible imports.
-#
-# 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:
-"""Implementation of reader for NeuroML 2 models.
-
+#
+"""Implementation of reader for NeuroML 2 models.
TODO: handle morphologies of more than one segment...
-
"""
-
from __future__ import print_function, division, absolute_import
try:
@@ -64,13 +27,19 @@ import math
import numpy as np
import neuroml as nml
from pyneuroml import pynml
-
import moose
import moose.utils as mu
from .units import SI
from . import hhfit
+def _unique( ls ):
+ res = [ ]
+ for l in ls:
+ if l not in res:
+ res.append( l )
+ return res
+
def sarea(comp):
"""
Return the surface area of compartment from length and
@@ -116,15 +85,13 @@ def setEk(comp, erev):
def getSegments(nmlcell, component, sg_to_segments):
"""Get the list of segments the `component` is applied to"""
sg = component.segment_groups
- #seg = component.segment
if sg is None:
- segments = nmlcell.morphology.segments
- elif sg == 'all': # Special case
+ segments = nmlcell.morphology.segments
+ elif sg == 'all':
segments = [seg for seglist in sg_to_segments.values() for seg in seglist]
else:
segments = sg_to_segments[sg]
-
- return list(set(segments))
+ return _unique(segments)
class NML2Reader(object):
@@ -148,10 +115,10 @@ class NML2Reader(object):
self.filename = None
self.nml_to_moose = {} # NeuroML object to MOOSE object
self.moose_to_nml = {} # Moose object to NeuroML object
- self.proto_cells = {} # map id to prototype cell in moose
- self.proto_chans = {} # map id to prototype channels in moose
- self.proto_pools = {} # map id to prototype pools (Ca2+, Mg2+)
- self.includes = {} # Included files mapped to other readers
+ self.proto_cells = {} # map id to prototype cell in moose
+ self.proto_chans = {} # map id to prototype channels in moose
+ self.proto_pools = {} # map id to prototype pools (Ca2+, Mg2+)
+ self.includes = {} # Included files mapped to other readers
self.lib = moose.Neutral('/library')
self.id_to_ionChannel = {}
self._cell_to_sg = {} # nml cell to dict - the dict maps segment groups to segments
@@ -198,11 +165,14 @@ class NML2Reader(object):
return self.cells_in_populations[pop_id][index]
def getComp(self, pop_id, cellIndex, segId):
- return moose.element('%s/%s/%s/%s' % (self.lib.path, pop_id, cellIndex, self.seg_id_to_comp_name[self.pop_to_cell_type[pop_id]][segId]))
+ return moose.element('%s/%s/%s/%s' % ( self.lib.path, pop_id, cellIndex
+ , self.seg_id_to_comp_name[self.pop_to_cell_type[pop_id]][segId])
+ )
def createPopulations(self):
for pop in self.network.populations:
- mpop = moose.Neutral('%s/%s' % (self.lib.path, pop.id))
+ ep = '%s/%s' % (self.lib.path, pop.id)
+ mpop = moose.element(ep) if moose.exists(ep) else moose.Neutral(ep)
self.cells_in_populations[pop.id] ={}
for i in range(pop.size):
mu.info("Creating %s/%s instances of %s under %s"%(i,pop.size,pop.component, mpop))
@@ -228,12 +198,15 @@ class NML2Reader(object):
for il in self.network.input_lists:
for ii in il.input:
input = self.getInput(il.component)
- moose.connect(input, 'output', self.getComp(il.populations,ii.get_target_cell_id(),ii.get_segment_id()), 'injectMsg')
+ moose.connect(input, 'output'
+ , self.getComp(il.populations,ii.get_target_cell__hash(),ii.get_segment__hash())
+ , 'injectMsg')
def createCellPrototype(self, cell, symmetric=True):
"""To be completed - create the morphology, channels in prototype"""
- nrn = moose.Neuron('%s/%s' % (self.lib.path, cell.id))
+ ep = '%s/%s' % (self.lib.path, cell.id)
+ nrn = moose.element(ep) if moose.exists(ep) else moose.Neuron(ep)
self.proto_cells[cell.id] = nrn
self.nml_to_moose[cell] = nrn
self.moose_to_nml[nrn] = cell
@@ -262,11 +235,13 @@ class NML2Reader(object):
self.seg_id_to_comp_name[nmlcell.id]={}
for seg in segments:
if seg.name is not None:
- id_to_comp[seg.id] = compclass('%s/%s' % (cellpath, seg.name))
+ ep = '%s/%s' % (cellpath, seg.name)
+ id_to_comp[seg.id] = moose.element(ep) if moose.exists(ep) else compclass(ep)
self.seg_id_to_comp_name[nmlcell.id][seg.id] = seg.name
else:
name = 'comp_%s'%seg.id
- id_to_comp[seg.id] = compclass('%s/%s' % (cellpath, name))
+ ep = '%s/%s' % (cellpath, name)
+ id_to_comp[seg.id] = moose.element(ep) if moose.exists(ep) else compclass(ep)
self.seg_id_to_comp_name[nmlcell.id][seg.id] = name
# Now assign the positions and connect up axial resistance
if not symmetric:
@@ -286,13 +261,15 @@ class NML2Reader(object):
if parent:
p0 = parent.distal
else:
- raise Exception('No proximal point and no parent segment for segment: name=%s, id=%s' % (segment.name, segment.id))
- comp.x0, comp.y0, comp.z0 = (x * self.lunit for x in map(float, (p0.x, p0.y, p0.z)))
+ raise Exception(
+ 'No proximal point and no parent segment for segment:'+\
+ 'name=%s, id=%s' % (segment.name, segment.id)
+ )
+ comp.x0, comp.y0, comp.z0 = (x*self.lunit for x in map(float, (p0.x, p0.y, p0.z)))
p1 = segment.distal
- comp.x, comp.y, comp.z = (x * self.lunit for x in map(float, (p1.x, p1.y, p1.z)))
- comp.length = np.sqrt((comp.x - comp.x0)**2
- + (comp.y - comp.y0)**2
- + (comp.z - comp.z0)**2)
+ comp.x, comp.y, comp.z = (x*self.lunit for x in map(float, (p1.x, p1.y, p1.z)))
+ comp.length = np.sqrt((comp.x-comp.x0)**2+(comp.y-comp.y0)**2+(comp.z-comp.z0)**2)
+
# This can pose problem with moose where both ends of
# compartment have same diameter. We are averaging the two
# - may be splitting the compartment into two is better?
@@ -377,10 +354,15 @@ class NML2Reader(object):
proto_pool = innerReader.proto_pools[species.concentrationModel]
break
if not proto_pool:
- raise Exception('No prototype pool for %s referred to by %s' % (species.concentration_model, species.id))
+ raise Exception('No prototype pool for %s referred to by %s' % (
+ species.concentration_model, species.id)
+ )
pool_id = moose.copy(proto_pool, comp, species.id)
pool = moose.element(pool_id)
- pool.B = pool.B / (np.pi * compartment.length * (0.5 * compartment.diameter + pool.thickness) * (0.5 * compartment.diameter - pool.thickness))
+ pool.B = pool.B / (np.pi * compartment.length * (
+ 0.5 * compartment.diameter + pool.thickness) *
+ (0.5 * compartment.diameter - pool.thickness)
+ )
return pool
def importAxialResistance(self, nmlcell, intracellularProperties):
@@ -417,8 +399,11 @@ class NML2Reader(object):
mu.info("Using %s to evaluate rate"%ct.name)
rate = []
for v in tab:
- vals = pynml.evaluate_component(ct,req_variables={'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()})
- '''mu.info vals'''
+ vals = pynml.evaluate_component(ct
+ , req_variables =
+ {'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()}
+ )
+ # mu.info vals
if 'x' in vals:
rate.append(vals['x'])
if 't' in vals:
@@ -483,32 +468,6 @@ class NML2Reader(object):
moose.connect(chan, 'channel', comp, 'channel')
return chan
- '''
- def importIncludes(self, doc):
- for include in doc.include:
- if self.verbose:
- mu.info(self.filename, 'Loading include', include)
- error = None
- inner = NML2Reader(self.verbose)
- paths = [include.href, os.path.join(os.path.dirname(self.filename), include.href)]
- for path in paths:
- try:
- inner.read(path)
- if self.verbose:
- mu.info(self.filename, 'Loaded', path, '... OK')
- except IOError as e:
- error = e
- else:
- self.includes[include.href] = inner
- self.id_to_ionChannel.update(inner.id_to_ionChannel)
- self.nml_to_moose.update(inner.nml_to_moose)
- self.moose_to_nml.update(inner.moose_to_nml)
- error = None
- break
- if error:
- mu.info(self.filename, 'Last exception:', error)
- raise IOError('Could not read any of the locations: %s' % (paths))'''
-
def _is_standard_nml_rate(self, rate):
return rate.type=='HHExpLinearRate' \
or rate.type=='HHExpRate' or \
@@ -517,7 +476,7 @@ class NML2Reader(object):
def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000):
mchan = moose.HHChannel('%s/%s' % (self.lib.path, chan.id))
- mgates = map(moose.element, (mchan.gateX, mchan.gateY, mchan.gateZ))
+ mgates = [moose.element(x) for x in [mchan.gateX, mchan.gateY, mchan.gateZ]]
assert(len(chan.gate_hh_rates) <= 3) # We handle only up to 3 gates in HHCHannel
if self.verbose:
@@ -551,19 +510,25 @@ class NML2Reader(object):
if ngate.q10_settings.type == 'q10Fixed':
q10_scale= float(ngate.q10_settings.fixed_q10)
elif ngate.q10_settings.type == 'q10ExpTemp':
- q10_scale = math.pow(float(ngate.q10_settings.q10_factor),(self._getTemperature()- SI(ngate.q10_settings.experimental_temp))/10)
- #mu.info('Q10: %s; %s; %s; %s'%(ngate.q10_settings.q10_factor, self._getTemperature(),SI(ngate.q10_settings.experimental_temp),q10_scale))
+ q10_scale = math.pow( float(ngate.q10_settings.q10_factor)
+ , (self._getTemperature()- SI(ngate.q10_settings.experimental_temp))/10
+ )
else:
- raise Exception('Unknown Q10 scaling type %s: %s'%(ngate.q10_settings.type,ngate.q10_settings))
+ raise Exception('Unknown Q10 scaling type %s: %s'%(
+ ngate.q10_settings.type,ngate.q10_settings)
+ )
if self.verbose:
- mu.info(' === Gate: %s; %s; %s; %s; %s; scale=%s'%(ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale))
+ mu.info(' === Gate: %s; %s; %s; %s; %s; scale=%s'% (
+ ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale)
+ )
if (fwd is not None) and (rev is not None):
alpha = self.calculateRateFn(fwd, vmin, vmax, vdivs)
beta = self.calculateRateFn(rev, vmin, vmax, vdivs)
mgate.tableA = q10_scale * (alpha)
mgate.tableB = q10_scale * (alpha + beta)
+
# Assuming the meaning of the elements in GateHHTauInf ...
if hasattr(ngate,'time_course') and hasattr(ngate,'steady_state') \
and (ngate.time_course is not None) and (ngate.steady_state is not None):
@@ -587,16 +552,25 @@ class NML2Reader(object):
return mchan
def createPassiveChannel(self, chan):
- mchan = moose.Leakage('%s/%s' % (self.lib.path, chan.id))
+ epath = '%s/%s' % (self.lib.path, chan.id)
+ if moose.exists( epath ):
+ mchan = moose.element(epath)
+ else:
+ mchan = moose.Leakage(epath)
if self.verbose:
mu.info(self.filename, 'Created', mchan.path, 'for', chan.id)
return mchan
def importInputs(self, doc):
- minputs = moose.Neutral('%s/inputs' % (self.lib.path))
- for pg_nml in doc.pulse_generators:
+ epath = '%s/inputs' % (self.lib.path)
+ if moose.exists( epath ):
+ minputs = moose.element( epath )
+ else:
+ minputs = moose.Neutral( epath )
- pg = moose.PulseGen('%s/%s' % (minputs.path, pg_nml.id))
+ for pg_nml in doc.pulse_generators:
+ epath = '%s/%s' % (minputs.path, pg_nml.id)
+ pg = moose.element(epath) if moose.exists(epath) else moose.PulseGen(epath)
pg.firstDelay = SI(pg_nml.delay)
pg.firstWidth = SI(pg_nml.duration)
pg.firstLevel = SI(pg_nml.amplitude)
@@ -639,14 +613,10 @@ class NML2Reader(object):
ca.CaBasal = SI(concModel.restingConc)
ca.tau = SI(concModel.decayConstant)
ca.thick = SI(concModel.shellThickness)
- ca.B = 5.2e-6 # B = 5.2e-6/(Ad) where A is the area of the shell and d is thickness - must divide by shell volume when copying
+ ca.B = 5.2e-6 # B = 5.2e-6/(Ad) where A is the area of the
+ # shell and d is thickness - must divide by
+ # shell volume when copying
self.proto_pools[concModel.id] = ca
self.nml_to_moose[concModel.id] = ca
self.moose_to_nml[ca] = concModel
logger.debug('Created moose element: %s for nml conc %s' % (ca.path, concModel.id))
-
-
-
-
-#
-# reader.py ends here
diff --git a/moose-core/python/setup.cmake.py b/moose-core/python/setup.cmake.py
index e6bbd1a7..c7254b40 100644
--- a/moose-core/python/setup.cmake.py
+++ b/moose-core/python/setup.cmake.py
@@ -32,11 +32,12 @@ with open( os.path.join( script_dir, 'VERSION'), 'r' ) as f:
version = f.read( )
print( 'Got %s from VERSION file' % version )
-# importlib is available only for python3.
+# importlib is available only for python3. Since we build wheels, prefer .so
+# extension. This way a wheel built by any python3.x will work with any python3.
suffix = '.so'
try:
import importlib.machinery
- suffix = importlib.machinery.EXTENSION_SUFFIXES[0]
+ suffix = importlib.machinery.EXTENSION_SUFFIXES[-1]
except Exception as e:
print( '[WARN] Failed to determine importlib suffix' )
suffix = '.so'
@@ -54,9 +55,9 @@ setup(
packages = [ 'rdesigneur', 'moose'
, 'moose.SBML', 'moose.genesis'
, 'moose.neuroml', 'moose.neuroml2'
- , 'moose.chemUtil', 'moose.chemMerge'
+ , 'moose.chemUtil', 'moose.chemMerge'
],
- install_requires = [ 'python-libsbml', 'numpy' ],
+ install_requires = [ 'numpy' ],
package_dir = { 'moose' : 'moose', 'rdesigneur' : 'rdesigneur' },
package_data = { 'moose' : ['_moose' + suffix, 'neuroml2/schema/NeuroMLCoreDimensions.xml'] },
)
diff --git a/moose-core/tests/python/test_Xreacs2.py b/moose-core/tests/python/test_Xreacs2.py
index 4089501c..2fb36cfb 100644
--- a/moose-core/tests/python/test_Xreacs2.py
+++ b/moose-core/tests/python/test_Xreacs2.py
@@ -7,7 +7,8 @@ import moose.fixXreacs as fixXreacs
def countCrossings( plot, thresh ):
vec = moose.element( plot ).vector
- #print (vec[:-1] < thresh)
+ print( vec )
+ # print (vec[:-1] <= thresh)
return sum( (vec[:-1] < thresh) * (vec[1:] >= thresh ) )
def main( standalone = False ):
@@ -43,7 +44,8 @@ def main( standalone = False ):
moose.reinit()
moose.start( runtime )
# I don't have an analytic way to assess oscillations
- assert( countCrossings( '/model/graphs/conc2/M.Co', 0.001 ) == 4 )
+ nCrossings = countCrossings( '/model/graphs/conc2/M.Co', 0.001 )
+ assert( nCrossings == 4 ), "Expected 4, got %d" % nCrossings
moose.delete( '/model' )
# Run the 'main' if this script is executed standalone.
diff --git a/moose-core/tests/python/test_dose_response.py b/moose-core/tests/python/test_dose_response.py
new file mode 100644
index 00000000..4a897301
--- /dev/null
+++ b/moose-core/tests/python/test_dose_response.py
@@ -0,0 +1,116 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function, division
+
+## 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__ ) )
+vals_ = [ ]
+
+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.
+ a = moose.element( '/model/kinetics/a' ).conc
+ vals_.append( a )
+ 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)
+
+ expected = (0.001411, 0.00092559)
+ got = ( np.mean(vals_), np.std(vals_) )
+ assert np.isclose(got, expected, atol=1e-4).all(), "Got %s, expected %s" % (got, expected)
+ print( '[TEST1] Passed. Concentration of a is same' )
+
+ 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-core/tests/python/test_neuroml2.py b/moose-core/tests/python/test_neuroml2.py
index 8460c092..541d7580 100644
--- a/moose-core/tests/python/test_neuroml2.py
+++ b/moose-core/tests/python/test_neuroml2.py
@@ -14,6 +14,18 @@ import numpy as np
SCRIPT_DIR = os.path.dirname( os.path.realpath( __file__ ) )
+# check if neuroml working properly.
+# NOTE: This script does not work with python3
+# See https://github.com/NeuroML/NeuroML2/issues/116 . If this bug is fixed then
+# remove this code block.
+import neuroml as nml
+a = nml.nml.nml.IonChannel()
+try:
+ b = {a : 1 }
+except TypeError as e:
+ print( 'Failed due to https://github.com/NeuroML/NeuroML2/issues/116' )
+ quit( 0 )
+
def run( nogui = True ):
global SCRIPT_DIR
filename = os.path.join(SCRIPT_DIR, 'test_files/passiveCell.nml' )
diff --git a/moose-core/tests/python/test_steady_state_solver.py b/moose-core/tests/python/test_steady_state_solver.py
new file mode 100644
index 00000000..fd3e13e1
--- /dev/null
+++ b/moose-core/tests/python/test_steady_state_solver.py
@@ -0,0 +1,152 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function, division
+
+# This program is part of 'MOOSE', the
+# Messaging Object Oriented Simulation Environment.
+# Copyright (C) 2013 Upinder S. Bhalla. and NCBS
+# It is made available under the terms of the
+# GNU Lesser General Public License version 2.1
+# See the file COPYING.LIB for the full notice.
+# Monday 17 September 2018 01:49:30 PM IST
+# - Converted to a test script
+
+import math
+import numpy as np
+import moose
+print( "[INFO ] Using moose from %s" % moose.__file__ )
+
+def main():
+ compartment = makeModel()
+ ksolve = moose.Ksolve( '/model/compartment/ksolve' )
+ stoich = moose.Stoich( '/model/compartment/stoich' )
+ stoich.compartment = compartment
+ stoich.ksolve = ksolve
+ stoich.path = "/model/compartment/##"
+ state = moose.SteadyState( '/model/compartment/state' )
+
+ moose.reinit()
+ state.stoich = stoich
+ state.convergenceCriterion = 1e-6
+ moose.seed( 111 ) # Used when generating the samples in state space
+
+ b = moose.element( '/model/compartment/b' )
+ a = moose.element( '/model/compartment/a' )
+ c = moose.element( '/model/compartment/c' )
+ a.concInit = 0.1
+ deltaA = 0.002
+ num = 150
+ avec = []
+ bvec = []
+ moose.reinit()
+
+ # Now go up.
+ for i in range( 0, num ):
+ moose.start( 1.0 ) # Run the model for 1 seconds.
+ state.settle() # This function finds the steady states.
+ avec.append( a.conc )
+ bvec.append( b.conc )
+ a.concInit += deltaA
+
+ aa, bb = avec, bvec
+ got = np.mean(aa), np.std(aa)
+ expected = 0.24899, 0.08660
+ assert np.isclose(got, expected, atol = 1e-4).all(), "Got %s, expected %s" % (got, expected)
+ print( "[INFO ] Test 1 PASSED" )
+
+
+ # Now go down.
+ avec = []
+ bvec = []
+ for i in range( 0, num ):
+ moose.start( 1.0 ) # Run the model for 1 seconds.
+ state.settle() # This function finds the steady states.
+ avec.append( a.conc )
+ bvec.append( b.conc )
+ a.concInit -= deltaA
+
+ aa, bb = avec, bvec
+ got = np.mean(aa), np.std(aa)
+ expected = 0.251, 0.0866
+ assert np.isclose(got, expected, atol = 1e-4).all(), "Got %s, expected %s" % (got, expected)
+ print( "[INFO ] Test 2 PASSED" )
+
+ # Now aim for the middle. We do this by judiciously choosing a
+ # start point that should be closer to the unstable fixed point.
+ avec = []
+ bvec = []
+ a.concInit = 0.28
+ b.conc = 0.15
+ for i in range( 0, 65 ):
+ moose.start( 1.0 ) # Run the model for 1 seconds.
+ state.settle() # This function finds the steady states.
+ avec.append( a.conc )
+ bvec.append( b.conc )
+ a.concInit -= deltaA
+
+ aa, bb = avec, bvec
+ got = np.mean(aa), np.std(aa)
+ expected = 0.216, 0.03752
+ assert np.isclose(got, expected, atol = 1e-4).all(), "Got %s, expected %s" % (got, expected)
+ print( "[INFO ] Test 3 PASSED" )
+ quit()
+
+def makeModel():
+ """ This function creates a bistable reaction system using explicit
+ MOOSE calls rather than load from a file.
+ The reaction is::
+
+ a ---b---> 2b # b catalyzes a to form more of b.
+ 2b ---c---> a # c catalyzes b to form a.
+ a <======> 2b # a interconverts to b.
+
+ """
+ # create container for model
+ model = moose.Neutral( 'model' )
+ compartment = moose.CubeMesh( '/model/compartment' )
+ compartment.volume = 1e-15
+ # the mesh is created automatically by the compartment
+ mesh = moose.element( '/model/compartment/mesh' )
+
+ # create molecules and reactions
+ a = moose.BufPool( '/model/compartment/a' )
+ b = moose.Pool( '/model/compartment/b' )
+ c = moose.Pool( '/model/compartment/c' )
+ enz1 = moose.Enz( '/model/compartment/b/enz1' )
+ enz2 = moose.Enz( '/model/compartment/c/enz2' )
+ cplx1 = moose.Pool( '/model/compartment/b/enz1/cplx' )
+ cplx2 = moose.Pool( '/model/compartment/c/enz2/cplx' )
+ reac = moose.Reac( '/model/compartment/reac' )
+
+ # connect them up for reactions
+ moose.connect( enz1, 'sub', a, 'reac' )
+ moose.connect( enz1, 'prd', b, 'reac' )
+ moose.connect( enz1, 'prd', b, 'reac' ) # Note 2 molecules of b.
+ moose.connect( enz1, 'enz', b, 'reac' )
+ moose.connect( enz1, 'cplx', cplx1, 'reac' )
+
+ moose.connect( enz2, 'sub', b, 'reac' )
+ moose.connect( enz2, 'sub', b, 'reac' ) # Note 2 molecules of b.
+ moose.connect( enz2, 'prd', a, 'reac' )
+ moose.connect( enz2, 'enz', c, 'reac' )
+ moose.connect( enz2, 'cplx', cplx2, 'reac' )
+
+ moose.connect( reac, 'sub', a, 'reac' )
+ moose.connect( reac, 'prd', b, 'reac' )
+ moose.connect( reac, 'prd', b, 'reac' ) # Note 2 order in b.
+
+ # Assign parameters
+ a.concInit = 1
+ b.concInit = 0
+ c.concInit = 0.01
+ enz1.kcat = 0.4
+ enz1.Km = 4
+ enz2.kcat = 0.6
+ enz2.Km = 0.01
+ reac.Kf = 0.001
+ reac.Kb = 0.01
+
+ return compartment
+
+# Run the 'main' if this script is executed standalone.
+if __name__ == '__main__':
+ main()
diff --git a/scripts/pull_subtree.sh b/scripts/pull_subtree.sh
index 71d70c29..4836384b 100755
--- a/scripts/pull_subtree.sh
+++ b/scripts/pull_subtree.sh
@@ -24,11 +24,11 @@ set -x
if [[ `pwd` == *"/moose" ]]; then
git subtree pull --prefix moose-core \
- https://github.com/BhallaLab/moose-core chhennapoda --squash
+ https://github.com/BhallaLab/moose-core master --squash
git subtree pull --prefix moose-examples \
- https://github.com/BhallaLab/moose-examples chhennapoda --squash
+ https://github.com/BhallaLab/moose-examples master --squash
git subtree pull --prefix moose-gui \
- https://github.com/BhallaLab/moose-gui chhennapoda --squash
+ https://github.com/BhallaLab/moose-gui master --squash
else
echo "Run this script from top-level git directory."
--
GitLab