From ec5baf0e268089f2cc1e7f1f088214b5ceb7685f Mon Sep 17 00:00:00 2001
From: Dilawar Singh <dilawars@ncbs.res.in>
Date: Fri, 25 May 2018 10:38:05 +0530
Subject: [PATCH] Squashed 'moose-core/' changes from 0cceb296d0..cddc3bd16d
cddc3bd16d Merge pull request #265 from upibhalla/master
5856ac1fa7 Fixes for travis.yml
219f474a2e Fix to HSolve to prevent segv during destruction. Update to rdesigneur to not create HSolve when turnOffElec flag is set
9ca16e8246 Merge branch 'master' of https://github.com/BhallaLab/moose-core
58e630b4dc Updates to the Voltage clamp modules. The key thing to note is to use the standalone VClamp object in preference to individually building the circuit with DiffAmp. RC and PID.
12cf83e17a Updates to Hsolve to be able to delete and recreate it. Still problematic for HHChannels. Also added Vclamp to the rdesigneurProtos.
3dc99b53e6 Changes made to object name which had number in starting
3ec52cd9e4 Small bugfix to rdesigneur
241c3d59cb Use sys.executable to execute test. It breaks on python3.
f527490ff0 Fixes to rdesigneur to use the new system for cross-compartment reactions
501ff45939 Major cleanup to get rid of legacy cross-compartment code
f4b4772f0a Numerical improvment for Dsolve::calcJunction done by splitting diffusion and channel calculations into two dt/2 segments bracketing the mass-exchange calculation of calcJnXfer.
35db4b5a31 Major changes to handling of cross-compartment reactions and to the flow of control in Ksolve, during chemical reac-diff simulations. Not yet optimized. Several tests cleared.
1b3b11504b Install pymoose in debug mode as well.
6d911de9fa Test install target with python2 on travis.
ba0bab2206 Fixes to installation prefix (reported by Harsha over email). Chamcham branch patch was not correctly applied.
1a57a54c5b Fixes for Debian.
fab88be5d8 Fix to ConcChan lookup logic
efcb412ff3 Fixes so that CubeMesh can be used with diffusion solver, needed now that Dsolve is required for multi-compartment reactions.
f54acebcb9 PyMOOSE installation directory can work.
fe3061b0a5 temp changes. [skip ci]
b4313cdcca Merge commit 'a039fc3afb0948e09befcb74d8d9c2b057416ee7'
1dc7a63356 Added numProducts field to EnzBase
0f1ecebcb0 Major update so that the diffusion system handles channels across membranes, which are implemented as ConcChans.
53b0e60dc3 Merge branch 'chennapoda' of github.com:BhallaLab/moose
db52562e0b Merge branch 'chennapoda'
f208e2b6fb Merge commit '33bf735fb61ae5a94a2a531fa14a03632ee7ca42'
ecce024d42 Merge commit 'a1f7f5ea38d23caa2be93d3882736a04691ce482'
643226e2cf Install in debian layout in packages. Added option in setup.cfg file Its delicate solution. Added a note.
c2a574289b Finally. Needs to write setup.cfg file to make /usr/local and /usr prefix consistent.
ac7e7e5cf4 bdist_egg is more consistent on both RPM and DEB based systems.
c866f9f2d9 Merge commit 'bcf5c52d54409a89f9af1cf08e5d369d9e9e0d7f' into chennapoda
aacb0300ee Fixed typo in command.
d773b400d0 Merge commit 'a1a6faa7ffb7c176d85d8843a812ef006bf8afc8' into chennapoda
1e23289141 multiple fix to build system. - dont use python setup.py install . It does not work anymore if --prefix directory is not in PYTHONPATH. - Use python setup.py bdist -d INSTALL_DIR and let cmake unarchive and copy it to desired location. This is much better solution.
f2537df577 Use bdist instead of install for creating binay distrution. During installation, untar this to destination.
43891c6900 Tweaks to setup.cmake.py file.
fa4a545df0 Test on travis.
16eea973b1 Merge commit '6d651b6013e1093a1cdb2541bb8615bbc3a66021'
ba1065cdb9 Added EndoMesh and updated other meshes to match. Updated diffusion code to do cross-compartment reactions.
545ee72491 Merge commit 'e0b8959ccefd77195b7eac2edc1db7cdcaa40fe6'
ef650c8ec0 On debian/ubuntu install-layout option is set to deb. [skip ci]
2a02c27dc4 Merge commit 'f80ae40d5b1e1b63eb9849b9699d76f33e4f2572'
f11b189779 This fixes moose but a better solution to install pymoose in moose-core must be found.
635dedf798 Merge commit '6f7697fb3256626a6d57bda24ae0a7bb5c0e90cc'
7687f3c446 Merge commit '627f1ce1cc73b0e1f6b4fe9051d3597e3a9113a5'
419d65dd2a Merge commit '7b3c355a04d22f279a88bfae4097e9cf3f62439b'
5cb77cc213 Squashed 'moose-core/' changes from d4417dcffd..2657e3345e (#222)
59bf6585ff Fixed the syntax causing failure.
448fb46d87 Handling MOOSE_VERSION is better way.
cc09546ecf Merge commit 'fba178451c55c2bf9f5d01b4a473f1819eefc856'
9bd9c55d3c Merge commit '750df4139010ce4edb9f34cc3896970894e9db1d'
c85e8fc03b Force --install-layout=deb when on debian/ubuntu.
cf20c6c086 Merge commit '99b934e97af8b699f81716254e6966788da102f2'
2bc6d025ed Merge commit '2f038f21bc94f16351c030dce68fb7bbafd1e62a'
8e8f75f08b Merge commit 'd0e1cee1d4e62b656772c77fc42e0cd7afad8f66'
9fd3829997 Merge commit '0014d8268ef8184db15745884208636c1f05e76a'
c39f533796 Follow up fixes to calcium-hsolve (#221)
c3204e010e Merge branch 'master' into asiaszmek
62b315f08a Merge commit 'e6e83b761e93923f3853512e1814a412296927a7' as 'moose-core'
3e88642a5f Calcium hsolve (#218)
a33619039f Merge commit '7fe1c2e7fe0b0d1cbe153cefb10618f62fe107cb'
4547a9f1d2 Merge commit 'af3c29a4cff12f8572527cd79095d50d16c5da84'
fbdaafb332 Fixes to allow moose module to be installed and imported (#220)
42fa11cac2 All changes from previous repository. No new file is added here.
9170767ac7 Merge commit '0f6f5d2769ab429a3bbb940c9add228981d9ab98'
d9819894f7 Added missing Revision macro file from chamcham branch.
5627ad198f Merge commit 'b4d20b7c10a5e94e4c10f8cce5f8da428cf1c7b6'
5ec49526f0 Fixed revision issue.
b3418608f7 Cmake changes from chamcham branch. Build for launchpad.
485c150c85 Do not copy VERSION file for python. Lets find a better and more portable solution.
40bfbf631e Updated default version in setup.py file.
baab7ef505 Merge commit '0ff1b77a6fa7152903e2a3f8e0fbe658fb314e5c'
79e1a65c0d VERSION string is fixed to today's date if not given.
f1adfd3ff9 Fixed for snappy builds.
68e138a565 Let cmake find path of Python.h file.
cf8faf32d8 Merge commit '4489d15994d2ab9ab2448afd937ba13c65a9907d'
91a4b30ab1 More changes to snap. Needs to add Python.h path [skip ci]
567c8013ce Removed <<<<, === and >>> from improper merge.
435188aa4b Added snapcraft file.
3131b29c76 VERSION conflict file is added
a3a81ace57 Merge commit '628138cb9426ec4c1dcf98d3ef8a5108d0ba783b'
4bbbf6ddff Merge commit '7c71305865364216fd1a00b4a2498da44f0cbb0b'
a205390e4e Script testing DifShells
0fcafe248d Merge commit 'a5122cfb6dd82f1a831f64d1fb1a29a7d1392334'
8b1c110945 Fixed docstring on example.
f87a802bdb Merge commit '5bceacfe2bebf128999e0565521b7e49ff2f6f25'
f6f5dd58cd Merge branch 'difshell' of https://github.com/asiaszmek/moose into asiaszmek-difshell
d9d0f6e382 Change order of operations to slightly improve accuracy
d7f378074f Small changes to make the code more alike genesis and other moose objects
4ce676215f Difshells and difbuffers initialize close to equilibrium
ebab1aee87 Add forgotten setting of helper variables to zero after calculations
04d7c5227e Add state in t(i-1)
6779a6e0a4 Fix wrong definitions
10f549bbfd Fix errors in volume calculation
f66bf285a8 Indentation
25feaf08f7 Add MMPump to scheduling and makefiles
c410524b90 Add MMPump to DifShell
2f31fa17cd Implementation of MMPump
f58668f93a Add ticks for DifBufferBase and DifShellBase
372b34f71a Fix bugs: DifShell now compiles and seems to give correct results
efd6ca0bd9 Fix a typo in zombify
4fa3f001bb Rewrite DifShell. Inherit from DifShellBase
82284b9b2c Write base class for DifShell
a3d45d71bd Indentation
e32896701b Implement exponential Euler for DifShell
8d0e9c3ace Move source messages to class. Correctly initialize Faraday's constant
0b2b708d21 Move source messages to class
6c2234bf82 Make DifBuffer look more like other moose classes and fix ticks
0de0b1ba06 Add DifBuffer to clocks
e1586d63a9 Add DifBuffer to compilation
bea1f0a18d Fix no ticks when initializing difshells warning
c986fd8172 Remove useless variable
65e67a0bd5 Reintroduce an old version of DifBuffer from moose
cede2afd61 Moved to cmake folder.
fff821e0d0 Merge commit '0d12ebb03bf7c3efce4b08bcbe40ed82137fa202'
ea023b6ae7 master branch.
340c310fea Fixes to BhallaLab/moose#204.
0b083f83d6 Merge commit '04170cb4d5a69b48041a12229181d7e4d2bf5afd'
bbd455a47f Updated.
ff2aec915b Fixed, python2 incompatible print function use.
725d2c9c47 Some updates to fix static hdf5 libraries.
b6d80ccd7e This project build pymoose, gui and its examples. TODO: ctest must check each example as well.
44b458500f Builds whole moose project.
bbedb768be Added some sanity tests for both moose and moogli.
32e2226897 Merge commit '1680bac6852381ac9e46c60bd501b80ad58a6fbf'
20b1010d24 Changes to cmake file so that we can find the correct static library.
527bc4d489 Fix a development warning in cmake.
30dda99116 Removing shebangs from the python scripts. These scripts were never called as standalone executables.
20641912cd A test commit to package hdf5-devel BhallaLab/moose#201.
077989b4d4 Fixed another typo in FindGSL. Library list was converted to string.
0ddc8bd15f Fixed a typo in FindGSL file.
ebe9a39888 FindGSL also looks for GSL_VERSION now.
07d432d580 Fixes to travis failure.
2c91dddfcc Updated CMake and FindGSL.
11c933e476 Fixes to extension bug.
fff7ce57ab Some fixes to static library dependencies [skip ci]
05856011a9 Cleanup in cmake files. Using bash scripts to generate subprojects. This way I can handle environment variables cleanly.
493c7c49b8 New hdf5 version generates static library with different name. Fix it.
c9f6b7767b This should pass on travis.
d3ac4b814d Merge commit 'd0a6060b975540dedc0f455f49ea4357da479eec'
4d37569e6c Merge commit 'b5a7c2e2a2d7a1a5250ad20f8c43ac18910b7382'
a39a84acb1 Merge commit '273124ec2f3f25e6af896ac0b3759bee2ee2e77a'
4c79628c60 Merge commit '8bfef2e62f38d0660e71e43dc30081fc9004f811'
49c322da4c Merge commit '6c8f20121ff0b38b0c47eb37aacef75ba5cb9acc'
f9e56dd7e3 Merge commit '96580978fd046d82b2545fdac8b72ea52bd4956c'
633b9e0503 A test case test_sbml_support.py which inturn calls to write moosewriteSBML which needs networkx which is added as a dependency for travis
f21e02bf7f Oops moment. Fixed the typo in compiler flag for apple-clang
3da3dec81b updated apple-clang flags in compiler. Fix to missing symbols from stl.
70d1a5d90a Merge commit '2af833a02ab1ad4d1e9caf12efe7705e73fb2c32'
24a742712e Moved docs to moose. Removed temp directory.
26f20a818e Merge commit '6b27e69e490e08191cf1f54bea6bf822a063e753' as 'moose-core'
git-subtree-dir: moose-core
git-subtree-split: cddc3bd16d3a8e9e3668bdeae524986b0f68b6d4
---
.travis.yml | 8 +-
.travis/travis_build_linux.sh | 2 +
CMakeLists.txt | 23 +-
biophysics/VClamp.cpp | 22 +-
device/DiffAmp.cpp | 1 +
device/RC.cpp | 25 +-
device/RC.h | 2 +-
diffusion/ConcChanInfo.h | 31 ++
diffusion/DiffJunction.h | 21 ++
diffusion/DiffPoolVec.cpp | 13 +-
diffusion/DiffPoolVec.h | 3 +
diffusion/Dsolve.cpp | 366 +++++++++++++++++---
diffusion/Dsolve.h | 29 +-
hsolve/HSolve.cpp | 37 +-
hsolve/HSolve.h | 2 +
hsolve/HSolveActive.cpp | 4 +-
hsolve/HSolveUtils.cpp | 5 +-
kinetics/CMakeLists.txt | 1 +
kinetics/ConcChan.cpp | 206 ++++++++++++
kinetics/ConcChan.h | 52 +++
kinetics/EnzBase.cpp | 15 +
kinetics/EnzBase.h | 1 +
kinetics/ReadKkit.cpp | 67 ++--
ksolve/Gsolve.cpp | 126 +------
ksolve/Gsolve.h | 1 -
ksolve/GssaVoxelPools.cpp | 28 --
ksolve/GssaVoxelPools.h | 11 -
ksolve/Ksolve.cpp | 152 +--------
ksolve/Ksolve.h | 19 --
ksolve/Stoich.cpp | 133 +-------
ksolve/Stoich.h | 30 --
ksolve/VoxelPools.cpp | 53 ---
ksolve/VoxelPools.h | 5 -
ksolve/ZombiePoolInterface.cpp | 302 +----------------
ksolve/ZombiePoolInterface.h | 38 +--
mesh/CMakeLists.txt | 1 +
mesh/ChemCompt.cpp | 26 +-
mesh/ChemCompt.h | 12 +
mesh/CubeMesh.cpp | 64 +++-
mesh/CylMesh.cpp | 61 ++--
mesh/EndoMesh.cpp | 464 ++++++++++++++++++++++++++
mesh/EndoMesh.h | 150 +++++++++
mesh/MeshEntry.h | 3 +-
mesh/NeuroMesh.cpp | 8 +
python/moose/moose_test.py | 5 +-
python/rdesigneur/fixXreacs.py | 193 +++++++++++
python/rdesigneur/rdesigneur.py | 25 +-
python/rdesigneur/rdesigneurProtos.py | 78 +++++
python/setup.cmake.py | 3 +
scheduling/Clock.cpp | 2 +
50 files changed, 1884 insertions(+), 1045 deletions(-)
create mode 100644 diffusion/ConcChanInfo.h
create mode 100644 kinetics/ConcChan.cpp
create mode 100644 kinetics/ConcChan.h
create mode 100644 mesh/EndoMesh.cpp
create mode 100644 mesh/EndoMesh.h
create mode 100644 python/rdesigneur/fixXreacs.py
diff --git a/.travis.yml b/.travis.yml
index 294e0896..b5f3f688 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -34,9 +34,7 @@ script:
- if type python3 > /dev/null; then python3 -m compileall -q . ; fi
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then ./.travis/travis_build_osx.sh; fi
- if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then ./.travis/travis_build_linux.sh; fi
-- set +e
-deploy:
- on: tags
- provider: script
- script: ./.travis/deploy_pypi.sh
+addons:
+ apt:
+ update: true
diff --git a/.travis/travis_build_linux.sh b/.travis/travis_build_linux.sh
index e475b2e0..3009af21 100755
--- a/.travis/travis_build_linux.sh
+++ b/.travis/travis_build_linux.sh
@@ -45,6 +45,8 @@ 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
+ $PYTHON2 -c 'import moose;print(moose.__file__);print(moose.version())'
)
(
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 349d580c..bd6aaf84 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -401,7 +401,7 @@ set(EXTRA_ARGS "--prefix ${CMAKE_INSTALL_PREFIX} ${DISTUTILS_EXTRA_ARGS}")
# On Debian/Ubuntu install using debian layout.
# NOTE: Also create setup.cfg file which setup prefix and install-layout
# suitable for DEBIAN systems.
-if(${_platform_desc} MATCHES ".*(Ubuntu|Debian).*")
+if(${_platform_desc} MATCHES ".*(Ubuntu|debian).*")
list(APPEND EXTRA_ARGS "--install-layout=deb")
file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/python/setup.cfg
"[install]\nprefix=/usr\ninstall-layout=deb"
@@ -415,7 +415,11 @@ find_package(PythonInterp REQUIRED)
if(NOT PYMOOSE_BDIST_DIR)
set(PYMOOSE_BDIST_DIR ${CMAKE_BINARY_DIR}/bdist)
endif( )
-set(PYMOOSE_BDIST_INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/pymoose_bdist_install)
+
+# If no one hsa set the PYMOOSE bdist installation directory, use default.
+if(NOT PYMOOSE_BDIST_INSTALL_DIR)
+ set(PYMOOSE_BDIST_INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/pymoose_bdist_install)
+endif()
file(MAKE_DIRECTORY ${PYMOOSE_BDIST_INSTALL_DIR})
# We need a custom name for bdist; same on all platform.
@@ -426,14 +430,14 @@ add_custom_target(bdist ALL DEPENDS ${PYMOOSE_BDIST_FILE} )
# Any command using setup.cmake.py must run in the same directory.
add_custom_command( OUTPUT ${PYMOOSE_BDIST_FILE}
- COMMAND ${PYTHON_EXECUTABLE} setup.cmake.py bdist_dumb -p ${_platform}
+ COMMAND ${PYTHON_EXECUTABLE} setup.cmake.py bdist_dumb -p ${_platform}
-d ${PYMOOSE_BDIST_DIR}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/python
COMMENT "bdist is saved to ${PYMOOSE_BDIST_DIR}"
VERBATIM
)
add_custom_command(TARGET bdist POST_BUILD
- COMMAND ${CMAKE_COMMAND} -E chdir ${PYMOOSE_BDIST_INSTALL_DIR} tar xvf ${PYMOOSE_BDIST_FILE}
+ COMMAND ${CMAKE_COMMAND} -E chdir ${PYMOOSE_BDIST_INSTALL_DIR} tar xf ${PYMOOSE_BDIST_FILE}
COMMENT "Unarchiving bdist file ${PYMOOSE_BDIST_FILE} to ${PYMOOSE_BDIST_INSTALL_DIR}"
VERBATIM
)
@@ -445,7 +449,7 @@ add_custom_target(copy_pymoose DEPENDS ${PYTHON_SRCS}
COMMAND ${CMAKE_COMMAND} -E copy_directory
${CMAKE_CURRENT_SOURCE_DIR}/python ${CMAKE_CURRENT_BINARY_DIR}/python
COMMENT "Copying required python files and other files to build directory"
- VERBATIM
+ VERBATIM
)
add_dependencies( bdist copy_pymoose )
@@ -455,14 +459,11 @@ add_dependencies( copy_pymoose _moose )
install(TARGETS moose.bin DESTINATION bin CONFIGURATIONS Debug)
install(TARGETS libmoose DESTINATION lib CONFIGURATIONS Debug)
-install(PROGRAMS ${CMAKE_CURRENT_SOURCE_DIR}/scripts/moose
- DESTINATION bin CONFIGURATIONS Debug
- )
-# install pymoose bdist.
-install(DIRECTORY ${PYMOOSE_BDIST_INSTALL_DIR}/
+# install pymoose bdist. The bdist comes with predefined /usr; remove it.
+install(DIRECTORY ${PYMOOSE_BDIST_INSTALL_DIR}/usr/
DESTINATION ${CMAKE_INSTALL_PREFIX}
- CONFIGURATIONS Release
+ CONFIGURATIONS Release Debug
)
# Print message to start build process
diff --git a/biophysics/VClamp.cpp b/biophysics/VClamp.cpp
index 9ed2be86..b28218d5 100644
--- a/biophysics/VClamp.cpp
+++ b/biophysics/VClamp.cpp
@@ -82,8 +82,9 @@ const Cinfo * VClamp::initCinfo()
"Shared message to receive Process messages from the scheduler",
processShared, sizeof(processShared) / sizeof(Finfo*));
- static ReadOnlyValueFinfo< VClamp, double> command("command",
+ static ValueFinfo< VClamp, double> command("command",
"Command input received by the clamp circuit.",
+ &VClamp::setCommand,
&VClamp::getCommand);
static ValueFinfo< VClamp, unsigned int> mode("mode",
"Working mode of the PID controller.\n"
@@ -210,14 +211,12 @@ VClamp::~VClamp()
void VClamp::setCommand(double value)
{
- // e2_ = 0;
- // e1_ = 0;
cmdIn_ = value;
}
double VClamp::getCommand() const
{
- return command_;
+ return cmdIn_;
}
void VClamp::setTi(double value)
@@ -334,13 +333,16 @@ void VClamp::reinit(const Eref& e, ProcPtr p)
tauByDt_ = tau_ / p->dt;
dtByTi_ = p->dt/ti_;
tdByDt_ = td_ / p->dt;
- if (Kp_ == 0){
- vector<Id> compartments;
- unsigned int numComp = e.element()->getNeighbors(compartments, currentOut());
- if (numComp > 0){
- double Cm = Field<double>::get(compartments[0], "Cm");
+
+ vector<Id> compartments;
+ unsigned int numComp = e.element()->getNeighbors(compartments, currentOut());
+ if (numComp > 0){
+ double Cm = Field<double>::get(compartments[0], "Cm");
+ if (Kp_ == 0){
Kp_ = Cm / p->dt;
- }
+ }
+ command_ = cmdIn_ = oldCmdIn_ =
+ Field<double>::get(compartments[0], "initVm");
}
}
diff --git a/device/DiffAmp.cpp b/device/DiffAmp.cpp
index 60bbb422..7f7059a4 100644
--- a/device/DiffAmp.cpp
+++ b/device/DiffAmp.cpp
@@ -206,6 +206,7 @@ void DiffAmp::reinit(const Eref& e, ProcPtr p)
output_ = 0.0;
plus_ = 0.0;
minus_ = 0.0;
+ outputOut()->send(e, output_);
}
//
diff --git a/device/RC.cpp b/device/RC.cpp
index 6ce69f23..b8b07f36 100644
--- a/device/RC.cpp
+++ b/device/RC.cpp
@@ -129,7 +129,7 @@ RC::RC():
state_(0),
inject_(0),
msg_inject_(0.0),
- exp_(0.0),
+ expTau_(0.0),
dt_tau_(0.0)
{
; // Do nothing
@@ -203,13 +203,30 @@ void RC::setInjectMsg( double inject )
void RC::process(const Eref& e, const ProcPtr proc )
{
+ /*
double sum_inject_prev = inject_ + msg_inject_;
double sum_inject = inject_ + msg_inject_;
double dVin = (sum_inject - sum_inject_prev) * resistance_;
double Vin = sum_inject * resistance_;
state_ = Vin + dVin - dVin / dt_tau_ +
- (state_ - Vin + dVin / dt_tau_) * exp_;
+ (state_ - Vin + dVin / dt_tau_) * expTau_;
sum_inject_prev = sum_inject;
+ msg_inject_ = 0.0;
+ outputOut()->send(e, state_);
+ */
+
+ ////////////////////////////////////////////////////////////////
+ // double A = inject + msgInject_;
+ // double B = 1.0/resistance_;
+ // double x = exp( -B * proc->dt / capacitance_ );
+ // state_ = state_ * x + (A/B) * (1.0-x);
+ // double x = exp( -dt_tau_ );
+ state_ = state_ * expTau_ +
+ resistance_*(inject_ + msg_inject_) * (1.0-expTau_);
+
+ /// OR, if we use simple Euler: ///
+ // state_ += (inject_ + msgInject_ - state_/resistance_ ) * proc->dt / capacitance_;
+
msg_inject_ = 0.0;
outputOut()->send(e, state_);
}
@@ -220,9 +237,9 @@ void RC::reinit(const Eref& e, const ProcPtr proc)
dt_tau_ = proc->dt / (resistance_ * capacitance_);
state_ = v0_;
if (dt_tau_ > 1e-15){
- exp_ = exp(-dt_tau_);
+ expTau_ = exp(-dt_tau_);
} else {// use approximation
- exp_ = 1 - dt_tau_;
+ expTau_ = 1 - dt_tau_;
}
msg_inject_ = 0.0;
outputOut()->send(e, state_);
diff --git a/device/RC.h b/device/RC.h
index 3023ddbb..abca9dbf 100644
--- a/device/RC.h
+++ b/device/RC.h
@@ -66,7 +66,7 @@ class RC{
double state_;
double inject_;
double msg_inject_;
- double exp_;
+ double expTau_;
double dt_tau_;
};
diff --git a/diffusion/ConcChanInfo.h b/diffusion/ConcChanInfo.h
new file mode 100644
index 00000000..cd16266b
--- /dev/null
+++ b/diffusion/ConcChanInfo.h
@@ -0,0 +1,31 @@
+/**********************************************************************
+** This program is part of 'MOOSE', the
+** Messaging Object Oriented Simulation Environment.
+** Copyright (C) 2003-2012 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.
+**********************************************************************/
+
+#ifndef _CONC_CHAN_INFO_H
+#define _CONC_CHAN_INFO_H
+
+/**
+ */
+class ConcChanInfo
+{
+ public:
+ ConcChanInfo();
+ ConcChanInfo( unsigned int my, unsigned int other,
+ unsigned int chan, double perm )
+ : myPool( my ), otherPool( other ), chanPool( chan ),
+ permeability( perm )
+ {;}
+
+ unsigned int myPool;
+ unsigned int otherPool;
+ unsigned int chanPool;
+ double permeability;
+};
+
+#endif // _CONC_CHAN_INFO_H
diff --git a/diffusion/DiffJunction.h b/diffusion/DiffJunction.h
index c2ac3070..a4e44d23 100644
--- a/diffusion/DiffJunction.h
+++ b/diffusion/DiffJunction.h
@@ -15,6 +15,17 @@
*
* The data includes all the pools that diffuse, which are identified
* by their indices in the respective Dsolves.
+ *
+ * It includes the pools that are transferred directly, due to cross-
+ * compartment reactions
+ *
+ * It includes the pools that go down a conc gradient due to a
+ * channel, and also specifies the channel. Equation is the same as for
+ * diffusion, but scaled by the channel conductance and amount. Channel
+ * may also do some local calculations for Nernst potentials if applicable.
+ * Note that this will be flawed unless current is somehow coupled to memb
+ * potential.
+ *
* It also includes the vector of VoxelJunctions, which specifies
* matching voxel indices in the two compartments, and the diffusion scale
* factor (based on geometry) for this junction.
@@ -25,5 +36,15 @@ class DiffJunction
unsigned int otherDsolve;
vector< unsigned int > myPools;
vector< unsigned int > otherPools;
+
+ vector< unsigned int > myXferSrc;
+ vector< unsigned int > otherXferDest;
+
+ vector< unsigned int > myXferDest;
+ vector< unsigned int > otherXferSrc;
+
+ vector< unsigned int > myChannels;
+ vector< unsigned int > otherChannels;
+
vector< VoxelJunction > vj;
};
diff --git a/diffusion/DiffPoolVec.cpp b/diffusion/DiffPoolVec.cpp
index e6faaa63..cd1da649 100644
--- a/diffusion/DiffPoolVec.cpp
+++ b/diffusion/DiffPoolVec.cpp
@@ -51,6 +51,12 @@ void DiffPoolVec::setN( unsigned int voxel, double v )
n_[ voxel ] = v;
}
+double DiffPoolVec::getPrev( unsigned int voxel ) const
+{
+ assert( voxel < n_.size() );
+ return prev_[ voxel ];
+}
+
const vector< double >& DiffPoolVec::getNvec() const
{
return n_;
@@ -71,6 +77,11 @@ void DiffPoolVec::setNvec( unsigned int start, unsigned int num,
*p++ = *q++;
}
+void DiffPoolVec::setPrevVec()
+{
+ prev_ = n_;
+}
+
double DiffPoolVec::getDiffConst() const
{
return diffConst_;
@@ -142,5 +153,5 @@ void DiffPoolVec::advance( double dt )
void DiffPoolVec::reinit() // Not called by the clock, but by parent.
{
assert( n_.size() == nInit_.size() );
- n_ = nInit_;
+ prev_ = n_ = nInit_;
}
diff --git a/diffusion/DiffPoolVec.h b/diffusion/DiffPoolVec.h
index 345403f6..587533eb 100644
--- a/diffusion/DiffPoolVec.h
+++ b/diffusion/DiffPoolVec.h
@@ -26,6 +26,7 @@ class DiffPoolVec
void setNinit( unsigned int vox, double value );
double getN( unsigned int vox ) const;
void setN( unsigned int vox, double value );
+ double getPrev( unsigned int vox ) const;
double getDiffConst() const;
void setDiffConst( double value );
@@ -46,6 +47,7 @@ class DiffPoolVec
void setNvec( const vector< double >& n );
void setNvec( unsigned int start, unsigned int num,
vector< double >::const_iterator q );
+ void setPrevVec(); /// Assigns prev_ = n_
void setOps( const vector< Triplet< double > >& ops_,
const vector< double >& diagVal_ ); /// Assign operations.
@@ -53,6 +55,7 @@ class DiffPoolVec
private:
unsigned int id_; /// Integer conversion of Id of pool handled.
vector< double > n_; /// Number of molecules of pool in each voxel
+ vector< double > prev_; /// # molecules of pool on previous timestep
vector< double > nInit_; /// Boundary condition: Initial 'n'.
double diffConst_; /// Diffusion const, assumed uniform
double motorConst_; /// Motor const, ie, transport rate.
diff --git a/diffusion/Dsolve.cpp b/diffusion/Dsolve.cpp
index 9185ab65..61ffb1bc 100644
--- a/diffusion/Dsolve.cpp
+++ b/diffusion/Dsolve.cpp
@@ -15,7 +15,9 @@
#include "../mesh/VoxelJunction.h"
#include "XferInfo.h"
#include "ZombiePoolInterface.h"
+#include "../kinetics/ConcChan.h"
#include "DiffPoolVec.h"
+#include "ConcChanInfo.h"
#include "FastMatrixElim.h"
#include "../mesh/VoxelJunction.h"
#include "DiffJunction.h"
@@ -304,22 +306,9 @@ static double integ( double myN, double rf, double rb, double dt )
return myN;
}
-/**
- * Computes flux through a junction between diffusion solvers.
- * Most used at junctions on spines and PSDs, but can also be used
- * when a given diff solver is decomposed. At present the lookups
- * on the other diffusion solver assume that the data is on the local
- * node. Once this works well I can figure out how to do across nodes.
- */
-void Dsolve::calcJunction( const DiffJunction& jn, double dt )
+void Dsolve::calcJnDiff( const DiffJunction& jn, Dsolve* other, double dt)
{
- const double EPSILON = 1e-15;
- Id oid( jn.otherDsolve );
- assert ( oid != Id() );
- assert ( oid.element()->cinfo()->isA( "Dsolve" ) );
-
- Dsolve* other = reinterpret_cast< Dsolve* >( oid.eref().data() );
-
+ const double EPSILON = 1e-16;
assert( jn.otherPools.size() == jn.myPools.size() );
for ( unsigned int i = 0; i < jn.myPools.size(); ++i ) {
DiffPoolVec& myDv = pools_[ jn.myPools[i] ];
@@ -357,17 +346,143 @@ void Dsolve::calcJunction( const DiffJunction& jn, double dt )
}
}
+void Dsolve::calcJnXfer( const DiffJunction& jn,
+ const vector< unsigned int >& srcXfer,
+ const vector< unsigned int >& destXfer,
+ Dsolve* srcDsolve, Dsolve* destDsolve )
+{
+ assert( destXfer.size() == srcXfer.size() );
+ for ( unsigned int i = 0; i < srcXfer.size(); ++i ) {
+ DiffPoolVec& srcDv = srcDsolve->pools_[ srcXfer[i] ];
+ DiffPoolVec& destDv = destDsolve->pools_[ destXfer[i] ];
+ for ( vector< VoxelJunction >::const_iterator
+ j = jn.vj.begin(); j != jn.vj.end(); ++j ) {
+ double prevSrc = srcDv.getPrev( j->first );
+ double prevDest = destDv.getPrev( j->second );
+ double srcN = srcDv.getN( j->first );
+ double destN = destDv.getN( j->second );
+ // Consider delta as sum of local dN, and reference as prevDest
+ // newN = (srcN - prevSrc + destN - prevDest) + prevDest
+ double newN = srcN + destN - prevSrc;
+ srcDv.setN( j->first, newN );
+ destDv.setN( j->second, newN );
+ }
+ }
+}
+
+void Dsolve::calcJnChan( const DiffJunction& jn, Dsolve* other, double dt )
+{
+ // Each jn has some channels
+ // Each channel has a chanPool, an intPool and an extPool.
+ // intPool is on self, extPool is on other, but we have a problem
+ // because the chanPool could be a third compartment, such as the memb
+ // If we stipulate it is is on self, that is easy but not general.
+ // Other alternative is to have a message to update the N of the chan,
+ // so it isn't in the domain of the solver at all except for here.
+ // In which case we will want to point to the Moose object for it.
+ //
+ //
+
+ for ( unsigned int i = 0; i < jn.myChannels.size(); ++i ) {
+ ConcChanInfo& myChan = channels_[ jn.myChannels[i] ];
+ DiffPoolVec& myDv = pools_[ myChan.myPool ];
+ DiffPoolVec& otherDv = other->pools_[ myChan.otherPool ];
+ DiffPoolVec& chanDv = pools_[ myChan.chanPool ];
+ /*
+ DiffPoolVec& myDv = pools_[ jn.myPools[myChan.myPool] ];
+ DiffPoolVec& otherDv =
+ other->pools_[ jn.otherPools[myChan.otherPool] ];
+ DiffPoolVec& chanDv = pools_[ jn.myPools[myChan.chanPool] ];
+ */
+ for ( vector< VoxelJunction >::const_iterator
+ j = jn.vj.begin(); j != jn.vj.end(); ++j ) {
+
+ double myN = myDv.getN( j->first );
+ double lastN = myN;
+ double otherN = otherDv.getN( j->second );
+ double perm = myChan.permeability * chanDv.getN( j->first );
+ myN = integ( myN, perm * myN/j->firstVol,
+ perm * otherN/j->secondVol, dt );
+ otherN += lastN - myN; // Mass consv
+ if ( otherN < 0.0 ) { // Avoid negatives
+ myN += otherN;
+ otherN = 0.0;
+ }
+ myDv.setN( j->first, myN );
+ otherDv.setN( j->second, otherN );
+ }
+ }
+}
+
+// Same as above, but now go through channels on other Dsolve.
+void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt )
+{
+ for ( unsigned int i = 0; i < jn.otherChannels.size(); ++i ) {
+ ConcChanInfo& otherChan = other->channels_[ jn.otherChannels[i] ];
+ // This is the DiffPoolVec for the pools on the other Dsolve,
+ // the one with the channel.
+ // DiffPoolVec& otherDv = other->pools_[ jn.otherPools[otherChan.myPool] ];
+ DiffPoolVec& otherDv = other->pools_[ otherChan.myPool ];
+ // Local diffPoolVec.
+ // DiffPoolVec& myDv = pools_[ jn.myPools[otherChan.otherPool] ];
+ DiffPoolVec& myDv = pools_[ otherChan.otherPool ];
+ DiffPoolVec& chanDv = other->pools_[ otherChan.chanPool ];
+ for ( vector< VoxelJunction >::const_iterator
+ j = jn.vj.begin(); j != jn.vj.end(); ++j ) {
+
+ double myN = myDv.getN( j->first );
+ double lastN = myN;
+ double otherN = otherDv.getN( j->second );
+ double perm = otherChan.permeability * chanDv.getN(j->second);
+ myN = integ( myN, perm * myN/j->firstVol,
+ perm * otherN/j->secondVol, dt );
+ otherN += lastN - myN; // Mass consv
+ if ( otherN < 0.0 ) { // Avoid negatives
+ myN += otherN;
+ otherN = 0.0;
+ }
+ myDv.setN( j->first, myN );
+ otherDv.setN( j->second, otherN );
+ }
+ }
+}
+
+/**
+ * Computes flux through a junction between diffusion solvers.
+ * Most used at junctions on spines and PSDs, but can also be used
+ * when a given diff solver is decomposed. At present the lookups
+ * on the other diffusion solver assume that the data is on the local
+ * node. Once this works well I can figure out how to do across nodes.
+ * Note that we split the diffusion and channel calculations before and
+ * after then calcJnXfer calculations. This improves accuracy by 5x.
+ */
+void Dsolve::calcJunction( const DiffJunction& jn, double dt )
+{
+ Id oid( jn.otherDsolve );
+ assert ( oid != Id() );
+ assert ( oid.element()->cinfo()->isA( "Dsolve" ) );
+
+ Dsolve* other = reinterpret_cast< Dsolve* >( oid.eref().data() );
+ calcJnDiff( jn, other, dt/2.0 );
+
+ calcJnChan( jn, other, dt/2.0 );
+ calcOtherJnChan( jn, other, dt/2.0 );
+
+ calcJnXfer( jn, jn.myXferSrc, jn.otherXferDest, this, other );
+ calcJnXfer( jn, jn.otherXferSrc, jn.myXferDest, other, this );
+
+ calcJnDiff( jn, other, dt/2.0 );
+
+ calcJnChan( jn, other, dt/2.0 );
+ calcOtherJnChan( jn, other, dt/2.0 );
+}
+
void Dsolve::process( const Eref& e, ProcPtr p )
{
for ( vector< DiffPoolVec >::iterator
i = pools_.begin(); i != pools_.end(); ++i ) {
i->advance( p->dt );
}
-
- for ( vector< DiffJunction >::const_iterator
- i = junctions_.begin(); i != junctions_.end(); ++i ) {
- calcJunction( *i, p->dt );
- }
}
void Dsolve::reinit( const Eref& e, ProcPtr p )
@@ -378,6 +493,14 @@ void Dsolve::reinit( const Eref& e, ProcPtr p )
i->reinit();
}
}
+
+void Dsolve::updateJunctions( double dt )
+{
+ for ( vector< DiffJunction >::const_iterator
+ i = junctions_.begin(); i != junctions_.end(); ++i ) {
+ calcJunction( *i, dt );
+ }
+}
//////////////////////////////////////////////////////////////
// Solver coordination and setup functions
//////////////////////////////////////////////////////////////
@@ -417,6 +540,54 @@ void Dsolve::setStoich( Id id )
*/
}
}
+ string chanpath = path_ + "[ISA=ConcChan]";
+ vector< ObjId > chans;
+ wildcardFind( chanpath, chans );
+ fillConcChans( chans );
+}
+
+void Dsolve::fillConcChans( const vector< ObjId >& chans )
+{
+ static const Cinfo* ccc = Cinfo::find( "ConcChan" );
+ static const Finfo* inPoolFinfo = ccc->findFinfo( "inPool" );
+ static const Finfo* outPoolFinfo = ccc->findFinfo( "outPool" );
+ static const Finfo* chanPoolFinfo = ccc->findFinfo( "setNumChan" );
+ FuncId fin = static_cast< const DestFinfo* >( inPoolFinfo )->getFid();
+ FuncId fout = static_cast< const DestFinfo* >(outPoolFinfo )->getFid();
+ FuncId fchan =
+ static_cast< const DestFinfo* >(chanPoolFinfo )->getFid();
+
+ // Find the in pools and the chan pools on the current compt.
+ // Save the Id of the outPool as an integer.
+ // Use these and the permeability to create the ConcChanInfo.
+ for ( auto i = chans.begin(); i != chans.end(); ++i ) {
+ vector< Id > ret;
+ if (i->element()->getNeighbors( ret, inPoolFinfo ) == 0 ) return;
+ ObjId inPool( ret[0] );
+ ret.clear();
+ if (i->element()->getNeighbors( ret, outPoolFinfo ) == 0 ) return;
+ ObjId outPool( ret[0] );
+ ret.clear();
+ if (i->element()->getNeighbors( ret, chanPoolFinfo ) == 0 ) return;
+ ObjId chanPool( ret[0] );
+ ret.clear();
+ /*
+ ObjId inPool = i->element()->findCaller( fin );
+ ObjId chanPool = i->element()->findCaller( fchan );
+ ObjId outPool = i->element()->findCaller( fout );
+ */
+ if ( !( inPool.bad() or chanPool.bad() ) ) {
+ unsigned int inPoolIndex = convertIdToPoolIndex( inPool.id );
+ unsigned int chanPoolIndex = convertIdToPoolIndex(chanPool.id);
+ if ( inPoolIndex != ~0U && chanPoolIndex != ~0U ) {
+ ConcChanInfo cci( inPoolIndex, outPool.id.value(),
+ chanPoolIndex,
+ Field< double >::get( *i, "permeability" )
+ );
+ channels_.push_back( cci );
+ }
+ }
+ }
}
Id Dsolve::getStoich() const
@@ -431,18 +602,18 @@ void Dsolve::setDsolve( Id dsolve )
void Dsolve::setCompartment( Id id )
{
const Cinfo* c = id.element()->cinfo();
- if ( c->isA( "NeuroMesh" ) || c->isA( "SpineMesh" ) ||
- c->isA( "PsdMesh" ) || c->isA( "CylMesh" ) ) {
- compartment_ = id;
- numVoxels_ = Field< unsigned int >::get( id, "numMesh" );
- /*
- const MeshCompt* m = reinterpret_cast< const MeshCompt* >(
- id.eref().data() );
- numVoxels_ = m->getStencil().nRows();
- */
- } else {
- cout << "Warning: Dsolve::setCompartment:: compartment must be "
- "NeuroMesh or CylMesh, you tried :" << c->name() << endl;
+ compartment_ = id;
+ numVoxels_ = Field< unsigned int >::get( id, "numMesh" );
+ if ( c->isA( "CubeMesh" ) ) { // we do only linear diffusion for now
+ unsigned int nx = Field< unsigned int >::get( id, "nx" );
+ unsigned int ny = Field< unsigned int >::get( id, "nx" );
+ unsigned int nz = Field< unsigned int >::get( id, "nx" );
+ if ( !( nx*ny == 1 || nx*nz == 1 || ny*nz == 1 ) ) {
+ cout << "Warning: Dsolve::setCompartment:: Cube mesh: " <<
+ c->name() << " found with >1 dimension of voxels. " <<
+ "Only 1-D diffusion supported for now.\n";
+ return;
+ }
}
}
@@ -474,6 +645,8 @@ void Dsolve::makePoolMapFromElist( const vector< ObjId >& elist,
stoich_ = Id();
poolMapStart_ = minId;
poolMap_.resize( 1 + maxId - minId );
+ for ( auto i = poolMap_.begin(); i != poolMap_.end(); ++i )
+ *i = ~0U;
for ( unsigned int i = 0; i < temp.size(); ++i ) {
unsigned int idValue = temp[i].value();
assert( idValue >= minId );
@@ -610,8 +783,8 @@ void Dsolve::buildNeuroMeshJunctions( const Eref& e, Id spineD, Id psdD )
return;
}
- innerBuildMeshJunctions( spineD, e.id() );
- innerBuildMeshJunctions( psdD, spineD );
+ innerBuildMeshJunctions( spineD, e.id(), false );
+ innerBuildMeshJunctions( psdD, spineD, false );
}
void Dsolve::buildMeshJunctions( const Eref& e, Id other )
@@ -621,8 +794,10 @@ void Dsolve::buildMeshJunctions( const Eref& e, Id other )
otherMesh = Field< Id >::get( other, "compartment" );
if ( compartment_.element()->cinfo()->isA( "ChemCompt" ) &&
otherMesh.element()->cinfo()->isA( "ChemCompt" ) ) {
- innerBuildMeshJunctions( e.id(), other );
- return;
+ bool isMembraneBound =
+ Field< bool >::get( compartment_, "isMembraneBound" );
+ innerBuildMeshJunctions( e.id(), other, isMembraneBound );
+ return;
}
}
cout << "Warning: Dsolve::buildMeshJunctions: one of '" <<
@@ -644,14 +819,11 @@ void printJunction( Id self, Id other, const DiffJunction& jn )
}
}
-// Static utility func for building junctions
-void Dsolve::innerBuildMeshJunctions( Id self, Id other )
+void Dsolve::mapDiffPoolsBetweenDsolves( DiffJunction& jn,
+ Id self, Id other)
{
- DiffJunction jn; // This is based on the Spine Dsolver.
- jn.otherDsolve = other.value();
- // Map pools between Dsolves
Dsolve* mySolve = reinterpret_cast< Dsolve* >( self.eref().data() );
- map< string, unsigned int > myPools;
+ unordered_map< string, unsigned int > myPools;
for ( unsigned int i = 0; i < mySolve->pools_.size(); ++i ) {
Id pool( mySolve->pools_[i].getId() );
assert( pool != Id() );
@@ -662,15 +834,83 @@ void Dsolve::innerBuildMeshJunctions( Id self, Id other )
other.eref().data() );
for ( unsigned int i = 0; i < otherSolve->pools_.size(); ++i ) {
Id otherPool( otherSolve->pools_[i].getId() );
- map< string, unsigned int >::iterator p =
+ unordered_map< string, unsigned int >::iterator p =
myPools.find( otherPool.element()->getName() );
if ( p != myPools.end() ) {
jn.otherPools.push_back( i );
jn.myPools.push_back( p->second );
}
}
+}
+
+/**
+ * static void mapXfersBetweenDsolves(...)
+ * Build a list of all the molecules that should transfer instantaneously
+ * to another compartment, for cross-compartment reactions.
+ * Here we look for src names of the form <name>_xfer_<destComptName>
+ * For example, if we had an enzyme in compartment 'src', whose product
+ * should go to pool 'bar' in compartment 'dest',
+ * the name of the enzyme product in compartment src would be
+ * bar_xfer_dest
+ */
+void Dsolve::mapXfersBetweenDsolves(
+ vector< unsigned int >& srcPools, vector< unsigned int >& destPools,
+ Id src, Id dest )
+{
+ Id destMesh = Field< Id >::get( dest, "compartment" );
+ string xferPost( string( "_xfer_" ) + destMesh.element()->getName() );
+ size_t xlen = xferPost.length();
+
+ Dsolve* srcSolve = reinterpret_cast< Dsolve* >( src.eref().data() );
+ unordered_map< string, unsigned int > srcMap;
+ for ( unsigned int i = 0; i < srcSolve->pools_.size(); ++i ) {
+ Id pool( srcSolve->pools_[i].getId() );
+ assert( pool != Id() );
+ string poolName = pool.element()->getName();
+ size_t prefixLen = poolName.length() - xlen;
+ if ( poolName.rfind( xferPost ) == prefixLen )
+ srcMap[ poolName.substr( 0, prefixLen) ] = i;
+ }
+
+ const Dsolve* destSolve = reinterpret_cast< const Dsolve* >(
+ dest.eref().data() );
+ for ( unsigned int i = 0; i < destSolve->pools_.size(); ++i ) {
+ Id destPool( destSolve->pools_[i].getId() );
+ unordered_map< string, unsigned int >::iterator p =
+ srcMap.find( destPool.element()->getName() );
+ if ( p != srcMap.end() ) {
+ destPools.push_back( i );
+ srcPools.push_back( p->second );
+ }
+ }
+}
+
+void Dsolve::mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other)
+{
+ Dsolve* otherSolve = reinterpret_cast< Dsolve* >(
+ other.eref().data() );
+ vector< ConcChanInfo >& ch = channels_;
+ unsigned int outIndex;
+ for ( unsigned int i = 0; i < ch.size(); ++i ) {
+ outIndex = otherSolve->convertIdToPoolIndex( ch[i].otherPool );
+ if ( outIndex != ~0U ) {
+ jn.myChannels.push_back(i);
+ ch[i].otherPool = outIndex; // replace the Id with the index.
+ }
+ }
+ // Now set up the other Dsolve.
+ vector< ConcChanInfo >& ch2 = otherSolve->channels_;
+ for ( unsigned int i = 0; i < ch2.size(); ++i ) {
+ outIndex = convertIdToPoolIndex( ch2[i].otherPool );
+ if ( outIndex != ~0U ) {
+ jn.otherChannels.push_back(i);
+ ch2[i].otherPool = outIndex; // replace the Id with the index
+ }
+ }
+}
- // Map voxels between meshes.
+static void mapVoxelsBetweenMeshes( DiffJunction& jn, Id self, Id other)
+{
Id myMesh = Field< Id >::get( self, "compartment" );
Id otherMesh = Field< Id >::get( other, "compartment" );
@@ -686,9 +926,25 @@ void Dsolve::innerBuildMeshJunctions( Id self, Id other )
i->firstVol = myVols[i->first];
i->secondVol = otherVols[i->second];
}
+}
+
+// Static utility func for building junctions
+void Dsolve::innerBuildMeshJunctions( Id self, Id other, bool selfIsMembraneBound )
+{
+ DiffJunction jn; // This is based on the Spine Dsolver.
+ jn.otherDsolve = other.value();
+ if ( selfIsMembraneBound ) {
+ mapChansBetweenDsolves( jn, self, other );
+ } else {
+ mapDiffPoolsBetweenDsolves( jn, self, other );
+ }
+ mapXfersBetweenDsolves( jn.myXferSrc, jn.otherXferDest, self, other );
+ mapXfersBetweenDsolves( jn.otherXferSrc, jn.myXferDest, other, self );
+
+ mapVoxelsBetweenMeshes( jn, self, other );
// printJunction( self, other, jn );
- mySolve->junctions_.push_back( jn );
+ junctions_.push_back( jn );
}
/////////////////////////////////////////////////////////////
@@ -712,18 +968,23 @@ void Dsolve::setNumAllVoxels( unsigned int num )
pools_[i].setNumVoxels( numVoxels_ );
}
-unsigned int Dsolve::convertIdToPoolIndex( const Eref& e ) const
+unsigned int Dsolve::convertIdToPoolIndex( const Id id ) const
{
- unsigned int i = e.id().value() - poolMapStart_;
+ unsigned int i = id.value() - poolMapStart_;
if ( i < poolMap_.size() ) {
return poolMap_[i];
}
cout << "Warning: Dsolve::convertIdToPoollndex: Id out of range, (" <<
- poolMapStart_ << ", " << e.id() << ", " <<
+ poolMapStart_ << ", " << id << ", " << id.path() << ", " <<
poolMap_.size() + poolMapStart_ << "\n";
return 0;
}
+unsigned int Dsolve::convertIdToPoolIndex( const Eref& e ) const
+{
+ return convertIdToPoolIndex( e.id() );
+}
+
void Dsolve::setN( const Eref& e, double v )
{
unsigned int pid = convertIdToPoolIndex( e );
@@ -848,6 +1109,15 @@ void Dsolve::getBlock( vector< double >& values ) const
}
}
+// Inefficient but easy to set up. Optimize later.
+void Dsolve::setPrev()
+{
+ for ( auto i = pools_.begin(); i != pools_.end(); ++i ) {
+ // if (i->getDiffConst() > 0.0 )
+ i->setPrevVec();
+ }
+}
+
void Dsolve::setBlock( const vector< double >& values )
{
unsigned int startVoxel = values[0];
diff --git a/diffusion/Dsolve.h b/diffusion/Dsolve.h
index 0a16e772..d6028e62 100644
--- a/diffusion/Dsolve.h
+++ b/diffusion/Dsolve.h
@@ -67,6 +67,9 @@ class Dsolve: public ZombiePoolInterface
void process( const Eref& e, ProcPtr p );
void reinit( const Eref& e, ProcPtr p );
+ //////////////////////////////////////////////////////////////////
+ void updateJunctions( double dt );
+
/**
* Builds junctions between Dsolves handling NeuroMesh, SpineMesh,
* and PsdMesh. Must only be called from the one handling the
@@ -93,7 +96,18 @@ class Dsolve: public ZombiePoolInterface
* the junction between any specified pair of Dsolves.
* Note that it builds the junction on the 'self' Dsolve.
*/
- static void innerBuildMeshJunctions( Id self, Id other );
+ void innerBuildMeshJunctions( Id self, Id other,
+ bool isMembraneBound );
+
+ /// Sets up map of matching pools for diffusion.
+ static void mapDiffPoolsBetweenDsolves( DiffJunction& jn,
+ Id self, Id other);
+
+ static void mapXfersBetweenDsolves(
+ vector< unsigned int >& srcPools,
+ vector< unsigned int >& destPools,
+ Id src, Id dest );
+ void mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other);
/**
* Computes flux through a junction between diffusion solvers.
@@ -123,6 +137,7 @@ class Dsolve: public ZombiePoolInterface
void getBlock( vector< double >& values ) const;
void setBlock( const vector< double >& values );
+ void setPrev();
// This one isn't used in Dsolve, but is defined as a dummy.
void setupCrossSolverReacs(
@@ -137,6 +152,7 @@ class Dsolve: public ZombiePoolInterface
//////////////////////////////////////////////////////////////////
// Model traversal and building functions
//////////////////////////////////////////////////////////////////
+ unsigned int convertIdToPoolIndex( Id id ) const;
unsigned int convertIdToPoolIndex( const Eref& e ) const;
unsigned int getPoolIndex( const Eref& e ) const;
@@ -156,6 +172,15 @@ class Dsolve: public ZombiePoolInterface
*/
void build( double dt );
void rebuildPools();
+ void calcJnDiff( const DiffJunction& jn, Dsolve* other, double dt );
+ void calcJnXfer( const DiffJunction& jn,
+ const vector< unsigned int >& srcXfer,
+ const vector< unsigned int >& destXfer,
+ Dsolve* srcDsolve, Dsolve* destDsolve );
+ void calcJnChan( const DiffJunction& jn, Dsolve* other, double dt );
+ void calcOtherJnChan( const DiffJunction& jn, Dsolve* other,
+ double dt );
+ void fillConcChans( const vector< ObjId >& chans );
/**
* Utility func for debugging: Prints N_ matrix
@@ -179,6 +204,8 @@ class Dsolve: public ZombiePoolInterface
/// Internal vector, one for each pool species managed by Dsolve.
vector< DiffPoolVec > pools_;
+ /// Internal vector, one for each ConcChan managed by Dsolve.
+ vector< ConcChanInfo > channels_;
/// smallest Id value for poolMap_
unsigned int poolMapStart_;
diff --git a/hsolve/HSolve.cpp b/hsolve/HSolve.cpp
index 05d5ff8e..e4b9b0e7 100644
--- a/hsolve/HSolve.cpp
+++ b/hsolve/HSolve.cpp
@@ -23,6 +23,7 @@
#include "../biophysics/ChanBase.h"
#include "../biophysics/ChanCommon.h"
#include "../biophysics/HHChannel.h"
+#include "../biophysics/CaConc.h"
#include "ZombieHHChannel.h"
#include "../shell/Shell.h"
@@ -192,6 +193,11 @@ HSolve::HSolve()
;
}
+HSolve::~HSolve()
+{
+ unzombify();
+}
+
///////////////////////////////////////////////////
// Dest function definitions
@@ -215,23 +221,48 @@ void HSolve::zombify( Eref hsolve ) const
for ( i = compartmentId_.begin(); i != compartmentId_.end(); ++i )
temp.push_back( ObjId( *i, 0 ) );
- for ( i = compartmentId_.begin(); i != compartmentId_.end(); ++i )
+ for ( i = compartmentId_.begin(); i != compartmentId_.end(); ++i ) {
CompartmentBase::zombify( i->eref().element(),
ZombieCompartment::initCinfo(), hsolve.id() );
+ }
temp.clear();
for ( i = caConcId_.begin(); i != caConcId_.end(); ++i )
temp.push_back( ObjId( *i, 0 ) );
// Shell::dropClockMsgs( temp, "process" );
- for ( i = caConcId_.begin(); i != caConcId_.end(); ++i )
+ for ( i = caConcId_.begin(); i != caConcId_.end(); ++i ) {
CaConcBase::zombify( i->eref().element(), ZombieCaConc::initCinfo(), hsolve.id() );
+ }
temp.clear();
for ( i = channelId_.begin(); i != channelId_.end(); ++i )
temp.push_back( ObjId( *i, 0 ) );
- for ( i = channelId_.begin(); i != channelId_.end(); ++i )
+ for ( i = channelId_.begin(); i != channelId_.end(); ++i ) {
HHChannelBase::zombify( i->eref().element(),
ZombieHHChannel::initCinfo(), hsolve.id() );
+ }
+}
+
+void HSolve::unzombify() const
+{
+ vector< Id >::const_iterator i;
+
+ for ( i = compartmentId_.begin(); i != compartmentId_.end(); ++i )
+ if ( i->element() ) {
+ CompartmentBase::zombify( i->eref().element(),
+ Compartment::initCinfo(), Id() );
+ }
+
+ for ( i = caConcId_.begin(); i != caConcId_.end(); ++i )
+ if ( i->element() ) {
+ CaConcBase::zombify( i->eref().element(), CaConc::initCinfo(), Id() );
+ }
+
+ for ( i = channelId_.begin(); i != channelId_.end(); ++i )
+ if ( i->element() ) {
+ HHChannelBase::zombify( i->eref().element(),
+ HHChannel::initCinfo(), Id() );
+ }
}
void HSolve::setup( Eref hsolve )
diff --git a/hsolve/HSolve.h b/hsolve/HSolve.h
index 9022a1f0..56c15d2d 100644
--- a/hsolve/HSolve.h
+++ b/hsolve/HSolve.h
@@ -17,6 +17,7 @@ class HSolve: public HSolveActive
{
public:
HSolve();
+ ~HSolve();
void process( const Eref& hsolve, ProcPtr p );
void reinit( const Eref& hsolve, ProcPtr p );
@@ -157,6 +158,7 @@ private:
void setup( Eref hsolve );
void zombify( Eref hsolve ) const;
+ void unzombify() const;
// Mapping global Id to local index. Defined in HSolveInterface.cpp.
void mapIds();
diff --git a/hsolve/HSolveActive.cpp b/hsolve/HSolveActive.cpp
index 9d34c638..ebc62b34 100644
--- a/hsolve/HSolveActive.cpp
+++ b/hsolve/HSolveActive.cpp
@@ -356,13 +356,13 @@ void HSolveActive::sendValues( ProcPtr info )
{
vector< unsigned int >::iterator i;
- for ( i = outVm_.begin(); i != outVm_.end(); ++i )
+ for ( i = outVm_.begin(); i != outVm_.end(); ++i ) {
Compartment::VmOut()->send(
//~ ZombieCompartment::VmOut()->send(
compartmentId_[ *i ].eref(),
V_[ *i ]
);
-
+ }
for ( i = outIk_.begin(); i != outIk_.end(); ++i ){
diff --git a/hsolve/HSolveUtils.cpp b/hsolve/HSolveUtils.cpp
index 77e0cc10..5038c6e4 100644
--- a/hsolve/HSolveUtils.cpp
+++ b/hsolve/HSolveUtils.cpp
@@ -320,9 +320,9 @@ int HSolveUtils::targets(
e->getNeighbors( all, f );
vector< Id >::iterator ia;
- if ( filter.empty() )
+ if ( filter.empty() ) {
target.insert( target.end(), all.begin(), all.end() );
- else
+ } else {
for ( ia = all.begin(); ia != all.end(); ++ia ) {
string className = (*ia).element()->cinfo()->name();
bool hit =
@@ -335,6 +335,7 @@ int HSolveUtils::targets(
if ( ( hit && include ) || ( !hit && !include ) )
target.push_back( *ia );
}
+ }
return target.size() - oldSize;
}
diff --git a/kinetics/CMakeLists.txt b/kinetics/CMakeLists.txt
index 403f3603..b836312c 100644
--- a/kinetics/CMakeLists.txt
+++ b/kinetics/CMakeLists.txt
@@ -11,6 +11,7 @@ add_library(kinetics
Enz.cpp
MMenz.cpp
Species.cpp
+ ConcChan.cpp
ReadKkit.cpp
WriteKkit.cpp
ReadCspace.cpp
diff --git a/kinetics/ConcChan.cpp b/kinetics/ConcChan.cpp
new file mode 100644
index 00000000..62170845
--- /dev/null
+++ b/kinetics/ConcChan.cpp
@@ -0,0 +1,206 @@
+/**********************************************************************
+** 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 "header.h"
+#include "ConcChan.h"
+
+static SrcFinfo2< double, double > *inPoolOut() {
+ static SrcFinfo2< double, double > inPoolOut(
+ "inPoolOut",
+ "Sends out increment to molecules on inside of membrane"
+ );
+ return &inPoolOut;
+}
+
+static SrcFinfo2< double, double > *outPoolOut() {
+ static SrcFinfo2< double, double > outPoolOut(
+ "outPoolOut",
+ "Sends out increment to molecules on outside of membrane"
+ );
+ return &outPoolOut;
+}
+
+const Cinfo* ConcChan::initCinfo()
+{
+ //////////////////////////////////////////////////////////////
+ // Field Definitions
+ //////////////////////////////////////////////////////////////
+ static ValueFinfo< ConcChan, double > permeability(
+ "permeability",
+ "Permability, in units of vol/(#s) i.e., 1/(numconc.s) "
+ "Flux (#/s) = permeability * N * (#out/vol_out - #in/vol_in)",
+ &ConcChan::setPermeability,
+ &ConcChan::getPermeability
+ );
+ static ValueFinfo< ConcChan, double > numChan(
+ "numChan",
+ "numChan is the number of molecules of the channel.",
+ &ConcChan::setNumChan,
+ &ConcChan::getNumChan
+ );
+ static ReadOnlyValueFinfo< ConcChan, double > flux(
+ "flux",
+ "Flux of molecules through channel, in units of #/s ",
+ &ConcChan::getFlux
+ );
+
+ //////////////////////////////////////////////////////////////
+ // MsgDest Definitions
+ //////////////////////////////////////////////////////////////
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< ConcChan >( &ConcChan::process ) );
+ static DestFinfo reinit( "reinit",
+ "Handles reinit call",
+ new ProcOpFunc< ConcChan >( &ConcChan::reinit ) );
+
+ //////////////////////////////////////////////////////////////
+ // Shared Msg Definitions
+ //////////////////////////////////////////////////////////////
+
+ static DestFinfo inPool( "inPool",
+ "Handles # of molecules of pool inside membrane",
+ new OpFunc1< ConcChan, double >( &ConcChan::inPool ) );
+ static DestFinfo outPool( "outPool",
+ "Handles # of molecules of pool outside membrane",
+ new OpFunc1< ConcChan, double >( &ConcChan::outPool ) );
+ static Finfo* inShared[] = {
+ inPoolOut(), &inPool
+ };
+ static Finfo* outShared[] = {
+ outPoolOut(), &outPool
+ };
+ static SharedFinfo in( "in",
+ "Connects to pool on inside",
+ inShared, sizeof( inShared ) / sizeof( const Finfo* )
+ );
+ static SharedFinfo out( "out",
+ "Connects to pool on outside",
+ outShared, sizeof( outShared ) / sizeof( const Finfo* )
+ );
+ static Finfo* procShared[] = {
+ &process, &reinit
+ };
+ static SharedFinfo proc( "proc",
+ "Shared message for process and reinit",
+ procShared, sizeof( procShared ) / sizeof( const Finfo* )
+ );
+
+
+ static Finfo* concChanFinfos[] = {
+ &permeability, // Value
+ &numChan, // Value
+ &flux, // ReadOnlyValue
+ &in, // SharedFinfo
+ &out, // SharedFinfo
+ &proc, // SharedFinfo
+ };
+
+ static string doc[] =
+ {
+ "Name", "ConcChan",
+ "Author", "Upinder S. Bhalla, 2018, NCBS",
+ "Description", "Channel for molecular flow between cellular "
+ "compartments. Need not be ions, and the flux is not a current, "
+ "but number of molecules/sec. ",
+ };
+
+ static Dinfo< ConcChan > dinfo;
+ static Cinfo concChanCinfo (
+ "ConcChan",
+ Neutral::initCinfo(),
+ concChanFinfos,
+ sizeof( concChanFinfos ) / sizeof ( Finfo* ),
+ &dinfo,
+ doc,
+ sizeof( doc ) / sizeof( string )
+ );
+
+ return &concChanCinfo;
+}
+
+ static const Cinfo* concChanCinfo = ConcChan::initCinfo();
+
+//////////////////////////////////////////////////////////////
+// ConcChan internal functions
+//////////////////////////////////////////////////////////////
+
+
+ConcChan::ConcChan( )
+ : permeability_( 0.1 ),
+ n_( 0.0 ),
+ flux_( 0.0 ),
+ charge_( 0.0 ),
+ Vm_( 0.0 )
+{
+ ;
+}
+
+ConcChan::~ConcChan( )
+{;}
+
+//////////////////////////////////////////////////////////////
+// MsgDest Definitions
+//////////////////////////////////////////////////////////////
+
+
+void ConcChan::process( const Eref& e, ProcPtr p )
+{
+}
+
+void ConcChan::reinit( const Eref& e, ProcPtr p )
+{
+}
+
+void ConcChan::inPool( double conc )
+{
+}
+
+void ConcChan::outPool( double conc )
+{
+}
+
+//////////////////////////////////////////////////////////////
+// Field Definitions
+//////////////////////////////////////////////////////////////
+
+void ConcChan::setPermeability( double v )
+{
+ permeability_ = v;
+}
+
+double ConcChan::getPermeability() const
+{
+ return permeability_;
+}
+
+void ConcChan::setNumChan( double v )
+{
+ n_ = v;
+}
+
+double ConcChan::getNumChan() const
+{
+ return n_;
+}
+
+double ConcChan::getFlux() const
+{
+ return flux_;
+}
+
+//////////////////////////////////////////////////////////////
+// Utility function
+//////////////////////////////////////////////////////////////
+double ConcChan::flux( double inConc, double outConc, double dt )
+{
+ return permeability_ * n_ * ( outConc - inConc ) * dt;
+ // Later add stuff for voltage dependence.
+ // Also use exp Euler
+}
diff --git a/kinetics/ConcChan.h b/kinetics/ConcChan.h
new file mode 100644
index 00000000..6ebf48b4
--- /dev/null
+++ b/kinetics/ConcChan.h
@@ -0,0 +1,52 @@
+/**********************************************************************
+** This program is part of 'MOOSE', the
+** Messaging Object Oriented Simulation Environment.
+** Copyright (C) 2018 Upinder S. Bhalla. and NCBS
+** It is made available under the terms of the
+** GNU Public License version 3 or later.
+** See the file COPYING.LIB for the full notice.
+**********************************************************************/
+
+#ifndef _CONC_CHAN_H
+#define _CONC_CHAN_H
+
+class ConcChan
+{
+ public:
+ ConcChan();
+ ~ConcChan();
+
+
+ //////////////////////////////////////////////////////////////////
+ // Field assignment stuff
+ //////////////////////////////////////////////////////////////////
+
+ void setPermeability( double v );
+ double getPermeability() const;
+ void setNumChan( double v );
+ double getNumChan() const;
+ double getFlux() const;
+
+ //////////////////////////////////////////////////////////////////
+ // Dest funcs
+ //////////////////////////////////////////////////////////////////
+
+ void process( const Eref& e, ProcPtr p );
+ void reinit( const Eref& e, ProcPtr p );
+ void inPool( double conc );
+ void outPool( double conc );
+ //////////////////////////////////////////////////////////////////
+ // Utility function
+ //////////////////////////////////////////////////////////////////
+ double flux( double inConc, double outConc, double dt );
+
+ static const Cinfo* initCinfo();
+ protected:
+ double permeability_; /// permeability in vol/(#.s) units.
+ double n_; /// Number of molecules of channel.
+ double flux_; /// Flux of molecule through channel
+ double charge_; /// Later: for including Nernst calculations
+ double Vm_; /// Later: for including Nernst calculations
+};
+
+#endif // CONC_CHAN_H
diff --git a/kinetics/EnzBase.cpp b/kinetics/EnzBase.cpp
index 0c4c136c..0812374c 100644
--- a/kinetics/EnzBase.cpp
+++ b/kinetics/EnzBase.cpp
@@ -63,6 +63,12 @@ const Cinfo* EnzBase::initCinfo()
&EnzBase::getNumSub
);
+ static ReadOnlyElementValueFinfo< EnzBase, unsigned int > numPrd(
+ "numProducts",
+ "Number of products in this MM reaction. Usually 1.",
+ &EnzBase::getNumPrd
+ );
+
//////////////////////////////////////////////////////////////
// Shared Msg Definitions
@@ -125,6 +131,7 @@ const Cinfo* EnzBase::initCinfo()
&numKm, // ElementValue
&kcat, // Value
&numSub, // ReadOnlyElementValue
+ &numPrd, // ReadOnlyElementValue
&enzDest, // DestFinfo
&sub, // SharedFinfo
&prd, // SharedFinfo
@@ -259,6 +266,14 @@ unsigned int EnzBase::getNumSub( const Eref& e ) const
return ( mfb->size() );
}
+unsigned int EnzBase::getNumPrd( const Eref& e ) const
+{
+ const vector< MsgFuncBinding >* mfb =
+ e.element()->getMsgAndFunc( prdOut()->getBindIndex() );
+ assert( mfb );
+ return ( mfb->size() );
+}
+
////////////////////////////////////////////////////////////////////////
// Zombie conversion routine to convert between Enz subclasses.
////////////////////////////////////////////////////////////////////////
diff --git a/kinetics/EnzBase.h b/kinetics/EnzBase.h
index 791e0532..bd5f37af 100644
--- a/kinetics/EnzBase.h
+++ b/kinetics/EnzBase.h
@@ -33,6 +33,7 @@ class EnzBase
// This doesn't need a virtual form, depends on standard msgs.
unsigned int getNumSub( const Eref& e ) const;
+ unsigned int getNumPrd( const Eref& e ) const;
//////////////////////////////////////////////////////////////////
// Virtual funcs for field assignment stuff
diff --git a/kinetics/ReadKkit.cpp b/kinetics/ReadKkit.cpp
index 2529919a..d3b5bd41 100644
--- a/kinetics/ReadKkit.cpp
+++ b/kinetics/ReadKkit.cpp
@@ -195,6 +195,7 @@ void makeSolverOnCompt( Shell* s, const vector< ObjId >& compts,
Field< Id >::set( stoich, "ksolve", ksolve );
Field< string >::set( stoich, "path", simpath );
}
+ /* Not needed now that we use xfer pools to cross compartments.
if ( stoichVec.size() == 2 ) {
SetGet1< Id >::set( stoichVec[1], "buildXreacs", stoichVec[0] );
}
@@ -206,6 +207,7 @@ void makeSolverOnCompt( Shell* s, const vector< ObjId >& compts,
i = stoichVec.begin(); i != stoichVec.end(); ++i ) {
SetGet0::set( *i, "filterXreacs" );
}
+ */
}
void setMethod( Shell* s, Id mgr, double simdt, double plotdt,
@@ -286,31 +288,9 @@ Id ReadKkit::read(
assignEnzCompartments();
assignMMenzCompartments();
- /*
- if ( moveOntoCompartment_ ) {
- assignPoolCompartments();
- assignReacCompartments();
- assignEnzCompartments();
- assignMMenzCompartments();
- }
- */
-
convertParametersToConcUnits();
- // s->doSetClock( 8, plotdt_ );
-
- // string plotpath = basePath_ + "/graphs/##[TYPE=Table]," + basePath_ + "/moregraphs/##[TYPE=Table]";
- // s->doUseClock( plotpath, "process", 8 );
-
setMethod( s, mgr, simdt_, plotdt_, method );
- /*
- if ( !moveOntoCompartment_ ) {
- assignPoolCompartments();
- assignReacCompartments();
- assignEnzCompartments();
- assignMMenzCompartments();
- }
- */
//Harsha: Storing solver and runtime at model level rather than model level
Id kinetics( basePath_+"/kinetics");
@@ -332,15 +312,6 @@ void ReadKkit::run()
shell_->doSetClock( 16, plotdt_ );
shell_->doSetClock( 17, plotdt_ );
shell_->doSetClock( 18, plotdt_ );
- /*
- string poolpath = basePath_ + "/kinetics/##[ISA=Pool]";
- string reacpath = basePath_ + "/kinetics/##[ISA!=Pool]";
- string plotpath = basePath_ + "/graphs/##[TYPE=Table]," +
- basePath_ + "/moregraphs/##[TYPE=Table]";
- shell_->doUseClock( reacpath, "process", 4 );
- shell_->doUseClock( poolpath, "process", 5 );
- shell_->doUseClock( plotpath, "process", 8 );
- */
shell_->doReinit();
if ( useVariableDt_ ) {
shell_->doSetClock( 11, fastdt_ );
@@ -540,16 +511,33 @@ string ReadKkit::cleanPath( const string& path ) const
// which later created a problem as the same name exist in moose when minus
// was replaced with underscore.
//So replacing minus with _minus_ like I do in SBML
- string ret = path;
+ size_t Iindex = 0;
+ string cleanDigit="/";
+ while(true)
+ {
+ size_t sindex = path.find('/',Iindex+1);
+ if (sindex == string::npos)
+ { if (isdigit((path.substr(Iindex+1,sindex-Iindex-1))[0]) )
+ cleanDigit += '_' ;
+ cleanDigit += path.substr(Iindex+1,sindex-Iindex-1);
+ break;
+ }
+ if (isdigit((path.substr(Iindex+1,sindex-Iindex-1))[0]))
+ cleanDigit+='_';
+ cleanDigit += path.substr(Iindex+1,sindex-Iindex);
+ Iindex = sindex;
+ }
+ string ret = cleanDigit;
+ //string ret = path;
string cleanString;
- for ( unsigned int i = 0; i < path.length(); ++i ) {
+ for ( unsigned int i = 0; i < cleanDigit.length(); ++i ) {
char c = ret[i];
if ( c == '*' )
cleanString += "_p";
else if ( c == '[' || c == ']' || c == '@' || c == ' ')
cleanString += '_';
else if (c == '-')
- cleanString += "_dash_";
+ cleanString += "_";
else
cleanString += c;
}
@@ -563,8 +551,7 @@ void assignArgs( map< string, int >& argConv, const vector< string >& args )
}
void ReadKkit::objdump( const vector< string >& args)
-{
- if ( args[1] == "kpool" )
+{ if ( args[1] == "kpool" )
assignArgs( poolMap_, args );
else if ( args[1] == "kreac" )
assignArgs( reacMap_, args );
@@ -619,8 +606,7 @@ void ReadKkit::textload( const vector< string >& args)
}
void ReadKkit::undump( const vector< string >& args)
-{
- if ( args[1] == "kpool" )
+{ if ( args[1] == "kpool" )
buildPool( args );
else if ( args[1] == "kreac" )
buildReac( args );
@@ -973,7 +959,6 @@ Id ReadKkit::buildInfo( Id parent,
map< string, int >& m, const vector< string >& args )
{
Id info = shell_->doCreate( "Annotator", parent, "info", 1 );
- //cout << "parent " << parent << " " << parent.path() << " info " << info.path();
assert( info != Id() );
double x = atof( args[ m[ "x" ] ].c_str() );
@@ -991,7 +976,6 @@ Id ReadKkit::buildGroup( const vector< string >& args )
{
string head;
string tail = pathTail( cleanPath( args[2] ), head );
-
Id pa = shell_->doFind( head ).id;
assert( pa != Id() );
Id group = shell_->doCreate( "Neutral", pa, tail, 1 );
@@ -1012,7 +996,7 @@ Id ReadKkit::buildGroup( const vector< string >& args )
* traffic are originally framed in terms of number of receptors, not conc.
*/
Id ReadKkit::buildPool( const vector< string >& args )
-{
+{
string head;
string clean = cleanPath( args[2] );
string tail = pathTail( clean, head );
@@ -1060,7 +1044,6 @@ Id ReadKkit::buildPool( const vector< string >& args )
poolVols_[pool] = vol;
Id info = buildInfo( pool, poolMap_, args );
-
/*
cout << setw( 20 ) << head << setw( 15 ) << tail << " " <<
setw( 12 ) << nInit << " " <<
diff --git a/ksolve/Gsolve.cpp b/ksolve/Gsolve.cpp
index 42bee6f7..25a5c060 100644
--- a/ksolve/Gsolve.cpp
+++ b/ksolve/Gsolve.cpp
@@ -33,23 +33,6 @@
const unsigned int OFFNODE = ~0;
-// static function
-SrcFinfo2< Id, vector< double > >* Gsolve::xComptOut()
-{
- static SrcFinfo2< Id, vector< double > > xComptOut( "xComptOut",
- "Sends 'n' of all molecules participating in cross-compartment "
- "reactions between any juxtaposed voxels between current compt "
- "and another compartment. This includes molecules local to this "
- "compartment, as well as proxy molecules belonging elsewhere. "
- "A(t+1) = (Alocal(t+1) + AremoteProxy(t+1)) - Alocal(t) "
- "A(t+1) = (Aremote(t+1) + Aproxy(t+1)) - Aproxy(t) "
- "Then we update A on the respective solvers with: "
- "Alocal(t+1) = Aproxy(t+1) = A(t+1) "
- "This is equivalent to sending dA over on each timestep. "
- );
- return &xComptOut;
-}
-
const Cinfo* Gsolve::initCinfo()
{
///////////////////////////////////////////////////////
@@ -161,11 +144,6 @@ const Cinfo* Gsolve::initCinfo()
"Handles initReinit call from Clock",
new ProcOpFunc< Gsolve >( &Gsolve::initReinit ) );
- static DestFinfo xComptIn( "xComptIn",
- "Handles arriving pool 'n' values used in cross-compartment "
- "reactions.",
- new EpFunc2< Gsolve, Id, vector< double > >( &Gsolve::xComptIn )
- );
///////////////////////////////////////////////////////
// Shared definitions
///////////////////////////////////////////////////////
@@ -188,16 +166,6 @@ const Cinfo* Gsolve::initCinfo()
initShared, sizeof( initShared ) / sizeof( const Finfo* )
);
- static Finfo* xComptShared[] =
- {
- xComptOut(), &xComptIn
- };
- static SharedFinfo xCompt( "xCompt",
- "Shared message for pool exchange for cross-compartment "
- "reactions. Exchanges latest values of all pools that "
- "participate in such reactions.",
- xComptShared, sizeof( xComptShared ) / sizeof( const Finfo* )
- );
///////////////////////////////////////////////////////
static Finfo* gsolveFinfos[] =
@@ -210,7 +178,6 @@ const Cinfo* Gsolve::initCinfo()
&voxelVol, // DestFinfo
&proc, // SharedFinfo
&init, // SharedFinfo
- &xCompt, // SharedFinfo
// Here we put new fields that were not there in the Ksolve.
&useRandInit, // Value
&useClockedUpdate, // Value
@@ -438,6 +405,7 @@ void Gsolve::process( const Eref& e, ProcPtr p )
dvalues[2] = 0;
dvalues[3] = stoichPtr_->getNumVarPools();
dsolvePtr_->getBlock( dvalues );
+ dsolvePtr_->setPrev();
// Here we need to convert to integers, just in case. Normally
// one would use a stochastic (integral) diffusion method with
@@ -459,31 +427,10 @@ void Gsolve::process( const Eref& e, ProcPtr p )
}
setBlock( dvalues );
}
- // Second, take the arrived xCompt reac values and update S with them.
- // Here the roundoff issues are handled by the GssaVoxelPools functions
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- // cout << xfer_.size() << " " << xf.xferVoxel.size() << endl;
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferIn( xf, j, &sys_ );
- }
- }
- // Third, record the current value of pools as the reference for the
- // next cycle.
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- }
- }
- // Fourth: Fix the rates if we have had any diffusion or xreacs
+ // Third: Fix the rates if we have had any diffusion or xreacs
// happening. This is very inefficient at this point, need to fix.
- if ( dsolvePtr_ || xfer_.size() > 0 )
+ if ( dsolvePtr_ )
{
for ( vector< GssaVoxelPools >::iterator
i = pools_.begin(); i != pools_.end(); ++i )
@@ -491,7 +438,7 @@ void Gsolve::process( const Eref& e, ProcPtr p )
i->refreshAtot( &sys_ );
}
}
- // Fifth, update the mol #s.
+ // Fourth, update the mol #s.
// First we advance the simulation.
size_t nvPools = pools_.size( );
@@ -537,6 +484,11 @@ void Gsolve::process( const Eref& e, ProcPtr p )
kvalues[3] = stoichPtr_->getNumVarPools();
getBlock( kvalues );
dsolvePtr_->setBlock( kvalues );
+
+ // Now use the values in the Dsolve to update junction fluxes
+ // for diffusion, channels, and xreacs
+ dsolvePtr_->updateJunctions( p->dt );
+ // Here the Gsolve may need to do something to convert to integers
}
}
@@ -552,30 +504,7 @@ void Gsolve::reinit( const Eref& e, ProcPtr p )
{
i->reinit( &sys_ );
}
-
- // Second, take the arrived xCompt reac values and update S with them.
- // Here the roundoff issues are handled by the GssaVoxelPools functions
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- const XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferInOnlyProxies(
- xf.xferPoolIdx, xf.values,
- stoichPtr_->getNumProxyPools(), j );
- }
- }
- // Third, record the current value of pools as the reference for the
- // next cycle.
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- }
- }
- // Fourth, update the atots.
+ // Second, update the atots.
for ( vector< GssaVoxelPools >::iterator
i = pools_.begin(); i != pools_.end(); ++i )
{
@@ -593,24 +522,7 @@ void Gsolve::reinit( const Eref& e, ProcPtr p )
// init operations.
//////////////////////////////////////////////////////////////
void Gsolve::initProc( const Eref& e, ProcPtr p )
-{
- if ( !stoichPtr_ )
- return;
- // vector< vector< double > > values( xfer_.size() );
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
- // values[i].resize( size, 0.0 );
- vector< double > values( size, 0.0 );
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- unsigned int vox = xf.xferVoxel[j];
- pools_[vox].xferOut( j, values, xf.xferPoolIdx );
- }
- xComptOut()->sendTo( e, xf.ksolve, e.id(), values );
- }
-}
+{;}
void Gsolve::initReinit( const Eref& e, ProcPtr p )
{
@@ -620,20 +532,6 @@ void Gsolve::initReinit( const Eref& e, ProcPtr p )
{
pools_[i].reinit( &sys_ );
}
- // vector< vector< double > > values( xfer_.size() );
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
- xf.lastValues.assign( size, 0.0 );
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- unsigned int vox = xf.xferVoxel[j];
- pools_[ vox ].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- // values[i] = xf.lastValues;
- }
- xComptOut()->sendTo( e, xf.ksolve, e.id(), xf.lastValues );
- }
}
//////////////////////////////////////////////////////////////
// Solver setup
@@ -1118,7 +1016,6 @@ void Gsolve::updateVoxelVol( vector< double > vols )
{
pools_[i].setVolumeAndDependencies( vols[i] );
}
- stoichPtr_->setupCrossSolverReacVols();
updateRateTerms( ~0U );
}
}
@@ -1130,7 +1027,6 @@ void Gsolve::updateRateTerms( unsigned int index )
// unsigned int numCrossRates = stoichPtr_->getNumRates() - stoichPtr_->getNumCoreRates();
for ( unsigned int i = 0 ; i < pools_.size(); ++i )
{
- // pools_[i].resetXreacScale( numCrossRates );
pools_[i].updateAllRateTerms( stoichPtr_->getRateTerms(),
stoichPtr_->getNumCoreRates() );
}
diff --git a/ksolve/Gsolve.h b/ksolve/Gsolve.h
index 9cf59a00..78f572ab 100644
--- a/ksolve/Gsolve.h
+++ b/ksolve/Gsolve.h
@@ -129,7 +129,6 @@ public:
#endif
//////////////////////////////////////////////////////////////////
- static SrcFinfo2< Id, vector< double > >* xComptOut();
static const Cinfo* initCinfo();
private:
diff --git a/ksolve/GssaVoxelPools.cpp b/ksolve/GssaVoxelPools.cpp
index 438652d8..27396430 100644
--- a/ksolve/GssaVoxelPools.cpp
+++ b/ksolve/GssaVoxelPools.cpp
@@ -343,7 +343,6 @@ void GssaVoxelPools::setStoich( const Stoich* stoichPtr )
void GssaVoxelPools::setVolumeAndDependencies( double vol )
{
VoxelPoolsBase::setVolumeAndDependencies( vol );
- stoichPtr_->setupCrossSolverReacVols();
updateAllRateTerms( stoichPtr_->getRateTerms(),
stoichPtr_->getNumCoreRates() );
}
@@ -415,30 +414,3 @@ void GssaVoxelPools::xferIn( XferInfo& xf,
// Does this fix the problem of negative concs?
refreshAtot( g );
}
-
-void GssaVoxelPools::xferInOnlyProxies(
- const vector< unsigned int >& poolIndex,
- const vector< double >& values,
- unsigned int numProxyPools,
- unsigned int voxelIndex )
-{
- unsigned int offset = voxelIndex * poolIndex.size();
- vector< double >::const_iterator i = values.begin() + offset;
- unsigned int proxyEndIndex = stoichPtr_->getNumVarPools() +
- stoichPtr_->getNumProxyPools();
- for ( vector< unsigned int >::const_iterator
- k = poolIndex.begin(); k != poolIndex.end(); ++k )
- {
- // if ( *k >= size() - numProxyPools )
- if ( *k >= stoichPtr_->getNumVarPools() && *k < proxyEndIndex )
- {
- double base = floor( *i );
- if ( rng_.uniform() >= (*i - base) )
- varSinit()[*k] = (varS()[*k] += base );
- else
- varSinit()[*k] = (varS()[*k] += base + 1.0 );
- // varSinit()[*k] = (varS()[*k] += round( *i ));
- }
- i++;
- }
-}
diff --git a/ksolve/GssaVoxelPools.h b/ksolve/GssaVoxelPools.h
index 1d3a3dba..716698a4 100644
--- a/ksolve/GssaVoxelPools.h
+++ b/ksolve/GssaVoxelPools.h
@@ -70,17 +70,6 @@ public:
void xferIn( XferInfo& xf,
unsigned int voxelIndex, const GssaSystem* g );
- /**
- * Used during initialization: Takes only the proxy pool values
- * from the incoming transfer data, and assigns it to the proxy
- * pools on current solver
- */
- void xferInOnlyProxies(
- const vector< unsigned int >& poolIndex,
- const vector< double >& values,
- unsigned int numProxyPools,
- unsigned int voxelIndex );
-
void setStoich( const Stoich* stoichPtr );
private:
diff --git a/ksolve/Ksolve.cpp b/ksolve/Ksolve.cpp
index a7ee4c7a..8b7b5ac9 100644
--- a/ksolve/Ksolve.cpp
+++ b/ksolve/Ksolve.cpp
@@ -17,7 +17,6 @@
#include "VoxelPoolsBase.h"
#include "VoxelPools.h"
#include "../mesh/VoxelJunction.h"
-#include "XferInfo.h"
#include "ZombiePoolInterface.h"
#include "RateTerm.h"
@@ -38,23 +37,6 @@
const unsigned int OFFNODE = ~0;
-// static function
-SrcFinfo2< Id, vector< double > >* Ksolve::xComptOut()
-{
- static SrcFinfo2< Id, vector< double > > xComptOut( "xComptOut",
- "Sends 'n' of all molecules participating in cross-compartment "
- "reactions between any juxtaposed voxels between current compt "
- "and another compartment. This includes molecules local to this "
- "compartment, as well as proxy molecules belonging elsewhere. "
- "A(t+1) = (Alocal(t+1) + AremoteProxy(t+1)) - Alocal(t) "
- "A(t+1) = (Aremote(t+1) + Aproxy(t+1)) - Aproxy(t) "
- "Then we update A on the respective solvers with: "
- "Alocal(t+1) = Aproxy(t+1) = A(t+1) "
- "This is equivalent to sending dA over on each timestep. "
- );
- return &xComptOut;
-}
-
const Cinfo* Ksolve::initCinfo()
{
///////////////////////////////////////////////////////
@@ -193,22 +175,6 @@ const Cinfo* Ksolve::initCinfo()
initShared, sizeof( initShared ) / sizeof( const Finfo* )
);
- static DestFinfo xComptIn( "xComptIn",
- "Handles arriving pool 'n' values used in cross-compartment "
- "reactions.",
- new EpFunc2< Ksolve, Id, vector< double > >( &Ksolve::xComptIn )
- );
- static Finfo* xComptShared[] =
- {
- xComptOut(), &xComptIn
- };
- static SharedFinfo xCompt( "xCompt",
- "Shared message for pool exchange for cross-compartment "
- "reactions. Exchanges latest values of all pools that "
- "participate in such reactions.",
- xComptShared, sizeof( xComptShared ) / sizeof( const Finfo* )
- );
-
static Finfo* ksolveFinfos[] =
{
&method, // Value
@@ -225,7 +191,6 @@ const Cinfo* Ksolve::initCinfo()
&estimatedDt, // ReadOnlyValue
&stoich, // ReadOnlyValue
&voxelVol, // DestFinfo
- &xCompt, // SharedFinfo
&proc, // SharedFinfo
&init, // SharedFinfo
};
@@ -542,47 +507,17 @@ void Ksolve::process( const Eref& e, ProcPtr p )
dvalues[1] = getNumLocalVoxels();
dvalues[2] = 0;
dvalues[3] = stoichPtr_->getNumVarPools();
- dsolvePtr_->getBlock( dvalues );
-
- /*
- vector< double >::iterator i = dvalues.begin() + 4;
- for ( ; i != dvalues.end(); ++i )
- cout << *i << " " << round( *i ) << endl;
- getBlock( kvalues );
- vector< double >::iterator d = dvalues.begin() + 4;
- for ( vector< double >::iterator
- k = kvalues.begin() + 4; k != kvalues.end(); ++k )
- *k++ = ( *k + *d )/2.0
- setBlock( kvalues );
- */
- setBlock( dvalues );
- }
-
- // Second, take the arrived xCompt reac values and update S with them.
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- const XferInfo& xf = xfer_[i];
- // cout << xfer_.size() << " " << xf.xferVoxel.size() << endl;
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferIn(
- xf.xferPoolIdx, xf.values, xf.lastValues, j );
- }
- }
- // Third, record the current value of pools as the reference for the
- // next cycle.
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
+ dsolvePtr_->getBlock( dvalues );
+ // Second, set the prev_ value in DiffPoolVec
+ dsolvePtr_->setPrev();
+ setBlock( dvalues );
}
size_t nvPools = pools_.size( );
#ifdef PARALLELIZE_KSOLVE_WITH_CPP11_ASYNC
- // Fourth, do the numerical integration for all reactions.
+ // Third, do the numerical integration for all reactions.
size_t grainSize = min( nvPools, 1 + (nvPools / numThreads_ ) );
size_t nWorkers = nvPools / grainSize;
@@ -615,7 +550,7 @@ void Ksolve::process( const Eref& e, ProcPtr p )
#endif
- // Finally, assemble and send the integrated values off for the Dsolve.
+ // Assemble and send the integrated values off for the Dsolve.
if ( dsolvePtr_ )
{
vector< double > kvalues( 4 );
@@ -625,6 +560,10 @@ void Ksolve::process( const Eref& e, ProcPtr p )
kvalues[3] = stoichPtr_->getNumVarPools();
getBlock( kvalues );
dsolvePtr_->setBlock( kvalues );
+
+ // Now use the values in the Dsolve to update junction fluxes
+ // for diffusion, channels, and xreacs
+ dsolvePtr_->updateJunctions( p->dt );
}
}
@@ -674,28 +613,6 @@ void Ksolve::reinit( const Eref& e, ProcPtr p )
return;
}
- // cout << "************************* path = " << e.id().path() << endl;
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- const XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferInOnlyProxies(
- xf.xferPoolIdx, xf.values,
- stoichPtr_->getNumProxyPools(),
- j );
- }
- }
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- pools_[xf.xferVoxel[j]].xferOut(
- j, xf.lastValues, xf.xferPoolIdx );
- }
- }
-
#if PARALLELIZE_KSOLVE_WITH_CPP11_ASYNC
if( 1 < getNumThreads( ) )
cout << "Debug: User wants Ksolve with " << numThreads_ << " threads" << endl;
@@ -708,40 +625,12 @@ void Ksolve::reinit( const Eref& e, ProcPtr p )
//////////////////////////////////////////////////////////////
void Ksolve::initProc( const Eref& e, ProcPtr p )
{
- // vector< vector< double > > values( xfer_.size() );
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
- // values[i].resize( size, 0.0 );
- vector< double > values( size, 0.0 );
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- unsigned int vox = xf.xferVoxel[j];
- pools_[vox].xferOut( j, values, xf.xferPoolIdx );
- }
- xComptOut()->sendTo( e, xf.ksolve, e.id(), values );
- }
- // xComptOut()->sendVec( e, values );
}
void Ksolve::initReinit( const Eref& e, ProcPtr p )
{
for ( unsigned int i = 0 ; i < pools_.size(); ++i )
pools_[i].reinit( p->dt );
-
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- XferInfo& xf = xfer_[i];
- unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
- xf.lastValues.assign( size, 0.0 );
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
- {
- unsigned int vox = xf.xferVoxel[j];
- pools_[ vox ].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- }
- xComptOut()->sendTo( e, xf.ksolve, e.id(), xf.lastValues );
- }
}
/**
@@ -919,7 +808,6 @@ void Ksolve::updateVoxelVol( vector< double > vols )
{
pools_[i].setVolumeAndDependencies( vols[i] );
}
- stoichPtr_->setupCrossSolverReacVols();
updateRateTerms( ~0U );
}
}
@@ -944,25 +832,5 @@ void Ksolve::print() const
cout << "method = " << method_ << ", stoich=" << stoich_.path() <<endl;
cout << "dsolve = " << dsolve_.path() << endl;
cout << "compartment = " << compartment_.path() << endl;
- cout << "xfer summary: numxfer = " << xfer_.size() << "\n";
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- cout << "xfer_[" << i << "] numValues=" <<
- xfer_[i].values.size() <<
- ", xferPoolIdx.size = " << xfer_[i].xferPoolIdx.size() <<
- ", xferVoxel.size = " << xfer_[i].xferVoxel.size() << endl;
- }
- cout << "xfer details:\n";
- for ( unsigned int i = 0; i < xfer_.size(); ++i )
- {
- cout << "xfer_[" << i << "] xferPoolIdx=\n";
- const vector< unsigned int>& xi = xfer_[i].xferPoolIdx;
- for ( unsigned int j = 0; j << xi.size(); ++j )
- cout << " " << xi[j];
- cout << "\nxfer_[" << i << "] xferVoxel=\n";
- const vector< unsigned int>& xv = xfer_[i].xferVoxel;
- for ( unsigned int j = 0; j << xv.size(); ++j )
- cout << " " << xv[j];
- }
}
diff --git a/ksolve/Ksolve.h b/ksolve/Ksolve.h
index 5e2153be..c7f4af8a 100644
--- a/ksolve/Ksolve.h
+++ b/ksolve/Ksolve.h
@@ -126,30 +126,11 @@ public:
*/
void updateRateTerms( unsigned int index );
- //////////////////////////////////////////////////////////////////
- // Functions for cross-compartment transfer
- //////////////////////////////////////////////////////////////////
- void setupXfer( Id myKsolve, Id otherKsolve,
- unsigned int numProxyMols,
- const vector< VoxelJunction >& vj );
-
- void assignXferIndex( unsigned int numProxyMols,
- unsigned int xferCompt,
- const vector< vector< unsigned int > >& voxy );
-
- void assignXferVoxels( unsigned int xferCompt );
-
- unsigned int assignProxyPools( const map< Id, vector< Id > >& xr,
- Id myKsolve, Id otherKsolve, Id otherComptId );
-
- void buildCrossReacVolScaling( Id otherKsolve,
- const vector< VoxelJunction >& vj );
//////////////////////////////////////////////////////////////////
// for debugging
void print() const;
//////////////////////////////////////////////////////////////////
- static SrcFinfo2< Id, vector< double > >* xComptOut();
static const Cinfo* initCinfo();
private:
diff --git a/ksolve/Stoich.cpp b/ksolve/Stoich.cpp
index a94dfff9..45815140 100644
--- a/ksolve/Stoich.cpp
+++ b/ksolve/Stoich.cpp
@@ -141,7 +141,7 @@ const Cinfo* Stoich::initCinfo()
"matrixEntry",
"The non-zero matrix entries in the sparse matrix. Their"
"column indices are in a separate vector and the row"
- "informatino in a third",
+ "information in a third",
&Stoich::getMatrixEntry
);
// Stuff here for getting Stoichiometry matrix to manipulate in
@@ -195,22 +195,6 @@ const Cinfo* Stoich::initCinfo()
"Restore all zombies to their native state",
new OpFunc0< Stoich >( &Stoich::unZombifyModel )
);
- static DestFinfo buildXreacs( "buildXreacs",
- "Build cross-reaction terms between current stoich and "
- "argument. This function scans the voxels at which there are "
- "junctions between different compartments, and orchestrates "
- "set up of interfaces between the Ksolves that implement "
- "the X reacs at those junctions. ",
- new EpFunc1< Stoich, Id >( &Stoich::buildXreacs )
- );
-
- static DestFinfo filterXreacs( "filterXreacs",
- "Filter cross-reaction terms on current stoich"
- "This function clears out absent rate terms that would "
- "otherwise try to compute cross reactions where the "
- "junctions are not present. ",
- new OpFunc0< Stoich >( &Stoich::filterXreacs )
- );
static DestFinfo scaleBufsAndRates( "scaleBufsAndRates",
"Args: voxel#, volRatio\n"
"Handles requests for runtime volume changes in the specified "
@@ -246,8 +230,6 @@ const Cinfo* Stoich::initCinfo()
&proxyPools, // ReadOnlyLookupValue
&status, // ReadOnlyLookupValue
&unzombify, // DestFinfo
- &buildXreacs, // DestFinfo
- &filterXreacs, // DestFinfo
&scaleBufsAndRates, // DestFinfo
};
@@ -405,7 +387,6 @@ void Stoich::setElist( const Eref& e, const vector< ObjId >& elist )
status_ = 16;
return;
}
- locateOffSolverReacs( compartment_, temp );
// allocateObjMap( temp );
allocateModel( temp );
@@ -426,55 +407,10 @@ void Stoich::setElist( const Eref& e, const vector< ObjId >& elist )
if ( kinterface_ )
{
kinterface_->setDsolve( dsolve_ );
- // kinterface_->setupCrossSolverReacs( offSolverPoolVec_ );
- kinterface_->setupCrossSolverReacVols( subComptVec_, prdComptVec_);
kinterface_->updateRateTerms();
}
}
-void Stoich::setupCrossSolverReacVols() const
-{
- kinterface_->setupCrossSolverReacVols( subComptVec_, prdComptVec_);
-}
-
-/*
-// Utility function for getting inputs to a single Function object.
-static vector< ObjId > inputsToFunc( Id func )
-{
- static const Finfo* inputFinfo =
- Cinfo::find( "Variable" )->findFinfo("input");
- static const DestFinfo* df = dynamic_cast< const DestFinfo* >( inputFinfo );
- assert( df );
- static FuncId fid = df->funcId();
- assert( func.element()->cinfo()->isA( "Function" );
- unsigned int numInputs = Field< unsigned int >::get( func, "numVars" );
- Id varId( func.value() );
- vector< ObjId > caller;
- varId.element()->getInputMsgs( caller, fid );
- cout << "NUM Inputs = " << numInputs << ", numObjId = " << caller.size() << endl;
- return caller;
-}
-
-// routine to get all inputs (pool Ids) to Pool Functions
-void Stoich::inputsToPoolFuncs(
- vector< pair< Id, vector< unsigned int> >& ret ) const
-{
- ret.clear();
- vector< Id >::iterator i;
- vector< ObjId >::iterator j;
- for ( i = poolFuncVec_.begin(); i != poolFuncVec_.end(); ++i )
- vector< ObjId > ovec = inputsToFunc( *i );
- vector< unsigned int > poolIndex;
- poolIndex.reserve( ovec.size() );
- for ( j = ovec.begin(); j != ovec.end(); ++j ) {
- unsigned int poolIndex =
- poolIndex.push_back( convertIdToPoolIndex( j->id ) );
- }
- ret.push_back( pair( *i, poolIndex) );
- return ret;
-}
-*/
-
//////////////////////////////////////////////////////////////////////
string Stoich::getPath( const Eref& e ) const
@@ -1276,48 +1212,6 @@ const KinSparseMatrix& Stoich::getStoichiometryMatrix() const
return N_;
}
-void Stoich::buildXreacs( const Eref& e, Id otherStoich )
-{
- if ( status_ == 0 )
- kinterface_->setupCrossSolverReacs( offSolverPoolMap_,otherStoich);
-}
-
-void Stoich::filterXreacs()
-{
- if ( status_ == 0 ) {
- kinterface_->filterCrossRateTerms( offSolverReacVec_, offSolverReacCompts_ );
- kinterface_->filterCrossRateTerms( offSolverEnzVec_, offSolverEnzCompts_ );
- kinterface_->filterCrossRateTerms( offSolverMMenzVec_, offSolverMMenzCompts_ );
- }
-}
-
-/*
-void Stoich::comptsOnCrossReacTerms( vector< pair< Id, Id > >& xr ) const
-{
- xr.clear();
- assert( offSolverReacVec_.size() == offSolverReacCompts_.size() );
- assert( offSolverEnzVec_.size() == offSolverEnzCompts_.size() );
- assert( offSolverMMenzVec_.size() == offSolverMMenzCompts_.size() );
- unsigned int numOffSolverRates = getNumRates() - getNumCoreRates();
- if ( numOffSolverRates == 0 )
- return;
- unsigned int k = getNumCoreRates();
- unsigned int j = 0;
-
- for ( unsigned int i = 0; i < offSolverReacs_.size() ; ++i ) {
- if ( i+1 < offSolverReacs_.size() )
- j = convertIdToReacIndex( offSolverReacs_[i+1] );
- else
- j = rates_.size();
- while ( k < j ) {
- xr.push_back( offSolverReacCompts_[i] );
- k++;
- }
- }
- assert( xr.size() == numXrates );
-}
-*/
-
//////////////////////////////////////////////////////////////
// Model zombification functions
//////////////////////////////////////////////////////////////
@@ -1334,31 +1228,6 @@ static Id findFuncMsgSrc( Id pa, const string& msg )
return ret[0];
}
return Id(); // failure
-
-
- /*
- cout << "fundFuncMsgSrc:: convert to DestFinfo\n";
- const DestFinfo *df = dynamic_cast< const DestFinfo* >( finfo );
- if ( !df ) return Id();
- FuncId fid = df->getFid();
-
- cout << "fundFuncMsgSrc:: findCaller\n";
- ObjId caller = pa.element()->findCaller( fid );
- if ( caller.bad() )
- return Id();
- cout << "Found Func for " << pa.path() << " on " << msg << endl;
- cout << "Func was " << caller.path() << endl;
- return caller.id;
- */
-
- /*
- vector< Id > kids = Neutral::getChildren( pa.eref() );
- for ( vector< Id >::iterator i = kids.begin(); i != kids.end(); ++i )
- if ( i->element()->cinfo()->isA( "Function" ) )
- return( *i );
- return Id();
- */
-
}
Id Stoich::zombifyPoolFuncWithScaling( Id pool )
diff --git a/ksolve/Stoich.h b/ksolve/Stoich.h
index ea517709..bea3d067 100644
--- a/ksolve/Stoich.h
+++ b/ksolve/Stoich.h
@@ -198,39 +198,9 @@ class Stoich
*/
void convertRatesToStochasticForm();
- /**
- * Builds cross-reaction terms between current stoich and
- * otherStoich. The function scans the voxels at which there
- * are jucntions between different compartments, and orchestrates
- * setup of interfaces between the Ksolves that implement the
- * cross-reactions at these junctions.
- */
- void buildXreacs( const Eref& e, Id otherStoich );
-
- void filterXreacs();
-
/// Used to handle run-time size updates for spines.
void scaleBufsAndRates( unsigned int index, double volScale );
- /**
- * Expands out list of compartment mappings of proxy reactions to
- * the appropriate entries on the rates_vector.
- void comptsOnCrossReacTerms( vector< pair< Id, Id > >& xr ) const;
- */
-
- /**
- * Used to set up and update all cross solver reac terms and
- * their rates, if there has been a change in volumes. Should
- * also work if there has been a change in voxelization.
- */
- void setupCrossSolverReacVols() const;
-
- /**
- * Finds all the input molecules contributing to any of the
- * Function cases: poolFunc, incrementFunc or reacFunc
- void inputsToPoolFuncs(
- vector< pair< Id, vector< unsigned int > > >& ret ) const;
- */
//////////////////////////////////////////////////////////////////
// Zombification functions.
//////////////////////////////////////////////////////////////////
diff --git a/ksolve/VoxelPools.cpp b/ksolve/VoxelPools.cpp
index 88405e72..64b07be2 100644
--- a/ksolve/VoxelPools.cpp
+++ b/ksolve/VoxelPools.cpp
@@ -371,60 +371,7 @@ void VoxelPools::print() const
void VoxelPools::setVolumeAndDependencies( double vol )
{
VoxelPoolsBase::setVolumeAndDependencies( vol );
- stoichPtr_->setupCrossSolverReacVols();
updateAllRateTerms( stoichPtr_->getRateTerms(),
stoichPtr_->getNumCoreRates() );
}
-
-////////////////////////////////////////////////////////////
-#if 0
-/**
- * Zeroes out rate terms that are involved in cross-reactions that
- * are not present on current voxel.
- */
-void VoxelPools::filterCrossRateTerms(
- const vector< pair< Id, Id > >&
- offSolverReacCompts )
-{
- /*
-From VoxelPoolsBase:proxyPoolVoxels[comptIndex][#] we know
-if specified compt has local proxies.
- Note that compt is identified by an index, and actually looks up
- the Ksolve.
-From Ksolve::compartment_ we know which compartment a given ksolve belongs
- in
-From Ksolve::xfer_[otherKsolveIndex].ksolve we have the id of the other
- Ksolves.
-From Stoich::offSolverReacCompts_ which is pair< Id, Id > we have the
- ids of the _compartments_ feeding into the specified rateTerms.
-
-Somewhere I need to make a map of compts to comptIndex.
-
-The ordering of the xfer vector is simply by the order of the script call
-for buildXfer.
-
-This has become too ugly
-Skip the proxyPoolVoxels info, or use the comptIndex here itself to
-build the table.
-comptIndex looks up xfer which holds the Ksolve Id. From that we can
-get the compt id. All this relies on this mapping being correct.
-Or I should pass in the compt when I build it.
-
-OK, now we have VoxelPoolsBase::proxyPoolCompts_ vector to match the
-comptIndex.
-
-*/
- unsigned int numCoreRates = stoichPtr_->getNumCoreRates();
- for ( unsigned int i = 0; i < offSolverReacCompts.size(); ++i ) {
- const pair< Id, Id >& p = offSolverReacCompts[i];
- if ( !isVoxelJunctionPresent( p.first, p.second) ) {
- unsigned int k = i + numCoreRates;
- assert( k < rates_.size() );
- if ( rates_[k] )
- delete rates_[k];
- rates_[k] = new ExternReac;
- }
- }
-}
-#endif
diff --git a/ksolve/VoxelPools.h b/ksolve/VoxelPools.h
index 05129e36..bd1e5ec5 100644
--- a/ksolve/VoxelPools.h
+++ b/ksolve/VoxelPools.h
@@ -95,11 +95,6 @@ public:
void updateReacVelocities(
const double* s, vector< double >& v ) const;
- /**
- * Changes cross rate terms to zero if there is no junction
- void filterCrossRateTerms( const vector< pair< Id, Id > >& vec );
- */
-
/// Used for debugging.
void print() const;
private:
diff --git a/ksolve/ZombiePoolInterface.cpp b/ksolve/ZombiePoolInterface.cpp
index 20fb4895..b75254d2 100644
--- a/ksolve/ZombiePoolInterface.cpp
+++ b/ksolve/ZombiePoolInterface.cpp
@@ -34,306 +34,12 @@ ZombiePoolInterface::ZombiePoolInterface()
isBuilt_( false )
{;}
-//////////////////////////////////////////////////////////////////////////
-// cross-compartment reaction stuff.
-//////////////////////////////////////////////////////////////////////////
-// void ZombiePoolInterface::xComptIn( const Eref& e, const ObjId& src,
-// vector< double > values )
-void ZombiePoolInterface::xComptIn( const Eref& e, Id srcZombiePoolInterface,
- vector< double > values )
-{
- // Identify the xfer_ that maps to the srcZombiePoolInterface. Assume only a small
- // number of them, otherwise we should use a map.
- unsigned int comptIdx ;
- for ( comptIdx = 0 ; comptIdx < xfer_.size(); ++comptIdx ) {
- if ( xfer_[comptIdx].ksolve == srcZombiePoolInterface ) break;
- }
- assert( comptIdx != xfer_.size() );
- XferInfo& xf = xfer_[comptIdx];
- // assert( values.size() == xf.values.size() );
- xf.values = values;
-// xfer_[comptIdx].lastValues = values;
-}
+void ZombiePoolInterface::updateJunctions( double dt )
+{;}
+void ZombiePoolInterface::setPrev()
+{;}
/////////////////////////////////////////////////////////////////////
-// Functions for setup of cross-compartment transfer.
-/////////////////////////////////////////////////////////////////////
-/**
- * Figures out which voxels are involved in cross-compt reactions. Stores
- * in the appropriate xfer_ entry.
- */
-void ZombiePoolInterface::assignXferVoxels( unsigned int xferCompt )
-{
- assert( xferCompt < xfer_.size() );
- XferInfo& xf = xfer_[xferCompt];
- for ( unsigned int i = 0; i < getNumLocalVoxels(); ++i ) {
- if ( pools(i)->hasXfer( xferCompt ) )
- xf.xferVoxel.push_back( i );
- }
- xf.values.resize( xf.xferVoxel.size() * xf.xferPoolIdx.size(), 0 );
- xf.lastValues.resize( xf.xferVoxel.size() * xf.xferPoolIdx.size(), 0 );
- xf.subzero.resize( xf.xferVoxel.size() * xf.xferPoolIdx.size(), 0 );
-}
-
-/**
- * Figures out indexing of the array of transferred pool n's used to fill
- * out proxies on each timestep.
- */
-void ZombiePoolInterface::assignXferIndex( unsigned int numProxyMols,
- unsigned int xferCompt,
- const vector< vector< unsigned int > >& voxy )
-{
- unsigned int idx = 0;
- for ( unsigned int i = 0; i < voxy.size(); ++i ) {
- const vector< unsigned int >& rpv = voxy[i];
- if ( rpv.size() > 0) { // There would be a transfer here
- for ( vector< unsigned int >::const_iterator
- j = rpv.begin(); j != rpv.end(); ++j ) {
- pools(*j)->addProxyTransferIndex( xferCompt, idx );
- }
- idx += numProxyMols;
- }
- }
-}
-
-/**
- * This function sets up the information about the pool transfer for
- * cross-compartment reactions. It consolidates the transfer into a
- * distinct vector for each direction of the transfer between each coupled
- * pair of ZombiePoolInterfaces.
- * This one call sets up the information about transfer on both sides
- * of the junction(s) between current ZombiePoolInterface and otherZombiePoolInterface.
- *
- * VoxelJunction refers to a specific junction between compartments.
- * It has the following relevant fields:
- * first, second: VoxelIndex for the first and second compartments.
- * firstVol, secondVol: VoxelVolume for the first and second compartments.
- */
-void ZombiePoolInterface::setupXfer( Id myZombiePoolInterface, Id otherZombiePoolInterface,
- unsigned int numProxyMols, const vector< VoxelJunction >& vj )
-{
- const ChemCompt *myCompt = reinterpret_cast< const ChemCompt* >(
- compartment_.eref().data() );
- ZombiePoolInterface* otherZombiePoolInterfacePtr =
- reinterpret_cast< ZombiePoolInterface* >(
- otherZombiePoolInterface.eref().data() );
- const ChemCompt *otherCompt = reinterpret_cast< const ChemCompt* >(
- otherZombiePoolInterfacePtr->compartment_.eref().data() );
- // Use this so we can figure out what the other side will send.
- vector< vector< unsigned int > > proxyVoxy( myCompt->getNumEntries() );
- vector< vector< unsigned int > > reverseProxyVoxy( otherCompt->getNumEntries() );
- assert( xfer_.size() > 0 && otherZombiePoolInterfacePtr->xfer_.size() > 0 );
- unsigned int myZombiePoolInterfaceIndex = xfer_.size() -1;
- unsigned int otherZombiePoolInterfaceIndex = otherZombiePoolInterfacePtr->xfer_.size() -1;
- for ( unsigned int i = 0; i < vj.size(); ++i ) {
- unsigned int j = vj[i].first;
- assert( j < getNumLocalVoxels() ); // Check voxel indices.
- proxyVoxy[j].push_back( vj[i].second );
- pools(j)->addProxyVoxy( myZombiePoolInterfaceIndex,
- otherZombiePoolInterfacePtr->compartment_, vj[i].second);
- unsigned int k = vj[i].second;
- assert( k < otherCompt->getNumEntries() );
- reverseProxyVoxy[k].push_back( vj[i].first );
- otherZombiePoolInterfacePtr->pools(k)->addProxyVoxy(
- otherZombiePoolInterfaceIndex, compartment_, vj[i].first );
- }
-
- // Build the indexing for the data values to transfer on each timestep
- assignXferIndex( numProxyMols, myZombiePoolInterfaceIndex, reverseProxyVoxy );
- otherZombiePoolInterfacePtr->assignXferIndex(
- numProxyMols, otherZombiePoolInterfaceIndex, proxyVoxy );
- // Figure out which voxels participate in data transfer.
- assignXferVoxels( myZombiePoolInterfaceIndex );
- otherZombiePoolInterfacePtr->assignXferVoxels( otherZombiePoolInterfaceIndex );
-}
-
-
-/**
- * Builds up the list of proxy pools on either side of the junction,
- * and assigns to the XferInfo data structures for use during runtime.
- */
-unsigned int ZombiePoolInterface::assignProxyPools( const map< Id, vector< Id > >& xr,
- Id myZombiePoolInterface, Id otherZombiePoolInterface, Id otherComptId )
-{
- map< Id, vector< Id > >::const_iterator i = xr.find( otherComptId );
- vector< Id > proxyMols;
- if ( i != xr.end() )
- proxyMols = i->second;
- ZombiePoolInterface* otherZombiePoolInterfacePtr =
- reinterpret_cast< ZombiePoolInterface* >(
- otherZombiePoolInterface.eref().data() );
-
- vector< Id > otherProxies = LookupField< Id, vector< Id > >::get(
- otherZombiePoolInterfacePtr->stoich_, "proxyPools", stoich_ );
-
- proxyMols.insert( proxyMols.end(),
- otherProxies.begin(), otherProxies.end() );
- // if ( proxyMols.size() == 0 )
- // return 0;
- sort( proxyMols.begin(), proxyMols.end() );
- xfer_.push_back( XferInfo( otherZombiePoolInterface ) );
-
- otherZombiePoolInterfacePtr->xfer_.push_back( XferInfo( myZombiePoolInterface ) );
- vector< unsigned int >& xfi = xfer_.back().xferPoolIdx;
- vector< unsigned int >& oxfi = otherZombiePoolInterfacePtr->xfer_.back().xferPoolIdx;
- xfi.resize( proxyMols.size() );
- oxfi.resize( proxyMols.size() );
- for ( unsigned int i = 0; i < xfi.size(); ++i ) {
- xfi[i] = getPoolIndex( proxyMols[i].eref() );
- oxfi[i] = otherZombiePoolInterfacePtr->getPoolIndex(
- proxyMols[i].eref() );
- }
- return proxyMols.size();
-}
-
-
-// This function cleans out the RateTerms of cross reactions that
-// don't have anything to connect to.
-// It should be called after all cross reacs have been assigned.
-void ZombiePoolInterface::filterCrossRateTerms( const vector< Id >& xreacs,
- const vector< pair< Id, Id > >& xrt )
-{
- for (unsigned int i = 0; i < getNumLocalVoxels(); ++i ) {
- pools(i)->filterCrossRateTerms( xreacs, xrt );
- }
-}
-
-/**
- * This function builds cross-solver reaction calculations. For the
- * specified pair of stoichs (this->stoich_, otherStoich) it identifies
- * interacting molecules, finds where the junctions are, sets up the
- * info to build the data transfer vector, and sets up the transfer
- * itself.
- */
-void ZombiePoolInterface::setupCrossSolverReacs( const map< Id, vector< Id > >& xr,
- Id otherStoich )
-{
- const ChemCompt *myCompt = reinterpret_cast< const ChemCompt* >(
- compartment_.eref().data() );
- Id otherComptId = Field< Id >::get( otherStoich, "compartment" );
- Id myZombiePoolInterface = Field< Id >::get( stoich_, "ksolve" );
- if ( myZombiePoolInterface == Id() )
- return;
- Id otherZombiePoolInterface = Field< Id >::get( otherStoich, "ksolve" );
- if ( otherZombiePoolInterface == Id() )
- return;
-
- // Establish which molecules will be exchanged.
- unsigned int numPools = assignProxyPools( xr, myZombiePoolInterface, otherZombiePoolInterface,
- otherComptId );
- if ( numPools == 0 ) return;
-
- // Then, figure out which voxels do the exchange.
- // Note that vj has a list of pairs of voxels on either side of a
- // junction. If one voxel on self touches 5 voxels on other, then
- // there will be five entries in vj for this contact.
- // If one voxel on self touches two different compartments, then
- // a distinct vj vector must be built for those contacts.
- const ChemCompt *otherCompt = reinterpret_cast< const ChemCompt* >(
- otherComptId.eref().data() );
- vector< VoxelJunction > vj;
- myCompt->matchMeshEntries( otherCompt, vj );
- if ( vj.size() == 0 )
- return;
-
- // This function sets up the information about the pool transfer on
- // both sides.
- setupXfer( myZombiePoolInterface, otherZombiePoolInterface, numPools, vj );
-
- /// This sets up the volume scaling from cross reac terms
- // Deprecated. Handled by setupCrossSolverReacVols.
- // buildCrossReacVolScaling( otherZombiePoolInterface, vj );
-
- // Here we set up the messaging.
- Shell *shell = reinterpret_cast< Shell* >( Id().eref().data() );
- shell->doAddMsg( "Single", myZombiePoolInterface, "xCompt", otherZombiePoolInterface, "xCompt" );
-}
-
-/**
- * This fills the vols vector with the volume of the abutting
- * voxel on compt. If there are no abutting voxels on a given
- * voxel then that entry of the vols vector is filled with a zero.
- * There is exactly one vols entry for each voxel of the local compt.
- */
-void ZombiePoolInterface::matchJunctionVols( vector< double >& vols, Id otherComptId )
- const
-{
- vols.resize( getNumLocalVoxels() );
- for ( unsigned int i = 0; i < vols.size(); ++i )
- vols[i] = volume(i);
- if ( otherComptId == compartment_ ) {
- // This may legitimately happen if the substrate or product is
- // on the local compartment.
- // cout << "Warning: ZombiePoolInterface::matchJunctionVols: self compt.\n";
- return;
- }
- const ChemCompt *myCompt = reinterpret_cast< const ChemCompt* >(
- compartment_.eref().data() );
- const ChemCompt *otherCompt = reinterpret_cast< const ChemCompt* >(
- otherComptId.eref().data() );
- vector< VoxelJunction > vj;
- myCompt->matchMeshEntries( otherCompt, vj );
- if ( vj.size() == 0 )
- return;
- for ( vector< VoxelJunction >::const_iterator
- i = vj.begin(); i != vj.end(); ++i ) {
- assert( i->first < vols.size() );
- /*
- if ( !doubleEq( vols[ i->first ], 0.0 ) )
- cout << "Warning: ZombiePoolInterface::matchJuntionVols: repeated voxel\n";
- */
- vols[ i->first ] = i->secondVol;
- }
-}
-
-/**
- * This function builds cross-solver reaction volume scaling.
- */
-void ZombiePoolInterface::setupCrossSolverReacVols(
- const vector< vector< Id > >& subCompts,
- const vector< vector< Id > >& prdCompts )
-{
- map< Id, vector< double > > comptVolMap;
- const Stoich* stoichPtr = reinterpret_cast< const Stoich* >(
- stoich_.eref().data() );
- unsigned int numOffSolverReacs =
- stoichPtr->getNumRates() - stoichPtr->getNumCoreRates();
- assert( subCompts.size() == numOffSolverReacs );
- assert( prdCompts.size() == numOffSolverReacs );
- for ( unsigned int i = 0 ; i < getNumLocalVoxels(); ++i )
- pools(i)->resetXreacScale( numOffSolverReacs );
- for( unsigned int i = 0; i < numOffSolverReacs; ++i ) {
- for ( unsigned int j = 0; j < subCompts[i].size(); ++j ) {
- map< Id, vector< double > >::iterator cvi;
- vector< double > vols;
- cvi = comptVolMap.find( subCompts[i][j] );
- if ( cvi == comptVolMap.end() ) {
- matchJunctionVols( vols, subCompts[i][j] );
- comptVolMap[subCompts[i][j]] = vols;
- } else {
- vols = cvi->second;
- }
- assert( vols.size() == getNumLocalVoxels() );
- for ( unsigned int k = 0; k < vols.size(); ++k )
- pools(k)->forwardReacVolumeFactor( i, vols[k] );
- }
-
- for ( unsigned int j = 0; j < prdCompts[i].size(); ++j ) {
- map< Id, vector< double > >::iterator cvi;
- vector< double > vols;
- cvi = comptVolMap.find( prdCompts[i][j] );
- if ( cvi == comptVolMap.end() ) {
- matchJunctionVols( vols, prdCompts[i][j] );
- comptVolMap[prdCompts[i][j]] = vols;
- } else {
- vols = cvi->second;
- }
- assert( vols.size() == getNumLocalVoxels() );
- for ( unsigned int k = 0; k < vols.size(); ++k )
- pools(k)->backwardReacVolumeFactor( i, vols[k] );
- }
- }
-}
Id ZombiePoolInterface::getCompartment() const
{
diff --git a/ksolve/ZombiePoolInterface.h b/ksolve/ZombiePoolInterface.h
index 10104fd2..fe5bb1d8 100644
--- a/ksolve/ZombiePoolInterface.h
+++ b/ksolve/ZombiePoolInterface.h
@@ -96,16 +96,10 @@ class ZombiePoolInterface
virtual void setCompartment( Id compartment );
Id getCompartment() const;
- /// Sets up cross-solver reactions.
- void setupCrossSolverReacs(
- const map< Id, vector< Id > >& xr,
- Id otherStoich );
- void setupCrossSolverReacVols(
- const vector< vector< Id > >& subCompts,
- const vector< vector< Id > >& prdCompts );
-
- void filterCrossRateTerms( const vector< Id >& xreacs,
- const vector< pair< Id, Id > >& xrt );
+ /// Used for telling Dsolver to handle all ops across Junctions
+ virtual void updateJunctions( double dt );
+ /// Used to tell Dsolver to assign 'prev' values.
+ virtual void setPrev();
/**
* Informs the solver that the rate terms or volumes have changed
* and that the parameters must be updated.
@@ -116,24 +110,6 @@ class ZombiePoolInterface
/// Return pool index, using Stoich ptr to do lookup.
virtual unsigned int getPoolIndex( const Eref& er ) const = 0;
- //////////////////////////////////////////////////////////////
- // Utility functions for Cross-compt reaction setup.
- //////////////////////////////////////////////////////////////
- void xComptIn( const Eref& e, Id srcZombiePoolInterface,
- vector< double > values );
- // void xComptOut( const Eref& e );
- void assignXferVoxels( unsigned int xferCompt );
- void assignXferIndex( unsigned int numProxyMols,
- unsigned int xferCompt,
- const vector< vector< unsigned int > >& voxy );
- void setupXfer( Id myZombiePoolInterface,
- Id otherZombiePoolInterface,
- unsigned int numProxyMols, const vector< VoxelJunction >& vj );
- unsigned int assignProxyPools( const map< Id, vector< Id > >& xr,
- Id myZombiePoolInterface, Id otherZombiePoolInterface,
- Id otherComptId );
- void matchJunctionVols( vector< double >& vols, Id otherComptId )
- const;
//////////////////////////////////////////////////////////////
protected:
@@ -146,12 +122,6 @@ class ZombiePoolInterface
/// Id of Chem compartment used to figure out volumes of voxels.
Id compartment_;
- /**
- * All the data transfer information from current to other solvers.
- * xfer_[otherKsolveIndex]
- */
- vector< XferInfo > xfer_;
-
/// Flag: True when solver setup has been completed.
bool isBuilt_;
};
diff --git a/mesh/CMakeLists.txt b/mesh/CMakeLists.txt
index 6046fc56..a04e31dd 100644
--- a/mesh/CMakeLists.txt
+++ b/mesh/CMakeLists.txt
@@ -13,5 +13,6 @@ add_library(mesh
SpineEntry.cpp
SpineMesh.cpp
PsdMesh.cpp
+ EndoMesh.cpp
testMesh.cpp
)
diff --git a/mesh/ChemCompt.cpp b/mesh/ChemCompt.cpp
index 32ae58e2..1d2b3f7f 100644
--- a/mesh/ChemCompt.cpp
+++ b/mesh/ChemCompt.cpp
@@ -89,6 +89,18 @@ const Cinfo* ChemCompt::initCinfo()
&ChemCompt::getStencilIndex
);
+ static ValueFinfo< ChemCompt, bool > isMembraneBound(
+ "isMembraneBound",
+ "Flag, set to True for meshes where each voxel is membrane "
+ "bound. \n"
+ "NeuroMesh and SpineMesh are false. \n"
+ "CubeMesh, CylMesh, and EndoMesh can be either. If they are "
+ "membrane bound they can still interact via channels and "
+ "cross-compartment reactions. ",
+ &ChemCompt::setIsMembraneBound,
+ &ChemCompt::getIsMembraneBound
+ );
+
//////////////////////////////////////////////////////////////
// MsgDest Definitions
//////////////////////////////////////////////////////////////
@@ -148,6 +160,7 @@ const Cinfo* ChemCompt::initCinfo()
&numDimensions, // ReadOnlyValue
&stencilRate, // ReadOnlyLookupValue
&stencilIndex, // ReadOnlyLookupValue
+ &isMembraneBound, // Value
voxelVolOut(), // SrcFinfo
&buildDefaultMesh, // DestFinfo
&setVolumeNotRates, // DestFinfo
@@ -184,7 +197,8 @@ static const Cinfo* chemMeshCinfo = ChemCompt::initCinfo();
ChemCompt::ChemCompt()
:
- entry_( this )
+ entry_( this ),
+ isMembraneBound_( false )
{
;
}
@@ -338,6 +352,16 @@ vector< unsigned int > ChemCompt::getStencilIndex( unsigned int row ) const
return this->getNeighbors( row );
}
+bool ChemCompt::getIsMembraneBound() const
+{
+ return isMembraneBound_;
+}
+
+void ChemCompt::setIsMembraneBound( bool v )
+{
+ isMembraneBound_ = v;
+}
+
//////////////////////////////////////////////////////////////
// Dest Definitions
diff --git a/mesh/ChemCompt.h b/mesh/ChemCompt.h
index eb77b965..ba909d2e 100644
--- a/mesh/ChemCompt.h
+++ b/mesh/ChemCompt.h
@@ -69,6 +69,9 @@ class ChemCompt
void setMethod( string method );
string getMethod() const;
+ bool getIsMembraneBound() const;
+ void setIsMembraneBound( bool v );
+
/**
* Function to return the stencil values used in the
* diffusion calculations for voxelized compartments.
@@ -349,6 +352,15 @@ class ChemCompt
* function.
*/
string method_;
+
+ /**
+ * Flag, set to True for meshes where each voxel is membrane bound.
+ * NeuroMesh and SpineMesh are false
+ * CubeMesh, CylMesh, and EndoMesh can be either. If they are
+ * membrane bound they can still interact via channels and
+ * cross-compartment reactions.
+ */
+ bool isMembraneBound_;
};
extern SrcFinfo5<
diff --git a/mesh/CubeMesh.cpp b/mesh/CubeMesh.cpp
index da90d72b..728c69ba 100644
--- a/mesh/CubeMesh.cpp
+++ b/mesh/CubeMesh.cpp
@@ -17,6 +17,7 @@
#include "ChemCompt.h"
#include "MeshCompt.h"
#include "CubeMesh.h"
+#include "EndoMesh.h"
const unsigned int CubeMesh::EMPTY = ~0;
const unsigned int CubeMesh::SURFACE = ~1;
@@ -223,13 +224,26 @@ const Cinfo* CubeMesh::initCinfo()
&surface, // Value
};
+ static string doc[] =
+ {
+ "Name", "CubeMesh",
+ "Author", "Upi Bhalla",
+ "Description", "Chemical compartment with cuboid grid. "
+ "Defaults to a cube of size 10 microns, with mesh size "
+ "also 10 microns, so that there is just 1 cubic voxel. "
+ "These defaults are similar to that of a typical cell. "
+ "Can be configured to have different x,y,z dimensions and "
+ "also different dx,dy,dz voxel sizes. ",
+ };
static Dinfo< CubeMesh > dinfo;
static Cinfo cubeMeshCinfo (
"CubeMesh",
ChemCompt::initCinfo(),
cubeMeshFinfos,
sizeof( cubeMeshFinfos ) / sizeof ( Finfo* ),
- &dinfo
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string)
);
return &cubeMeshCinfo;
@@ -252,12 +266,12 @@ CubeMesh::CubeMesh()
x0_( 0.0 ),
y0_( 0.0 ),
z0_( 0.0 ),
- x1_( 1.0 ),
- y1_( 1.0 ),
- z1_( 1.0 ),
- dx_( 1.0 ),
- dy_( 1.0 ),
- dz_( 1.0 ),
+ x1_( 10e-6 ),
+ y1_( 10e-6 ),
+ z1_( 10e-6 ),
+ dx_( 10e-6 ),
+ dy_( 10e-6 ),
+ dz_( 10e-6 ),
nx_( 1 ),
ny_( 1 ),
nz_( 1 ),
@@ -887,9 +901,16 @@ void CubeMesh::innerSetNumEntries( unsigned int n )
cout << "Warning: CubeMesh::innerSetNumEntries is readonly.\n";
}
+// We assume this is linear diffusion for now. Fails for 2 or 3-D diffusion
vector< unsigned int > CubeMesh::getParentVoxel() const
{
- static vector< unsigned int > ret;
+ unsigned int numEntries = innerGetNumEntries();
+ vector< unsigned int > ret( numEntries );
+ if ( numEntries > 0 )
+ ret[0] = static_cast< unsigned int >( -1 );
+ for (unsigned int i = 1; i < numEntries; ++i )
+ ret[i] = i-1;
+
return ret;
}
@@ -921,14 +942,27 @@ const vector< double >& CubeMesh::vGetVoxelMidpoint() const
const vector< double >& CubeMesh::getVoxelArea() const
{
static vector< double > area;
- assert( 0 ); // Not yet operational
+ if ( nx_ * ny_ == 1 )
+ area.resize( nz_, dx_ * dy_ );
+ else if ( nx_ * nz_ == 1 )
+ area.resize( ny_, dx_ * dz_ );
+ else if ( ny_ * nz_ == 1 )
+ area.resize( nx_, dy_ * dz_ );
+ else
+ area.resize( nx_, dy_ * dz_ ); // Just put in a number.
+ assert( area.size() == nx_ * ny_ * nz_ );
return area;
}
const vector< double >& CubeMesh::getVoxelLength() const
{
static vector< double > length;
- assert( 0 ); // Not yet operational
+ if ( dx_ > dy_ && dx_ > dz_ )
+ length.assign( innerGetNumEntries(), dx_ );
+ else if ( dy_ > dz_ )
+ length.assign( innerGetNumEntries(), dy_ );
+ else
+ length.assign( innerGetNumEntries(), dz_ );
return length;
}
@@ -1100,9 +1134,15 @@ void CubeMesh::matchMeshEntries( const ChemCompt* other,
cout << "Warning:CubeMesh::matchMeshEntries: cannot yet handle unequal meshes\n";
}
*/
- } else {
- cout << "Warning:CubeMesh::matchMeshEntries: cannot yet handle Neuro or Cyl meshes.\n";
+ return;
+ }
+ const EndoMesh* em = dynamic_cast< const EndoMesh* >( other );
+ if ( em ) {
+ em->matchMeshEntries( this, ret );
+ flipRet( ret );
+ return;
}
+ cout << "Warning:CubeMesh::matchMeshEntries: cannot yet handle Neuro or Cyl meshes.\n";
}
unsigned int CubeMesh::numDims() const
diff --git a/mesh/CylMesh.cpp b/mesh/CylMesh.cpp
index ed2c6c43..82cd895f 100644
--- a/mesh/CylMesh.cpp
+++ b/mesh/CylMesh.cpp
@@ -22,6 +22,7 @@
// #include "NeuroStencil.h"
#include "NeuroMesh.h"
#include "CylMesh.h"
+#include "EndoMesh.h"
#include "../utility/numutil.h"
const Cinfo* CylMesh::initCinfo()
{
@@ -131,13 +132,26 @@ const Cinfo* CylMesh::initCinfo()
&totLength, // ReadOnlyValue
};
+ static string doc[] =
+ {
+ "Name", "CylMesh",
+ "Author", "Upi Bhalla",
+ "Description", "Chemical compartment with cylindrical geometry. "
+ "Defaults to a uniform cylinder of radius 1 micron, "
+ "length 100 microns, and voxel length 1 micron so there "
+ "are 100 voxels in the cylinder. "
+ "The cylinder can be given a linear taper, by assigning "
+ "different radii r0 and r1 to the two ends. ",
+ };
static Dinfo< CylMesh > dinfo;
static Cinfo cylMeshCinfo (
"CylMesh",
ChemCompt::initCinfo(),
cylMeshFinfos,
sizeof( cylMeshFinfos ) / sizeof ( Finfo* ),
- &dinfo
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string)
);
return &cylMeshCinfo;
@@ -154,20 +168,20 @@ static const Cinfo* cylMeshCinfo = CylMesh::initCinfo();
//////////////////////////////////////////////////////////////////
CylMesh::CylMesh()
:
- numEntries_( 1 ),
+ numEntries_( 100 ),
useCaps_( 0 ),
isToroid_( false ),
x0_( 0.0 ),
y0_( 0.0 ),
z0_( 0.0 ),
- x1_( 1.0 ),
+ x1_( 100.0e-6 ),
y1_( 0.0 ),
z1_( 0.0 ),
- r0_( 1.0 ),
- r1_( 1.0 ),
- diffLength_( 1.0 ),
+ r0_( 1.0e-6 ),
+ r1_( 1.0e-6 ),
+ diffLength_( 1.0e-6 ),
surfaceGranularity_( 0.1 ),
- totLen_( 1.0 ),
+ totLen_( 100.0e-6 ),
rSlope_( 0.0 ),
lenSlope_( 0.0 )
{
@@ -818,20 +832,27 @@ void CylMesh::matchMeshEntries( const ChemCompt* other,
const CylMesh* cyl = dynamic_cast< const CylMesh* >( other );
if ( cyl ) {
matchCylMeshEntries( cyl, ret );
- } else {
- const CubeMesh* cube = dynamic_cast< const CubeMesh* >( other );
- if ( cube ) {
- matchCubeMeshEntries( cube, ret );
- } else {
- const NeuroMesh* nm = dynamic_cast< const NeuroMesh* >( other );
- if ( nm ) {
- matchNeuroMeshEntries( nm, ret );
- } else {
- cout << "Warning:CylMesh::matchMeshEntries: " <<
- " unknown mesh type\n";
- }
- }
+ return;
+ }
+ const EndoMesh* em = dynamic_cast< const EndoMesh* >( other );
+ if ( em )
+ {
+ em->matchMeshEntries( this, ret );
+ flipRet( ret );
+ return;
+ }
+ const CubeMesh* cube = dynamic_cast< const CubeMesh* >( other );
+ if ( cube ) {
+ matchCubeMeshEntries( cube, ret );
+ return;
+ }
+ const NeuroMesh* nm = dynamic_cast< const NeuroMesh* >( other );
+ if ( nm ) {
+ matchNeuroMeshEntries( nm, ret );
+ return;
}
+ cout << "Warning:CylMesh::matchMeshEntries: " <<
+ " unknown mesh type\n";
}
// Look for end-to-end diffusion, not sideways for now.
diff --git a/mesh/EndoMesh.cpp b/mesh/EndoMesh.cpp
new file mode 100644
index 00000000..33c13bdb
--- /dev/null
+++ b/mesh/EndoMesh.cpp
@@ -0,0 +1,464 @@
+/**********************************************************************
+** This program is part of 'MOOSE', the
+** Messaging Object Oriented Simulation Environment.
+** Copyright (C) 2011 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 "header.h"
+#include "SparseMatrix.h"
+#include "ElementValueFinfo.h"
+#include "../utility/Vec.h"
+#include "Boundary.h"
+#include "MeshEntry.h"
+// #include "Stencil.h"
+#include "ChemCompt.h"
+#include "MeshCompt.h"
+#include "CubeMesh.h"
+#include "CylBase.h"
+#include "NeuroNode.h"
+// #include "NeuroStencil.h"
+#include "NeuroMesh.h"
+#include "EndoMesh.h"
+#include "../utility/numutil.h"
+const Cinfo* EndoMesh::initCinfo()
+{
+ //////////////////////////////////////////////////////////////
+ // Field Definitions
+ //////////////////////////////////////////////////////////////
+ static ElementValueFinfo< EndoMesh, double > rPower(
+ "rPower",
+ "Used in rEndo = rScale * pow(surroundVol, rPower)."
+ "Used to obtain radius of each endo voxel from matching "
+ "surround voxel. Defaults to 1/3",
+ &EndoMesh::setRpower,
+ &EndoMesh::getRpower
+ );
+ static ElementValueFinfo< EndoMesh, double > rScale(
+ "rScale",
+ "Used in rEndo = rScale * pow(surroundVol, rPower)."
+ "Used to obtain radius of each endo voxel from matching "
+ "surround voxel. Defaults to 0.5",
+ &EndoMesh::setRscale,
+ &EndoMesh::getRscale
+ );
+ static ElementValueFinfo< EndoMesh, double > aPower(
+ "aPower",
+ "Used in rEndo = aScale * pow(surroundVol, aPower)."
+ "Used to obtain area of each endo voxel from matching "
+ "surround voxel. Defaults to 1/2",
+ &EndoMesh::setApower,
+ &EndoMesh::getApower
+ );
+ static ElementValueFinfo< EndoMesh, double > aScale(
+ "aScale",
+ "Used in rEndo = aScale * pow(surroundVol, aPower)."
+ "Used to obtain area of each endo voxel from matching "
+ "surround voxel. Defaults to 0.25",
+ &EndoMesh::setAscale,
+ &EndoMesh::getAscale
+ );
+ static ElementValueFinfo< EndoMesh, double > vPower(
+ "vPower",
+ "Used in rEndo = vScale * pow(surroundVol, vPower)."
+ "Used to obtain volume of each endo voxel from matching "
+ "surround voxel. Defaults to 1.",
+ &EndoMesh::setVpower,
+ &EndoMesh::getVpower
+ );
+ static ElementValueFinfo< EndoMesh, double > vScale(
+ "vScale",
+ "Used in rEndo = vScale * pow(surroundVol, vPower)."
+ "Used to obtain volume of each endo voxel from matching "
+ "surround voxel. Defaults to 0.125",
+ &EndoMesh::setVscale,
+ &EndoMesh::getVscale
+ );
+
+ static ElementValueFinfo< EndoMesh, ObjId > surround(
+ "surround",
+ "ObjId of compartment surrounding current EndoMesh",
+ &EndoMesh::setSurround,
+ &EndoMesh::getSurround
+ );
+
+ static ElementValueFinfo< EndoMesh, bool > doAxialDiffusion(
+ "doAxialDiffusion",
+ "Flag to determine if endoMesh should undertake axial "
+ "diffusion. Defaults to 'false'. "
+ "Should only be used within NeuroMesh and CylMesh. "
+ "Must be assigned before Dsolver is built.",
+ &EndoMesh::setDoAxialDiffusion,
+ &EndoMesh::getDoAxialDiffusion
+ );
+
+
+ //////////////////////////////////////////////////////////////
+ // MsgDest Definitions
+ //////////////////////////////////////////////////////////////
+
+ //////////////////////////////////////////////////////////////
+ // Field Elements
+ //////////////////////////////////////////////////////////////
+
+ static Finfo* endoMeshFinfos[] = {
+ &rPower, // Value
+ &rScale, // Value
+ &aPower, // Value
+ &aScale, // Value
+ &vPower, // Value
+ &vScale, // Value
+ &surround, // Value
+ &doAxialDiffusion, // Value
+ };
+
+ static Dinfo< EndoMesh > dinfo;
+ static Cinfo endoMeshCinfo (
+ "EndoMesh",
+ ChemCompt::initCinfo(),
+ endoMeshFinfos,
+ sizeof( endoMeshFinfos ) / sizeof ( Finfo* ),
+ &dinfo
+ );
+
+ return &endoMeshCinfo;
+}
+
+//////////////////////////////////////////////////////////////
+// Basic class Definitions
+//////////////////////////////////////////////////////////////
+
+static const Cinfo* endoMeshCinfo = EndoMesh::initCinfo();
+
+//////////////////////////////////////////////////////////////////
+// Class stuff.
+//////////////////////////////////////////////////////////////////
+EndoMesh::EndoMesh()
+ :
+ parent_( 0 ),
+ rPower_( 1.0 / 3.0 ),
+ rScale_( 0.5 ),
+ aPower_( 0.5 ),
+ aScale_( 0.25 ),
+ vPower_( 1.0 ),
+ vScale_( 0.125 ),
+ doAxialDiffusion_( false )
+{
+ ;
+}
+
+EndoMesh::~EndoMesh()
+{
+ ;
+}
+
+//////////////////////////////////////////////////////////////////
+// Field assignment stuff
+//////////////////////////////////////////////////////////////////
+//
+void EndoMesh::setRpower( const Eref& e, double v )
+{
+ rPower_ = v;
+}
+
+double EndoMesh::getRpower( const Eref& e ) const
+{
+ return rPower_;
+}
+
+void EndoMesh::setRscale( const Eref& e, double v )
+{
+ rScale_ = v;
+}
+
+double EndoMesh::getRscale( const Eref& e ) const
+{
+ return rScale_;
+}
+
+void EndoMesh::setApower( const Eref& e, double v )
+{
+ aPower_ = v;
+}
+
+double EndoMesh::getApower( const Eref& e ) const
+{
+ return aPower_;
+}
+
+void EndoMesh::setAscale( const Eref& e, double v )
+{
+ aScale_ = v;
+}
+
+double EndoMesh::getAscale( const Eref& e ) const
+{
+ return aScale_;
+}
+
+void EndoMesh::setVpower( const Eref& e, double v )
+{
+ vPower_ = v;
+}
+
+double EndoMesh::getVpower( const Eref& e ) const
+{
+ return vPower_;
+}
+
+void EndoMesh::setVscale( const Eref& e, double v )
+{
+ vScale_ = v;
+}
+
+double EndoMesh::getVscale( const Eref& e ) const
+{
+ return vScale_;
+}
+
+void EndoMesh::setSurround( const Eref& e, ObjId v )
+{
+ if ( !v.element()->cinfo()->isA( "ChemCompt" ) ) {
+ cout << "Warning: 'surround' may only be set to an object of class 'ChemCompt'\n";
+ cout << v.path() << " is of class " << v.element()->cinfo()->name() << endl;
+ return;
+ }
+ surround_ = v;
+ parent_ = reinterpret_cast< const MeshCompt* >( v.data() );
+}
+
+ObjId EndoMesh::getSurround( const Eref& e ) const
+{
+ return surround_;
+}
+
+void EndoMesh::setDoAxialDiffusion( const Eref& e, bool v )
+{
+ doAxialDiffusion_ = v;
+}
+
+bool EndoMesh::getDoAxialDiffusion( const Eref& e ) const
+{
+ return doAxialDiffusion_;
+}
+
+//////////////////////////////////////////////////////////////////
+// FieldElement assignment stuff for MeshEntries
+//////////////////////////////////////////////////////////////////
+
+/// Virtual function to return MeshType of specified entry.
+unsigned int EndoMesh::getMeshType( unsigned int fid ) const
+{
+ return ENDO;
+}
+
+/// Virtual function to return dimensions of specified entry.
+unsigned int EndoMesh::getMeshDimensions( unsigned int fid ) const
+{
+ return 3;
+}
+
+/// Virtual function to return # of spatial dimensions of mesh
+unsigned int EndoMesh::innerGetDimensions() const
+{
+ return 1;
+}
+/// Virtual function to return volume of mesh Entry.
+double EndoMesh::getMeshEntryVolume( unsigned int fid ) const
+{
+ double vpa = parent_->getMeshEntryVolume( fid );
+ return vScale_ * pow( vpa, vPower_ );
+}
+
+/// Virtual function to return coords of mesh Entry.
+/// For Endo mesh, coords are just middle of parent.
+vector< double > EndoMesh::getCoordinates( unsigned int fid ) const
+{
+ vector< double > temp = parent_->getCoordinates( fid );
+ vector< double > ret;
+ if ( temp.size() > 6 ) {
+ ret.resize( 4 );
+ ret[0] = 0.5 * (temp[0] + temp[3] );
+ ret[1] = 0.5 * (temp[1] + temp[4] );
+ ret[2] = 0.5 * (temp[2] + temp[5] );
+ ret[3] = 0;
+ }
+ return ret;
+}
+
+/// Virtual function to return diffusion X-section area for each neighbor
+/// This applies if we have endo mesh voxels diffusing with each other
+vector< double > EndoMesh::getDiffusionArea( unsigned int fid ) const
+{
+ return vector< double >( parent_->getNumEntries(), 0.0 );
+}
+
+/// Virtual function to return scale factor for diffusion. 1 here.
+/// This applies if we have endo mesh voxels diffusing with each other
+vector< double > EndoMesh::getDiffusionScaling( unsigned int fid ) const
+{
+ return vector< double >( parent_->getNumEntries(), 0.0 );
+}
+
+/// Virtual function to return volume of mesh Entry, including
+/// for diffusively coupled voxels from other solvers.
+double EndoMesh::extendedMeshEntryVolume( unsigned int fid ) const
+{
+ double vpa = parent_->extendedMeshEntryVolume( fid );
+ return vScale_ * pow( vpa, vPower_ );
+}
+
+//////////////////////////////////////////////////////////////////
+// Dest funcs
+//////////////////////////////////////////////////////////////////
+
+/// More inherited virtual funcs: request comes in for mesh stats
+void EndoMesh::innerHandleRequestMeshStats( const Eref& e,
+ const SrcFinfo2< unsigned int, vector< double > >* meshStatsFinfo
+ )
+{
+ vector< double > ret( vGetEntireVolume() / parent_->getNumEntries(), 1);
+ meshStatsFinfo->send( e, 1, ret );
+}
+
+void EndoMesh::innerHandleNodeInfo(
+ const Eref& e,
+ unsigned int numNodes, unsigned int numThreads )
+{
+}
+//////////////////////////////////////////////////////////////////
+
+/**
+ * Inherited virtual func. Returns number of MeshEntry in array
+ */
+unsigned int EndoMesh::innerGetNumEntries() const
+{
+ return parent_->innerGetNumEntries();
+}
+
+/**
+ * Inherited virtual func. Dummy, because the EndoMesh just does what its
+ * parent does.
+ */
+void EndoMesh::innerSetNumEntries( unsigned int n )
+{
+}
+
+
+void EndoMesh::innerBuildDefaultMesh( const Eref& e,
+ double volume, unsigned int numEntries )
+{
+}
+
+/**
+ * This means the sibling voxel immediately adjacent to it, not the
+ * voxel surrounding it. Should we wish to do diffusion we would need
+ * to use the parent voxel of the surround voxel. Otherewise use
+ * EMTPY_VOXEL == -1U.
+ */
+vector< unsigned int > EndoMesh::getParentVoxel() const
+{
+ if ( doAxialDiffusion_ )
+ return parent_->getParentVoxel();
+
+ vector< unsigned int > ret( parent_->innerGetNumEntries(), -1U );
+ return ret;
+}
+
+const vector< double >& EndoMesh::vGetVoxelVolume() const
+{
+ static vector< double > ret;
+ ret = parent_->vGetVoxelVolume();
+ for ( auto i = ret.begin(); i != ret.end(); ++i )
+ *i = vScale_ * pow( *i, vPower_ );
+ return ret;
+}
+
+const vector< double >& EndoMesh::vGetVoxelMidpoint() const
+{
+ return parent_->vGetVoxelMidpoint();
+}
+
+const vector< double >& EndoMesh::getVoxelArea() const
+{
+ static vector< double > ret;
+ ret = parent_->vGetVoxelVolume();
+ for ( auto i = ret.begin(); i != ret.end(); ++i )
+ *i = aScale_ * pow( *i, aPower_ );
+ return ret;
+}
+
+const vector< double >& EndoMesh::getVoxelLength() const
+{
+ static vector< double > ret;
+ ret = parent_->vGetVoxelVolume();
+ for ( auto i = ret.begin(); i != ret.end(); ++i )
+ *i = vScale_ * pow( *i, vPower_ );
+ return ret;
+}
+
+double EndoMesh::vGetEntireVolume() const
+{
+ double vol = 0.0;
+ auto vec = vGetVoxelVolume();
+ for ( auto i = vec.begin(); i != vec.end(); ++i )
+ vol += *i;
+ return vol;
+}
+
+bool EndoMesh::vSetVolumeNotRates( double volume )
+{
+ return true; // maybe should be false? Unsure
+}
+
+//////////////////////////////////////////////////////////////////
+// Utility function for junctions
+//////////////////////////////////////////////////////////////////
+
+void EndoMesh::matchMeshEntries( const ChemCompt* other,
+ vector< VoxelJunction >& ret ) const
+{
+ /**
+ * Consider concentric EndoMeshes in an outer cyl/sphere of radius r_0.
+ * Radii, from out to in, are r_0, r_1, r_2.... r_innermost.
+ * Consider that we are on EndoMesh n.
+ * The diffusion length in each case is 1/2(r_n-1 - r_n+1).
+ * For the innermost we use r_n+1 = 0
+ * The diffusion XA is the area of the EndoMesh.
+ */
+ vector< double > vol = other->vGetVoxelVolume();
+ const MeshCompt* mc = static_cast< const MeshCompt* >( other );
+ vector< double > len = mc->getVoxelLength();
+ assert( vol.size() == len.size() );
+ ret.resize( vol.size() );
+ vector< double > myVol = vGetVoxelVolume();
+ assert( myVol.size() == vol.size() );
+ vector< double > myArea = getVoxelArea();
+ assert( myArea.size() == vol.size() );
+
+ for ( unsigned int i = 0; i < vol.size(); ++i ) {
+ double rSurround = sqrt( vol[i] / (len[i] * PI ) );
+ ret[i].first = ret[i].second = i;
+ ret[i].firstVol = myVol[i];
+ ret[i].secondVol = vol[i];
+
+ // For now we don't consider internal EndoMeshes.
+ ret[i].diffScale = 2.0 * myArea[i] / rSurround;
+ }
+}
+
+// This function returns the index.
+double EndoMesh::nearest( double x, double y, double z,
+ unsigned int& index ) const
+{
+ return parent_->nearest( x, y, z, index );
+}
+
+// This function returns the index.
+void EndoMesh::indexToSpace( unsigned int index,
+ double &x, double &y, double &z ) const
+{
+ parent_->indexToSpace( index, x, y, z );
+}
diff --git a/mesh/EndoMesh.h b/mesh/EndoMesh.h
new file mode 100644
index 00000000..5567bda7
--- /dev/null
+++ b/mesh/EndoMesh.h
@@ -0,0 +1,150 @@
+/**********************************************************************
+** This program is part of 'MOOSE', the
+** Messaging Object Oriented Simulation Environment.
+** Copyright (C) 2011 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.
+**********************************************************************/
+
+#ifndef _ENDO_MESH_H
+#define _ENDO_MESH_H
+
+/**
+ * The EndoMesh is a cellular compartment entirely contained within another
+ * compartment. It supports diffusion only with its immediate surround,
+ * of which there can only be one, and with its immediate interior, which
+ * can be one or more
+ * EndoMeshes. Each voxel in the EndoMesh is well mixed, and does not
+ * exchange with any other voxels.
+ * Typically used in modelling endosomes, ER, mitochondria, and other
+ * organelles.
+ */
+class EndoMesh: public MeshCompt
+{
+ public:
+ EndoMesh();
+ ~EndoMesh();
+ //////////////////////////////////////////////////////////////////
+ // Field assignment stuff
+ //////////////////////////////////////////////////////////////////
+ void setRpower( const Eref& e, double v );
+ double getRpower( const Eref& e ) const;
+ void setRscale( const Eref& e, double v );
+ double getRscale( const Eref& e ) const;
+ void setApower( const Eref& e, double v );
+ double getApower( const Eref& e ) const;
+ void setAscale( const Eref& e, double v );
+ double getAscale( const Eref& e ) const;
+ void setVpower( const Eref& e, double v );
+ double getVpower( const Eref& e ) const;
+ void setVscale( const Eref& e, double v );
+ double getVscale( const Eref& e ) const;
+ void setSurround( const Eref& e, ObjId v );
+ ObjId getSurround( const Eref& e ) const;
+ void setDoAxialDiffusion( const Eref& e, bool v );
+ bool getDoAxialDiffusion( const Eref& e ) const;
+
+ //////////////////////////////////////////////////////////////////
+ // FieldElement assignment stuff for MeshEntries
+ //////////////////////////////////////////////////////////////////
+
+ /// Virtual function to return MeshType of specified entry.
+ unsigned int getMeshType( unsigned int fid ) const;
+ /// Virtual function to return dimensions of specified entry.
+ unsigned int getMeshDimensions( unsigned int fid ) const;
+ unsigned int innerGetDimensions() const;
+ /// Virtual function to return volume of mesh Entry.
+ double getMeshEntryVolume( unsigned int fid ) const;
+ /// Virtual function to return coords of mesh Entry.
+ vector< double > getCoordinates( unsigned int fid ) const;
+ /// Virtual function to return diffusion X-section area
+ vector< double > getDiffusionArea( unsigned int fid ) const;
+ /// Virtual function to return scale factor for diffusion. 1 here.
+ vector< double > getDiffusionScaling( unsigned int fid ) const;
+ /// Volume of mesh Entry including abutting diff-coupled voxels
+ double extendedMeshEntryVolume( unsigned int fid ) const;
+
+ //////////////////////////////////////////////////////////////////
+
+ /**
+ * Inherited virtual func. Returns number of MeshEntry in array
+ */
+ unsigned int innerGetNumEntries() const;
+ /// Inherited virtual func.
+ void innerSetNumEntries( unsigned int n );
+
+ /// Inherited virtual, do nothing for now.
+ vector< unsigned int > getParentVoxel() const;
+ const vector< double >& vGetVoxelVolume() const;
+ const vector< double >& vGetVoxelMidpoint() const;
+ const vector< double >& getVoxelArea() const;
+ const vector< double >& getVoxelLength() const;
+
+ /// Inherited virtual. Returns entire volume of compartment.
+ double vGetEntireVolume() const;
+
+ /// Inherited virtual. Resizes len and dia of each voxel.
+ bool vSetVolumeNotRates( double volume );
+ //////////////////////////////////////////////////////////////////
+ // Dest funcs
+ //////////////////////////////////////////////////////////////////
+
+ /// Virtual func to make a mesh with specified Volume and numEntries
+ void innerBuildDefaultMesh( const Eref& e,
+ double volume, unsigned int numEntries );
+
+ void innerHandleRequestMeshStats(
+ const Eref& e,
+ const SrcFinfo2< unsigned int, vector< double > >*
+ meshStatsFinfo
+ );
+
+ void innerHandleNodeInfo(
+ const Eref& e,
+ unsigned int numNodes, unsigned int numThreads );
+
+ //////////////////////////////////////////////////////////////////
+ // inherited virtual funcs for Boundary
+ //////////////////////////////////////////////////////////////////
+
+ void matchMeshEntries( const ChemCompt* other,
+ vector< VoxelJunction > & ret ) const;
+
+ double nearest( double x, double y, double z,
+ unsigned int& index ) const;
+
+ void indexToSpace( unsigned int index,
+ double& x, double& y, double& z ) const;
+
+ //////////////////////////////////////////////////////////////////
+
+ static const Cinfo* initCinfo();
+
+ private:
+ /**
+ * Surrounding mesh. Use this for calculating
+ * all volume and diffusion terms, don't maintain any local
+ * variables.
+ */
+ ObjId surround_;
+ const MeshCompt* parent_;
+
+ Id insideMeshes_; /// EndoMeshes inside. Used to update.
+
+ /**
+ * The reason why the powers don't have to be 1/3, 1/2 and 1 is
+ * because some organelles are not a simple linear vol relationship
+ * with their surround. I want to keep flexibility.
+ */
+ double rPower_; /// rEndo = rScale * pow( surroundVol, rPower_)
+ double rScale_; /// rEndo = rScale * pow( surroundVol, rPower_)
+ double aPower_; /// aEndo = aScale * pow( surroundVol, aPower_)
+ double aScale_; /// aEndo = aScale * pow( surroundVol, aPower_)
+ double vPower_; /// vEndo = vScale * pow( surroundVol, vPower_)
+ double vScale_; /// vEndo = vScale * pow( surroundVol, vPower_)
+
+ bool doAxialDiffusion_; /// Flag for axial diffusion
+};
+
+#endif // _ENDO_MESH_H
diff --git a/mesh/MeshEntry.h b/mesh/MeshEntry.h
index 9737ead2..85521f17 100644
--- a/mesh/MeshEntry.h
+++ b/mesh/MeshEntry.h
@@ -19,7 +19,7 @@ enum MeshType {
CUBOID,
CYL, CYL_SHELL, CYL_SHELL_SEG,
SPHERE, SPHERE_SHELL, SPHERE_SHELL_SEG,
- TETRAHEDRON, DISK
+ TETRAHEDRON, DISK, ENDO
};
class ChemCompt;
@@ -58,6 +58,7 @@ class MeshEntry
* 6: spherical shell
* 7: spherical shell segment
* 8: Tetrahedral
+ * 9: EndoMesh
*/
unsigned int getMeshType( const Eref& e ) const;
diff --git a/mesh/NeuroMesh.cpp b/mesh/NeuroMesh.cpp
index 2cf28411..cb0ccc56 100644
--- a/mesh/NeuroMesh.cpp
+++ b/mesh/NeuroMesh.cpp
@@ -23,6 +23,7 @@
#include "NeuroMesh.h"
#include "SpineEntry.h"
#include "SpineMesh.h"
+#include "EndoMesh.h"
#include "../utility/numutil.h"
#include "../utility/strutil.h"
#include "../shell/Wildcard.h"
@@ -1256,6 +1257,13 @@ void NeuroMesh::matchMeshEntries( const ChemCompt* other,
matchCubeMeshEntries( other, ret );
return;
}
+ const EndoMesh* em = dynamic_cast< const EndoMesh* >( other );
+ if ( em )
+ {
+ em->matchMeshEntries( this, ret );
+ flipRet( ret );
+ return;
+ }
const SpineMesh* sm = dynamic_cast< const SpineMesh* >( other );
if ( sm )
{
diff --git a/python/moose/moose_test.py b/python/moose/moose_test.py
index 9cda2b80..eeab9f94 100644
--- a/python/moose/moose_test.py
+++ b/python/moose/moose_test.py
@@ -87,6 +87,9 @@ class Command(object):
self.process.terminate()
thread.join()
+ if self.process.stderr is not None:
+ _logger.warn( '%s' % self.process.stderr.read() )
+
return self.process.returncode
def init_test_dir( ):
@@ -116,7 +119,7 @@ def run_test( index, testfile, timeout, **kwargs):
"""
global test_status_
global total_
- pyExec = os.environ.get( 'PYTHON_EXECUTABLE', '/usr/bin/python' )
+ pyExec = os.environ.get( 'PYTHON_EXECUTABLE', sys.executable )
cmd = Command( [ pyExec, testfile ] )
ti = time.time( )
diff --git a/python/rdesigneur/fixXreacs.py b/python/rdesigneur/fixXreacs.py
new file mode 100644
index 00000000..19379f5f
--- /dev/null
+++ b/python/rdesigneur/fixXreacs.py
@@ -0,0 +1,193 @@
+####################################################################
+# fixXreacs.py
+# The program is meant to take a model and reconfigure any cross-compartment
+# reactions so that they are split into portions on either side coupled
+# either by diffusion, or by a reaction-driven translocation process.
+#
+# Copyright (C) Upinder S. Bhalla 2018
+# This program is free software. It is licensed under the terms of
+# GPL version 3 or later at your discretion.
+# This program carries no warranty whatsoever.
+####################################################################
+
+
+import sys
+import moose
+
+msgSeparator = "_xMsg_"
+
+def findCompt( elm ):
+ elm = moose.element( elm )
+ pa = elm.parent
+ while pa.path != '/':
+ if moose.Neutral(pa).isA[ 'ChemCompt' ]:
+ return pa.path
+ pa = pa.parent
+ print( 'Error: No compartment parent found for ' + elm.path )
+ return '/'
+
+# http://stackoverflow.com/q/3844948/
+def checkEqual(lst):
+ return not lst or lst.count(lst[0]) == len(lst)
+
+def findXreacs( basepath, reacType ):
+ reacs = moose.wildcardFind( basepath + '/##[ISA=' + reacType + 'Base]' )
+ ret = []
+ for i in reacs:
+ reacc = findCompt( i )
+ subs = i.neighbors['subOut']
+ prds = i.neighbors['prdOut']
+ subc = [findCompt(j) for j in subs]
+ prdc = [findCompt(j) for j in prds]
+
+ enzc = []
+ if reacType == 'Enz':
+ enzc = [reacc]
+ if not checkEqual( subc + prdc + enzc ):
+ ret.append( [i, reacc, subs, subc, prds, prdc ])
+ return ret
+
+def removeEnzFromPool( pool ):
+ kids = moose.wildcardFind( pool.path + "/#" )
+ for i in kids:
+ if i.isA[ 'EnzBase' ]:
+ moose.delete( i )
+ elif i.isA[ 'Function' ]:
+ moose.delete( i )
+
+# If a pool is not in the same compt as reac, make a proxy in the reac
+# compt, connect it up, and disconnect the one in the old compt.
+def proxify( reac, reacc, direction, pool, poolc ):
+ reacc_elm = moose.element( reacc )
+ reac_elm = moose.element( reac )
+ # Preserve the rates which were set up for the x-compt reacn
+ #moose.showfield( reac )
+ dupname = pool.name + '_xfer_' + moose.element(poolc).name
+ #print "#############", pool, dupname, poolc
+ if moose.exists( reacc + '/' + dupname ):
+ duppool = moose.element( reacc + '/' + dupname )
+ else:
+ # This also deals with cases where the duppool is buffered.
+ duppool = moose.copy(pool, reacc_elm, dupname )
+ duppool.diffConst = 0 # diffusion only happens in original compt
+ removeEnzFromPool( duppool )
+ disconnectReactant( reac, pool, duppool )
+ moose.connect( reac, direction, duppool, 'reac' )
+ #moose.showfield( reac )
+ #moose.showmsg( reac )
+
+def enzProxify( enz, enzc, direction, pool, poolc ):
+ if enzc == poolc:
+ return
+ enze = moose.element( enz )
+ # kcat and k2 are indept of volume, just time^-1
+ km = enze.numKm
+ proxify( enz, enzc, direction, pool, poolc )
+ enze.numKm = km
+
+def reacProxify( reac, reacc, direction, pool, poolc ):
+ if reacc == poolc:
+ return
+ reac_elm = moose.element( reac )
+ kf = reac_elm.numKf
+ kb = reac_elm.numKb
+ proxify( reac, reacc, direction, pool, poolc )
+ reac_elm.numKf = kf
+ reac_elm.numKb = kb
+
+def identifyMsg( src, srcOut, dest ):
+ if src.isA[ 'ReacBase' ] or src.isA[ 'EnzBase' ]:
+ if srcOut == 'subOut':
+ return msgSeparator + src.path + ' sub ' + dest.path + ' reac'
+ if srcOut == 'prdOut':
+ return msgSeparator + src.path + ' prd ' + dest.path + ' reac'
+ return ''
+
+def disconnectReactant( reacOrEnz, reactant, duppool ):
+ outMsgs = reacOrEnz.msgOut
+ infoPath = duppool.path + '/info'
+ if moose.exists( infoPath ):
+ info = moose.element( infoPath )
+ else:
+ info = moose.Annotator( infoPath )
+
+ #moose.le( reactant )
+ notes = ""
+ #moose.showmsg( reacOrEnz )
+ for i in outMsgs:
+ #print "killing msg from {} to {}\nfor {} and {}".format( reacOrEnz.path, reactant.path, i.srcFieldsOnE1[0], i.srcFieldsOnE2[0] )
+ if i.e1 == reactant:
+ msgStr = identifyMsg( i.e2, i.e2.srcFieldsOnE2[0], i.e1 )
+ if len( msgStr ) > 0:
+ notes += msgStr
+ moose.delete( i )
+ elif i.e2 == reactant:
+ msgStr = identifyMsg( i.e1[0], i.srcFieldsOnE1[0], i.e2[0] )
+ if len( msgStr ) > 0:
+ notes += msgStr
+ moose.delete( i )
+ #print "MSGS to rebuild:", notes
+ info.notes += notes
+
+def fixXreacs( basepath ):
+ xr = findXreacs( basepath, 'Reac' )
+ xe = findXreacs( basepath, 'Enz' )
+
+ for i in (xr):
+ reac, reacc, subs, subc, prds, prdc = i
+ for j in range( len( subs )):
+ reacProxify( reac, reacc, 'sub', subs[j], subc[j] )
+ for j in range( len( prds )):
+ reacProxify( reac, reacc, 'prd', prds[j], prdc[j] )
+
+ for i in (xe):
+ reac, reacc, subs, subc, prds, prdc = i
+ for j in range( len( subs )):
+ enzProxify( reac, reacc, 'sub', subs[j], subc[j] )
+ for j in range( len( prds )):
+ enzProxify( reac, reacc, 'prd', prds[j], prdc[j] )
+
+#####################################################################
+
+def getOldRates( msgs ):
+ if len( msgs ) > 1 :
+ m1 = msgs[1].split( msgSeparator )[0]
+ elm = moose.element( m1.split( ' ' )[0] )
+ if elm.isA[ 'ReacBase' ]:
+ return [elm.numKf, elm.numKb]
+ elif elm.isA[ 'EnzBase' ]:
+ return [elm.numKm,]
+ print( "Warning: getOldRates did not have any messages" )
+ return [0,]
+
+def restoreOldRates( oldRates, msgs ):
+ #print oldRates, msgs
+ if len( msgs ) > 1 :
+ m1 = msgs[1].split( msgSeparator )[0]
+ elm = moose.element( m1.split( ' ' )[0] )
+ if elm.isA[ 'ReacBase' ]:
+ elm.numKf = oldRates[0]
+ elm.numKb = oldRates[1]
+ elif elm.isA[ 'enzBase' ]:
+ elm.numKm = oldRates[0]
+
+
+
+def restoreXreacs( basepath ):
+ proxyInfo = moose.wildcardFind( basepath + "/##/#_xfer_#/info" )
+ for i in proxyInfo:
+ msgs = i.notes.split( msgSeparator )
+ oldRates = getOldRates( msgs )
+ #print( "Deleting {}".format( i.parent.path ) )
+ #print msgs
+ moose.delete( i.parent )
+ for j in msgs:
+ if len( j ) > 0:
+ args = j.split( ' ' )
+ assert( len( args ) == 4 )
+ #moose.showfield( args[0] )
+ moose.connect( args[0], args[1], args[2], args[3] )
+ #print( "Reconnecting {}".format( args ) )
+ #moose.showfield( args[0] )
+ restoreOldRates( oldRates, msgs )
+
diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py
index 30906a25..0ef8d5c4 100644
--- a/python/rdesigneur/rdesigneur.py
+++ b/python/rdesigneur/rdesigneur.py
@@ -28,6 +28,7 @@ import sys
import time
import rdesigneur.rmoogli as rmoogli
+import rdesigneur.fixXreacs as fixXreacs
from rdesigneur.rdesigneurProtos import *
from moose.neuroml.NeuroML import NeuroML
@@ -196,6 +197,7 @@ class rdesigneur:
self._buildPlots()
self._buildMoogli()
self._buildStims()
+ self._configureHSolve()
self._configureClocks()
if self.verbose:
self._printModelStats()
@@ -896,7 +898,7 @@ class rdesigneur:
'inject':('CompartmentBase', 'setInject'),
'Ca':('CaConcBase', 'getCa'),
'n':('PoolBase', 'setN'),
- 'conc':('PoolBase''setConc')
+ 'conc':('PoolBase', 'setConc')
}
stims = moose.Neutral( self.modelPath + '/stims' )
k = 0
@@ -918,7 +920,13 @@ class rdesigneur:
moose.connect( func, 'valueOut', q, stimField )
################################################################
- # Utility function for setting up clocks.
+ def _configureHSolve( self ):
+ if not self.turnOffElec:
+ hsolve = moose.HSolve( self.elecid.path + '/hsolve' )
+ hsolve.dt = self.elecDt
+ hsolve.target = self.soma.path
+
+# Utility function for setting up clocks.
def _configureClocks( self ):
if self.turnOffElec:
elecDt = 1e6
@@ -937,9 +945,6 @@ class rdesigneur:
if not self.turnOffElec: # Assign the Function clock
moose.setClock( 12, self.funcDt )
moose.setClock( 18, self.chemPlotDt )
- hsolve = moose.HSolve( self.elecid.path + '/hsolve' )
- hsolve.dt = elecDt
- hsolve.target = self.soma.path
sys.stdout.flush()
################################################################
################################################################
@@ -1184,6 +1189,7 @@ class rdesigneur:
return
if not hasattr( self, 'dendCompt' ):
raise BuildError( "configureSolvers: no chem meshes defined." )
+ fixXreacs.fixXreacs( self.modelPath )
dmksolve = moose.Ksolve( self.dendCompt.path + '/ksolve' )
dmdsolve = moose.Dsolve( self.dendCompt.path + '/dsolve' )
dmstoich = moose.Stoich( self.dendCompt.path + '/stoich' )
@@ -1216,15 +1222,10 @@ class rdesigneur:
pmstoich.dsolve = pmdsolve
pmstoich.path = self.psdCompt.path + "/##"
+ # Here we should test what kind of geom we have to use
# Put in cross-compartment diffusion between ksolvers
dmdsolve.buildNeuroMeshJunctions( smdsolve, pmdsolve )
- # Put in cross-compartment reactions between ksolvers
- smstoich.buildXreacs( pmstoich )
- #pmstoich.buildXreacs( smstoich )
- smstoich.buildXreacs( dmstoich )
- dmstoich.filterXreacs()
- smstoich.filterXreacs()
- pmstoich.filterXreacs()
+ # All the Xreac stuff in the solvers is now eliminated.
# set up the connections so that the spine volume scaling can happen
self.elecid.setSpineAndPsdMesh( self.spineCompt, self.psdCompt)
diff --git a/python/rdesigneur/rdesigneurProtos.py b/python/rdesigneur/rdesigneurProtos.py
index 43213b09..1c22ea0c 100644
--- a/python/rdesigneur/rdesigneurProtos.py
+++ b/python/rdesigneur/rdesigneurProtos.py
@@ -325,6 +325,84 @@ def make_LCa( name = 'LCa', parent = '/library' ):
ygate.tableA = yA
ygate.tableB = yB
return Ca
+######################################################################
+
+# Derived from : squid/electronics.py
+# Description:
+# Author: Subhasis Ray
+# Maintainer:
+# Created: Wed Feb 22 00:53:38 2012 (+0530)
+# Version:
+# Last-Updated: Fri May 04 16:35:40 2018 (+0530)
+# By: Upi
+# Update #: 221
+
+# Change log:
+#
+# 2012-02-22 23:22:30 (+0530) Subha - the circuitry put in a class.
+# 2018-05-04 23:22:30 (+0530) Upi - Adapted for Rdesigneur
+#
+
+# Code:
+
+class ClampCircuit(moose.Neutral):
+ """Container for a Voltage-Clamp/Current clamp circuit."""
+ defaults = {
+ 'level1': 25.0e-3,
+ 'width1': 50.0e-3,
+ 'delay1': 2.0e-3,
+ 'delay2': 1e3,
+ 'trigMode': 0,
+ 'delay3': 1e6
+ }
+ def __init__(self, path ):
+ moose.Neutral.__init__(self, path)
+ '''
+ self.pulsegen = moose.PulseGen(path+"/pulse") # holding voltage/current generator
+ self.pulsegen.count = 2
+ self.pulsegen.baseLevel = -65.0e-3
+ self.pulsegen.firstLevel = -40.0e-3
+ self.pulsegen.firstWidth = 50.0e-3
+ self.pulsegen.firstDelay = 2.0e-3
+ self.pulsegen.secondDelay = 0.0
+ self.pulsegen.trigMode = 2
+ self.gate = moose.PulseGen(path+"/gate") # holding voltage/current generator
+ self.gate.level[0] = 1.0
+ self.gate.delay[0] = 0.0
+ self.gate.width[0] = 1e3
+ moose.connect(self.gate, 'output', self.pulsegen, 'input')
+ '''
+ self.lowpass = moose.RC(path+"/lowpass") # lowpass filter
+ self.lowpass.R = 1.0
+ self.lowpass.C = 0.03
+ self.vclamp = moose.DiffAmp(path+"/vclamp")
+ self.vclamp.gain = 1.0
+ self.vclamp.saturation = 1e10
+ self.iclamp = moose.DiffAmp(path+"/iclamp")
+ self.iclamp.gain = 0.0
+ self.iclamp.saturation = 1e10
+ self.pid = moose.PIDController(path+"/pid")
+ self.pid.gain = 0.5
+ self.pid.tauI = 0.02e-3
+ self.pid.tauD = 0.005e-3
+ self.pid.saturation = 1e7
+ # Connect voltage clamp circuitry
+ #moose.connect(self.pulsegen, "output", self.lowpass, "injectIn")
+ moose.connect(self.lowpass, "output", self.vclamp, "plusIn")
+ moose.connect(self.vclamp, "output", self.pid, "commandIn")
+ #moose.connect(compartment, "VmOut", self.pid, "sensedIn")
+ #moose.connect(self.pid, "output", compartment, "injectMsg")
+ addmsg1 = moose.Mstring( path + '/addmsg1' )
+ addmsg1.value = './pid output .. injectMsg'
+ addmsg2 = moose.Mstring( path + '/addmsg2' )
+ addmsg2.value = '.. VmOut ./pid sensedIn'
+
+def make_vclamp( name = 'Vclamp', parent = '/library' ):
+ if moose.exists( '/library/' + name ):
+ return
+ vclamp = ClampCircuit( parent + '/' + name )
+
+######################################################################
################################################################
# API function for building spine prototypes. Here we put in the
diff --git a/python/setup.cmake.py b/python/setup.cmake.py
index a3eaff22..d74a5242 100644
--- a/python/setup.cmake.py
+++ b/python/setup.cmake.py
@@ -21,6 +21,9 @@ import sys
# guidelines. This caused havoc on our OBS build.
from distutils.core import setup
+
+# Read version from VERSION created by cmake file. This file must be present for
+# setup.cmake.py to work perfectly.
script_dir = os.path.dirname( os.path.abspath( __file__ ) )
# Version file must be available. It MUST be written by cmake. Or create
diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp
index bd1e6ca8..dc45989f 100644
--- a/scheduling/Clock.cpp
+++ b/scheduling/Clock.cpp
@@ -958,10 +958,12 @@ void Clock::buildDefaultTick()
defaultTick_["ChemCompt"] = ~0U;
defaultTick_["Cinfo"] = ~0U;
defaultTick_["Clock"] = ~0U;
+ defaultTick_["ConcChan"] = ~0U;
defaultTick_["CubeMesh"] = ~0U;
defaultTick_["CylMesh"] = ~0U;
defaultTick_["DiagonalMsg"] = ~0U;
defaultTick_["Double"] = ~0U;
+ defaultTick_["EndoMesh"] = ~0U;
defaultTick_["Finfo"] = ~0U;
defaultTick_["Group"] = ~0U;
defaultTick_["HHGate"] = ~0U;
--
GitLab