From 5bb85aeccaeb2977a39bbfc0e5c9695292f8b84a Mon Sep 17 00:00:00 2001 From: bhalla Date: Tue, 9 Aug 2022 10:55:35 +0530 Subject: [PATCH 1/3] Added readout of per-voxel rate terms for reactions using Ksolve::rateVec --- diffusion/Dsolve.cpp | 8 ++++++ diffusion/Dsolve.h | 1 + kinetics/Reac.cpp | 29 +++++++++++++++------ ksolve/FuncRateTerm.h | 9 ++++++- ksolve/Gsolve.cpp | 8 ++++++ ksolve/Gsolve.h | 1 + ksolve/Ksolve.cpp | 53 +++++++++++++++++++++++++++++++++++++++ ksolve/Ksolve.h | 5 ++++ ksolve/KsolveBase.h | 4 +++ ksolve/Stoich.cpp | 10 ++++++++ ksolve/Stoich.h | 11 +++++--- ksolve/VoxelPoolsBase.cpp | 6 +++++ ksolve/VoxelPoolsBase.h | 1 + 13 files changed, 135 insertions(+), 11 deletions(-) diff --git a/diffusion/Dsolve.cpp b/diffusion/Dsolve.cpp index 68256762c6..6f8141f3a4 100644 --- a/diffusion/Dsolve.cpp +++ b/diffusion/Dsolve.cpp @@ -1158,6 +1158,14 @@ double Dsolve::getN( const Eref& e ) const return 0.0; } +double Dsolve::getR1( unsigned int reacIdx, const Eref& e ) const +{ + // Should not look for a reaction rate in the Dsolve. + return 0.0; +} + + + void Dsolve::setConcInit( const Eref& e, double v ) { unsigned int pid = convertIdToPoolIndex( e ); diff --git a/diffusion/Dsolve.h b/diffusion/Dsolve.h index 52d50a8489..368389ac36 100644 --- a/diffusion/Dsolve.h +++ b/diffusion/Dsolve.h @@ -130,6 +130,7 @@ class Dsolve: public KsolveBase void setConcInit( const Eref& e, double value ); double getN( const Eref& e ) const; void setN( const Eref& e, double value ); + double getR1( unsigned int reacIdx, const Eref& e ) const; double getVolumeOfPool( const Eref& e ) const; double getDiffConst( const Eref& e ) const; void setDiffConst( const Eref& e, double value ); diff --git a/kinetics/Reac.cpp b/kinetics/Reac.cpp index 7c861cece1..97cc74c261 100644 --- a/kinetics/Reac.cpp +++ b/kinetics/Reac.cpp @@ -208,12 +208,19 @@ void Reac::setNumKf( const Eref& e, double v ) double Reac::getNumKf( const Eref& e ) const { - // Return value for voxel 0. Conceivably I might want to use the - // DataId part to specify which voxel to use, but that isn't in the - // current definition for Reacs as being a single entity for the entire - // compartment. - double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 ); - return concKf_ / volScale; + // Updated Aug 2022: Return voxel-wise numKf. Really only needed when + // handling func adjustments to rate, earlier was just returning the + // local value of Kf. + if ( stoich_ ) { + return stoich_->getReacNumKf( e ); + } else { + // Return value for voxel 0. Conceivably I might want to use the + // DataId part to specify which voxel to use, but that isn't in the + // current definition for Reacs as being a single entity for the + // entire compartment. + double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 ); + return concKf_ / volScale; + } } void Reac::setNumKb( const Eref& e, double v ) @@ -239,7 +246,15 @@ void Reac::setConcKf( const Eref& e, double v ) double Reac::getConcKf( const Eref& e ) const { - return concKf_; + // Updated Aug 2022: Return voxel-wise numKf. Really only needed when + // handling func adjustments to rate, earlier was just returning the + // local value of Kf. + if ( stoich_ ) { + double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 ); + return stoich_->getReacNumKf( e ) * volScale; + } else { + return concKf_; + } } void Reac::setConcKb( const Eref& e, double v ) diff --git a/ksolve/FuncRateTerm.h b/ksolve/FuncRateTerm.h index 720026fe19..404da6d0ed 100644 --- a/ksolve/FuncRateTerm.h +++ b/ksolve/FuncRateTerm.h @@ -32,6 +32,7 @@ class FuncRate: public ExternReac double operator() ( const double* S ) const { double t = Field< double >::get( Id(1), "currentTime" ); auto v = (*func_)( S, t ); // get rate from func calculation. + *(const_cast< double * >( &k_ ) ) = v; assert(! std::isnan(v)); return v; } @@ -76,6 +77,10 @@ class FuncRate: public ExternReac return ret; } + double getR1() const { + return k_; + } + protected: double k_; shared_ptr func_; @@ -108,7 +113,9 @@ class FuncReac: public FuncRate double operator() ( const double* S ) const { // double ret = k_ * func_( S, 0.0 ); // get rate from func calculation. - double ret = (*func_)( S, 0.0 ); // get rate from func calculation. + double t = Field< double >::get( Id(1), "currentTime" ); + double ret = (*func_)( S, t ); //get rate from func calculation. + *(const_cast< double * >( &k_ ) ) = ret; vector< unsigned int >::const_iterator i; for ( i = v_.begin(); i != v_.end(); i++) { assert( !std::isnan( S[ *i ] ) ); diff --git a/ksolve/Gsolve.cpp b/ksolve/Gsolve.cpp index b58d10755c..e820180a30 100644 --- a/ksolve/Gsolve.cpp +++ b/ksolve/Gsolve.cpp @@ -969,6 +969,14 @@ double Gsolve::getN( const Eref& e ) const return 0.0; } +double Gsolve::getR1( unsigned int reacIdx, const Eref& e ) const +{ + unsigned int vox = getVoxelIndex( e ); + if ( vox != OFFNODE ) + return pools_[vox].getR1( reacIdx ); + return 0.0; +} + void Gsolve::setConcInit( const Eref& e, double v ) { unsigned int vox = getVoxelIndex( e ); diff --git a/ksolve/Gsolve.h b/ksolve/Gsolve.h index 90500e2d85..77c65cf5ce 100644 --- a/ksolve/Gsolve.h +++ b/ksolve/Gsolve.h @@ -98,6 +98,7 @@ class Gsolve: public KsolveBase void setN( const Eref& e, double v ); double getN( const Eref& e ) const; + double getR1( unsigned int reacIdx, const Eref& e ) const; void setConcInit( const Eref& e, double v ); double getConcInit( const Eref& e ) const; double getVolumeOfPool( const Eref& e ) const; diff --git a/ksolve/Ksolve.cpp b/ksolve/Ksolve.cpp index ef9b04381f..444c1a8e10 100644 --- a/ksolve/Ksolve.cpp +++ b/ksolve/Ksolve.cpp @@ -102,6 +102,12 @@ const Cinfo* Ksolve::initCinfo() &Ksolve::getNvec ); + static ReadOnlyLookupValueFinfo< Ksolve, string, vector< double > > rateVec( + "rateVec", + "vector of forward rate consts of specified reaction.", + &Ksolve::getRateVecFromPath + ); + static ValueFinfo< Ksolve, unsigned int > numAllVoxels( "numAllVoxels", "Number of voxels in the entire reac-diff system, " @@ -201,6 +207,7 @@ const Cinfo* Ksolve::initCinfo() &compartment, // Value &numLocalVoxels, // ReadOnlyValue &nVec, // LookupValue + &rateVec, // ReadOnlyLookupValue &numAllVoxels, // ReadOnlyValue &numPools, // Value &estimatedDt, // ReadOnlyValue @@ -490,6 +497,44 @@ void Ksolve::setNvec( unsigned int voxel, vector< double > nVec ) } } +/// Unlike getNvec, this returns vector of R1 for this reac across voxels +vector< double > Ksolve::getR1vec( unsigned int reacIdx ) const +{ + vector< double > ret( pools_.size(), 0.0 ); + for ( unsigned int ii = 0; ii < pools_.size(); ++ii ) { + ret[ii] = pools_[ii].getR1( reacIdx ); + } + return ret; +} + +/// Unlike getNvec, this returns vector of R1 for this reac across voxels +vector< double > Ksolve::getRateVecFromId( Id reacId ) const +{ + if ( reacId != Id() ) { + unsigned int idx = stoichPtr_->convertIdToReacIndex( reacId ); + if ( idx != ~0U ) + return getR1vec( idx ); + } + return vector< double >( pools_.size(), 0.0 ); +} + +/// Unlike getNvec, this returns vector of R1 for this reac across voxels +vector< double > Ksolve::getRateVecFromPath( string reacPath ) const +{ + Id reacId( reacPath ); + if ( reacId == Id() ) { + cout << "Error: object not found on " << reacPath << endl; + return vector< double >( pools_.size(), 0.0 ); + } + + if ( reacId != Id() ) { + unsigned int idx = stoichPtr_->convertIdToReacIndex( reacId ); + if ( idx != ~0U ) + return getR1vec( idx ); + } + return vector< double >( pools_.size(), 0.0 ); +} + double Ksolve::getEstimatedDt() const { @@ -714,6 +759,14 @@ double Ksolve::getN( const Eref& e ) const return 0.0; } +double Ksolve::getR1( unsigned int reacIdx, const Eref& e ) const +{ + unsigned int vox = getVoxelIndex( e ); + if ( vox != OFFNODE ) + return pools_[vox].getR1( reacIdx ); + return 0.0; +} + void Ksolve::setConcInit( const Eref& e, double v ) { unsigned int vox = getVoxelIndex( e ); diff --git a/ksolve/Ksolve.h b/ksolve/Ksolve.h index 53d63726b4..6c5d573f16 100644 --- a/ksolve/Ksolve.h +++ b/ksolve/Ksolve.h @@ -75,6 +75,10 @@ class Ksolve: public KsolveBase */ double getEstimatedDt() const; + vector getRateVecFromId( Id reacId ) const; //field func + vector getRateVecFromPath( string reacPath ) const; //field func + vector getR1vec( unsigned int reacIdx ) const; // Utility func + ////////////////////////////////////////////////////////////////// // Dest Finfos ////////////////////////////////////////////////////////////////// @@ -97,6 +101,7 @@ class Ksolve: public KsolveBase // KsolveBase inherited functions void setN( const Eref& e, double v ); double getN( const Eref& e ) const; + double getR1( unsigned int reacIdx, const Eref& e ) const; void setConcInit( const Eref& e, double v ); double getConcInit( const Eref& e ) const; diff --git a/ksolve/KsolveBase.h b/ksolve/KsolveBase.h index 7cb849fefd..0f9664ca21 100644 --- a/ksolve/KsolveBase.h +++ b/ksolve/KsolveBase.h @@ -30,6 +30,10 @@ class KsolveBase /// Get # of molecules in given pool and voxel. Varies with time. virtual double getN( const Eref& e ) const = 0; + /// Get rate const in a given reac and voxel. Usually fixed but + /// may vary if the reac is controlled by a function. + virtual double getR1( unsigned int reacIdx, const Eref& e ) const = 0; + /// Set buffer status of pool. Changes class between Pool and BufPool virtual void setIsBuffered( const Eref& e, bool val ) {;} diff --git a/ksolve/Stoich.cpp b/ksolve/Stoich.cpp index 248959a2f3..00ebd1ebd4 100644 --- a/ksolve/Stoich.cpp +++ b/ksolve/Stoich.cpp @@ -1682,6 +1682,16 @@ void Stoich::setReacKf(const Eref& e, double v) const } } +double Stoich::getReacNumKf(const Eref& e ) const +{ + unsigned int i = convertIdToReacIndex(e.id()); + if(i != ~0U) { + return kinterface_->getR1( i, e ); + // return rates_[i]->getR1(e); + } + return 0.0; +} + /** * For now assume a single rate term. */ diff --git a/ksolve/Stoich.h b/ksolve/Stoich.h index ccd5eb3611..5987294720 100644 --- a/ksolve/Stoich.h +++ b/ksolve/Stoich.h @@ -334,7 +334,7 @@ class Stoich { void setSpecies(unsigned int poolIndex, unsigned int s); /** - * Sets the forward rate v (given in millimoloar concentration units) + * Sets the forward rate v (given in millimolor concentration units) * for the specified reaction throughout the compartment in which the * reaction lives. Internally the stoich uses #/voxel units so this * involves querying the volume subsystem about volumes for each @@ -342,8 +342,13 @@ class Stoich { */ void setReacKf(const Eref& e, double v) const; + /** + * Return forward reaction rate in #/second units. + */ + double getReacNumKf(const Eref& e ) const; + /** - * Sets the reverse rate v (given in millimoloar concentration units) + * Sets the reverse rate v (given in millimolor concentration units) * for the specified reaction throughout the compartment in which the * reaction lives. Internally the stoich uses #/voxel units so this * involves querying the volume subsystem about volumes for each @@ -367,7 +372,7 @@ class Stoich { double getMMenzKcat(const Eref& e) const; /** - * Sets the rate v (given in millimoloar concentration units) + * Sets the rate v (given in millimolor concentration units) * for the forward enzyme reaction of binding substrate to enzyme. * Does this throughout the compartment in which the * enzyme lives. Internally the stoich uses #/voxel units so this diff --git a/ksolve/VoxelPoolsBase.cpp b/ksolve/VoxelPoolsBase.cpp index 62351fa171..4dd50f4a55 100644 --- a/ksolve/VoxelPoolsBase.cpp +++ b/ksolve/VoxelPoolsBase.cpp @@ -155,6 +155,12 @@ double VoxelPoolsBase::getN( unsigned int i ) const return S_[i]; } +double VoxelPoolsBase::getR1( unsigned int i ) const +{ + return rates_[i]->getR1(); +} + + void VoxelPoolsBase::setConcInit( unsigned int i, double v ) { Cinit_[i] = v; diff --git a/ksolve/VoxelPoolsBase.h b/ksolve/VoxelPoolsBase.h index b85d8643d1..374765da5d 100644 --- a/ksolve/VoxelPoolsBase.h +++ b/ksolve/VoxelPoolsBase.h @@ -94,6 +94,7 @@ class VoxelPoolsBase void setN( unsigned int i, double v ); double getN( unsigned int ) const; + double getR1( unsigned int ) const; void setConcInit( unsigned int, double v ); double getConcInit( unsigned int ) const; void setDiffConst( unsigned int, double v ); From 0cb99b90e89cfc96fa182c650390bffa2bc69a79 Mon Sep 17 00:00:00 2001 From: bhalla Date: Tue, 9 Aug 2022 11:19:06 +0530 Subject: [PATCH 2/3] Reverting to old getConcKf and getNumKf code for Reac, but the hooks for individual voxel access remain. --- kinetics/Reac.cpp | 29 +++++++---------------------- 1 file changed, 7 insertions(+), 22 deletions(-) diff --git a/kinetics/Reac.cpp b/kinetics/Reac.cpp index 97cc74c261..263cebb2f5 100644 --- a/kinetics/Reac.cpp +++ b/kinetics/Reac.cpp @@ -208,19 +208,12 @@ void Reac::setNumKf( const Eref& e, double v ) double Reac::getNumKf( const Eref& e ) const { - // Updated Aug 2022: Return voxel-wise numKf. Really only needed when - // handling func adjustments to rate, earlier was just returning the - // local value of Kf. - if ( stoich_ ) { - return stoich_->getReacNumKf( e ); - } else { - // Return value for voxel 0. Conceivably I might want to use the - // DataId part to specify which voxel to use, but that isn't in the - // current definition for Reacs as being a single entity for the - // entire compartment. - double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 ); - return concKf_ / volScale; - } + // Return value for voxel 0. Conceivably I might want to use the + // DataId part to specify which voxel to use, but that isn't in the + // current definition for Reacs as being a single entity for the + // entire compartment. + double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 ); + return concKf_ / volScale; } void Reac::setNumKb( const Eref& e, double v ) @@ -246,15 +239,7 @@ void Reac::setConcKf( const Eref& e, double v ) double Reac::getConcKf( const Eref& e ) const { - // Updated Aug 2022: Return voxel-wise numKf. Really only needed when - // handling func adjustments to rate, earlier was just returning the - // local value of Kf. - if ( stoich_ ) { - double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 ); - return stoich_->getReacNumKf( e ) * volScale; - } else { - return concKf_; - } + return concKf_; } void Reac::setConcKb( const Eref& e, double v ) From 66feaa1974da331b4e55894e5ebde9c26dc9085a Mon Sep 17 00:00:00 2001 From: bhalla Date: Tue, 30 Aug 2022 12:19:24 +0530 Subject: [PATCH 3/3] Fixes for LSODA to work with diffusion. --- ksolve/Ksolve.cpp | 4 +++- ksolve/VoxelPools.cpp | 22 ++++++++-------------- ksolve/VoxelPools.h | 2 +- ksolve/VoxelPoolsBase.cpp | 5 +++++ ksolve/VoxelPoolsBase.h | 7 +++++++ python/rdesigneur/rdesigneur.py | 6 ++++-- 6 files changed, 28 insertions(+), 18 deletions(-) diff --git a/ksolve/Ksolve.cpp b/ksolve/Ksolve.cpp index 444c1a8e10..8706e7e136 100644 --- a/ksolve/Ksolve.cpp +++ b/ksolve/Ksolve.cpp @@ -659,8 +659,10 @@ void Ksolve::reinit( const Eref& e, ProcPtr p ) if ( isBuilt_ ) { - for ( unsigned int i = 0 ; i < pools_.size(); ++i ) + for ( unsigned int i = 0 ; i < pools_.size(); ++i ) { + pools_[i].setNumVoxels( pools_.size() ); pools_[i].reinit( p->dt ); + } } else { diff --git a/ksolve/VoxelPools.cpp b/ksolve/VoxelPools.cpp index 45e6d5c261..68d24f9e25 100644 --- a/ksolve/VoxelPools.cpp +++ b/ksolve/VoxelPools.cpp @@ -36,6 +36,7 @@ using namespace boost::numeric; VoxelPools::VoxelPools() : pLSODA(nullptr) { + lsodaState_ = 1; #ifdef USE_GSL driver_ = 0; #endif @@ -57,6 +58,7 @@ VoxelPools::~VoxelPools() void VoxelPools::reinit( double dt ) { VoxelPoolsBase::reinit(); + lsodaState_ = 1; #ifdef USE_GSL if ( !driver_ ) return; @@ -106,27 +108,19 @@ 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( stoichPtr_->getKsolve().eref().data() ); if( getMethod() == "lsoda" ) { + // True if first step or restart, or if diffusion. Tells LSODA to + // recalculate using new pool n values, which slows it down a bit. + if ( p->isStart() || (numVoxels_ > 1) ) + lsodaState_ = 1; size_t totVar = stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools(); vector yout(size()+1); pLSODA->lsoda_update( &VoxelPools::lsodaSys, size() , Svec(), yout , &t - , p->currTime, &lsodaState, this + , p->currTime, &lsodaState_, this ); // Now update the y from yout. This is different thant normal GSL or @@ -135,7 +129,7 @@ void VoxelPools::advance( const ProcInfo* p ) for (size_t i = 0; i < totVar; i++) varS()[i] = yout[i+1]; - if( lsodaState == 0 ) + if( lsodaState_ == 0 ) { cerr << "Error: VoxelPools::advance: LSODA integration error at time " << t << "\n"; diff --git a/ksolve/VoxelPools.h b/ksolve/VoxelPools.h index 07cf4b8d98..ba034d9068 100644 --- a/ksolve/VoxelPools.h +++ b/ksolve/VoxelPools.h @@ -101,7 +101,7 @@ class VoxelPools: public VoxelPoolsBase std::shared_ptr pLSODA; LSODA_ODE_SYSTEM_TYPE lsodaSystem; - int lsodaState = 1; + int lsodaState_; #ifdef USE_GSL gsl_odeiv2_driver* driver_; diff --git a/ksolve/VoxelPoolsBase.cpp b/ksolve/VoxelPoolsBase.cpp index 4dd50f4a55..d9f00dd281 100644 --- a/ksolve/VoxelPoolsBase.cpp +++ b/ksolve/VoxelPoolsBase.cpp @@ -139,6 +139,11 @@ void VoxelPoolsBase::scaleVolsBufsRates(double ratio, const Stoich* stoichPtr) } } +void VoxelPoolsBase::setNumVoxels( unsigned int n ) +{ + numVoxels_ = n; +} + ////////////////////////////////////////////////////////////// // Zombie Pool Access functions ////////////////////////////////////////////////////////////// diff --git a/ksolve/VoxelPoolsBase.h b/ksolve/VoxelPoolsBase.h index 374765da5d..7cec73f7f2 100644 --- a/ksolve/VoxelPoolsBase.h +++ b/ksolve/VoxelPoolsBase.h @@ -205,12 +205,19 @@ class VoxelPoolsBase void scaleVolsBufsRates( double ratio, const Stoich* stoichPtr ); + void setNumVoxels( unsigned int ); + /// Debugging utility void print() const; protected: const Stoich* stoichPtr_; vector< RateTerm* > rates_; + /** + * Number of voxels. If > 1, set flag for LSODA to handle molecule + * flux on each timestep during diffusion, which slows it down. + */ + unsigned int numVoxels_; private: /** diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py index e750280719..7f93dd58be 100644 --- a/python/rdesigneur/rdesigneur.py +++ b/python/rdesigneur/rdesigneur.py @@ -133,7 +133,7 @@ def __init__(self, moogList = [], outputFileList = [], # List of all file save specifications. modelFileNameList = [], # List of any files used to build. - ode_method = "gsl", # gsl, lsoda, gssa, gillespie + ode_method = "lsoda", # gsl, lsoda, gssa, gillespie isLegacyMethod = False, params = None ): @@ -1588,10 +1588,12 @@ def _buildNeuroMesh( self ): self.chemid.name = 'temp_chem' newChemid = moose.Neutral( self.model.path + '/chem' ) comptlist = self._assignComptNamesFromKkit_SBML() - if len( comptlist ) == 1 and comptlist[0].name == 'kinetics': + #if len( comptlist ) == 1 and comptlist[0].name == 'kinetics': + if len( comptlist ) == 1: comptlist[0].name = 'dend' comptdict = { i.name:i for i in comptlist } if len(comptdict) == 1 or 'dend' in comptdict: + #print( "COMPTDICT = ", comptdict ) self.dendCompt = moose.NeuroMesh( newChemid.path + '/dend' ) self.dendCompt.geometryPolicy = 'cylinder' self.dendCompt.separateSpines = 0