upi commit dec 2019 (#380)

* Fixed Ca channel proto, thanks to Harshavardan

* Fix to vclamp in rdesigneur. Fix to ConcChan handling involving many files

* Fix to HSolve so that assignment of passive properties Em, Cm and Rm reflects in the solver calculations without having to rebuild the entire solver.

* Fix conversion for ConcChan in ReadKkit, corresponding fix to calculations in Dsolve

* Fixed test_Xchan1.py to work with fixed ConcChan permeability units

Author: upibhalla

diffusion/ConcChanInfo.h

@@ -14,22 +14,26 @@
  */
 class ConcChanInfo
 {
-public:
-    ConcChanInfo();
-    ConcChanInfo( unsigned int my, unsigned int other,
-                  unsigned int chan, double perm, bool isSwapped )
-        : myPool( my ), otherPool( other ), chanPool( chan ),
-          swapped( isSwapped ),
-          permeability( perm )
-    {;}
+	public:
+		ConcChanInfo();
+		ConcChanInfo( unsigned int my, unsigned int other, 
+						unsigned int chan, double perm, bool isSwapped, bool _isLocal )
+				: myPool( my ), otherPool( other ), chanPool( chan ), 
+				swapped( isSwapped ),
+				isLocal( _isLocal ),
+				permeability( perm )
+		{;}
 
-    unsigned int myPool;
-    unsigned int otherPool;
-    unsigned int chanPool;
-    bool swapped;	// Flag records whether myPool/otherPool are
-    // swapped wrt inPool/outPool. System expects
-    // that inPool should be in same compt as chanPool.
-    double permeability;
+		unsigned int myPool; /// This is not Id, it is internal PoolIndex
+		/// Obtained using convertIdToPoolIndex.
+		unsigned int otherPool; /// Internal PoolIndex
+		unsigned int chanPool; /// Internal PoolIndex 
+		bool swapped;	// Flag records whether myPool/otherPool are 
+						// swapped wrt inPool/outPool. System expects
+						// that inPool should be in same compt as chanPool.
+		bool isLocal;	// Flag, records odd cases where the channel is
+						// entirely local to one compartment.
+		double permeability;
 };
 
 #endif // _CONC_CHAN_INFO_H

diffusion/Dsolve.cpp

@@ -412,6 +412,8 @@ void Dsolve::calcJnChan( const DiffJunction& jn, Dsolve* other, double dt )
 
     for ( unsigned int i = 0; i < jn.myChannels.size(); ++i )
     {
+		if ( channels_[ jn.myChannels[i] ].isLocal )
+			continue;
         ConcChanInfo& myChan = channels_[ jn.myChannels[i] ];
         DiffPoolVec& myDv = pools_[ myChan.myPool ];
         DiffPoolVec& otherDv = other->pools_[ myChan.otherPool ];
@@ -426,7 +428,7 @@ void Dsolve::calcJnChan( const DiffJunction& jn, Dsolve* other, double dt )
             double chanN = chanDv.getN( j->first );
 			// Stick in a conversion factor for the myN and otherN into
 			// concentrations. Note that SI is millimolar.
-            double perm = myChan.permeability * chanN * 1000.0 / NA;
+            double perm = myChan.permeability * chanN / NA;
             myN = integ( myN, perm * myN/j->firstVol,
                          perm * otherN/j->secondVol, dt );
             otherN += lastN - myN;	// Mass consv
@@ -446,6 +448,8 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt )
 {
     for ( unsigned int i = 0; i < jn.otherChannels.size(); ++i )
     {
+		if ( other->channels_[ jn.otherChannels[i] ].isLocal )
+			continue;
         ConcChanInfo& otherChan = other->channels_[ jn.otherChannels[i] ];
         // This is the DiffPoolVec for the pools on the other Dsolve,
         // the one with the channel.
@@ -465,7 +469,7 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt )
             double chanN = chanDv.getN( j->second );
 			// Stick in a conversion factor for the myN and otherN into
 			// concentrations. Note that SI is millimolar.
-            double perm = otherChan.permeability * chanN * 1000.0 / NA;
+            double perm = otherChan.permeability * chanN / NA;
             myN = integ( myN, perm * myN/j->firstVol,
                          perm * otherN/j->secondVol, dt );
             otherN += lastN - myN;	// Mass consv
@@ -480,6 +484,46 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt )
     }
 }
 
+void Dsolve::calcLocalChan( double dt )
+{
+	// There may be some chans which connect within the compartment.
+	// This is odd biologically, but happens for modelling convenience.
+	// The channels have a flag indicating if this is so. It is a simple
+	// calculation, but lacks numerical accuracy as it is outside the 
+	// solver framework. 
+	
+    ChemCompt* cc = reinterpret_cast< ChemCompt* >( compartment_.eref().data() );
+	const vector< double >& vols = cc->vGetVoxelVolume();
+	for (auto ch = channels_.begin(); ch != channels_.end(); ++ch ) {
+		if ( ch->isLocal ) {
+        	DiffPoolVec& myDv = pools_[ ch->myPool ];
+        	DiffPoolVec& otherDv = pools_[ ch->otherPool ];
+        	DiffPoolVec& chanDv = pools_[ ch->chanPool ];
+			assert( myDv.getNumVoxels() == numVoxels_ );
+			for (unsigned int i = 0; i < numVoxels_; ++i ) {
+            	double myN = myDv.getN( i );
+            	double lastN = myN;
+            	double otherN = otherDv.getN( i );
+            	double chanN = chanDv.getN( i );
+				double vol = vols[i];
+				// Stick in a conversion factor for the myN and otherN into
+				// concentrations. Note that SI is millimolar.
+            	double perm = ch->permeability * chanN / NA;
+            	myN = integ( myN, perm * myN/vol, perm * otherN/vol, dt );
+            	otherN += lastN - myN;	// Mass consv
+            	if ( otherN < 0.0 )   // Avoid negatives
+            	{
+                	myN += otherN;
+                	otherN = 0.0;
+            	}
+            	myDv.setN( i, myN );
+            	otherDv.setN( i, otherN );
+			}
+		}
+	}
+}
+
+
 /**
  * Computes flux through a junction between diffusion solvers.
  * Most used at junctions on spines and PSDs, but can also be used
@@ -525,6 +569,7 @@ void Dsolve::reinit( const Eref& e, ProcPtr p )
 
 void Dsolve::updateJunctions( double dt )
 {
+	calcLocalChan( dt );
     for (auto i = junctions_.begin(); i != junctions_.end(); ++i )
         calcJunction( *i, dt );
 }
@@ -614,7 +659,8 @@ void Dsolve::fillConcChans( const vector< ObjId >& chans )
         ret.clear();
         unsigned int outPoolValue = outPool.id.value();
         bool swapped = false;
-        if ( !( inPool.bad() or chanPool.bad() ) )
+        bool isLocal = false;
+        if ( !( inPool.bad() || chanPool.bad() ) )
         {
             unsigned int inPoolIndex = convertIdToPoolIndex( inPool.id );
             unsigned int chanPoolIndex = convertIdToPoolIndex(chanPool.id);
@@ -623,14 +669,22 @@ void Dsolve::fillConcChans( const vector< ObjId >& chans )
                 inPoolIndex = convertIdToPoolIndex( outPool.id );
                 outPoolValue = inPool.id.value();
                 swapped = true;
-            }
+            } else {
+            	unsigned int outPoolIndex = convertIdToPoolIndex( outPool.id );
+				// edge case where both in and out are on same compartment.
+				// Not much sense in bio, but a common convenience hack.
+				if ( outPoolIndex != ~0U ) {
+					outPoolValue = outPoolIndex;
+					isLocal = true;
+				}
+
+			}
             if ( ( inPoolIndex != ~0U) && (chanPoolIndex != ~0U ) )
             {
                 ConcChanInfo cci(
                     inPoolIndex, outPoolValue, chanPoolIndex,
                     Field< double >::get( *i, "permeability" ),
-                    ////// Fix it below ////////
-                    swapped
+                    swapped, isLocal
                 );
                 channels_.push_back( cci );
             }
@@ -972,27 +1026,29 @@ void Dsolve::mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other)
     unsigned int outIndex;
     for ( unsigned int i = 0; i < ch.size(); ++i )
     {
-        unsigned int chanIndex = ch[i].chanPool;
-        outIndex = otherSolve->convertIdToPoolIndex( ch[i].otherPool );
-        if ( (outIndex != ~0U) && (chanIndex != ~0U ) )
-        {
+		if ( ch[i].isLocal ) {
             jn.myChannels.push_back(i);
-            ch[i].otherPool = outIndex;	// replace the Id with the index.
-            ch[i].chanPool = chanIndex; //chanIndex may be on either Dsolve
+		} else {
+        	outIndex = otherSolve->convertIdToPoolIndex( ch[i].otherPool );
+        	if (outIndex != ~0U) {
+            	ch[i].otherPool = outIndex;	// replace Id with the index.
+            	jn.myChannels.push_back(i);
+			}
         }
     }
     // Now set up the other Dsolve.
     vector< ConcChanInfo >& ch2 = otherSolve->channels_;
     for ( unsigned int i = 0; i < ch2.size(); ++i )
     {
-        unsigned int chanIndex = ch2[i].chanPool;
-        outIndex = selfSolve->convertIdToPoolIndex( ch2[i].otherPool );
-        if ( (outIndex != ~0U) && (chanIndex != ~0U) )
-        {
+		if ( ch2[i].isLocal ) {
             jn.otherChannels.push_back(i);
-            ch2[i].otherPool = outIndex;  // replace the Id with the index
-            ch2[i].chanPool = chanIndex;  //chanIndex may be on either Dsolve
-        }
+		} else {
+        	outIndex = selfSolve->convertIdToPoolIndex( ch2[i].otherPool );
+        	if ( outIndex != ~0U ) {
+            	ch2[i].otherPool = outIndex;  // replace Id with the index
+            	jn.otherChannels.push_back(i);
+			}
+		}
     }
 }
 
@@ -1063,15 +1119,18 @@ void Dsolve::setNumAllVoxels( unsigned int num )
 
 unsigned int Dsolve::convertIdToPoolIndex( const Id id ) const
 {
+	bool verbose = false;
     unsigned int i  = id.value() - poolMapStart_;
     if ( i < poolMap_.size() )
     {
         return poolMap_[i];
     }
-    cout << "Warning: Dsolve::convertIdToPoollndex: Id out of range, (" <<
+	if ( verbose ) {
+    	cout << "Warning: Dsolve::convertIdToPoollndex: Id out of range, (" <<
          poolMapStart_ << ", " << id << ", " << id.path() << ", " <<
          poolMap_.size() + poolMapStart_ << "\n";
-    return 0;
+	}
+    return ~0U;
 }
 
 unsigned int Dsolve::convertIdToPoolIndex( const Eref& e ) const
@@ -1083,7 +1142,7 @@ void Dsolve::setN( const Eref& e, double v )
 {
     unsigned int pid = convertIdToPoolIndex( e );
     // Ignore silently, as this may be a valid pid for the ksolve to use.
-    if ( pid >= pools_.size() )
+    if ( pid == ~0U || pid >= pools_.size() )
         return;
     unsigned int vox = e.dataIndex();
     if ( vox < numVoxels_ )
@@ -1098,7 +1157,7 @@ void Dsolve::setN( const Eref& e, double v )
 double Dsolve::getN( const Eref& e ) const
 {
     unsigned int pid = convertIdToPoolIndex( e );
-    if ( pid >= pools_.size() ) return 0.0; // ignore silently
+    if ( pid == ~0U || pid >= pools_.size() ) return 0.0; //ignore silently
     unsigned int vox = e.dataIndex();
     if ( vox <  numVoxels_ )
     {
@@ -1112,7 +1171,7 @@ double Dsolve::getN( const Eref& e ) const
 void Dsolve::setNinit( const Eref& e, double v )
 {
     unsigned int pid = convertIdToPoolIndex( e );
-    if ( pid >= pools_.size() )  // Ignore silently
+    if ( pid == ~0U || pid >= pools_.size() )  // Ignore silently
         return;
     unsigned int vox = e.dataIndex();
     if ( vox < numVoxels_ )
@@ -1127,7 +1186,7 @@ void Dsolve::setNinit( const Eref& e, double v )
 double Dsolve::getNinit( const Eref& e ) const
 {
     unsigned int pid = convertIdToPoolIndex( e );
-    if ( pid >= pools_.size() ) return 0.0; // ignore silently
+    if ( pid == ~0U || pid >= pools_.size() ) return 0.0; //ignore silently
     unsigned int vox = e.dataIndex();
     if ( vox < numVoxels_ )
     {
@@ -1141,36 +1200,52 @@ double Dsolve::getNinit( const Eref& e ) const
 void Dsolve::setDiffConst( const Eref& e, double v )
 {
     unsigned int pid = convertIdToPoolIndex( e );
-    if ( pid >= pools_.size() )   // Ignore silently, out of range.
+    if ( pid == ~0U || pid >= pools_.size() )   // Ignore silently, out of range.
         return;
-    pools_[ convertIdToPoolIndex( e ) ].setDiffConst( v );
+    pools_[ pid ].setDiffConst( v );
 }
 
 double Dsolve::getDiffConst( const Eref& e ) const
 {
     unsigned int pid = convertIdToPoolIndex( e );
-    if ( pid >= pools_.size() )   // Ignore silently, out of range.
+    if ( pid == ~0U || pid >= pools_.size() )   // Ignore silently, out of range.
         return 0.0;
-    return pools_[ convertIdToPoolIndex( e ) ].getDiffConst();
+    return pools_[ pid ].getDiffConst();
 }
 
 void Dsolve::setMotorConst( const Eref& e, double v )
 {
     unsigned int pid = convertIdToPoolIndex( e );
-    if ( pid >= pools_.size() )   // Ignore silently, out of range.
+    if ( pid == ~0U || pid >= pools_.size() )   // Ignore silently, out of range.
         return;
-    pools_[ convertIdToPoolIndex( e ) ].setMotorConst( v );
+    pools_[ pid ].setMotorConst( v );
 }
 
-void Dsolve::setNumPools( unsigned int numPoolSpecies )
+void Dsolve::setNumVarTotPools( unsigned int var, unsigned int tot )
 {
     // Decompose numPoolSpecies here, assigning some to each node.
-    numTotPools_ = numPoolSpecies;
-    numLocalPools_ = numPoolSpecies;
+    numTotPools_ = tot;
+    numLocalPools_ = var;
     poolStartIndex_ = 0;
 
-    pools_.resize( numLocalPools_ );
-    for ( unsigned int i = 0 ; i < numLocalPools_; ++i )
+    pools_.resize( numTotPools_ );
+    for ( unsigned int i = 0 ; i < numTotPools_; ++i )
+    {
+        pools_[i].setNumVoxels( numVoxels_ );
+        // pools_[i].setId( reversePoolMap_[i] );
+        // pools_[i].setParent( me );
+    }
+}
+
+void Dsolve::setNumPools( unsigned int numVarPoolSpecies )
+{
+    // Decompose numPoolSpecies here, assigning some to each node.
+    numTotPools_ = numVarPoolSpecies;
+    numLocalPools_ = numVarPoolSpecies;
+    poolStartIndex_ = 0;
+
+    pools_.resize( numTotPools_ );
+    for ( unsigned int i = 0 ; i < numTotPools_; ++i )
     {
         pools_[i].setNumVoxels( numVoxels_ );
         // pools_[i].setId( reversePoolMap_[i] );

diffusion/Dsolve.h

@@ -136,6 +136,8 @@ class Dsolve: public ZombiePoolInterface
 
     void setNumPools( unsigned int num );
     unsigned int getNumPools() const;
+    void setNumVarTotPools( unsigned int var, unsigned int tot );
+
     unsigned int getNumLocalVoxels() const;
     VoxelPoolsBase* pools( unsigned int i );
     double volume( unsigned int i ) const;
@@ -185,6 +187,7 @@ class Dsolve: public ZombiePoolInterface
     void calcJnChan( const DiffJunction& jn, Dsolve* other, double dt );
     void calcOtherJnChan( const DiffJunction& jn, Dsolve* other,
                           double dt );
+	void calcLocalChan( double dt );
     void fillConcChans( const vector< ObjId >& chans );
 
     /**

hsolve/HSolveInterface.cpp

@@ -89,8 +89,10 @@ void HSolve::setCm( Id id, double value )
 {
     unsigned int index = localIndex( id );
     assert( index < tree_.size() );
-    // Also update data structures used for calculations.
     tree_[ index ].Cm = value;
+    // Also update data structures used for calculations.
+	assert( tree_.size() == compartment_.size() );
+	compartment_[index].CmByDt = 2.0 * value / dt_;
 }
 
 double HSolve::getEm( Id id ) const
@@ -104,8 +106,10 @@ void HSolve::setEm( Id id, double value )
 {
     unsigned int index = localIndex( id );
     assert( index < tree_.size() );
-    // Also update data structures used for calculations.
     tree_[ index ].Em = value;
+    // Also update data structures used for calculations.
+	assert( tree_.size() == compartment_.size() );
+	compartment_[index].EmByRm = value / tree_[index].Rm;
 }
 
 double HSolve::getRm( Id id ) const
@@ -119,8 +123,10 @@ void HSolve::setRm( Id id, double value )
 {
     unsigned int index = localIndex( id );
     assert( index < tree_.size() );
-    // Also update data structures used for calculations.
     tree_[ index ].Rm = value;
+    // Also update data structures used for calculations.
+	assert( tree_.size() == compartment_.size() );
+	compartment_[index].EmByRm = tree_[index].Em / value;
 }
 
 double HSolve::getRa( Id id ) const

kinetics/ReadKkit.cpp

@@ -1192,7 +1192,8 @@ Id ReadKkit::buildChan( const vector< string >& args )
 	//
 	double permeability = atof( args[ chanMap_["perm"] ].c_str() );
 	Id chan = shell_->doCreate( "ConcChan", pa, tail, 1 );
-	Field< double >::set( chan, "permeability", permeability );
+	// Convert from perm in uM in GENESIS, to mM for MOOSE.
+	Field< double >::set( chan, "permeability", permeability *1000.0 );
 	assert( chan != Id() );
 	string chanPath = clean.substr( 10 );
 	chanIds_[ chanPath ] = chan;