SteadyStateBoost.cpp

Go to the documentation of this file.

/***
 *       Filename:  SteadyStateBoost.cpp
 *
 *    Description:  Steady state solver. Copy of ../ksolve/SteadyStateGsl.cpp
 *    but it uses boost solver.
 *
 *        Created:  2016-05-10
 *
 *         Author:  Dilawar Singh <dilawars@ncbs.res.in>
 *   Organization:  NCBS Bangalore
 *
 * This program works out a steady-state value for a reaction system.
 * It uses boost-ublas and lapack heavily.
 *
 * It finds the ss value closest to the initial conditions.
 *
 * If you want to find multiple stable states, it is best to do this
 * in Python as it gives a lot of flexibility in working out how to
 * find steady states.
 *
 * Likewise, if you want to carry out a dose-response calculation.
 */


#include "header.h"
#include "global.h"
#include "SparseMatrix.h"
#include "KinSparseMatrix.h"
#include "RateTerm.h"
#include "FuncTerm.h"
#include "VoxelPoolsBase.h"
#include "../mesh/VoxelJunction.h"
#include "XferInfo.h"
#include "ZombiePoolInterface.h"
#include "Stoich.h"
#include "../randnum/RNG.h"

/* Root finding algorithm is implemented here */
#include "NonlinearSystem.h"

/*
 * Bindings to lapack. These headers are not part of standard boost
 * distribution. They are available with moose distribution. See 'external'
 * folder.
 */
#include "boost/numeric/bindings/lapack/lapack.hpp"
#include "boost/numeric/bindings/lapack/geev.hpp"

#include "VoxelPoolsBase.h"
#include "OdeSystem.h"
#include "VoxelPools.h"
#include "SteadyStateBoost.h"

using namespace boost::numeric::bindings;
using namespace boost::numeric;

void ss_func( const vector_type& x, void* params, vector_type& f );
unsigned int rankUsingBoost( ublas::matrix<double>& U );

// Limit below which small numbers are treated as zero.
const double SteadyState::EPSILON = 1e-9;

// This fraction of molecules is used as an increment in computing the
// Jacobian
const double SteadyState::DELTA = 1e-6;
struct reac_info
{

    int rank;
    int num_reacs;
    size_t num_mols;
    int nIter;
    double convergenceCriterion;

    double* T;
    VoxelPools* pool;
    vector< double > nVec;

    ublas::matrix< double >* Nr;
    ublas::matrix< double >* gamma;
};

const Cinfo* SteadyState::initCinfo()
{
    // 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;
}

static const Cinfo* steadyStateCinfo = SteadyState::initCinfo();

// Class function definitions

SteadyState::SteadyState()
    :
    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()
{

}

// Field function definitions

Id SteadyState::getStoich() const
{
    return stoich_;
}

void SteadyState::setStoich( Id value )
{
    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;
}

bool SteadyState::badStoichiometry() const
{
    return badStoichiometry_;
}

bool SteadyState::isInitialized() const
{
    return isInitialized_;
}

unsigned int SteadyState::getNiter() const
{
    return nIter_;
}

string SteadyState::getStatus() const
{
    return status_;
}

unsigned int SteadyState::getMaxIter() const
{
    return maxIter_;
}

void SteadyState::setMaxIter( unsigned int value )
{
    maxIter_ = value;
}

unsigned int SteadyState::getRank() const
{
    return rank_;
}

unsigned int SteadyState::getNumVarPools() const
{
    return numVarPools_;
}

unsigned int SteadyState::getStateType() const
{
    return stateType_;
}

unsigned int SteadyState::getNnegEigenvalues() const
{
    return nNegEigenvalues_;
}

unsigned int SteadyState::getNposEigenvalues() const
{
    return nPosEigenvalues_;
}

unsigned int SteadyState::getSolutionStatus() const
{
    return solutionStatus_;
}

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

double SteadyState::getConvergenceCriterion() const
{
    return convergenceCriterion_;
}

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

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

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

// Dest function definitions

// Static func
void SteadyState::setupMatrix()
{
    setupSSmatrix();
}
void SteadyState::settleFunc()
{
    settle( 0 );
}

void SteadyState::resettleFunc()
{
    settle( 1 );
}

// Dummy function
void SteadyState::assignY( double* S )
{
}

void SteadyState::showMatrices()
{
    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;
}

void SteadyState::setupSSmatrix()
{

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

}

void SteadyState::classifyState( const double* T )
{
    /* 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
    }
}

static bool isSolutionValid( const vector< double >& x )
{
    for( size_t i = 0; i < x.size(); i++ )
    {
        double v = x[i];
        if ( std::isnan( v ) or std::isinf( v ) )
        {
            cout << "Warning: SteadyState iteration gave nan/inf concs\n";
            return false;
        }
        else if( v < 0.0 )
        {
            cout << "Warning: SteadyState iteration gave negative concs\n";
            return false;
        }
    }
    return true;
}

void SteadyState::settle( bool forceSetup )
{

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


/*-----------------------------------------------------------------------------
 *  These functions computes rank of a matrix.
 *-----------------------------------------------------------------------------*/

void swapRows( ublas::matrix< double >& mat, unsigned int r1, unsigned int r2)
{
    ublas::vector<value_type> temp( mat.size2() );
    for (size_t i = 0; i < mat.size2(); i++)
    {
        temp[i] = mat(r1, i );
        mat(r1, i ) = mat(r2, i );
    }

    for (size_t i = 0; i < mat.size2(); i++)
        mat(r2, i) = temp[i];
}


int reorderRows( ublas::matrix< double >& U, int start, int leftCol )
{
    int leftMostRow = start;
    int numReacs = U.size2() - U.size1();
    int newLeftCol = numReacs;
    for ( size_t i = start; i < U.size1(); ++i )
    {
        for ( int j = leftCol; j < numReacs; ++j )
        {
            if ( fabs( U(i,j )) > SteadyState::EPSILON )
            {
                if ( j < newLeftCol )
                {
                    newLeftCol = j;
                    leftMostRow = i;
                }
                break;
            }
        }
    }

    if ( leftMostRow != start )   // swap them.
        swapRows( U, start, leftMostRow );

    return newLeftCol;
}

void eliminateRowsBelow( ublas::matrix< double >& U, int start, int leftCol )
{
    int numMols = U.size1();
    double pivot = U( start, leftCol );
    assert( fabs( pivot ) > SteadyState::EPSILON );
    for ( int i = start + 1; i < numMols; ++i )
    {
        double factor = U(i, leftCol);
        if( fabs ( factor ) > SteadyState::EPSILON )
        {
            factor = factor / pivot;
            for ( size_t j = leftCol + 1; j < U.size2(); ++j )
            {
                double x = U(i,j);
                double y = U( start, j );
                x -= y * factor;
                if ( fabs( x ) < SteadyState::EPSILON )
                    x = 0.0;
                U( i, j ) = x;
            }
        }
        U(i, leftCol) = 0.0;
    }
}

unsigned int rankUsingBoost( ublas::matrix<double>& U )
{
    int numMols = U.size1();
    int numReacs = U.size2() - numMols;
    int i;
    // Start out with a nonzero entry at 0,0
    int leftCol = reorderRows( U, 0, 0 );

    for ( i = 0; i < numMols - 1; ++i )
    {
        eliminateRowsBelow( U, i, leftCol );
        leftCol = reorderRows( U, i + 1, leftCol );
        if ( leftCol == numReacs )
            break;
    }
    return i + 1;
}


static bool checkAboveZero( const vector< double >& y )
{
    for ( vector< double >::const_iterator
            i = y.begin(); i != y.end(); ++i )
    {
        if ( *i < 0.0 )
            return false;
    }
    return true;
}

void recalcTotal( vector< double >& tot, ublas::matrix<double>& g, const double* S )
{
    assert( g.size1() == tot.size() );
    for ( size_t i = 0; i < g.size1(); ++i )
    {
        double t = 0.0;
        for ( unsigned int j = 0; j < g.size2(); ++j )
            t += g( i, j ) * S[j];
        tot[ i ] = t;
    }
}

void SteadyState::randomizeInitialCondition( const Eref& me )
{
    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 );
}

void SteadyState::fitConservationRules(
    ublas::matrix<double>& U, const vector< double >& eliminatedTotal,
    vector< double >&y
)
{
    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;
            }
        }
    }
}