From 1fd6b2a3cd15b6c4d864acb43e4557075e4a176c Mon Sep 17 00:00:00 2001 From: bhalla Date: Sun, 27 Mar 2022 12:45:33 +0530 Subject: [PATCH 1/3] Fixes to get the LSODA ODE engine to work with chem solvers. --- ksolve/VoxelPools.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/ksolve/VoxelPools.cpp b/ksolve/VoxelPools.cpp index e1c5e48e30..2fccf37e8c 100644 --- a/ksolve/VoxelPools.cpp +++ b/ksolve/VoxelPools.cpp @@ -110,6 +110,7 @@ void VoxelPools::advance( const ProcInfo* p ) if( getMethod() == "lsoda" ) { + size_t totVar = stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools(); vector yout(size()+1); pLSODA->lsoda_update( &VoxelPools::lsodaSys, size() , Svec(), yout , &t @@ -118,7 +119,7 @@ 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++) + for (size_t i = 0; i < totVar; i++) varS()[i] = yout[i+1]; if( lsodaState == 0 ) @@ -341,6 +342,12 @@ 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; ii < vp->size(); ii++ ) { + y[ii] = vp->Svec()[ii]; + } + vp->stoichPtr_->updateFuncs( y, t ); vp->updateRates( y, dydt ); } From 7709ae174f2c70e892e6f1674f737e8720105e10 Mon Sep 17 00:00:00 2001 From: bhalla Date: Wed, 6 Apr 2022 17:23:21 +0530 Subject: [PATCH 2/3] Fixes to VoxelPools to handle LSODA a bit better --- ksolve/VoxelPools.cpp | 4 ++++ python/rdesigneur/rdesigneur.py | 6 +++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/ksolve/VoxelPools.cpp b/ksolve/VoxelPools.cpp index 2fccf37e8c..a8198c8a90 100644 --- a/ksolve/VoxelPools.cpp +++ b/ksolve/VoxelPools.cpp @@ -106,6 +106,10 @@ const string VoxelPools::getMethod( ) void VoxelPools::advance( const ProcInfo* p ) { double t = p->currTime - p->dt; + if (t < p->dt ) + lsodaState = 1; + else + lsodaState = 3; Ksolve* k = reinterpret_cast( stoichPtr_->getKsolve().eref().data() ); if( getMethod() == "lsoda" ) diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py index 4026433b6b..e750280719 100644 --- a/python/rdesigneur/rdesigneur.py +++ b/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 From d5bdc1a163360652ee1fb2d7a44f74494bb94637 Mon Sep 17 00:00:00 2001 From: bhalla Date: Mon, 18 Jul 2022 11:04:20 +0530 Subject: [PATCH 3/3] LSODA bug-fixes. SeqSynHandler bug-fixes. --- basecode/ProcInfo.h | 13 +++++++++- ksolve/Stoich.cpp | 50 +++++++++++++++++++++++++++++++++++---- ksolve/Stoich.h | 11 +++++++++ ksolve/VoxelPools.cpp | 16 +++++++++++-- scheduling/Clock.cpp | 7 ++++++ synapse/SeqSynHandler.cpp | 46 +++++++++++++++++++++++++++++++---- 6 files changed, 132 insertions(+), 11 deletions(-) diff --git a/basecode/ProcInfo.h b/basecode/ProcInfo.h index 4a08e5ea75..af94b2c8a4 100644 --- a/basecode/ProcInfo.h +++ b/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; diff --git a/ksolve/Stoich.cpp b/ksolve/Stoich.cpp index 9bc4f7d57f..248959a2f3 100644 --- a/ksolve/Stoich.cpp +++ b/ksolve/Stoich.cpp @@ -92,9 +92,15 @@ const Cinfo* Stoich::initCinfo() static ReadOnlyValueFinfo numBufPools( "numBufPools", "Number of buffered pools to be computed by the " - "numerical engine. Includes pools controlled by functions.", + "numerical engine.", &Stoich::getNumBufPools); + static ReadOnlyValueFinfo numFuncPools( + "numFuncPools", + "Number of pools controlled by functions computed by the " + "numerical engine.", + &Stoich::getNumFuncPools); + static ReadOnlyValueFinfo 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& 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,7 +925,7 @@ void Stoich::allocateModel(const vector& 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::iterator i; @@ -898,6 +933,8 @@ void Stoich::buildPoolLookup() 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() diff --git a/ksolve/Stoich.h b/ksolve/Stoich.h index c4377879b0..ccd5eb3611 100644 --- a/ksolve/Stoich.h +++ b/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 funcTarget_; + /** + * vector of id of pools controlled by a func. + */ + vector funcTargetPoolVec_; + /** * Vector of funcs controlling pool increment, that is dN/dt * This is handled as a rateTerm. diff --git a/ksolve/VoxelPools.cpp b/ksolve/VoxelPools.cpp index a8198c8a90..45e6d5c261 100644 --- a/ksolve/VoxelPools.cpp +++ b/ksolve/VoxelPools.cpp @@ -106,10 +106,18 @@ 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 = 3; + 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( stoichPtr_->getKsolve().eref().data() ); if( getMethod() == "lsoda" ) @@ -123,6 +131,7 @@ void VoxelPools::advance( const ProcInfo* p ) // Now update the y from yout. This is different thant normal GSL or // BOOST based approach. + // totVar += stoichPtr_->getNumFuncPools(); for (size_t i = 0; i < totVar; i++) varS()[i] = yout[i+1]; @@ -348,11 +357,14 @@ void VoxelPools::lsodaSys( double t, double* y, double* dydt, void* 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; ii < vp->size(); ii++ ) { + 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 ); } diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp index 0bc3a8970b..1f5c2c83b0 100644 --- a/scheduling/Clock.cpp +++ b/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; } diff --git a/synapse/SeqSynHandler.cpp b/synapse/SeqSynHandler.cpp index c56101cf21..14cab13192 100644 --- a/synapse/SeqSynHandler.cpp +++ b/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,6 +579,9 @@ 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 ) @@ -576,13 +589,26 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p ) 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; @@ -590,6 +616,7 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p ) 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 );