#include <SteadyStateBoost.h>

Collaboration diagram for SteadyState:

Public Member Functions
bool	badStoichiometry () const

bool	badStoichiometry () const

void	classifyState (const double *T)

void	classifyState (const double *T)

void	fitConservationRules (boost::numeric::ublas::matrix< value_type_ > &U, const vector< value_type_ > &eliminatedTotal, vector< value_type_ > &yi)

double	getConvergenceCriterion () const

double	getConvergenceCriterion () const

double	getEigenvalue (const unsigned int i) const

double	getEigenvalue (const unsigned int i) const

unsigned int	getMaxIter () const

unsigned int	getMaxIter () const

unsigned int	getNiter () const

unsigned int	getNiter () const

unsigned int	getNnegEigenvalues () const

unsigned int	getNnegEigenvalues () const

unsigned int	getNposEigenvalues () const

unsigned int	getNposEigenvalues () const

unsigned int	getNumVarPools () const

unsigned int	getNumVarPools () const

unsigned int	getRank () const

unsigned int	getRank () const

unsigned int	getSolutionStatus () const

unsigned int	getSolutionStatus () const

unsigned int	getStateType () const

unsigned int	getStateType () const

string	getStatus () const

string	getStatus () const

Id	getStoich () const

Id	getStoich () const

double	getTotal (const unsigned int i) const

double	getTotal (const unsigned int i) const

bool	isInitialized () const

bool	isInitialized () const

void	randomizeInitialCondition (const Eref &e)

void	randomizeInitialCondition (const Eref &e)

void	resettleFunc ()

void	resettleFunc ()

void	setConvergenceCriterion (double value)

void	setConvergenceCriterion (double value)

void	setEigenvalue (double val, const unsigned int i)

void	setEigenvalue (double val, const unsigned int i)

void	setMaxIter (unsigned int value)

void	setMaxIter (unsigned int value)

void	setStoich (Id s)

void	setStoich (Id s)

void	settle (bool forceSetup)

void	settle (bool forceSetup)

void	settleFunc ()

void	settleFunc ()

void	setTotal (const unsigned int i, double val)

void	setTotal (const unsigned int i, double val)

void	setupMatrix ()

void	setupMatrix ()

void	showMatrices ()

void	showMatrices ()

void	showMatricesFunc ()

void	showMatricesFunc ()

	SteadyState ()

	SteadyState ()

	~SteadyState ()

	~SteadyState ()

Static Public Member Functions
static void	assignY (double *S)

static void	assignY (double *S)

static const Cinfo *	initCinfo ()

static const Cinfo *	initCinfo ()

static void	setMolN (double y, unsigned int i)

static void	setMolN (double y, unsigned int i)

Static Public Attributes
static const double	DELTA = 1e-6

static const double	EPSILON = 1e-9

Private Member Functions
void	setupSSmatrix ()

void	setupSSmatrix ()

Private Attributes
bool	badStoichiometry_

double	convergenceCriterion_

vector< double >	eigenvalues_

boost::numeric::ublas::matrix < value_type_ >	gamma_

bool	isInitialized_

bool	isSetup_

boost::numeric::ublas::matrix < value_type_ >	LU_

unsigned int	maxIter_

unsigned int	nIter_

unsigned int	nNegEigenvalues_

unsigned int	nPosEigenvalues_

boost::numeric::ublas::matrix < value_type_ >	Nr_

unsigned int	nReacs_

unsigned int	numFailed_

unsigned int	numVarPools_

size_t	numVarPools_

VoxelPools	pool_

unsigned int	rank_

bool	reassignTotal_

unsigned int	solutionStatus_

unsigned int	stateType_

string	status_

Id	stoich_

vector< double >	total_

Detailed Description

Definition at line 20 of file SteadyStateBoost.h.

Constructor & Destructor Documentation

SteadyState::SteadyState ( )

Definition at line 322 of file SteadyStateBoost.cpp.

     :
     nIter_( 0 ),
     maxIter_( 100 ),
     badStoichiometry_( 0 ),
     status_( "OK" ),
     isInitialized_( 0 ),
     isSetup_( 0 ),
     convergenceCriterion_( 1e-7 ),
     stoich_(),
     numVarPools_( 0 ),
     nReacs_( 0 ),
     rank_( 0 ),
     reassignTotal_( 0 ),
     nNegEigenvalues_( 0 ),
     nPosEigenvalues_( 0 ),
     stateType_( 0 ),
     solutionStatus_( 0 ),
     numFailed_( 0 )
 {
 }

SteadyState::~SteadyState ( )

Definition at line 344 of file SteadyStateBoost.cpp.

 {
 
 }

SteadyState::SteadyState ( )

SteadyState::~SteadyState ( )

Member Function Documentation

static void SteadyState::assignY ( double * S )

static

void SteadyState::assignY ( double * S )

static

Definition at line 505 of file SteadyStateBoost.cpp.

506 {

507 }

bool SteadyState::badStoichiometry ( ) const

Definition at line 380 of file SteadyStateBoost.cpp.

References badStoichiometry_.

Referenced by initCinfo().

 {
     return badStoichiometry_;
 }

Here is the caller graph for this function:

void SteadyState::classifyState ( const double * T )

This does the iteration, using the specified method. First try gsl_multiroot_fsolver_hybrids If that doesn't work try gsl_multiroot_fsolver_dnewton Returns the gsl status.

Definition at line 609 of file SteadyStateBoost.cpp.

References Eref::data(), DELTA, eigenvalues_, EPSILON, Id::eref(), Stoich::getKsolve(), nNegEigenvalues_, nPosEigenvalues_, numVarPools_, pool_, rank_, solutionStatus_, stateType_, stoich_, and VoxelPools::updateRates().

 {
     /* column_major trait is needed for fortran */
     ublas::matrix<double, ublas::column_major> J(numVarPools_, numVarPools_);
 
     double tot = 0.0;
     Stoich* s = reinterpret_cast< Stoich* >( stoich_.eref().data() );
     vector< double > nVec = LookupField< unsigned int, vector< double > >::get(
                                 s->getKsolve(), "nVec", 0 );
     for ( unsigned int i = 0; i < numVarPools_; ++i )
         tot += nVec[i];
 
     tot *= DELTA;
 
     vector< double > yprime( nVec.size(), 0.0 );
     // Fill up Jacobian
     for ( unsigned int i = 0; i < numVarPools_; ++i )
     {
         double orig = nVec[i];
         if ( std::isnan( orig ) or std::isinf( orig ) )
         {
             cout << "Warning: SteadyState::classifyState: orig=nan\n";
             solutionStatus_ = 2; // Steady state OK, eig failed
             J.clear();
             return;
         }
         if ( std::isnan( tot ) or std::isinf( tot ))
         {
             cout << "Warning: SteadyState::classifyState: tot=nan\n";
             solutionStatus_ = 2; // Steady state OK, eig failed
             J.clear();
             return;
         }
         nVec[i] = orig + tot;
 
         pool_.updateRates( &nVec[0], &yprime[0] );
         nVec[i] = orig;
 
         // Assign the rates for each mol.
         for ( unsigned int j = 0; j < numVarPools_; ++j )
         {
             if( std::isnan( yprime[j] ) or std::isinf( yprime[j] ) )
             {
                 cout << "Warning: Overflow/underflow. Can't continue " << endl;
                 solutionStatus_ = 2;
                 return;
             }
             J(i, j) = yprime[j];
         }
     }
 
     // Jacobian is now ready. Find eigenvalues.
     ublas::vector< std::complex< double > > eigenVec ( J.size1() );
 
     ublas::matrix< std::complex<double>, ublas::column_major >* vl, *vr;
     vl = NULL; vr = NULL;
 
     /*-----------------------------------------------------------------------------
      *  INFO: Calling lapack routine geev to compute eigen vector of matrix J.
      *
      *  Argument 3 and 4 are left- and right-eigenvectors. Since we do not need
      *  them, they are set to NULL. Argument 2 holds eigen-vector and result is
      *  copied to it (output ).
      *-----------------------------------------------------------------------------*/
     int status = lapack::geev( J, eigenVec, vl, vr, lapack::optimal_workspace() );
 
     eigenvalues_.clear();
     eigenvalues_.resize( numVarPools_, 0.0 );
     if ( status != 0 )
     {
         cout << "Warning: SteadyState::classifyState failed to find eigenvalues. Status = " <<
              status << endl;
         solutionStatus_ = 2; // Steady state OK, eig classification failed
     }
     else     // Eigenvalues are ready. Classify state.
     {
         nNegEigenvalues_ = 0;
         nPosEigenvalues_ = 0;
         for ( unsigned int i = 0; i < numVarPools_; ++i )
         {
             std::complex<value_type> z = eigenVec[ i ];
             double r = z.real();
             nNegEigenvalues_ += ( r < -EPSILON );
             nPosEigenvalues_ += ( r > EPSILON );
             eigenvalues_[i] = r;
             // We have a problem here because numVarPools_ usually > rank
             // This means we have several zero eigenvalues.
         }
         if ( nNegEigenvalues_ == rank_ )
             stateType_ = 0; // Stable
         else if ( nPosEigenvalues_ == rank_ ) // Never see it.
             stateType_ = 1; // Unstable
         else  if (nPosEigenvalues_ == 1)
             stateType_ = 2; // Saddle
         else if ( nPosEigenvalues_ >= 2 )
             stateType_ = 3; // putative oscillatory
         else if ( nNegEigenvalues_ == ( rank_ - 1) && nPosEigenvalues_ == 0 )
             stateType_ = 4; // one zero or unclassified eigenvalue. Messy.
         else
             stateType_ = 5; // Other
     }
 }

Here is the call graph for this function:

void SteadyState::fitConservationRules	(	boost::numeric::ublas::matrix< value_type_ > &	U,
		const vector< value_type_ > &	eliminatedTotal,
		vector< value_type_ > &	yi
	)

This does the actual work of generating random numbers and making sure they fit.

Definition at line 1030 of file SteadyStateBoost.cpp.

References EPSILON, moose::mtrand(), numVarPools_, and total_.

Referenced by randomizeInitialCondition().

 {
     int numConsv = total_.size();
     int lastJ = numVarPools_;
     for ( int i = numConsv - 1; i >= 0; --i )
     {
         for ( unsigned int j = 0; j < numVarPools_; ++j )
         {
             double g = U( i, j );
             if ( fabs( g ) > EPSILON )
             {
                 // double ytot = calcTot( g, i, j, lastJ );
                 double ytot = 0.0;
                 for ( int k = j; k < lastJ; ++k )
                 {
                     y[k] = moose::mtrand();
                     ytot += y[k] * U( i, k );
                 }
                 assert( fabs( ytot ) > EPSILON );
                 double lastYtot = 0.0;
                 for ( unsigned int k = lastJ; k < numVarPools_; ++k )
                 {
                     lastYtot += y[k] * U( i, k );
                 }
                 double scale = ( eliminatedTotal[i] - lastYtot ) / ytot;
                 for ( int k = j; k < lastJ; ++k )
                 {
                     y[k] *= scale;
                 }
                 lastJ = j;
                 break;
             }
         }
     }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double SteadyState::getConvergenceCriterion ( ) const

Definition at line 450 of file SteadyStateBoost.cpp.

References convergenceCriterion_.

Referenced by initCinfo().

 {
     return convergenceCriterion_;
 }

Here is the caller graph for this function:

double SteadyState::getEigenvalue ( const unsigned int i ) const

Definition at line 476 of file SteadyStateBoost.cpp.

References eigenvalues_.

Referenced by initCinfo().

 {
     if ( i < eigenvalues_.size() )
         return eigenvalues_[i];
     cout << "Warning: SteadyState::getEigenvalue: index " << i <<
          " out of range " << eigenvalues_.size() << endl;
     return 0.0;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getMaxIter ( ) const

Definition at line 400 of file SteadyStateBoost.cpp.

References maxIter_.

Referenced by initCinfo().

 {
     return maxIter_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getNiter ( ) const

Definition at line 390 of file SteadyStateBoost.cpp.

References nIter_.

Referenced by initCinfo().

 {
     return nIter_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getNnegEigenvalues ( ) const

Definition at line 425 of file SteadyStateBoost.cpp.

References nNegEigenvalues_.

Referenced by initCinfo().

 {
     return nNegEigenvalues_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getNposEigenvalues ( ) const

Definition at line 430 of file SteadyStateBoost.cpp.

References nPosEigenvalues_.

Referenced by initCinfo().

 {
     return nPosEigenvalues_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getNumVarPools ( ) const

Definition at line 415 of file SteadyStateBoost.cpp.

References numVarPools_.

Referenced by initCinfo().

 {
     return numVarPools_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getRank ( ) const

Definition at line 410 of file SteadyStateBoost.cpp.

References rank_.

Referenced by initCinfo().

 {
     return rank_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getSolutionStatus ( ) const

Definition at line 435 of file SteadyStateBoost.cpp.

References solutionStatus_.

Referenced by initCinfo().

 {
     return solutionStatus_;
 }

Here is the caller graph for this function:

unsigned int SteadyState::getStateType ( ) const

Definition at line 420 of file SteadyStateBoost.cpp.

References stateType_.

Referenced by initCinfo().

 {
     return stateType_;
 }

Here is the caller graph for this function:

string SteadyState::getStatus ( ) const

Definition at line 395 of file SteadyStateBoost.cpp.

References status_.

Referenced by initCinfo().

 {
     return status_;
 }

Here is the caller graph for this function:

Id SteadyState::getStoich ( ) const

Definition at line 353 of file SteadyStateBoost.cpp.

References stoich_.

Referenced by initCinfo().

 {
     return stoich_;
 }

Here is the caller graph for this function:

double SteadyState::getTotal ( const unsigned int i ) const

Definition at line 455 of file SteadyStateBoost.cpp.

References total_.

Referenced by initCinfo().

 {
     if ( i < total_.size() )
         return total_[i];
     cout << "Warning: SteadyState::getTotal: index " << i <<
          " out of range " << total_.size() << endl;
     return 0.0;
 }

Here is the caller graph for this function:

static const Cinfo* SteadyState::initCinfo ( )

static

const Cinfo * SteadyState::initCinfo ( )

static

This picks up the entire Stoich data structure static Finfo* gslShared[] = { new SrcFinfo( "reinitSrc", Ftype0() ), new DestFinfo( "assignStoich", Ftype1< void* >(), RFCAST( &SteadyState::assignStoichFunc ) ), new DestFinfo( "setMolN", Ftype2< double, unsigned int >(), RFCAST( &SteadyState::setMolN ) ), new SrcFinfo( "requestYsrc", Ftype0() ), new DestFinfo( "assignY", Ftype1< double* >(), RFCAST( &SteadyState::assignY ) ), };

These are the fields of the SteadyState class

Definition at line 86 of file SteadyStateBoost.cpp.

References badStoichiometry(), getConvergenceCriterion(), getEigenvalue(), getMaxIter(), getNiter(), getNnegEigenvalues(), getNposEigenvalues(), getNumVarPools(), getRank(), getSolutionStatus(), getStateType(), getStatus(), getStoich(), getTotal(), Neutral::initCinfo(), isInitialized(), randomizeInitialCondition(), resettleFunc(), setConvergenceCriterion(), setMaxIter(), setStoich(), settleFunc(), setTotal(), setupMatrix(), showMatrices(), and steadyStateCinfo.

 {
     // Field definitions
     static ValueFinfo< SteadyState, Id > stoich(
         "stoich",
         "Specify the Id of the stoichiometry system to use",
         &SteadyState::setStoich,
         &SteadyState::getStoich
     );
     static ReadOnlyValueFinfo< SteadyState, bool > badStoichiometry(
         "badStoichiometry",
         "Bool: True if there is a problem with the stoichiometry",
         &SteadyState::badStoichiometry
     );
     static ReadOnlyValueFinfo< SteadyState, bool > isInitialized(
         "isInitialized",
         "True if the model has been initialized successfully",
         &SteadyState::isInitialized
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int > nIter(
         "nIter",
         "Number of iterations done by steady state solver",
         &SteadyState::getNiter
     );
     static ReadOnlyValueFinfo< SteadyState, string > status(
         "status",
         "Status of solver",
         &SteadyState::getStatus
     );
     static ValueFinfo< SteadyState, unsigned int > maxIter(
         "maxIter",
         "Max permissible number of iterations to try before giving up",
         &SteadyState::setMaxIter,
         &SteadyState::getMaxIter
     );
     static ValueFinfo< SteadyState, double> convergenceCriterion(
         "convergenceCriterion",
         "Fractional accuracy required to accept convergence",
         &SteadyState::setConvergenceCriterion,
         &SteadyState::getConvergenceCriterion
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int > numVarPools(
         "numVarPools",
         "Number of variable molecules in reaction system.",
         &SteadyState::getNumVarPools
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int > rank(
         "rank",
         "Number of independent molecules in reaction system",
         &SteadyState::getRank
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int > stateType(
         "stateType",
         "0: stable; 1: unstable; 2: saddle; 3: osc?; 4: one near-zero eigenvalue; 5: other",
         &SteadyState::getStateType
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int >
     nNegEigenvalues (
         "nNegEigenvalues",
         "Number of negative eigenvalues: indicates type of solution",
         &SteadyState::getNnegEigenvalues
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int >
     nPosEigenvalues(
         "nPosEigenvalues",
         "Number of positive eigenvalues: indicates type of solution",
         &SteadyState::getNposEigenvalues
     );
     static ReadOnlyValueFinfo< SteadyState, unsigned int > solutionStatus(
         "solutionStatus",
         "0: Good; 1: Failed to find steady states; "
         "2: Failed to find eigenvalues",
         &SteadyState::getSolutionStatus
     );
     static LookupValueFinfo< SteadyState, unsigned int, double > total(
         "total",
         "Totals table for conservation laws. The exact mapping of"
         "this to various sums of molecules is given by the "
         "conservation matrix, and is currently a bit opaque."
         "The value of 'total' is set to initial conditions when"
         "the 'SteadyState::settle' function is called."
         "Assigning values to the total is a special operation:"
         "it rescales the concentrations of all the affected"
         "molecules so that they are at the specified total."
         "This happens the next time 'settle' is called.",
         &SteadyState::setTotal,
         &SteadyState::getTotal
     );
     static ReadOnlyLookupValueFinfo<
     SteadyState, unsigned int, double > eigenvalues(
         "eigenvalues",
         "Eigenvalues computed for steady state",
         &SteadyState::getEigenvalue
     );
     // MsgDest definitions
     static DestFinfo setupMatrix( "setupMatrix",
                                   "This function initializes and rebuilds the matrices used "
                                   "in the calculation.",
                                   new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
                                 );
 
     static DestFinfo settle( "settle",
                              "Finds the nearest steady state to the current initial "
                              "conditions. This function rebuilds the entire calculation "
                              "only if the object has not yet been initialized.",
                              new OpFunc0< SteadyState >( &SteadyState::settleFunc )
                            );
     static DestFinfo resettle( "resettle",
                                "Finds the nearest steady state to the current initial "
                                "conditions. This function rebuilds the entire calculation ",
                                new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
                              );
     static DestFinfo showMatrices( "showMatrices",
                                    "Utility function to show the matrices derived for the calculations on the reaction system. Shows the Nr, gamma, and total matrices",
                                    new OpFunc0< SteadyState >( &SteadyState::showMatrices )
                                  );
     static DestFinfo randomInit( "randomInit",
                                  "Generate random initial conditions consistent with the mass"
                                  "conservation rules. Typically invoked in order to scan"
                                  "states",
                                  new EpFunc0< SteadyState >(
                                      &SteadyState::randomizeInitialCondition )
                                );
     // Shared definitions
 
     static Finfo * steadyStateFinfos[] =
     {
         &stoich,                // Value
         &badStoichiometry,      // ReadOnlyValue
         &isInitialized,         // ReadOnlyValue
         &nIter,                 // ReadOnlyValue
         &status,                // ReadOnlyValue
         &maxIter,               // Value
         &convergenceCriterion,  // ReadOnlyValue
         &numVarPools,           // ReadOnlyValue
         &rank,                  // ReadOnlyValue
         &stateType,             // ReadOnlyValue
         &nNegEigenvalues,       // ReadOnlyValue
         &nPosEigenvalues,       // ReadOnlyValue
         &solutionStatus,        // ReadOnlyValue
         &total,                 // LookupValue
         &eigenvalues,           // ReadOnlyLookupValue
         &setupMatrix,           // DestFinfo
         &settle,                // DestFinfo
         &resettle,              // DestFinfo
         &showMatrices,          // DestFinfo
         &randomInit,            // DestFinfo
 
 
     };
 
     static string doc[] =
     {
         "Name", "SteadyState",
         "Author", "Upinder S. Bhalla, 2009, updated 2014, NCBS",
         "Description", "SteadyState: works out a steady-state value for "
         "a reaction system. "
         "This class uses the multidimensional root finder algorithms "
         "to find the fixed points closest to the "
         "current molecular concentrations. "
         "When it finds the fixed points, it figures out eigenvalues of "
         "the solution, as a way to help classify the fixed points. "
         "Note that the method finds unstable as well as stable fixed "
         "points.\n "
         "The SteadyState class also provides a utility function "
         "*randomInit()* to "
         "randomly initialize the concentrations, within the constraints "
         "of stoichiometry. This is useful if you are trying to find "
         "the major fixed points of the system. Note that this is "
         "probabilistic. If a fixed point is in a very narrow range of "
         "state space the probability of finding it is small and you "
         "will have to run many iterations with different initial "
         "conditions to find it.\n "
         "The numerical calculations used by the SteadyState solver are "
         "prone to failing on individual calculations. All is not lost, "
         "because the system reports the solutionStatus. "
         "It is recommended that you test this field after every "
         "calculation, so you can simply ignore "
         "cases where it failed and try again with different starting "
         "conditions.\n "
         "Another rule of thumb is that the SteadyState object is more "
         "likely to succeed in finding solutions from a new starting point "
         "if you numerically integrate the chemical system for a short "
         "time (typically under 1 second) before asking it to find the "
         "fixed point. "
     };
 
     static Dinfo< SteadyState > dinfo;
     static Cinfo steadyStateCinfo(
         "SteadyState",
         Neutral::initCinfo(),
         steadyStateFinfos,
         sizeof( steadyStateFinfos )/sizeof(Finfo *),
         &dinfo,
         doc, sizeof( doc ) / sizeof( string )
     );
 
     return &steadyStateCinfo;
 }

Here is the call graph for this function:

bool SteadyState::isInitialized ( ) const

Definition at line 385 of file SteadyStateBoost.cpp.

References isInitialized_.

Referenced by initCinfo().

 {
     return isInitialized_;
 }

Here is the caller graph for this function:

void SteadyState::randomizeInitialCondition ( const Eref & e )

void SteadyState::randomizeInitialCondition ( const Eref & me )

Generates a new set of values for the S vector that is a) random and b) obeys the conservation rules.

Definition at line 968 of file SteadyStateBoost.cpp.

References checkAboveZero(), EPSILON, fitConservationRules(), gamma_, Field< A >::get(), numVarPools_, rankUsingBoost(), recalcTotal(), stoich_, and total_.

Referenced by initCinfo().

 {
     Id ksolve = Field< Id >::get( stoich_, "ksolve" );
     vector< double > nVec =
         LookupField< unsigned int, vector< double > >::get(
                 ksolve,"nVec", 0 );
     int numConsv = total_.size();
     recalcTotal( total_, gamma_, &nVec[0] );
     // The reorderRows function likes to have an I matrix at the end of
     // numVarPools_, so we provide space for it, although only its first
     // column is used for the total vector.
     ublas::matrix<double> U(numConsv, numVarPools_ + numConsv, 0);
 
     for ( int i = 0; i < numConsv; ++i )
     {
         for ( unsigned int j = 0; j < numVarPools_; ++j )
             U(i, j) = gamma_(i, j);
 
         U(i, numVarPools_) = total_[i];
     }
     // Do the forward elimination
     int rank = rankUsingBoost( U );
     assert( rank = numConsv );
 
     vector< double > eliminatedTotal( numConsv, 0.0 );
     for ( int i = 0; i < numConsv; ++i )
         eliminatedTotal[i] = U( i, numVarPools_ );
 
     // Put Find a vector Y that fits the consv rules.
     vector< double > y( numVarPools_, 0.0 );
     do
     {
         fitConservationRules( U, eliminatedTotal, y );
     }
     while ( !checkAboveZero( y ) );
 
     // Sanity check. Try the new vector with the old gamma and tots
     for ( int i = 0; i < numConsv; ++i )
     {
         double tot = 0.0;
         for ( unsigned int j = 0; j < numVarPools_; ++j )
         {
             tot += y[j] * gamma_( i, j );
         }
         assert( fabs( tot - total_[i] ) / tot < EPSILON );
     }
 
     // Put the new values into S.
     // cout << endl;
     for ( unsigned int j = 0; j < numVarPools_; ++j )
     {
         nVec[j] = y[j];
         // cout << y[j] << " ";
     }
     LookupField< unsigned int, vector< double > >::set(
         ksolve,"nVec", 0, nVec );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::resettleFunc ( )

Definition at line 499 of file SteadyStateBoost.cpp.

References settle().

Referenced by initCinfo().

 {
     settle( 1 );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::setConvergenceCriterion ( double value )

Definition at line 440 of file SteadyStateBoost.cpp.

References convergenceCriterion_, and value.

Referenced by initCinfo().

 {
     if ( value > 1e-10 )
         convergenceCriterion_ = value;
     else
         cout << "Warning: Convergence criterion " << value <<
              " too small. Old value " <<
              convergenceCriterion_ << " retained\n";
 }

Here is the caller graph for this function:

void SteadyState::setEigenvalue	(	double	val,
		const unsigned int	i
	)

void SteadyState::setEigenvalue	(	double	val,
		const unsigned int	i
	)

void SteadyState::setMaxIter ( unsigned int value )

Definition at line 405 of file SteadyStateBoost.cpp.

References maxIter_, and value.

Referenced by initCinfo().

 {
     maxIter_ = value;
 }

Here is the caller graph for this function:

static void SteadyState::setMolN	(	double	y,
		unsigned int	i
	)

static

static void SteadyState::setMolN	(	double	y,
		unsigned int	i
	)

static

void SteadyState::setStoich ( Id s )

Definition at line 358 of file SteadyStateBoost.cpp.

References Element::cinfo(), Eref::data(), Id::element(), Id::eref(), Field< A >::get(), LookupField< L, A >::get(), Stoich::getCompartment(), Stoich::getNumCoreRates(), Stoich::getRateTerms(), Cinfo::isA(), isInitialized_, nReacs_, numVarPools_, pool_, VoxelPools::setStoich(), setupSSmatrix(), VoxelPoolsBase::setVolume(), stoich_, VoxelPools::updateAllRateTerms(), and value.

Referenced by initCinfo().

 {
     if ( !value.element()->cinfo()->isA( "Stoich" ) )
     {
         cout << "Error: SteadyState::setStoich: Must be of Stoich class\n";
         return;
     }
 
     stoich_ = value;
     Stoich* stoichPtr = reinterpret_cast< Stoich* >( value.eref().data());
     numVarPools_ = Field< unsigned int >::get( stoich_, "numVarPools" );
     nReacs_ = Field< unsigned int >::get( stoich_, "numRates" );
     setupSSmatrix();
     double vol = LookupField< unsigned int, double >::get(
                      stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
     pool_.setVolume( vol );
     pool_.setStoich( stoichPtr, 0 );
     pool_.updateAllRateTerms( stoichPtr->getRateTerms(),
                               stoichPtr->getNumCoreRates() );
     isInitialized_ = 1;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::settle ( bool forceSetup )

The settle function computes the steady state nearest the initial conditions.

Definition at line 735 of file SteadyStateBoost.cpp.

References convergenceCriterion_, gamma_, Field< A >::get(), init(), isInitialized_, isSetup_, isSolutionValid(), Nr_, nReacs_, numVarPools_, pool_, rank_, reassignTotal_, setupSSmatrix(), solutionStatus_, status_, stoich_, and total_.

Referenced by resettleFunc(), and settleFunc().

 {
 
     if ( !isInitialized_ )
     {
         cout << "Error: SteadyState object has not been initialized. No calculations done\n";
         return;
     }
 
     if ( forceSetup || isSetup_ == 0 )
         setupSSmatrix();
 
     // Setting up matrices and vectors for the calculation.
     unsigned int nConsv = numVarPools_ - rank_;
     double * T = (double *) calloc( nConsv, sizeof( double ) );
 
     // Setting up matrices and vectors for the calculation.
     Id ksolve = Field< Id >::get( stoich_, "ksolve" );
 
     ss = new NonlinearSystem( numVarPools_ );
 
     ss->ri.rank = rank_;
     ss->ri.num_reacs = nReacs_;
     ss->ri.num_mols = numVarPools_;
     ss->ri.T = T;
     ss->ri.Nr = Nr_;
     ss->ri.gamma = gamma_;
     ss->ri.pool = &pool_;
     ss->ri.nVec = LookupField< unsigned int, vector< double > >::get(
             ksolve,"nVec", 0
             );
     ss->ri.convergenceCriterion = convergenceCriterion_;
 
     // This gives the starting point for finding the solution.
     vector<double> init( numVarPools_ );
 
     // Instead of starting at sqrt( x ),
     for( size_t i = 0; i < numVarPools_; ++i )
         init[i] = max( 0.0, sqrt(ss->ri.nVec[i]) );
 
     ss->initialize<vector<double>>( init );
 
     // Fill up boundary condition values
     if ( reassignTotal_ )   // The user has defined new conservation values.
     {
         for (size_t i = 0; i < nConsv; ++i )
             T[i] = total_[i];
         reassignTotal_ = 0;
     }
     else
     {
         for ( size_t i = 0; i < nConsv; ++i )
             for ( size_t j = 0; j < numVarPools_; ++j )
                 T[i] += gamma_( i, j ) * ss->ri.nVec[ j ];
         total_.assign( T, T + nConsv );
     }
 
     vector< double > repair( numVarPools_, 0.0 );
 
     for ( unsigned int j = 0; j < numVarPools_; ++j )
         repair[j] = ss->ri.nVec[j];
 
 
     int status = 1;
 
     // Find roots . If successful, set status to 0.
     if( ss->find_roots_gnewton( ) )
         status = 0;
 
     if ( status == 0 && isSolutionValid( ss->ri.nVec ) )
     {
         solutionStatus_ = 0; // Good solution
 
         LookupField< unsigned int, vector< double > >::set(
             ksolve, "nVec", 0, ss->ri.nVec
             );
         // Check what we set
         vector<double> t = LookupField< unsigned int, vector< double > >::get(
              ksolve,"nVec", 0
              );
 
         //classifyState( T );
     }
     else
     {
         cout << "Warning: SteadyState iteration failed, status = " <<
              status_ << ", nIter = " << ss->ri.nIter << endl;
 
         for ( unsigned int j = 0; j < numVarPools_; j++ )
             ss->ri.nVec[j] = repair[j];
 
         solutionStatus_ = 1; // Steady state failed.
         LookupField< unsigned int, vector< double > >::set(
             ksolve, "nVec", 0, ss->ri.nVec
             );
 
     }
 
     // Clean up.
     free( T );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::settleFunc ( )

Definition at line 494 of file SteadyStateBoost.cpp.

References settle().

Referenced by initCinfo().

 {
     settle( 0 );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::setTotal	(	const unsigned int	i,
		double	val
	)

void SteadyState::setTotal	(	const unsigned int	i,
		double	val
	)

Definition at line 464 of file SteadyStateBoost.cpp.

References reassignTotal_, and total_.

Referenced by initCinfo().

 {
     if ( i < total_.size() )
     {
         total_[i] = val;
         reassignTotal_ = 1;
         return;
     }
     cout << "Warning: SteadyState::setTotal: index " << i <<
          " out of range " << total_.size() << endl;
 }

Here is the caller graph for this function:

void SteadyState::setupMatrix ( )

Definition at line 490 of file SteadyStateBoost.cpp.

References setupSSmatrix().

Referenced by initCinfo().

 {
     setupSSmatrix();
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::setupSSmatrix ( )

private

void SteadyState::setupSSmatrix ( )

private

Definition at line 525 of file SteadyStateBoost.cpp.

References gamma_, Field< A >::get(), isSetup_, LU_, Nr_, nReacs_, numVarPools_, rank_, rankUsingBoost(), stoich_, and total_.

Referenced by setStoich(), settle(), and setupMatrix().

 {
 
     if ( numVarPools_ == 0 || nReacs_ == 0 )
         return;
 
     int nTot = numVarPools_ + nReacs_;
 
     ublas::matrix<double> N(numVarPools_, nReacs_, 0.0);
 
     LU_ = ublas::matrix< double >( numVarPools_, nTot, 0.0);
 
     vector< int > entry = Field< vector< int > >::get(
                               stoich_, "matrixEntry" );
     vector< unsigned int > colIndex = Field< vector< unsigned int > >::get(
                                           stoich_, "columnIndex" );
     vector< unsigned int > rowStart = Field< vector< unsigned int > >::get(
                                           stoich_, "rowStart" );
 
     for ( unsigned int i = 0; i < numVarPools_; ++i )
     {
         LU_(i, i + nReacs_) = 1;
         unsigned int k = rowStart[i];
         for ( unsigned int j = 0; j < nReacs_; ++j )
         {
             double x = 0;
             if ( j == colIndex[k] && k < rowStart[i+1] )
             {
                 x = entry[k++];
             }
             N(i,j) = x;
             LU_(i,j) = x;
         }
     }
 
     // This function reorgranize LU_.
     rank_ = rankUsingBoost( LU_ );
     Nr_ = ublas::matrix< double >( rank_, nReacs_ );
     Nr_.assign( ublas::zero_matrix< double >( rank_, nReacs_ ) );
 
     unsigned int nConsv = numVarPools_ - rank_;
     if ( nConsv == 0 )
     {
         cout << "SteadyState::setupSSmatrix(): Number of conserved species == 0. Aborting\n";
         return;
     }
 
     // Fill up Nr.
     for ( unsigned int i = 0; i < rank_; i++)
         for ( unsigned int j = i; j < nReacs_; j++)
             Nr_(i,j) = LU_(i, j);
 
     gamma_ = ublas::matrix< double >( nConsv, numVarPools_, 0.0 );
 
     // Fill up gamma
     for ( unsigned int i = rank_; i < numVarPools_; ++i )
         for ( unsigned int j = 0; j < numVarPools_; ++j )
             gamma_(i-rank_, j) = LU_(i, j+nReacs_);
 
     // Fill up boundary condition values
     total_.resize( nConsv );
     total_.assign( nConsv, 0.0 );
 
     Id ksolve = Field< Id >::get( stoich_, "ksolve" );
     vector< double > nVec =
         LookupField< unsigned int, vector< double > >::get(
             ksolve,"nVec", 0 );
 
     if ( nVec.size() >= numVarPools_ )
     {
         for ( unsigned int i = 0; i < nConsv; ++i )
             for ( unsigned int j = 0; j < numVarPools_; ++j )
                 total_[i] += gamma_(i,j) * nVec[ j ];
         isSetup_ = 1;
     }
     else
     {
         cout << "Error: SteadyState::setupSSmatrix(): unable to get"
              "pool numbers from ksolve.\n";
         isSetup_ = 0;
     }
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SteadyState::showMatrices ( )

Definition at line 509 of file SteadyStateBoost.cpp.

References gamma_, isInitialized_, LU_, Nr_, numVarPools_, rank_, and total_.

Referenced by initCinfo().

 {
     if ( !isInitialized_ ) {
         cout << "SteadyState::showMatrices: Sorry, the system is not yet initialized.\n";
         return;
     }
     int numConsv = numVarPools_ - rank_;
     cout << "Totals:    ";
     for ( int i = 0; i < numConsv; ++i )
         cout << total_[i] << "  ";
     cout << endl;
     cout << "gamma " << gamma_ << endl;
     cout << "Nr " << Nr_ << endl;
     cout << "LU " << LU_ << endl;
 }