diff --git a/.travis.yml b/.travis.yml
index 294e0896243791f983832cc2272705c819a76702..b5f3f688e66c00d97989013b897183b478ca3c70 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 e475b2e0bf622b47ac46a08f0b14eaa308c08b03..3009af21b0f185f8dbe2760711cc4b5663ec6f5f 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 349d580c9a7357cc4c243538a8e2ee228a66006c..bd6aaf84449898d713efba60273018fa55b0bf0d 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 9ed2be8632130943b4c1d91a366ea6ed77e2dfe4..b28218d5261621b28f143addbc0ddbee50728389 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 60bbb422fdf41448ff274596f27d1e58734d8bb8..7f7059a41e492d5daae67cdf334ba4c71bac9a70 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 6ce69f23a741d469e316a31da7cb476034c5bdec..b8b07f36101cec335e171759afe8c88beda34186 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 3023ddbb1324e26277463a7b71e11e16c39b87d1..abca9dbfcc54b67cbb5948aa61826ce2e3b503fb 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 0000000000000000000000000000000000000000..cd16266b5e5ea895ce7a2fec5722bb18de339a50
--- /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 c2ac3070ef8fc4de657554604afbc91d77f343c4..a4e44d23c88d6e9898a38881cdba23af434784af 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 e6faaa63133205c2dc3ec6be1e9f8197eb300e3c..cd1da649567100c4a3de2008d96c0d2d0ebd24ff 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 345403f691cb6a12bb581a576a43a62bfa664a72..587533eb35d0ba4715f4c7cf54947abb4359edff 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 9185ab6537ffab493c4df27be085420e8bc18671..61ffb1bc612c01dde5f917fd15faabd54169e0be 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 0a16e7727432fe9aea86e117512e02a12d9c2caa..d6028e62a4c7fb8bdbbe267fa87f54f8174f2c16 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 05d5ff8ecbe3bb3e0b35d71f806faf49826e6922..e4b9b0e74091ae27a2763d93dc604edecd682100 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 9022a1f0dc579f796fcab9928d4c1b27e2477163..56c15d2dc712371a5f0ac3003a9b0cc0ff3a1bcb 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 9d34c6380024df603ba9e27e0bd34bb2263b046a..ebc62b340c7ca12538a0a108f805ff1e2352d09c 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 77e0cc108dc882c059c9577c9ec863fd95958a9b..5038c6e4357a7b9a0b23fb2630a463a2d687db30 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 403f360349e71b083dd0d10c3383c4eae0a96945..b836312c7e44d67aacfddbcdca578517b42a6f75 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 0000000000000000000000000000000000000000..6217084580e80802edc8b61af191ddb5bee74b6a
--- /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 0000000000000000000000000000000000000000..6ebf48b417f9b8711d0805edf2a03d9a2f888d97
--- /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 0c4c136ceae179987e213bdcf60aa89d1bc21198..0812374c50886a4fe10587186fda234f54eca9d3 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 791e05322e327fb7a1827d755d2dfec720e1323f..bd5f37af13e0979d5ab909accf9255bdee4e9322 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 2529919a8714fb64013324656d13a545bfda446d..d3b5bd4154067e43753f4a4e60d1453da2bb7b4d 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 42bee6f72bf4d4de08e0d8724388eb0ce7a56764..25a5c060e7fcc6cb14c6be0550a1121b5842d548 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 9cf59a008d9c9159b8470ce5f6ac859900d628c4..78f572abc3a32d03af85a67812f67bb34ee04398 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 438652d8af1f748a9226d4befabbd04a0de234ce..273964301a5e40d0a8b853fdee1b198ccd3c616f 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 1d3a3dbac7d18c341ddc3e64002c307753758121..716698a48037aab51363745fc19348b9536f415b 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 a7ee4c7adb45109227d2262613ec202426dcba6d..8b7b5ac9ddeaa261724ad8eda3954cfc496d2010 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 5e2153bee804a5fa10e812ff0c6bcfd4561b14f5..c7f4af8a378b0764fa17a26b9d731c9e12fb2daa 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 a94dfff9d6143a8e9b19d2a2080d8ef19fa7bd86..458151401e2c2a0bcb36a9f04671e1e3b3947891 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 ea5177093e7b18a99f858df0614f3e270cf98986..bea3d067f1b560863ac02f8240841274c1e75af4 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 88405e72994e7f1196152349a8904d2adc13e391..64b07be22f8ab7c16761881de85461b5f606433a 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 05129e361ddcecc1a6e3b4b1db56fc6699f2721b..bd1e5ec5003e23d33699723d04ed21340791458a 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 20fb4895bdc290fbc59187445218166e93a28688..b75254d23441f59160fc54e4ca62993adcf8bd95 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 10104fd2f5a54edde76928326d6e2215f6a23e58..fe5bb1d8d98e6b0dd93d8d1b888a4b6bf2c11b13 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 6046fc568dd566e42da372e50e3338ac999157b8..a04e31ddec7611a56febe8af1f51661d86ea6bcd 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 32ae58e2432ca26bcb67fb99636639a14ba3ca3a..1d2b3f7f289a3c1bd4399bf32115170612540198 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 eb77b965fc883bb585ff2f412ea1ae2646c80cb5..ba909d2e825285d4b1b89f031d29ab94f510680d 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 da90d72b5afef9776b522e50cad676be8bc2b8f5..728c69ba095b1176e9d731be2a9f70f69f50c6a6 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 ed2c6c4372832ed0ee1033ccaa2c03612f2c14a5..82cd895ffa2fcfdad6bc5aaee9318981c3c9e57b 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 0000000000000000000000000000000000000000..33c13bdbb4c6c7aa9ccb2d6d09e8eb5d05b0307d
--- /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 0000000000000000000000000000000000000000..5567bda72de18894cce29902d82d450711bc29b7
--- /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 9737ead2d72e9a5b6b9c5a382bd68ae232d88ad9..85521f17d6e03f99cc92659acdfd3afb8a1d6e12 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 2cf28411d06d8df1ed5d2b0edf922279e0eb2245..cb0ccc56a4ef31529ccafee18ddbcae3cc541c53 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 9cda2b80d81b4ffc20ff77bf0dae656c47763e7c..eeab9f9467b434bc71a2b6fdc2fe23f1e48cf342 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 0000000000000000000000000000000000000000..19379f5f189b60a02a4dcc3e17899509e5ce3e9c
--- /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 30906a2502bbe8143afb96d91bd8cc9ff3211af5..0ef8d5c45c81c01f2a230e9f6cdcce9ffc50cd9d 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 43213b09be3209e52442b68c79d5f886c96c25df..1c22ea0cbc68a695f7c653cf0f4a0bba3e023231 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 a3eaff221156e64e27cded558038f9fad2ac4e0a..d74a52426fa622a500bb935adeb8824c286b1734 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 bd1e6ca8b820e269d513958fcce55293002a502a..dc45989f3a781825554bb42f8958fc98a251f3c1 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;