Field in Ksolve to view forward rates across voxels

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 );

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 );

kinetics/Reac.cpp

@@ -210,8 +210,8 @@ 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.
+	// current definition for Reacs as being a single entity for the
+	// entire compartment.
 	double volScale = convertConcToNumRateUsingMesh( e, subOut(), 0 );
 	return concKf_ / volScale;
 }

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<FuncTerm> 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 ] ) );

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 );

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;

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 );

ksolve/Ksolve.h

@@ -75,6 +75,10 @@ class Ksolve: public KsolveBase
      */
     double getEstimatedDt() const;
 
+    vector<double> getRateVecFromId( Id reacId ) const; //field func
+    vector<double> getRateVecFromPath( string reacPath ) const; //field func
+    vector<double> 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;

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 ) {;}
 

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.
  */

ksolve/Stoich.h

@@ -334,16 +334,21 @@ 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
      * voxel, and scaling accordingly.
      */
     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

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;

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 );