#include <MarkovGslSolver.h>

Collaboration diagram for MarkovGslSolver:

Public Member Functions
double	getAbsoluteAccuracy () const

double	getInternalDt () const

bool	getIsInitialized () const

string	getMethod () const

double	getRelativeAccuracy () const

void	handleQ (vector< vector< double > >)

void	init (vector< double >)

	MarkovGslSolver ()

void	process (const Eref &e, ProcPtr info)

void	reinit (const Eref &e, ProcPtr info)

void	setAbsoluteAccuracy (double value)

void	setInternalDt (double value)

void	setMethod (string method)

void	setRelativeAccuracy (double value)

	~MarkovGslSolver ()

Static Public Member Functions
static const Cinfo *	initCinfo ()

Static Private Member Functions
static int	evalSystem (double, const double , double , void *)

Private Attributes
double	absAccuracy_

gsl_odeiv_control *	gslControl_

gsl_odeiv_evolve *	gslEvolve_

gsl_odeiv_step *	gslStep_

const gsl_odeiv_step_type *	gslStepType_

gsl_odeiv_system	gslSys_

vector< double >	initialState_

double	internalStepSize_

bool	isInitialized_

string	method_

unsigned int	nVars_

vector< vector< double > >	Q_

double	relAccuracy_

vector< double >	state_

double *	stateGsl_

Detailed Description

Definition at line 25 of file MarkovGslSolver.h.

Constructor & Destructor Documentation

MarkovGslSolver::MarkovGslSolver ( )

Definition at line 127 of file MarkovGslSolver.cpp.

References absAccuracy_, gslControl_, gslEvolve_, gslStep_, gslStepType_, internalStepSize_, isInitialized_, method_, nVars_, relAccuracy_, and stateGsl_.

 {
     isInitialized_ = 0;
     method_ = "rk5";
     gslStepType_ = gsl_odeiv_step_rkf45;
     gslStep_ = 0;
     nVars_ = 0;
     absAccuracy_ = 1.0e-8;
     relAccuracy_ = 1.0e-8;
     internalStepSize_ = 1.0e-6;
     stateGsl_ = 0;
     gslEvolve_ = NULL;
     gslControl_ = NULL;
 }

MarkovGslSolver::~MarkovGslSolver ( )

Definition at line 142 of file MarkovGslSolver.cpp.

References gslControl_, gslEvolve_, gslStep_, and stateGsl_.

 {
     if ( gslEvolve_ )
         gsl_odeiv_evolve_free( gslEvolve_ );
     if ( gslControl_ )
         gsl_odeiv_control_free( gslControl_ );
     if ( gslStep_ )
         gsl_odeiv_step_free( gslStep_ );
 
     if ( stateGsl_ )
         delete[] stateGsl_;
 }

Member Function Documentation

int MarkovGslSolver::evalSystem	(	double	t,
		const double *	state,
		double *	f,
		void *	params
	)

staticprivate

Definition at line 155 of file MarkovGslSolver.cpp.

Referenced by init().

 {
     vector< vector< double > >* Q = static_cast< vector< vector< double > >* >( params );
     unsigned int nVars = Q->size();
 
     //Matrix being accessed along columns, which is a very bad thing in terms of
     //cache optimality. Transposing the matrix during reinit() would be a good idea.
     for ( unsigned int i = 0; i < nVars; ++i)
     {
         f[i] = 0;
         for ( unsigned int j = 0; j < nVars; ++j)
             f[i] += state[j] * ((*Q)[j][i]);
     }
 
     return GSL_SUCCESS;
 }

Here is the caller graph for this function:

double MarkovGslSolver::getAbsoluteAccuracy ( ) const

Definition at line 229 of file MarkovGslSolver.cpp.

References absAccuracy_.

Referenced by initCinfo().

 {
     return absAccuracy_;
 }

Here is the caller graph for this function:

double MarkovGslSolver::getInternalDt ( ) const

Definition at line 238 of file MarkovGslSolver.cpp.

References internalStepSize_.

Referenced by initCinfo().

 {
     return internalStepSize_;
 }

Here is the caller graph for this function:

bool MarkovGslSolver::getIsInitialized ( ) const

Definition at line 176 of file MarkovGslSolver.cpp.

References isInitialized_.

Referenced by initCinfo().

 {
     return isInitialized_;
 }

Here is the caller graph for this function:

string MarkovGslSolver::getMethod ( ) const

Definition at line 181 of file MarkovGslSolver.cpp.

References method_.

Referenced by initCinfo().

 {
     return method_;
 }

Here is the caller graph for this function:

double MarkovGslSolver::getRelativeAccuracy ( ) const

Definition at line 219 of file MarkovGslSolver.cpp.

References relAccuracy_.

Referenced by initCinfo().

 {
     return relAccuracy_;
 }

Here is the caller graph for this function:

void MarkovGslSolver::handleQ ( vector< vector< double > > )

Definition at line 347 of file MarkovGslSolver.cpp.

References Q_.

Referenced by initCinfo().

 {
     Q_ = Q;
 }

Here is the caller graph for this function:

void MarkovGslSolver::init ( vector< double > initialState )

Definition at line 253 of file MarkovGslSolver.cpp.

References absAccuracy_, evalSystem(), gslControl_, gslEvolve_, gslStep_, gslStepType_, gslSys_, initialState_, isInitialized_, nVars_, Q_, relAccuracy_, resize(), state_, and stateGsl_.

Referenced by initCinfo().

 {
     nVars_ = initialState.size();
 
     if ( stateGsl_ == 0 )
         stateGsl_ = new double[ nVars_ ];
 
     state_ = initialState;
     initialState_ = initialState;
 
     Q_.resize( nVars_ );
 
     for ( unsigned int i = 0; i < nVars_; ++i )
         Q_[i].resize( nVars_, 0.0 );
 
     isInitialized_ = 1;
 
     assert( gslStepType_ != 0 );
     if ( gslStep_ )
         gsl_odeiv_step_free(gslStep_);
 
     gslStep_ = gsl_odeiv_step_alloc( gslStepType_, nVars_ );
     assert( gslStep_ != 0 );
 
     if ( !gslEvolve_ )
         gslEvolve_ = gsl_odeiv_evolve_alloc(nVars_);
     else
         gsl_odeiv_evolve_reset(gslEvolve_);
 
     assert(gslEvolve_ != 0);
 
     if ( !gslControl_ )
         gslControl_ = gsl_odeiv_control_y_new( absAccuracy_, relAccuracy_ );
     else
         gsl_odeiv_control_init(gslControl_,absAccuracy_, relAccuracy_, 1, 0);
 
     assert(gslControl_!= 0);
 
     gslSys_.function = &MarkovGslSolver::evalSystem;
     gslSys_.jacobian = 0;
     gslSys_.dimension = nVars_;
     gslSys_.params = static_cast< void * >( &Q_ );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

const Cinfo * MarkovGslSolver::initCinfo ( )

static

Definition at line 22 of file MarkovGslSolver.cpp.

References getAbsoluteAccuracy(), getInternalDt(), getIsInitialized(), getMethod(), getRelativeAccuracy(), handleQ(), init(), Neutral::initCinfo(), MarkovGslSolverCinfo, process(), reinit(), setAbsoluteAccuracy(), setInternalDt(), setMethod(), setRelativeAccuracy(), and stateOut().

 {
         // Field definitions
         static ReadOnlyValueFinfo< MarkovGslSolver, bool > isInitialized(
             "isInitialized",
             "True if the message has come in to set solver parameters.",
             &MarkovGslSolver::getIsInitialized
         );
         static ValueFinfo< MarkovGslSolver, string > method( "method",
             "Numerical method to use.",
             &MarkovGslSolver::setMethod,
             &MarkovGslSolver::getMethod
         );
         static ValueFinfo< MarkovGslSolver, double > relativeAccuracy(
             "relativeAccuracy",
             "Accuracy criterion",
             &MarkovGslSolver::setRelativeAccuracy,
             &MarkovGslSolver::getRelativeAccuracy
         );
         static ValueFinfo< MarkovGslSolver, double > absoluteAccuracy(
             "absoluteAccuracy",
             "Another accuracy criterion",
             &MarkovGslSolver::setAbsoluteAccuracy,
             &MarkovGslSolver::getAbsoluteAccuracy
         );
         static ValueFinfo< MarkovGslSolver, double > internalDt(
             "internalDt",
             "internal timestep to use.",
             &MarkovGslSolver::setInternalDt,
             &MarkovGslSolver::getInternalDt
         );
 
         // DestFinfo definitions
         static DestFinfo init( "init",
             "Initialize solver parameters.",
             new OpFunc1< MarkovGslSolver, vector< double > >
             ( &MarkovGslSolver::init )
         );
 
         static DestFinfo handleQ( "handleQ",
             "Handles information regarding the instantaneous rate matrix from "
             "the MarkovRateTable class.",
             new OpFunc1< MarkovGslSolver, vector< vector< double > > >( &MarkovGslSolver::handleQ) );
 
         static DestFinfo process( "process",
             "Handles process call",
             new ProcOpFunc< MarkovGslSolver >( &MarkovGslSolver::process ) );
         static DestFinfo reinit( "reinit",
             "Handles reinit call",
             new ProcOpFunc< MarkovGslSolver >( &MarkovGslSolver::reinit ) );
         // Shared definitions
         static Finfo* procShared[] = {
             &process, &reinit
         };
         static SharedFinfo proc( "proc",
             "Shared message for process and reinit",
             procShared, sizeof( procShared ) / sizeof( const Finfo* )
         );
 
     static Finfo* MarkovGslFinfos[] =
     {
         &isInitialized,         // ValueFinfo
         &method,                        // ValueFinfo
         &relativeAccuracy,  // ValueFinfo
         &absoluteAccuracy,  // ValueFinfo
         &internalDt,                // ValueFinfo
         &init,                          // DestFinfo
         &handleQ,                       // DestFinfo
         &proc,                          // SharedFinfo
         stateOut(),                 // SrcFinfo
     };
 
     static string doc[] =
     {
         "Name", "MarkovGslSolver",
         "Author", "Vishaka Datta S, 2011, NCBS",
         "Description", "Solver for Markov Channel."
     };
 
     static Dinfo< MarkovGslSolver > dinfo;
     static Cinfo MarkovGslSolverCinfo(
         "MarkovGslSolver",
         Neutral::initCinfo(),
         MarkovGslFinfos,
         sizeof(MarkovGslFinfos)/sizeof(Finfo *),
         &dinfo,
         doc,
         sizeof(doc) / sizeof(string)
     );
 
     return &MarkovGslSolverCinfo;
 }

Here is the call graph for this function:

void MarkovGslSolver::process	(	const Eref &	e,
		ProcPtr	info
	)

Definition at line 300 of file MarkovGslSolver.cpp.

References ProcInfo::currTime, ProcInfo::dt, gslControl_, gslEvolve_, gslStep_, gslSys_, internalStepSize_, nVars_, state_, stateGsl_, and stateOut().

Referenced by initCinfo().

 {
     double nextt = info->currTime + info->dt;
     double t = info->currTime;
     double sum = 0;
 
     for ( unsigned int i = 0; i < nVars_; ++i )
         stateGsl_[i] = state_[i];
 
     while ( t < nextt ) {
         int status = gsl_odeiv_evolve_apply (
             gslEvolve_, gslControl_, gslStep_, &gslSys_,
             &t, nextt,
             &internalStepSize_, stateGsl_);
 
         //Simple idea borrowed from Dieter Jaeger's implementation of a Markov
         //channel to deal with potential round-off error.
         sum = 0;
         for ( unsigned int i = 0; i < nVars_; i++ )
             sum += stateGsl_[i];
 
         for ( unsigned int i = 0; i < nVars_; i++ )
             stateGsl_[i] /= sum;
 
         if ( status != GSL_SUCCESS )
             break;
     }
 
     for ( unsigned int i = 0; i < nVars_; ++i )
         state_[i] = stateGsl_[i];
 
     stateOut()->send( e, state_ );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void MarkovGslSolver::reinit	(	const Eref &	e,
		ProcPtr	info
	)

Definition at line 334 of file MarkovGslSolver.cpp.

References initialState_, state_, and stateOut().

Referenced by initCinfo().

 {
     state_ = initialState_;
     if ( initialState_.empty() )
     {
         cerr << "MarkovGslSolver::reinit : "
                  "Initial state has not been set. Solver has not been initialized."
                  "Call init() before running.\n";
     }
 
     stateOut()->send( e, state_ );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void MarkovGslSolver::setAbsoluteAccuracy ( double value )

Definition at line 233 of file MarkovGslSolver.cpp.

References absAccuracy_, and value.

Referenced by initCinfo().

 {
     absAccuracy_ = value;
 }

Here is the caller graph for this function:

void MarkovGslSolver::setInternalDt ( double value )

Definition at line 243 of file MarkovGslSolver.cpp.

References internalStepSize_, and value.

Referenced by initCinfo().

 {
     internalStepSize_ = value;
 }

Here is the caller graph for this function:

void MarkovGslSolver::setMethod ( string method )

Definition at line 186 of file MarkovGslSolver.cpp.

References gslStepType_, and method_.

Referenced by initCinfo().

 {
     method_ = method;
     gslStepType_ = 0;
 
     if ( method == "rk2" ) {
         gslStepType_ = gsl_odeiv_step_rk2;
     } else if ( method == "rk4" ) {
         gslStepType_ = gsl_odeiv_step_rk4;
     } else if ( method == "rk5" ) {
         gslStepType_ = gsl_odeiv_step_rkf45;
     } else if ( method == "rkck" ) {
         gslStepType_ = gsl_odeiv_step_rkck;
     } else if ( method == "rk8pd" ) {
         gslStepType_ = gsl_odeiv_step_rk8pd;
     } else if ( method == "rk2imp" ) {
         gslStepType_ = gsl_odeiv_step_rk2imp;
     } else if ( method == "rk4imp" ) {
         gslStepType_ = gsl_odeiv_step_rk4imp;
     } else if ( method == "bsimp" ) {
         gslStepType_ = gsl_odeiv_step_rk4imp;
         cout << "Warning: implicit Bulirsch-Stoer method not yet implemented: needs Jacobian\n";
     } else if ( method == "gear1" ) {
         gslStepType_ = gsl_odeiv_step_gear1;
     } else if ( method == "gear2" ) {
         gslStepType_ = gsl_odeiv_step_gear2;
     } else {
         cout << "Warning: MarkovGslSolver::innerSetMethod: method '" <<
             method << "' not known, using rk5\n";
         gslStepType_ = gsl_odeiv_step_rkf45;
     }
 }