Zombie gone

basecode/ProcInfo.h

@@ -13,10 +13,21 @@ class ProcInfo
 {
     public:
         ProcInfo()
-            : dt( 1.0 ), currTime( 0.0 )
+            : dt( 1.0 ), currTime( 0.0 ), status_( 0 )
         {;}
         double dt;
         double currTime;
+		bool isFirstStep() const { return (status_==0x01); }; // first start
+		bool isStart() const { return (status_ & 0x03); }; // any start call
+		bool isContinue() const { return (status_ == 0x02); }; // Start call to continue from current nonzero time of simulation
+		bool isReinit() const { return (status_ == 0x04); };
+		void setFirstStep() { status_ = 0x01 ; };
+		void setContinue() { status_ = 0x02 ; };
+		void setRunning() { status_ = 0x0 ; };
+		void setReinit() { status_ = 0x04 ; };
+
+	private:
+		unsigned int status_;	// bit 0: firstStep. bit 1: continue Bit 2: reinit.
 };
 
 typedef const ProcInfo* ProcPtr;

ksolve/Stoich.cpp

@@ -92,9 +92,15 @@ const Cinfo* Stoich::initCinfo()
     static ReadOnlyValueFinfo<Stoich, unsigned int> numBufPools(
         "numBufPools",
         "Number of buffered pools to be computed by the "
-        "numerical engine. Includes pools controlled by functions.",
+        "numerical engine.",
         &Stoich::getNumBufPools);
 
+    static ReadOnlyValueFinfo<Stoich, unsigned int> numFuncPools(
+        "numFuncPools",
+        "Number of pools controlled by functions computed by the "
+        "numerical engine.",
+        &Stoich::getNumFuncPools);
+
     static ReadOnlyValueFinfo<Stoich, unsigned int> numAllPools(
         "numAllPools",
         "Total number of pools handled by the numerical engine. "
@@ -195,6 +201,7 @@ const Cinfo* Stoich::initCinfo()
         &allowNegative,      // Value
         &numVarPools,        // ReadOnlyValue
         &numBufPools,        // ReadOnlyValue
+        &numFuncPools,        // ReadOnlyValue
         &numAllPools,        // ReadOnlyValue
         &numProxyPools,      // ReadOnlyValue
         &poolIdMap,          // ReadOnlyValue
@@ -471,9 +478,14 @@ unsigned int Stoich::getNumBufPools() const
     return bufPoolVec_.size();
 }
 
+unsigned int Stoich::getNumFuncPools() const
+{
+    return funcTargetPoolVec_.size();
+}
+
 unsigned int Stoich::getNumAllPools() const
 {
-    return varPoolVec_.size() + offSolverPoolVec_.size() + bufPoolVec_.size();
+    return varPoolVec_.size() + offSolverPoolVec_.size() + funcTargetPoolVec_.size() + bufPoolVec_.size();
 }
 
 unsigned int Stoich::getNumProxyPools() const
@@ -789,6 +801,8 @@ void Stoich::allocateModelObject(Id id)
             poolFuncVec_.push_back(ei->id());
             // objMap_[ id.value() - objMapStart_ ] = numFunctions_;
             // ++numFunctions_;
+			assert( tgt.size() == 1 );
+			funcTargetPoolVec_.push_back( tgt[0] );
         }
     }
 }
@@ -801,9 +815,27 @@ void myUnique(vector<Id>& v)
     v.erase(last, v.end());
 }
 
+void Stoich::clearFuncTargetPools()
+{
+	vector< Id > temp;
+	for ( auto pid : varPoolVec_ ) {
+		if ( std::find( funcTargetPoolVec_.begin(), funcTargetPoolVec_.end(), pid ) == funcTargetPoolVec_.end() )
+			temp.push_back( pid );
+	}
+	varPoolVec_ = temp;
+
+	temp.clear();
+	for ( auto pid : bufPoolVec_ ) {
+		if ( std::find( funcTargetPoolVec_.begin(), funcTargetPoolVec_.end(), pid ) == funcTargetPoolVec_.end() )
+			temp.push_back( pid );
+	}
+	bufPoolVec_ = temp;
+}
+
 /// Using the computed array sizes, now allocate space for them.
 void Stoich::resizeArrays()
 {
+    myUnique(funcTargetPoolVec_);
     myUnique(varPoolVec_);
     myUnique(bufPoolVec_);
     myUnique(offSolverPoolVec_);
@@ -814,8 +846,10 @@ void Stoich::resizeArrays()
     myUnique(mmEnzVec_);
     myUnique(offSolverMMenzVec_);
 
+	clearFuncTargetPools();
+
     unsigned int totNumPools =
-        varPoolVec_.size() + bufPoolVec_.size() + +offSolverPoolVec_.size();
+        varPoolVec_.size() + bufPoolVec_.size() + funcTargetPoolVec_.size() + offSolverPoolVec_.size();
 
     species_.resize(totNumPools, 0);
 
@@ -841,6 +875,7 @@ void Stoich::deAllocateModel()
 {
     varPoolVec_.clear();
     bufPoolVec_.clear();
+    funcTargetPoolVec_.clear();
     reacVec_.clear();
     enzVec_.clear();
     mmEnzVec_.clear();
@@ -890,14 +925,16 @@ void Stoich::allocateModel(const vector<Id>& elist)
 ///////////////////////////////////////////////////////////////////
 void Stoich::buildPoolLookup()
 {
-    // The order of pools is: varPools, offSolverVarPools, bufPools.
+    // The order of pools is: varPools, offSolverVarPools, funcTargetPools, bufPools.
     poolLookup_.clear();
     int poolNum = 0;
     vector<Id>::iterator i;
     for(i = varPoolVec_.begin(); i != varPoolVec_.end(); ++i)
         poolLookup_[*i] = poolNum++;
     for(i = offSolverPoolVec_.begin(); i != offSolverPoolVec_.end(); ++i)
         poolLookup_[*i] = poolNum++;
+    for(i = funcTargetPoolVec_.begin(); i != funcTargetPoolVec_.end(); ++i)
+        poolLookup_[*i] = poolNum++;
     for(i = bufPoolVec_.begin(); i != bufPoolVec_.end(); ++i)
         poolLookup_[*i] = poolNum++;
 }
@@ -1228,6 +1265,11 @@ void Stoich::unZombifyPools()
         if(e && !e->isDoomed())
 			SetGet2< ObjId, ObjId >::set( *i, "setSolvers", root, root );
     }
+    for(auto i = funcTargetPoolVec_.begin(); i != funcTargetPoolVec_.end(); ++i) {
+        Element* e = i->element();
+        if(e && !e->isDoomed())
+			SetGet2< ObjId, ObjId >::set( *i, "setSolvers", root, root );
+    }
 }
 
 void Stoich::unZombifyModel()

ksolve/Stoich.h

@@ -70,6 +70,9 @@ class Stoich {
     /// Returns number of local buffered pools.
     unsigned int getNumBufPools() const;
 
+    /// Returns number of local func target pools.
+    unsigned int getNumFuncPools() const;
+
     /**
      *  Returns total number of pools. Includes the pools whose
      *  actual calculations happen on another solver, but are given a
@@ -191,6 +194,9 @@ class Stoich {
 	/// Clears out old model.
     void deAllocateModel();
 
+	/// Removes funcTargetPools from the list of varPools and bufPools
+	void clearFuncTargetPools();
+
     /// Functions to build the maps between Ids and internal indices
     void buildPoolLookup();
     void buildRateTermLookup();
@@ -625,6 +631,11 @@ class Stoich {
      */
     vector<unsigned int> funcTarget_;
 
+	/**
+	 * vector of id of pools controlled by a func.
+	 */
+    vector<Id> funcTargetPoolVec_;
+
     /**
      * Vector of funcs controlling pool increment, that is dN/dt
      * This is handled as a rateTerm.

ksolve/VoxelPools.cpp

@@ -106,10 +106,23 @@ const string VoxelPools::getMethod( )
 void VoxelPools::advance( const ProcInfo* p )
 {
     double t = p->currTime - p->dt;
+	if ( p->isStart() ) // True if first step or restart.
+		lsodaState = 1;
+	// lsodaState = 1;
+	/**
+	if (t < p->dt )
+		lsodaState = 1;
+	else
+		lsodaState = 1; // Earlier this was set to 3. This made it 
+		// impossible to alter non-buffered molecules from the script after
+		// the simulation had started, even if the control returned to the
+		// parser.
+	**/
     Ksolve* k = reinterpret_cast<Ksolve*>( stoichPtr_->getKsolve().eref().data() );
 
     if( getMethod() == "lsoda" )
     {
+    	size_t totVar = stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools();
         vector<double> yout(size()+1);
         pLSODA->lsoda_update( &VoxelPools::lsodaSys, size()
                 , Svec(), yout , &t
@@ -118,7 +131,8 @@ void VoxelPools::advance( const ProcInfo* p )
 
         // Now update the y from yout. This is different thant normal GSL or
         // BOOST based approach.
-        for (size_t i = 0; i < size(); i++)
+		// totVar += stoichPtr_->getNumFuncPools();
+        for (size_t i = 0; i < totVar; i++)
             varS()[i] = yout[i+1];
 
         if( lsodaState == 0 )
@@ -341,7 +355,16 @@ void VoxelPools::lsodaSys( double t, double* y, double* dydt, void* param)
 {
     VoxelPools* vp = reinterpret_cast< VoxelPools* >( param );
     // Fill in the values.
+	// Ensure that the buffered values are assigned to y.
+   	size_t totVar = vp->stoichPtr_->getNumVarPools() + vp->stoichPtr_->getNumProxyPools();
+	for( size_t ii = totVar + vp->stoichPtr_->getNumFuncPools(); ii < vp->size(); ii++ ) {
+		y[ii] = vp->Svec()[ii];
+	}
+
     vp->stoichPtr_->updateFuncs( y, t );
+	for( size_t ii = totVar; ii < totVar + vp->stoichPtr_->getNumFuncPools(); ii++ ) {
+		vp->Svec()[ii] = y[ii];
+	}
     vp->updateRates( y, dydt );
 }
 

python/rdesigneur/rdesigneur.py

@@ -868,7 +868,7 @@ def _parseComptField( self, comptList, plotSpec, knownFields ):
         if kf[0] in ['CaConcBase', 'ChanBase', 'NMDAChan', 'VClamp']:
             objList = self._collapseElistToPathAndClass( comptList, plotSpec.relpath, kf[0] )
             return objList, kf[1]
-        elif field in [ 'n', 'conc', 'volume', 'increment']:
+        elif field in [ 'n', 'conc', 'nInit', 'concInit', 'volume', 'increment']:
             path = plotSpec.relpath
             pos = path.find( '/' )
             if pos == -1:   # Assume it is in the dend compartment.
@@ -1534,6 +1534,10 @@ def _assignComptNamesFromKkit_SBML( self ):
         needs to talk both to PSD and to spine bulk.
         '''
         comptList = moose.wildcardFind( self.chemid.path + '/##[ISA=ChemCompt]' )
+        #if len( comptList ) == 0 and moose.exists( self.chemid.path + '/kinetics' ):
+        print( "LEN = ", len( comptList ) )
+        if len( comptList ) == 0:
+            print( "EMPTY comptlist, found kinetics" )
         oldNaming = len([i.name for i in comptList if (i.name.find( "compartment_") == 0)])
         if oldNaming == 0:
             return comptList

scheduling/Clock.cpp

@@ -712,6 +712,10 @@ void Clock::handleStep( const Eref& e, unsigned long numSteps )
     char now[80];
 
     buildTicks( e );
+	if (nSteps_ == 0 )
+		info_.setFirstStep();
+	else
+		info_.setContinue();
     assert( currentStep_ == nSteps_ );
     assert( activeTicks_.size() == activeTicksMap_.size() );
     nSteps_ += numSteps;
@@ -769,6 +773,7 @@ void Clock::handleStep( const Eref& e, unsigned long numSteps )
             ++k;
         }
 #endif
+		info_.setRunning();
 
         // When 10% of simulation is over, notify user when notify_ is set to
         // true.
@@ -810,6 +815,7 @@ void Clock::handleReinit( const Eref& e )
     doingReinit_ = true;
     // Curr time is end of current step.
     info_.currTime = 0.0;
+	info_.setReinit();
     vector< unsigned int >::const_iterator k = activeTicksMap_.begin();
     for ( vector< unsigned int>::iterator j =
                 activeTicks_.begin(); j != activeTicks_.end(); ++j )
@@ -819,6 +825,7 @@ void Clock::handleReinit( const Eref& e )
     }
 
     info_.dt = dt_;
+	info_.setRunning();
     doingReinit_ = false;
 }
 

synapse/SeqSynHandler.cpp

@@ -528,7 +528,9 @@ void SeqSynHandler::addSpike(unsigned int index, double time, double weight)
     // for the latestSpikes slice.
     //
     // Here we reorder the entries in latestSpikes by the synapse order.
-    latestSpikes_[ synapseOrder_[index] ] += weight;
+	// Don't need to do this here at all - just examine the events_ 
+	// at Process and populate latestSpikes_ there.
+    // latestSpikes_[ synapseOrder_[index] ] += weight;
 }
 
 double SeqSynHandler::getTopSpike( unsigned int index ) const
@@ -557,6 +559,14 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p )
 {
     // Here we look at the correlations and do something with them.
     int nh = numHistory();
+	vector< PreSynEvent > newEvents;
+    while( !events_.empty() && events_.top().time <= p->currTime )
+    {
+        newEvents.push_back( events_.top() );
+		latestSpikes_[ synapseOrder_[ events_.top().synIndex ] ] += 
+			events_.top().weight;
+        events_.pop();
+    }
 
     // Check if we need to do correlations at all.
     if ( nh > 0 && kernel_.size() > 0 )
@@ -569,27 +579,44 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p )
             history_.sumIntoRow( latestSpikes_, 0 );
             latestSpikes_.assign( vGetNumSynapses(), 0.0 );
 
+			// I don't understand why we iterate over nh here.
+			// We just want the current correlation.
+			/**
             // Build up the sum of correlations over time
             vector< double > correlVec( vGetNumSynapses(), 0.0 );
             for ( int i = 0; i < nh; ++i )
                 history_.correl( correlVec, kernel_[i], i );
             if ( sequenceScale_ > 0.0 )   // Sum all responses, send to chan
             {
                 seqActivation_ = 0.0;
-                for ( vector< double >::iterator y = correlVec.begin();
-                        y != correlVec.end(); ++y )
-                    seqActivation_ += pow( *y, sequencePower_ );
+                // for ( vector< double >::iterator y = correlVec.begin(); y != correlVec.end(); ++y )
+                //   seqActivation_ += pow( *y, sequencePower_ );
+				for (const auto y : correlVec )
+                    seqActivation_ += pow( y, sequencePower_ );
 
                 // We'll use the seqActivation_ to send a special msg.
                 seqActivation_ *= sequenceScale_;
             }
+			*/
+
+            seqActivation_ = 0.0;
+            if ( sequenceScale_ > 0.0 )   // Sum all responses, send to chan
+            {
+            	for ( int i = 0; i < nh; ++i )
+                	seqActivation_ += pow( history_.dotProduct( kernel_[i], i, 0 ), sequencePower_ );
+                seqActivation_ *= sequenceScale_;
+            }
+
+			/*
+			 * I think this is altering the correl vec for next cycle.
             if ( plasticityScale_ > 0.0 )   // Short term changes in individual wts
             {
                 weightScaleVec_ = correlVec;
                 for ( vector< double >::iterator y=weightScaleVec_.begin();
                         y != weightScaleVec_.end(); ++y )
                     *y *= plasticityScale_;
             }
+			*/
         }
     }
 
@@ -599,20 +626,31 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p )
     double activation = seqActivation_; // Start with seq activation
     if ( plasticityScale_ > 0.0 )
     {
+		for (const auto ee : newEvents ) {
+            activation += ee.weight * baseScale_ *
+                          weightScaleVec_[ ee.synIndex ] / p->dt;
+		}
+		/*
         while( !events_.empty() && events_.top().time <= p->currTime )
         {
             activation += events_.top().weight * baseScale_ *
                           weightScaleVec_[ events_.top().synIndex ] / p->dt;
             events_.pop();
         }
+		*/
     }
     else
     {
+		for (const auto ee : newEvents ) {
+            activation += baseScale_ * ee.weight / p->dt;
+		}
+		/*
         while( !events_.empty() && events_.top().time <= p->currTime )
         {
             activation += baseScale_ * events_.top().weight / p->dt;
             events_.pop();
         }
+		*/
     }
     if ( activation != 0.0 )
         SynHandlerBase::activationOut()->send( e, activation );