diff --git a/.travis.yml b/.travis.yml index d4ee1e7fc6..e93cf13311 100644 --- a/.travis.yml +++ b/.travis.yml @@ -32,8 +32,3 @@ before_script: script: - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then ./.ci/travis_build_osx.sh; fi - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then ./.ci/travis_build_linux.sh; fi - -after_success: - - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then\ - mkdir _sdist && cd _sdist && cmake .. && make sdist && make upload_test;\ - fi diff --git a/CMakeLists.txt b/CMakeLists.txt index c61efdfc86..30b6b79ade 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -102,6 +102,14 @@ if(WITH_ASAN AND "${CMAKE_BUILD_TYPE}" STREQUAL "Debug") set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -fsanitize=address") endif() +if(WITH_GSL) + set(WITH_BOOST OFF) + set(WITH_BOOST_ODE OFF) +elseif(WITH_BOOST) + set(WITH_BOOST_ODE ON) + set(WITH_GSL OFF) +endif() + ################################### TARGETS #################################### link_directories(${CMAKE_BINARY_DIR}) @@ -112,17 +120,6 @@ add_executable(moose.bin basecode/main.cpp) ################################### SETUP BUILD ################################ # default include paths. include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ) - -if(WITH_BOOST) - set(WITH_BOOST_ODE ON) - set(WITH_GSL OFF) -endif(WITH_BOOST) - -# If using BOOST ODE2 library to solve ODE system, then don't use GSL. -if(WITH_BOOST_ODE) - set(WITH_GSL OFF) -endif(WITH_BOOST_ODE) - set_target_properties(libmoose PROPERTIES COMPILE_DEFINITIONS "MOOSE_LIB") set_target_properties(libmoose PROPERTIES PREFIX "") @@ -235,6 +232,7 @@ add_subdirectory(benchmarks) add_subdirectory(kinetics) add_subdirectory(synapse) add_subdirectory(intfire) +add_subdirectory(external/libsoda) ###################################### LINKING ################################# list(APPEND MOOSE_LIBRARIES @@ -255,6 +253,7 @@ list(APPEND MOOSE_LIBRARIES signeur diffusion ksolve + lsoda device basecode ) diff --git a/basecode/global.cpp b/basecode/global.cpp index c2c365428e..c5035f45db 100644 --- a/basecode/global.cpp +++ b/basecode/global.cpp @@ -144,6 +144,22 @@ namespace moose { if( p.size() == 0 ) return true; +#ifdef USE_BOOST_FILESYSTEM + try + { + boost::filesystem::path pdirs( p ); + boost::filesystem::create_directories( pdirs ); + LOG( moose::info, "Created directory " << p ); + return true; + } + catch(const boost::filesystem::filesystem_error& e) + { + LOG( moose::warning, "create_directories(" << p << ") failed with " + << e.code().message() + ); + return false; + } +#else /* ----- not USE_BOOST_FILESYSTEM ----- */ string command( "mkdir -p "); command += p; int ret = system( command.c_str() ); @@ -163,6 +179,7 @@ namespace moose { LOG( moose::warning, p << " is no directory" ); return false; } +#endif /* ----- not USE_BOOST_FILESYSTEM ----- */ return true; } diff --git a/basecode/global.h b/basecode/global.h index 3beab75dc5..76ca1cdd75 100644 --- a/basecode/global.h +++ b/basecode/global.h @@ -17,7 +17,12 @@ #include #include -#include "../randnum/RNG.h" /* Use inbuilt rng */ + +#ifdef USE_BOOST_FILESYSTEM +#include +#endif + +#include "randnum/RNG.h" /* Use inbuilt rng */ #include "../utility/print_function.hpp" using namespace std; diff --git a/basecode/main.cpp b/basecode/main.cpp index dc194c0801..ec4a7c6d87 100644 --- a/basecode/main.cpp +++ b/basecode/main.cpp @@ -35,7 +35,6 @@ int testIndex = 0; f; \ cout << std::right << " [DONE]" << endl; \ -extern void testSync(); extern void testAsync(); extern void testSyncArray( unsigned int size, unsigned int method ); extern void testShell(); @@ -67,9 +66,9 @@ extern void testSigNeurProcess(); extern unsigned int initMsgManagers(); extern void destroyMsgManagers(); // void regressionTests(); -#endif -extern void speedTestMultiNodeIntFireNetwork( - unsigned int size, unsigned int runsteps ); +#endif // DO_UNIT_TESTS + +extern void speedTestMultiNodeIntFireNetwork(unsigned int size, unsigned int runsteps); #ifdef USE_SMOLDYN extern void testSmoldyn(); diff --git a/biophysics/CMakeLists.txt b/biophysics/CMakeLists.txt index 390117a92e..0c02172b7a 100644 --- a/biophysics/CMakeLists.txt +++ b/biophysics/CMakeLists.txt @@ -2,7 +2,11 @@ cmake_minimum_required(VERSION 2.8) include(${CMAKE_CURRENT_SOURCE_DIR}/../CheckCXXCompiler.cmake) if(WITH_GSL) + find_package(GSL 1.16) include_directories(${GSL_INCLUDE_DIRS}) +elseif(WITH_BOOST_ODE) + find_package(Boost) + include_directories(${Boost_INCLUDE_DIRS}) endif(WITH_GSL) set(BIOPHYSICS_SRCS @@ -32,7 +36,6 @@ set(BIOPHYSICS_SRCS ReadSwc.cpp SynChan.cpp NMDAChan.cpp - testBiophysics.cpp IzhikevichNrn.cpp DifShellBase.cpp DifShell.cpp @@ -48,10 +51,9 @@ set(BIOPHYSICS_SRCS MarkovSolver.cpp VClamp.cpp Spine.cpp + MarkovOdeSolver.cpp + testBiophysics.cpp ) -if(WITH_GSL) - list(APPEND BIOPHYSICS_SRCS MarkovGslSolver.cpp) -endif(WITH_GSL) +add_library(biophysics ${BIOPHYSICS_SRCS}) -add_library( biophysics ${BIOPHYSICS_SRCS} ) diff --git a/biophysics/MarkovChannel.cpp b/biophysics/MarkovChannel.cpp index a79b17b49b..bca0cc2541 100644 --- a/biophysics/MarkovChannel.cpp +++ b/biophysics/MarkovChannel.cpp @@ -21,207 +21,207 @@ const Cinfo* MarkovChannel::initCinfo() { - /////////////////////// - //Field information. - /////////////////////// - static ValueFinfo< MarkovChannel, double > ligandconc( "ligandConc", - "Ligand concentration.", - &MarkovChannel::setLigandConc, - &MarkovChannel::getLigandConc - ); - - static ValueFinfo< MarkovChannel, double > vm( "Vm", - "Membrane voltage.", - &MarkovChannel::setVm, - &MarkovChannel::getVm - ); - - static ValueFinfo< MarkovChannel, unsigned int > numstates( "numStates", - "The number of states that the channel can occupy.", - &MarkovChannel::setNumStates, - &MarkovChannel::getNumStates - ); - - - static ValueFinfo< MarkovChannel, unsigned int > numopenstates( "numOpenStates", - "The number of states which are open/conducting.", - &MarkovChannel::setNumOpenStates, - &MarkovChannel::getNumOpenStates - ); - - static ValueFinfo< MarkovChannel, vector< string > > labels("labels", - "Labels for each state.", - &MarkovChannel::setStateLabels, - &MarkovChannel::getStateLabels - ); - - static ReadOnlyValueFinfo< MarkovChannel, vector< double > > state( "state", - "This is a row vector that contains the probabilities of finding the channel in each state.", - &MarkovChannel::getState - ); - - static ValueFinfo< MarkovChannel, vector< double > > initialstate( "initialState", - "This is a row vector that contains the probabilities of finding the channel in each state at t = 0. The state of the channel is reset to this value during a call to reinit()", - &MarkovChannel::setInitialState, - &MarkovChannel::getInitialState - ); - - static ValueFinfo< MarkovChannel, vector< double > > gbar( "gbar", - "A row vector containing the conductance associated with each of the open/conducting states.", - &MarkovChannel::setGbars, - &MarkovChannel::getGbars - ); - - //MsgDest functions - static DestFinfo handleligandconc( "handleLigandConc", - "Deals with incoming messages containing information of ligand concentration", - new OpFunc1< MarkovChannel, double >(&MarkovChannel::handleLigandConc) ); - - static DestFinfo handlestate("handleState", - "Deals with incoming message from MarkovSolver object containing state information of the channel.\n", - new OpFunc1< MarkovChannel, vector< double > >(&MarkovChannel::handleState) ); - - /////////////////////////////////////////// - static Finfo* MarkovChannelFinfos[] = - { - &ligandconc, - &vm, - &numstates, - &numopenstates, - &state, - &initialstate, - &labels, - &gbar, - &handleligandconc, - &handlestate, - }; - - static string doc[] = - { - "Name", "MarkovChannel", - "Author", "Vishaka Datta S, 2011, NCBS", - "Description", "MarkovChannel : Multistate ion channel class." - "It deals with ion channels which can be found in one of multiple states, " - "some of which are conducting. This implementation assumes the occurence " - "of first order kinetics to calculate the probabilities of the channel " - "being found in all states. Further, the rates of transition between these " - "states can be constant, voltage-dependent or ligand dependent (only one " - "ligand species). The current flow obtained from the channel is calculated " - "in a deterministic method by solving the system of differential equations " - "obtained from the assumptions above." - }; - - static Dinfo< MarkovChannel > dinfo; - static Cinfo MarkovChannelCinfo( - "MarkovChannel", - ChanBase::initCinfo(), - MarkovChannelFinfos, - sizeof( MarkovChannelFinfos )/ sizeof( Finfo* ), - &dinfo, + /////////////////////// + //Field information. + /////////////////////// + static ValueFinfo< MarkovChannel, double > ligandconc( "ligandConc", + "Ligand concentration.", + &MarkovChannel::setLigandConc, + &MarkovChannel::getLigandConc + ); + + static ValueFinfo< MarkovChannel, double > vm( "Vm", + "Membrane voltage.", + &MarkovChannel::setVm, + &MarkovChannel::getVm + ); + + static ValueFinfo< MarkovChannel, unsigned int > numstates( "numStates", + "The number of states that the channel can occupy.", + &MarkovChannel::setNumStates, + &MarkovChannel::getNumStates + ); + + + static ValueFinfo< MarkovChannel, unsigned int > numopenstates( "numOpenStates", + "The number of states which are open/conducting.", + &MarkovChannel::setNumOpenStates, + &MarkovChannel::getNumOpenStates + ); + + static ValueFinfo< MarkovChannel, vector< string > > labels("labels", + "Labels for each state.", + &MarkovChannel::setStateLabels, + &MarkovChannel::getStateLabels + ); + + static ReadOnlyValueFinfo< MarkovChannel, vector< double > > state( "state", + "This is a row vector that contains the probabilities of finding the channel in each state.", + &MarkovChannel::getState + ); + + static ValueFinfo< MarkovChannel, vector< double > > initialstate( "initialState", + "This is a row vector that contains the probabilities of finding the channel in each state at t = 0. The state of the channel is reset to this value during a call to reinit()", + &MarkovChannel::setInitialState, + &MarkovChannel::getInitialState + ); + + static ValueFinfo< MarkovChannel, vector< double > > gbar( "gbar", + "A row vector containing the conductance associated with each of the open/conducting states.", + &MarkovChannel::setGbars, + &MarkovChannel::getGbars + ); + + //MsgDest functions + static DestFinfo handleligandconc( "handleLigandConc", + "Deals with incoming messages containing information of ligand concentration", + new OpFunc1< MarkovChannel, double >(&MarkovChannel::handleLigandConc) ); + + static DestFinfo handlestate("handleState", + "Deals with incoming message from MarkovSolver object containing state information of the channel.\n", + new OpFunc1< MarkovChannel, vector< double > >(&MarkovChannel::handleState) ); + + /////////////////////////////////////////// + static Finfo* MarkovChannelFinfos[] = + { + &ligandconc, + &vm, + &numstates, + &numopenstates, + &state, + &initialstate, + &labels, + &gbar, + &handleligandconc, + &handlestate, + }; + + static string doc[] = + { + "Name", "MarkovChannel", + "Author", "Vishaka Datta S, 2011, NCBS Bangalore", + "Description", "MarkovChannel.\n" + "It deals with ion channels which can be found in one of multiple states, " + "some of which are conducting. This implementation assumes the occurence " + "of first order kinetics to calculate the probabilities of the channel " + "being found in all states. Further, the rates of transition between these " + "states can be constant, voltage-dependent or ligand dependent (only one " + "ligand species). The current flow obtained from the channel is calculated " + "in a deterministic method by solving the system of differential equations " + "obtained from the assumptions above." + }; + + static Dinfo< MarkovChannel > dinfo; + static Cinfo MarkovChannelCinfo( + "MarkovChannel", + ChanBase::initCinfo(), + MarkovChannelFinfos, + sizeof( MarkovChannelFinfos )/ sizeof( Finfo* ), + &dinfo, doc, sizeof(doc) / sizeof(string) - ); + ); - return &MarkovChannelCinfo; + return &MarkovChannelCinfo; } static const Cinfo* markovChannelCinfo = MarkovChannel::initCinfo(); MarkovChannel::MarkovChannel() : - g_(0), - ligandConc_(0), - numStates_(0), - numOpenStates_(0) + g_(0), + ligandConc_(0), + numStates_(0), + numOpenStates_(0) { ; } MarkovChannel::MarkovChannel(unsigned int numStates, unsigned int numOpenStates) : - g_(0), ligandConc_(0), numStates_(numStates), numOpenStates_(numOpenStates) + g_(0), ligandConc_(0), numStates_(numStates), numOpenStates_(numOpenStates) { - stateLabels_.resize( numStates ); - state_.resize( numStates ); - initialState_.resize( numStates ); - Gbars_.resize( numOpenStates ) ; + stateLabels_.resize( numStates ); + state_.resize( numStates ); + initialState_.resize( numStates ); + Gbars_.resize( numOpenStates ) ; } MarkovChannel::~MarkovChannel( ) { - ; + ; } double MarkovChannel::getVm( ) const { - return Vm_; + return Vm_; } void MarkovChannel::setVm( double Vm ) { - Vm_ = Vm; + Vm_ = Vm; } double MarkovChannel::getLigandConc( ) const { - return ligandConc_; + return ligandConc_; } void MarkovChannel::setLigandConc( double ligandConc ) { - ligandConc_ = ligandConc; + ligandConc_ = ligandConc; } unsigned int MarkovChannel::getNumStates( ) const { - return numStates_; + return numStates_; } void MarkovChannel::setNumStates( unsigned int numStates ) { - numStates_ = numStates; + numStates_ = numStates; } unsigned int MarkovChannel::getNumOpenStates( ) const { - return numOpenStates_; + return numOpenStates_; } void MarkovChannel::setNumOpenStates( unsigned int numOpenStates ) { - numOpenStates_ = numOpenStates; + numOpenStates_ = numOpenStates; } vector< string > MarkovChannel::getStateLabels( ) const { - return stateLabels_; + return stateLabels_; } void MarkovChannel::setStateLabels( vector< string > stateLabels ) { - stateLabels_ = stateLabels; + stateLabels_ = stateLabels; } vector< double > MarkovChannel::getState ( ) const { - return state_; + return state_; } vector< double > MarkovChannel::getInitialState() const { - return initialState_; + return initialState_; } void MarkovChannel::setInitialState( vector< double > initialState ) { - initialState_ = initialState; - state_ = initialState; + initialState_ = initialState; + state_ = initialState; } vector< double > MarkovChannel::getGbars() const { - return Gbars_; + return Gbars_; } void MarkovChannel::setGbars( vector< double > Gbars ) { - Gbars_ = Gbars; + Gbars_ = Gbars; } ///////////////////////////// @@ -230,41 +230,41 @@ void MarkovChannel::setGbars( vector< double > Gbars ) void MarkovChannel::vProcess( const Eref& e, const ProcPtr p ) { - g_ = 0.0; + g_ = 0.0; - //Cannot use the Gbar_ variable of the ChanBase class. The conductance - //Gk_ calculated here is the "expected conductance" of the channel due to its - //stochastic nature. + //Cannot use the Gbar_ variable of the ChanBase class. The conductance + //Gk_ calculated here is the "expected conductance" of the channel due to its + //stochastic nature. - for( unsigned int i = 0; i < numOpenStates_; ++i ) - g_ += Gbars_[i] * state_[i]; + for( unsigned int i = 0; i < numOpenStates_; ++i ) + g_ += Gbars_[i] * state_[i]; - setGk( e, g_ ); - updateIk(); + setGk( e, g_ ); + updateIk(); - sendProcessMsgs( e, p ); + sendProcessMsgs( e, p ); } void MarkovChannel::vReinit( const Eref& e, const ProcPtr p ) { - g_ = 0.0; + g_ = 0.0; - if ( initialState_.empty() ) - { - cerr << "MarkovChannel::reinit : Initial state has not been set.!\n"; - return; - } - state_ = initialState_; + if ( initialState_.empty() ) + { + cerr << "MarkovChannel::reinit : Initial state has not been set.!\n"; + return; + } + state_ = initialState_; - sendReinitMsgs( e, p ); + sendReinitMsgs( e, p ); } void MarkovChannel::handleLigandConc( double ligandConc ) { - ligandConc_ = ligandConc; + ligandConc_ = ligandConc; } void MarkovChannel::handleState( vector< double > state ) { - state_ = state; + state_ = state; } diff --git a/biophysics/MarkovGslSolver.cpp b/biophysics/MarkovGslSolver.cpp deleted file mode 100644 index ca90e7a3a3..0000000000 --- a/biophysics/MarkovGslSolver.cpp +++ /dev/null @@ -1,350 +0,0 @@ -/********************************************************************** -** This program is part of 'MOOSE', the -** Messaging Object Oriented Simulation Environment. -** Copyright (C) 2003-2011 Upinder S. Bhalla. and NCBS -** It is made available under the terms of the -** GNU Lesser General Public License version 2.1 -** See the file COPYING.LIB for the full notice. -**********************************************************************/ - -#include "../basecode/header.h" -#include -#include -#include "MarkovGslSolver.h" - -static SrcFinfo1< vector >* stateOut() -{ - static SrcFinfo1< vector< double > > stateOut( "stateOut", - "Sends updated state to the MarkovChannel class." ); - return &stateOut; -} - -const Cinfo* MarkovGslSolver::initCinfo() -{ - /////////////////////////////////////////////////////// - // 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; -} - -static const Cinfo* MarkovGslSolverCinfo = MarkovGslSolver::initCinfo(); - -/////////////////////////////////////////////////// -// Basic class function definitions -/////////////////////////////////////////////////// - -MarkovGslSolver::MarkovGslSolver() -{ - 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() -{ - 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_; -} - -int MarkovGslSolver::evalSystem( double t, const double* state, double* f, void *params) -{ - 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; -} - -/////////////////////////////////////////////////// -// Field function definitions -/////////////////////////////////////////////////// - -bool MarkovGslSolver::getIsInitialized() const -{ - return isInitialized_; -} - -string MarkovGslSolver::getMethod() const -{ - return method_; -} - -void MarkovGslSolver::setMethod( string method ) -{ - 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; - } -} - -double MarkovGslSolver::getRelativeAccuracy() const -{ - return relAccuracy_; -} - -void MarkovGslSolver::setRelativeAccuracy( double value ) -{ - relAccuracy_ = value; -} - -double MarkovGslSolver::getAbsoluteAccuracy() const -{ - return absAccuracy_; -} -void MarkovGslSolver::setAbsoluteAccuracy( double value ) -{ - absAccuracy_ = value; -} - -double MarkovGslSolver::getInternalDt() const -{ - return internalStepSize_; -} - -void MarkovGslSolver::setInternalDt( double value ) -{ - internalStepSize_ = value; -} - -/////////////////////////////////////////////////// -// Dest function definitions -/////////////////////////////////////////////////// - -//Handles data from MarkovChannel class. -void MarkovGslSolver::init( vector< double > initialState ) -{ - 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_ ); -} - -////////////////////////// -//MsgDest functions. -///////////////////////// -void MarkovGslSolver::process( const Eref& e, ProcPtr info ) -{ - 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_ ); -} - -void MarkovGslSolver::reinit( const Eref& e, ProcPtr info ) -{ - 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_ ); -} - -void MarkovGslSolver::handleQ( vector< vector< double > > Q ) -{ - Q_ = Q; -} diff --git a/biophysics/MarkovGslSolver.h b/biophysics/MarkovGslSolver.h deleted file mode 100644 index a1ae31f125..0000000000 --- a/biophysics/MarkovGslSolver.h +++ /dev/null @@ -1,78 +0,0 @@ -/********************************************************************** -** This program is part of 'MOOSE', the -** Messaging Object Oriented Simulation Environment. -** Copyright (C) 2003-2011 Upinder S. Bhalla. and NCBS -** It is made available under the terms of the -** GNU Lesser General Public License version 2.1 -** See the file COPYING.LIB for the full notice. -**********************************************************************/ - -#ifndef _MARKOVGSLSOLVER_H -#define _MARKOVGSLSOLVER_H - -//////////////////////////////////////////////////// -// Class : MarkovGslSolver -// Author : Vishaka Datta S, 2011, NCBS. -// -// The GslIntegrator class in ksolve deals with a system whose coefficients stay -// constant. In the case of a Markov channel, the coefficients of the system -// vary with time. -// -// This makes it necessary for the system to keep track of changes in the system -// matrix, which is implemented by the message handler. -/////////////////////////////////////////////////// - -class MarkovGslSolver -{ - public: - MarkovGslSolver(); - ~MarkovGslSolver(); - -/////////////////////////////////////////////////// -// Field function definitions -/////////////////////////////////////////////////// - bool getIsInitialized() const; - string getMethod() const; - void setMethod( string method ); - double getRelativeAccuracy() const; - void setRelativeAccuracy( double value ); - double getAbsoluteAccuracy() const; - void setAbsoluteAccuracy( double value ); - double getInternalDt() const; - void setInternalDt( double value ); - -/////////////////////////////////////////////////// -// Dest function definitions -/////////////////////////////////////////////////// - - void process( const Eref& e, ProcPtr info ); - void reinit( const Eref& e, ProcPtr info ); - - void init( vector< double > ); - void handleQ( vector< vector< double > > ); - - static const Cinfo* initCinfo(); - private: - bool isInitialized_; - string method_; - double absAccuracy_; - double relAccuracy_; - double internalStepSize_; - double* stateGsl_; - - //The following four variables should be members of any solver class that - //will be implmented. - unsigned int nVars_; - vector< double > state_; - vector< double > initialState_; - vector< vector< double > > Q_; - - const gsl_odeiv_step_type* gslStepType_; - gsl_odeiv_step* gslStep_; - gsl_odeiv_control* gslControl_; - gsl_odeiv_evolve* gslEvolve_; - gsl_odeiv_system gslSys_; - - static int evalSystem( double, const double*, double*, void* ); -}; -#endif diff --git a/biophysics/MarkovOdeSolver.cpp b/biophysics/MarkovOdeSolver.cpp new file mode 100644 index 0000000000..9e951e7475 --- /dev/null +++ b/biophysics/MarkovOdeSolver.cpp @@ -0,0 +1,439 @@ +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment. +** Copyright (C) 2003-2011 Upinder S. Bhalla. and NCBS +** It is made available under the terms of the +** GNU Lesser General Public License version 2.1 +** See the file COPYING.LIB for the full notice. +**********************************************************************/ + +#include "../basecode/header.h" +#if USE_BOOSE_ODE +#include "../utility/boost_ode.h" +#endif +#include "MarkovOdeSolver.h" + +#ifdef USE_GSL +#include +#include +#endif /* ----- not USE_GSL ----- */ + +#ifdef USE_BOOST_ODE +#include +using namespace boost::numeric; +#endif /* ----- not USE_BOOST_ODE ----- */ + +static SrcFinfo1< vector >* stateOut() +{ + static SrcFinfo1< vector< double > > stateOut( "stateOut", + "Sends updated state to the MarkovChannel class." ); + return &stateOut; +} + +const Cinfo* MarkovOdeSolver::initCinfo() +{ + /////////////////////////////////////////////////////// + // Field definitions + /////////////////////////////////////////////////////// + static ReadOnlyValueFinfo< MarkovOdeSolver, bool > isInitialized( + "isInitialized", + "True if the message has come in to set solver parameters.", + &MarkovOdeSolver::getIsInitialized + ); + static ValueFinfo< MarkovOdeSolver, string > method( "method", + "Numerical method to use.", + &MarkovOdeSolver::setMethod, + &MarkovOdeSolver::getMethod + ); + static ValueFinfo< MarkovOdeSolver, double > relativeAccuracy( + "relativeAccuracy", + "Accuracy criterion", + &MarkovOdeSolver::setRelativeAccuracy, + &MarkovOdeSolver::getRelativeAccuracy + ); + static ValueFinfo< MarkovOdeSolver, double > absoluteAccuracy( + "absoluteAccuracy", + "Another accuracy criterion", + &MarkovOdeSolver::setAbsoluteAccuracy, + &MarkovOdeSolver::getAbsoluteAccuracy + ); + static ValueFinfo< MarkovOdeSolver, double > internalDt( + "internalDt", + "internal timestep to use.", + &MarkovOdeSolver::setInternalDt, + &MarkovOdeSolver::getInternalDt + ); + + /////////////////////////////////////////////////////// + // DestFinfo definitions + /////////////////////////////////////////////////////// + static DestFinfo init( "init", + "Initialize solver parameters.", + new OpFunc1< MarkovOdeSolver, vector< double > > + ( &MarkovOdeSolver::init ) + ); + + static DestFinfo handleQ( "handleQ", + "Handles information regarding the instantaneous rate matrix from " + "the MarkovRateTable class.", + new OpFunc1< MarkovOdeSolver, vector< vector< double > > >( &MarkovOdeSolver::handleQ) ); + + static DestFinfo process( "process", + "Handles process call", + new ProcOpFunc< MarkovOdeSolver >( &MarkovOdeSolver::process ) ); + static DestFinfo reinit( "reinit", + "Handles reinit call", + new ProcOpFunc< MarkovOdeSolver >( &MarkovOdeSolver::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* MarkovOdeFinfos[] = + { + &isInitialized, // ValueFinfo + &method, // ValueFinfo + &relativeAccuracy, // ValueFinfo + &absoluteAccuracy, // ValueFinfo + &internalDt, // ValueFinfo + &init, // DestFinfo + &handleQ, // DestFinfo + &proc, // SharedFinfo + stateOut(), // SrcFinfo + }; + + static string doc[] = + { + "Name", "MarkovOdeSolver", + "Author", "Vishaka Datta (c) 2011, Dilawar Singh (c) 2018", + "Description", "Solver for Markov Channel." + }; + + static Dinfo< MarkovOdeSolver > dinfo; + static Cinfo MarkovOdeSolverCinfo( + "MarkovOdeSolver", + Neutral::initCinfo(), + MarkovOdeFinfos, + sizeof(MarkovOdeFinfos)/sizeof(Finfo *), + &dinfo, + doc, + sizeof(doc) / sizeof(string) + ); + + return &MarkovOdeSolverCinfo; +} + +static const Cinfo* MarkovOdeSolverCinfo = MarkovOdeSolver::initCinfo(); + +/////////////////////////////////////////////////// +// Basic class function definitions +/////////////////////////////////////////////////// + +MarkovOdeSolver::MarkovOdeSolver() +{ + isInitialized_ = 0; + method_ = "rk5"; +#ifdef USE_GSL + gslStepType_ = gsl_odeiv_step_rkf45; + gslEvolve_ = NULL; + gslControl_ = NULL; + gslStep_ = 0; +#endif /* ----- not USE_GSL ----- */ + nVars_ = 0; + absAccuracy_ = 1e-8; + relAccuracy_ = 1e-8; + internalStepSize_ = 1.0e-6; +} + +MarkovOdeSolver::~MarkovOdeSolver() +{ +#ifdef USE_GSL + if ( gslEvolve_ ) + gsl_odeiv_evolve_free( gslEvolve_ ); + if ( gslControl_ ) + gsl_odeiv_control_free( gslControl_ ); + if ( gslStep_ ) + gsl_odeiv_step_free( gslStep_ ); +#endif /* ----- not USE_GSL ----- */ +} + +#ifdef USE_GSL +int MarkovOdeSolver::evalSystem( double t, const double* state, double* f, void *params) +{ + 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 0; +} +#endif /* ----- not USE_GSL ----- */ + +#ifdef USE_BOOST_ODE +/* --------------------------------------------------------------------------*/ +/** + * @Synopsis OdeSystem. Similar to MarkovOdeSolver::evalSystem + * + * @Param y + * @Param dydt + */ +/* ----------------------------------------------------------------------------*/ +void MarkovOdeSolver::OdeSystem( const vector& y, vector& dydt ) +{ + 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) + { + dydt[i] = 0; + for ( unsigned int j = 0; j < nVars; ++j) + dydt[i] += y[j] * Q_[j][i]; + } +} +#endif + +/////////////////////////////////////////////////// +// Field function definitions +/////////////////////////////////////////////////// + +bool MarkovOdeSolver::getIsInitialized() const +{ + return isInitialized_; +} + +string MarkovOdeSolver::getMethod() const +{ + return method_; +} + +void MarkovOdeSolver::setMethod( string method ) +{ + method_ = method; +#ifdef USE_GSL + 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: MarkovOdeSolver::innerSetMethod: method '" << + method << "' not known, using rk5\n"; + gslStepType_ = gsl_odeiv_step_rkf45; + } +#endif +} + +double MarkovOdeSolver::getRelativeAccuracy() const +{ + return relAccuracy_; +} + +void MarkovOdeSolver::setRelativeAccuracy( double value ) +{ + relAccuracy_ = value; +} + +double MarkovOdeSolver::getAbsoluteAccuracy() const +{ + return absAccuracy_; +} +void MarkovOdeSolver::setAbsoluteAccuracy( double value ) +{ + absAccuracy_ = value; +} + +double MarkovOdeSolver::getInternalDt() const +{ + return internalStepSize_; +} + +void MarkovOdeSolver::setInternalDt( double value ) +{ + internalStepSize_ = value; +} + +/////////////////////////////////////////////////// +// Dest function definitions +/////////////////////////////////////////////////// + +//Handles data from MarkovChannel class. +void MarkovOdeSolver::init( vector< double > initialState ) +{ + + nVars_ = initialState.size(); + state_ = initialState; + initialState_ = initialState; + Q_.resize( nVars_ ); + for ( unsigned int i = 0; i < nVars_; ++i ) + Q_[i].resize( nVars_, 0.0 ); + + stateOde_.resize( nVars_ ); + isInitialized_ = 1; + +#ifdef USE_GSL + 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 = &MarkovOdeSolver::evalSystem; + gslSys_.jacobian = 0; + gslSys_.dimension = nVars_; + gslSys_.params = static_cast< void * >( &Q_ ); +#endif /* ----- not USE_GSL ----- */ + +} + +////////////////////////// +//MsgDest functions. +///////////////////////// +void MarkovOdeSolver::process( const Eref& e, ProcPtr p ) +{ + double nextt = p->currTime + p->dt; + double t = p->currTime; + double sum = 0; + + for ( unsigned int i = 0; i < nVars_; ++i ) + stateOde_[i] = state_[i]; + + +#ifdef USE_GSL + while( t < nextt ) + { + int status = gsl_odeiv_evolve_apply ( gslEvolve_, gslControl_, gslStep_, &gslSys_, + &t, nextt, &internalStepSize_, &stateOde_[0] + ); + //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 += stateOde_[i]; + + for ( unsigned int i = 0; i < nVars_; i++ ) + stateOde_[i] /= sum; + + if( status != GSL_SUCCESS ) + break; + } +#endif + +#if USE_BOOST_ODE + auto sys = [this](const vector& dy, vector& dydt, const double t) + { + this->OdeSystem(dy, dydt); + }; + + // It works well in practice for this setup. + if( method_ == "rk5" || method_ == "gsl" || method_ == "boost" ) + odeint::integrate( sys, stateOde_, t, nextt, p->dt ); + else if( method_ == "rk5a" || method_ == "adaptive" ) + odeint::integrate_adaptive( + odeint::make_controlled( absAccuracy_, relAccuracy_ ) + , sys, stateOde_, t, nextt, p->dt + ); + else if( method_ == "rk4" ) + odeint::integrate_const( rk4_stepper_type_() + , sys, stateOde_, t, nextt, p->dt + ); + else if ("rk54" == method_ ) + odeint::integrate_const( rk_karp_stepper_type_() + , sys, stateOde_, t, nextt, p->dt + ); + else if ("rkck" == method_ || "rkcka" == method_) + odeint::integrate_adaptive( + odeint::make_controlled( absAccuracy_, relAccuracy_ ) + , sys, stateOde_, t, nextt, p->dt + ); + else if( method_ == "rk8a" || "rk8" == method_ ) + odeint::integrate_adaptive( rk_felhberg_stepper_type_() + , sys, stateOde_, t, nextt, p->dt + ); + else + odeint::integrate( sys, stateOde_, t, nextt, p->dt ); +#endif /* ----- not USE_GSL ----- */ + + + for ( unsigned int i = 0; i < nVars_; ++i ) + state_[i] = stateOde_[i]; + + stateOut()->send( e, state_ ); +} + +void MarkovOdeSolver::reinit( const Eref& e, ProcPtr info ) +{ + state_ = initialState_; + if ( initialState_.empty() ) + { + cerr << "MarkovOdeSolver::reinit : " + "Initial state has not been set. Solver has not been initialized." + "Call init() before running.\n"; + } + stateOut()->send( e, state_ ); +} + +void MarkovOdeSolver::handleQ( vector< vector< double > > Q ) +{ + Q_ = Q; +} diff --git a/biophysics/MarkovOdeSolver.h b/biophysics/MarkovOdeSolver.h new file mode 100644 index 0000000000..a12c1ea16f --- /dev/null +++ b/biophysics/MarkovOdeSolver.h @@ -0,0 +1,88 @@ +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment. +** Copyright (C) 2003-2011 Upinder S. Bhalla. and NCBS +** It is made available under the terms of the +** GNU Lesser General Public License version 2.1 +** See the file COPYING.LIB for the full notice. +**********************************************************************/ + +#ifndef _MARKOVGSLSOLVER_H +#define _MARKOVGSLSOLVER_H + +//////////////////////////////////////////////////// +// Class : MarkovOdeSolver +// Author : Vishaka Datta S, 2011, NCBS. +// +// The GslIntegrator class in ksolve deals with a system whose coefficients stay +// constant. In the case of a Markov channel, the coefficients of the system +// vary with time. +// +// This makes it necessary for the system to keep track of changes in the system +// matrix, which is implemented by the message handler. +/////////////////////////////////////////////////// +// +#ifdef USE_GSL +#include +#include +#endif /* ----- not USE_GSL ----- */ + +class MarkovOdeSolver +{ +public: + MarkovOdeSolver(); + ~MarkovOdeSolver(); + +/////////////////////////////////////////////////// +// Field function definitions +/////////////////////////////////////////////////// + bool getIsInitialized() const; + string getMethod() const; + void setMethod( string method ); + double getRelativeAccuracy() const; + void setRelativeAccuracy( double value ); + double getAbsoluteAccuracy() const; + void setAbsoluteAccuracy( double value ); + double getInternalDt() const; + void setInternalDt( double value ); + +/////////////////////////////////////////////////// +// Dest function definitions +/////////////////////////////////////////////////// + + void process( const Eref& e, ProcPtr info ); + void reinit( const Eref& e, ProcPtr info ); + + void init( vector< double > ); + void handleQ( vector< vector< double > > ); + + static const Cinfo* initCinfo(); + +private: + bool isInitialized_; + string method_; + double absAccuracy_; + double relAccuracy_; + double internalStepSize_; + vector stateOde_; + + //The following four variables should be members of any solver class that + //will be implmented. + unsigned int nVars_; + vector< double > state_; + vector< double > initialState_; + vector< vector< double > > Q_; + +#ifdef USE_GSL + const gsl_odeiv_step_type* gslStepType_; + gsl_odeiv_step* gslStep_; + gsl_odeiv_control* gslControl_; + gsl_odeiv_evolve* gslEvolve_; + gsl_odeiv_system gslSys_; + static int evalSystem( double, const double*, double*, void* ); +#else /* ----- not USE_GSL ----- */ + // Using Boost + void OdeSystem( const vector& y, vector& dydt ); +#endif /* ----- not USE_GSL ----- */ +}; +#endif diff --git a/biophysics/MarkovSolver.cpp b/biophysics/MarkovSolver.cpp index 3f1a0e006e..2d368010ea 100644 --- a/biophysics/MarkovSolver.cpp +++ b/biophysics/MarkovSolver.cpp @@ -7,7 +7,6 @@ ** See the file COPYING.LIB for the full notice. **********************************************************************/ -#include //Needed for DBL_MAX and DBL_MIN #include "../basecode/header.h" #include "MatrixOps.h" @@ -20,264 +19,264 @@ const Cinfo* MarkovSolver::initCinfo() { - ////////////////////// - //DestFinfos - ////////////////////// - - static DestFinfo process( "process", - "Handles process call", - new ProcOpFunc< MarkovSolver >( &MarkovSolver::process ) ); - - static DestFinfo reinit( "reinit", - "Handles reinit call", - new ProcOpFunc< MarkovSolver >( &MarkovSolver::reinit ) ); - - static Finfo* processShared[] = - { - &process, &reinit - }; - - static SharedFinfo proc( "proc", - "This is a shared message to receive Process message from the" - "scheduler. The first entry is a MsgDest for the Process " - "operation. It has a single argument, ProcInfo, which " - "holds lots of information about current time, thread, dt and" - "so on. The second entry is a MsgDest for the Reinit " - "operation. It also uses ProcInfo.", - processShared, sizeof( processShared ) / sizeof( Finfo* ) - ); - - - static Finfo* markovSolverFinfos[] = - { - &proc, //SharedFinfo - }; - - static Dinfo < MarkovSolver > dinfo; - static Cinfo markovSolverCinfo( - "MarkovSolver", - MarkovSolverBase::initCinfo(), - markovSolverFinfos, - sizeof( markovSolverFinfos ) / sizeof( Finfo* ), - &dinfo - ); - - return &markovSolverCinfo; + ////////////////////// + //DestFinfos + ////////////////////// + + static DestFinfo process( "process", + "Handles process call", + new ProcOpFunc< MarkovSolver >( &MarkovSolver::process ) ); + + static DestFinfo reinit( "reinit", + "Handles reinit call", + new ProcOpFunc< MarkovSolver >( &MarkovSolver::reinit ) ); + + static Finfo* processShared[] = + { + &process, &reinit + }; + + static SharedFinfo proc( "proc", + "This is a shared message to receive Process message from the" + "scheduler. The first entry is a MsgDest for the Process " + "operation. It has a single argument, ProcInfo, which " + "holds lots of information about current time, thread, dt and" + "so on. The second entry is a MsgDest for the Reinit " + "operation. It also uses ProcInfo.", + processShared, sizeof( processShared ) / sizeof( Finfo* ) + ); + + + static Finfo* markovSolverFinfos[] = + { + &proc, //SharedFinfo + }; + + static Dinfo < MarkovSolver > dinfo; + static Cinfo markovSolverCinfo( + "MarkovSolver", + MarkovSolverBase::initCinfo(), + markovSolverFinfos, + sizeof( markovSolverFinfos ) / sizeof( Finfo* ), + &dinfo + ); + + return &markovSolverCinfo; } static const Cinfo* markovSolverCinfo = MarkovSolver::initCinfo(); MarkovSolver::MarkovSolver() { - ; + ; } MarkovSolver::~MarkovSolver() { - ; + ; } Matrix* MarkovSolver::computePadeApproximant( Matrix* Q1, - unsigned int degreeIndex ) + unsigned int degreeIndex ) { - Matrix *expQ; - Matrix *U, *VplusU, *VminusU, *invVminusU, *Qpower; - vector< unsigned int >* swaps = new vector< unsigned int >; - unsigned int n = Q1->size(); - unsigned int degree = mCandidates[degreeIndex]; - double *padeCoeffs = NULL; - Matrix *V = matAlloc(n); - - //Vector of Matrix pointers. Each entry is an even power of Q. - vector< Matrix* > QevenPowers; - - //Selecting the right coefficient array. - switch (degree) - { - case 13: - padeCoeffs = b13; - break; - - case 9: - padeCoeffs = b9; - break; - - case 7: - padeCoeffs = b7; - break; - - case 5: - padeCoeffs = b5; - break; - - case 3: - padeCoeffs = b3; - break; - } - - - ///////// - //Q2 = Q^2 is computed for all degrees. - //Q4 = Q^4 = Q^2 * Q^2 is computed when degree = 5,7,9,13. - //Q6 = Q^6 = Q^4 * Q^2 is computed when degree = 7,9,13. - //Q8 = Q^8 = Q^4 * Q^4 is computed when degree = 7,9. - //Note that the formula for the 13th degree approximant used here - //is different from the one used for smaller degrees. - //////// - switch( degree ) - { - case 3 : - case 5 : - case 7 : - case 9 : - U = matAlloc( n ); - - QevenPowers.push_back( Q1 ); - - for( unsigned int i = 1; i < (degree + 1)/2 ; ++i ) - { - Qpower = QevenPowers.back(); - QevenPowers.push_back( matMatMul( Qpower, Qpower ) ); - } - - //Computation of U. - for ( int i = degree; i > 1; i -= 2 ) - matMatAdd( U, QevenPowers[i/2], 1.0, padeCoeffs[i], FIRST ); - - //Adding b0 * I . - matEyeAdd( U, padeCoeffs[1], 0 ); - matMatMul( Q1, U, SECOND ); - - //Computation of V. - for ( int i = degree - 1; i > 0; i -= 2 ) - matMatAdd( V, QevenPowers[i/2], 1.0, padeCoeffs[i], FIRST ); - - //Adding b1 * I - matEyeAdd( V, padeCoeffs[0], DUMMY ); - - while (!QevenPowers.empty()) - { - delete QevenPowers.back(); - QevenPowers.pop_back(); - } - break; - - case 13: - Matrix *Q2, *Q4, *Q6; - Matrix *temp; - - Q2 = matMatMul( Q1, Q1 ); - Q4 = matMatMul( Q2, Q2 ); - Q6 = matMatMul( Q4, Q2 ); - - //Long and rather messy expression for U and V. - //Refer paper mentioned in header for more details. - //Storing the result in temporaries is a better idea as it gives us - //control on how many temporaries are being created and also to - //help in memory deallocation. - - //Computation of U. - temp = matScalShift( Q6, b13[13], 0.0 ); - matMatAdd( temp, Q4, 1.0, b13[11], FIRST ); - matMatAdd( temp, Q2, 1.0, b13[9], FIRST ); - - matMatMul( Q6, temp, SECOND ); - - matMatAdd( temp, Q6, 1.0, b13[7], FIRST ); - matMatAdd( temp, Q4, 1.0, b13[5], FIRST ); - matMatAdd( temp, Q2, 1.0, b13[3], FIRST ); - matEyeAdd( temp, b13[1], DUMMY ); - U = matMatMul( Q1, temp ); - delete temp; - - //Computation of V - temp = matScalShift( Q6, b13[12], 0.0 ); - matMatAdd( temp, Q4, 1.0, b13[10], FIRST ); - matMatAdd( temp, Q2, 1.0, b13[8], FIRST ); - matMatMul( Q6, temp, SECOND ); - matMatAdd( temp, Q6, 1.0, b13[6], FIRST ); - matMatAdd( temp, Q4, 1.0, b13[4], FIRST ); - matMatAdd( temp, Q2, 1.0, b13[2], FIRST ); - delete( V ); - V = matEyeAdd( temp, b13[0] ); - delete temp; - - delete Q2; - delete Q4; - delete Q6; - break; - } - - VplusU = matMatAdd( U, V, 1.0, 1.0 ); - VminusU = matMatAdd( U, V, -1.0, 1.0 ); - - invVminusU = matAlloc( n ); - matInv( VminusU, swaps, invVminusU ); - expQ = matMatMul( invVminusU, VplusU ); - - - //Memory cleanup. - delete U; - delete V; - delete VplusU; - delete VminusU; - delete invVminusU; - delete swaps; - - return expQ; + Matrix *expQ; + Matrix *U, *VplusU, *VminusU, *invVminusU, *Qpower; + vector< unsigned int >* swaps = new vector< unsigned int >; + unsigned int n = Q1->size(); + unsigned int degree = mCandidates[degreeIndex]; + double *padeCoeffs = NULL; + Matrix *V = matAlloc(n); + + //Vector of Matrix pointers. Each entry is an even power of Q. + vector< Matrix* > QevenPowers; + + //Selecting the right coefficient array. + switch (degree) + { + case 13: + padeCoeffs = b13; + break; + + case 9: + padeCoeffs = b9; + break; + + case 7: + padeCoeffs = b7; + break; + + case 5: + padeCoeffs = b5; + break; + + case 3: + padeCoeffs = b3; + break; + } + + + ///////// + //Q2 = Q^2 is computed for all degrees. + //Q4 = Q^4 = Q^2 * Q^2 is computed when degree = 5,7,9,13. + //Q6 = Q^6 = Q^4 * Q^2 is computed when degree = 7,9,13. + //Q8 = Q^8 = Q^4 * Q^4 is computed when degree = 7,9. + //Note that the formula for the 13th degree approximant used here + //is different from the one used for smaller degrees. + //////// + switch( degree ) + { + case 3 : + case 5 : + case 7 : + case 9 : + U = matAlloc( n ); + + QevenPowers.push_back( Q1 ); + + for( unsigned int i = 1; i < (degree + 1)/2 ; ++i ) + { + Qpower = QevenPowers.back(); + QevenPowers.push_back( matMatMul( Qpower, Qpower ) ); + } + + //Computation of U. + for ( int i = degree; i > 1; i -= 2 ) + matMatAdd( U, QevenPowers[i/2], 1.0, padeCoeffs[i], FIRST ); + + //Adding b0 * I . + matEyeAdd( U, padeCoeffs[1], 0 ); + matMatMul( Q1, U, SECOND ); + + //Computation of V. + for ( int i = degree - 1; i > 0; i -= 2 ) + matMatAdd( V, QevenPowers[i/2], 1.0, padeCoeffs[i], FIRST ); + + //Adding b1 * I + matEyeAdd( V, padeCoeffs[0], DUMMY ); + + while (!QevenPowers.empty()) + { + delete QevenPowers.back(); + QevenPowers.pop_back(); + } + break; + + case 13: + Matrix *Q2, *Q4, *Q6; + Matrix *temp; + + Q2 = matMatMul( Q1, Q1 ); + Q4 = matMatMul( Q2, Q2 ); + Q6 = matMatMul( Q4, Q2 ); + + //Long and rather messy expression for U and V. + //Refer paper mentioned in header for more details. + //Storing the result in temporaries is a better idea as it gives us + //control on how many temporaries are being created and also to + //help in memory deallocation. + + //Computation of U. + temp = matScalShift( Q6, b13[13], 0.0 ); + matMatAdd( temp, Q4, 1.0, b13[11], FIRST ); + matMatAdd( temp, Q2, 1.0, b13[9], FIRST ); + + matMatMul( Q6, temp, SECOND ); + + matMatAdd( temp, Q6, 1.0, b13[7], FIRST ); + matMatAdd( temp, Q4, 1.0, b13[5], FIRST ); + matMatAdd( temp, Q2, 1.0, b13[3], FIRST ); + matEyeAdd( temp, b13[1], DUMMY ); + U = matMatMul( Q1, temp ); + delete temp; + + //Computation of V + temp = matScalShift( Q6, b13[12], 0.0 ); + matMatAdd( temp, Q4, 1.0, b13[10], FIRST ); + matMatAdd( temp, Q2, 1.0, b13[8], FIRST ); + matMatMul( Q6, temp, SECOND ); + matMatAdd( temp, Q6, 1.0, b13[6], FIRST ); + matMatAdd( temp, Q4, 1.0, b13[4], FIRST ); + matMatAdd( temp, Q2, 1.0, b13[2], FIRST ); + delete( V ); + V = matEyeAdd( temp, b13[0] ); + delete temp; + + delete Q2; + delete Q4; + delete Q6; + break; + } + + VplusU = matMatAdd( U, V, 1.0, 1.0 ); + VminusU = matMatAdd( U, V, -1.0, 1.0 ); + + invVminusU = matAlloc( n ); + matInv( VminusU, swaps, invVminusU ); + expQ = matMatMul( invVminusU, VplusU ); + + + //Memory cleanup. + delete U; + delete V; + delete VplusU; + delete VminusU; + delete invVminusU; + delete swaps; + + return expQ; } Matrix* MarkovSolver::computeMatrixExponential() { - double mu, norm; - unsigned int n = Q_->size(); - Matrix *expQ, *Q1; - - mu = matTrace( Q_ )/n; - - //Q1 <- Q - mu*I - //This reduces the norm of the matrix. The idea is that a lower - //order approximant will suffice if the norm is smaller. - Q1 = matEyeAdd( Q_, -mu ); - - //We cycle through the first four candidate values of m. The moment the norm - //satisfies the theta_M bound, we choose that m and compute the Pade' - //approximant to the exponential. We can then directly return the exponential. - norm = matColNorm( Q1 ); - for ( unsigned int i = 0; i < 4; ++i ) - { - if ( norm < thetaM[i] ) - { - expQ = computePadeApproximant( Q1, i ); - matScalShift( expQ, exp( mu ), 0, DUMMY ); - return expQ; - } - } - - //In case none of the candidates were satisfactory, we scale down the norm - //by dividing A by 2^s until ||A|| < 1. We then use a 13th degree - //Pade approximant. - double sf = ceil( log( norm / thetaM[4] ) / log( (double)2 ) ); - unsigned int s = 0; - - if ( sf > 0 ) - { - s = static_cast< unsigned int >( sf ); - matScalShift( Q1, 1.0/(2 << (s - 1)), 0, DUMMY ); - } - expQ = computePadeApproximant( Q1, 4 ); - - //Upto this point, the matrix stored in expQ is r13, the 13th degree - //Pade approximant corresponding to A/2^s, not A. - //Now we repeatedly square r13 's' times to get the exponential - //of A. - for ( unsigned int i = 0; i < s; ++i ) - matMatMul( expQ, expQ, FIRST ); - - matScalShift( expQ, exp( mu ), 0, DUMMY ); - - delete Q1; - return expQ; + double mu, norm; + unsigned int n = Q_->size(); + Matrix *expQ, *Q1; + + mu = matTrace( Q_ )/n; + + //Q1 <- Q - mu*I + //This reduces the norm of the matrix. The idea is that a lower + //order approximant will suffice if the norm is smaller. + Q1 = matEyeAdd( Q_, -mu ); + + //We cycle through the first four candidate values of m. The moment the norm + //satisfies the theta_M bound, we choose that m and compute the Pade' + //approximant to the exponential. We can then directly return the exponential. + norm = matColNorm( Q1 ); + for ( unsigned int i = 0; i < 4; ++i ) + { + if ( norm < thetaM[i] ) + { + expQ = computePadeApproximant( Q1, i ); + matScalShift( expQ, exp( mu ), 0, DUMMY ); + return expQ; + } + } + + //In case none of the candidates were satisfactory, we scale down the norm + //by dividing A by 2^s until ||A|| < 1. We then use a 13th degree + //Pade approximant. + double sf = ceil( log( norm / thetaM[4] ) / log( (double)2 ) ); + unsigned int s = 0; + + if ( sf > 0 ) + { + s = static_cast< unsigned int >( sf ); + matScalShift( Q1, 1.0/(2 << (s - 1)), 0, DUMMY ); + } + expQ = computePadeApproximant( Q1, 4 ); + + //Upto this point, the matrix stored in expQ is r13, the 13th degree + //Pade approximant corresponding to A/2^s, not A. + //Now we repeatedly square r13 's' times to get the exponential + //of A. + for ( unsigned int i = 0; i < s; ++i ) + matMatMul( expQ, expQ, FIRST ); + + matScalShift( expQ, exp( mu ), 0, DUMMY ); + + delete Q1; + return expQ; } /////////////// @@ -285,22 +284,22 @@ Matrix* MarkovSolver::computeMatrixExponential() ////////////// void MarkovSolver::reinit( const Eref& e, ProcPtr p ) { - MarkovSolverBase::reinit( e, p ); + MarkovSolverBase::reinit( e, p ); } void MarkovSolver::process( const Eref& e, ProcPtr p ) { - MarkovSolverBase::process( e, p ); + MarkovSolverBase::process( e, p ); } #ifdef DO_UNIT_TESTS void assignMat( Matrix* A, double testMat[3][3] ) { - for ( unsigned int i = 0; i < 3; ++i ) - { - for ( unsigned int j = 0; j < 3; ++j ) - (*A)[i][j] = testMat[i][j]; - } + for ( unsigned int i = 0; i < 3; ++i ) + { + for ( unsigned int j = 0; j < 3; ++j ) + (*A)[i][j] = testMat[i][j]; + } } #if 0 @@ -308,219 +307,227 @@ void assignMat( Matrix* A, double testMat[3][3] ) //we test out all degrees of the Pade approximant. void testMarkovSolver() { - MarkovSolver solver; - - Matrix *expQ; - - solver.Q_ = matAlloc( 3 ); - - double testMats[5][3][3] = { - { //Will require 3rd degree Pade approximant. - { 0.003859554140797, 0.003828667792972, 0.000567545354509 }, - { 0.000630452326710, 0.001941502594891, 0.001687045130505 }, - { 0.003882371127042, 0.003201121875555, 0.003662942100756 } - }, - { //Will require 5th degree Pade approximant. - { 0.009032772098686, 0.000447799046538, 0.009951718245937 }, - { 0.004122293240791, 0.005703676675533, 0.010337598714782 }, - { 0.012352886634899, 0.004960259942209, 0.002429343859207 } - }, - { //Will require 7th degree Pade approximant. - { 0.249033679156721, 0.026663192545146, 0.193727616177876 }, - { 0.019543882188296, 0.240474520213763, 0.204325805163358 }, - { 0.110669567443862, 0.001158556033517, 0.217173676340877 } - }, - { //Will require 9th degree Pade approximant. - { 0.708590392291234, 0.253289557366033, 0.083402066470341 }, - { 0.368148069351060, 0.675040384813246, 0.585189051240853 }, - { 0.366939478800014, 0.276935085840161, 0.292304127720940 }, - }, - { //Will require 13th degree Pade approximant., - { 5.723393958255834, 2.650265678621879, 2.455089725038183 }, - { 5.563819918184171, 5.681063207977340, 6.573010933999208 }, - { 4.510226911355842, 3.729779121596184, 6.131599680450886 } - } - }; - - double correctExps[5][3][3] = { - { - { 1.003869332885896, 0.003840707339656, 0.000572924780299 }, - { 0.000635569997632, 1.001947309951925, 0.001691961742250 }, - { 0.003898019965821, 0.003217566115560, 1.003673477658833 } - }, - { - { 1.009136553319587, 0.000475947724581, 0.010011555682222 }, - { 0.004217120071231, 1.005746704105642, 0.010400661735110 }, - { 0.012434554824033, 0.004983402186283, 1.002519822649746 } - }, - { - { 1.296879503336410, 0.034325768324765, 0.248960074106229 }, - { 0.039392602584525, 1.272463533413523, 0.260228538022164 }, - { 0.140263068698347, 0.003332731847933, 1.256323839764616 } - }, - { - { 2.174221102665550, 0.550846463313377, 0.279203836454105 }, - { 0.963674962388503, 2.222317715620410, 1.033020817699738 }, - { 0.733257221615105, 0.559435366776953, 1.508376826849517 } - }, - { - { 3.274163243250245e5, 2.405867301580962e5, 3.034390382218154e5 }, - { 5.886803379935408e5, 4.325844111569120e5, 5.456065763194024e5 }, - { 4.671930521670584e5, 3.433084310645007e5, 4.330101744194682e5 } - } - }; - - double correctColumnNorms[5] = { - 1.009005583407142, - 1.025788228214852, - 1.765512451893010, - 3.871153286669158, - 1.383289714485624e+06 - }; - - for ( unsigned int i = 0; i < 5; ++i ) - { - assignMat( solver.Q_, testMats[i] ); - expQ = solver.computeMatrixExponential(); - assert( doubleEq( matColNorm( expQ ), correctColumnNorms[i] ) ); - - //Comparing termwise just to be doubly sure. - for( unsigned int j = 0; j < 3; ++j ) - { - for( unsigned int k = 0; k < 3; ++k ) - assert( doubleEq( (*expQ)[j][k], correctExps[i][j][k] ) ); - } - - delete expQ; - } - - ///////////////// - //Testing state space interpolation. - //////////////// - const Cinfo* rateTableCinfo = MarkovRateTable::initCinfo(); - const Cinfo* interpol2dCinfo = Interpol2D::initCinfo(); - const Cinfo* vectorTableCinfo = VectorTable::initCinfo(); - const Cinfo* markovSolverCinfo = MarkovSolver::initCinfo(); - - vector< DimInfo > single; - - Id rateTable2dId = Id::nextId(); - Id rateTable1dId = Id::nextId(); - Id int2dId = Id::nextId(); - Id vecTableId = Id::nextId(); - Id solver2dId = Id::nextId(); - Id solver1dId = Id::nextId(); - - ObjId eRateTable2d = Element( rateTable2dId, rateTableCinfo, - "rateTable2d", single, 1, true ); - Element *eRateTable1d = new Element( rateTable1dId, rateTableCinfo, - "rateTable1d", single, 1, true ); - Element *eInt2d = new Element( int2dId, interpol2dCinfo, "int2d", single, 1 ); - Element *eVecTable = new Element( vecTableId, vectorTableCinfo, "vecTable", - single, 1, true ); - Element *eSolver2d = new Element( solver2dId, markovSolverCinfo, - "solver2d", single, 1, true ); - Element *eSolver1d = new Element( solver1dId, markovSolverCinfo, - "solver1d", single, 1, true ); - - Eref rateTable2dEref( eRateTable2d, 0 ); - Eref rateTable1dEref( eRateTable1d, 0 ); - Eref int2dEref( eInt2d, 0 ); - Eref vecTableEref( eVecTable, 0 ); - Eref solver2dEref( eSolver2d, 0 ); - Eref solver1dEref( eSolver1d, 0 ); - - vector< double > table1d; - vector< vector< double > > table2d; - double v, conc; - - MarkovRateTable* rateTable2d = reinterpret_cast< MarkovRateTable* > - ( rateTable2dEref.data() ); - MarkovRateTable* rateTable1d = reinterpret_cast< MarkovRateTable* > - ( rateTable1dEref.data() ); - VectorTable* vecTable = reinterpret_cast< VectorTable* > - ( vecTableEref.data() ); - Interpol2D* int2d = reinterpret_cast< Interpol2D* > - ( int2dEref.data() ); - MarkovSolver* markovSolver2d = reinterpret_cast< MarkovSolver* > - ( solver2dEref.data() ); - MarkovSolver* markovSolver1d = reinterpret_cast< MarkovSolver* > - ( solver1dEref.data() ); - - rateTable2d->init( 3 ); - rateTable1d->init( 3 ); - - vecTable->setMin( -0.10 ); - vecTable->setMax( 0.10 ); - vecTable->setDiv( 200 ); - - v = vecTable->getMin(); - for ( unsigned int i = 0; i < 201; ++i ) - { - table1d.push_back( 1e3 * exp( 9 * v - 0.45 ) ); - v += 0.001; - } - - vecTable->setTable( table1d ); - - rateTable2d->setVtChildTable( 1, 3, vecTableId, 0 ); - rateTable1d->setVtChildTable( 1, 3, vecTableId, 0 ); - - int2d->setXmin( -0.10 ); - int2d->setXmax( 0.10 ); - int2d->setYmin( 0 ); - int2d->setYmax( 50e-6 ); - - v = int2d->getXmin(); - table2d.resize( 201 ); - for ( unsigned int i = 0; i < 201; ++i ) - { - conc = int2d->getYmin(); - for ( unsigned int j = 0; j < 50; ++j ) - { - table2d[i].push_back( 1e3 * conc * exp( -45 * v + 0.65 ) ); - conc += 1e-6; - } - v += 0.001; - } - - int2d->setTableVector( table2d ); - - rateTable2d->setInt2dChildTable( 2, 3, int2dId ); - - rateTable2d->setConstantRate( 3, 1, 0.652 ); - rateTable2d->setConstantRate( 2, 1, 1.541 ); - - rateTable1d->setConstantRate( 3, 1, 0.652 ); - rateTable1d->setConstantRate( 2, 1, 1.541 ); - - markovSolver2d->init( rateTable2dId, 0.1 ); - markovSolver1d->init( rateTable1dId, 0.1 ); - - markovSolver2d->setVm( 0.0533 ); - markovSolver2d->setLigandConc( 3.41e-6 ); - - markovSolver1d->setVm( 0.0533 ); - - Vector initState; - - initState.push_back( 0.2 ); - initState.push_back( 0.4 ); - initState.push_back( 0.4 ); - - markovSolver2d->setInitialState( initState ); - markovSolver2d->computeState(); - - Vector interpState = markovSolver2d->getState(); - - rateTable2dId.destroy(); - rateTable1dId.destroy(); - int2dId.destroy(); - vecTableId.destroy(); - solver2dId.destroy(); - solver1dId.destroy(); - - cout << "." << flush; + MarkovSolver solver; + + Matrix *expQ; + + solver.Q_ = matAlloc( 3 ); + + double testMats[5][3][3] = + { + { + //Will require 3rd degree Pade approximant. + { 0.003859554140797, 0.003828667792972, 0.000567545354509 }, + { 0.000630452326710, 0.001941502594891, 0.001687045130505 }, + { 0.003882371127042, 0.003201121875555, 0.003662942100756 } + }, + { + //Will require 5th degree Pade approximant. + { 0.009032772098686, 0.000447799046538, 0.009951718245937 }, + { 0.004122293240791, 0.005703676675533, 0.010337598714782 }, + { 0.012352886634899, 0.004960259942209, 0.002429343859207 } + }, + { + //Will require 7th degree Pade approximant. + { 0.249033679156721, 0.026663192545146, 0.193727616177876 }, + { 0.019543882188296, 0.240474520213763, 0.204325805163358 }, + { 0.110669567443862, 0.001158556033517, 0.217173676340877 } + }, + { + //Will require 9th degree Pade approximant. + { 0.708590392291234, 0.253289557366033, 0.083402066470341 }, + { 0.368148069351060, 0.675040384813246, 0.585189051240853 }, + { 0.366939478800014, 0.276935085840161, 0.292304127720940 }, + }, + { + //Will require 13th degree Pade approximant., + { 5.723393958255834, 2.650265678621879, 2.455089725038183 }, + { 5.563819918184171, 5.681063207977340, 6.573010933999208 }, + { 4.510226911355842, 3.729779121596184, 6.131599680450886 } + } + }; + + double correctExps[5][3][3] = + { + { + { 1.003869332885896, 0.003840707339656, 0.000572924780299 }, + { 0.000635569997632, 1.001947309951925, 0.001691961742250 }, + { 0.003898019965821, 0.003217566115560, 1.003673477658833 } + }, + { + { 1.009136553319587, 0.000475947724581, 0.010011555682222 }, + { 0.004217120071231, 1.005746704105642, 0.010400661735110 }, + { 0.012434554824033, 0.004983402186283, 1.002519822649746 } + }, + { + { 1.296879503336410, 0.034325768324765, 0.248960074106229 }, + { 0.039392602584525, 1.272463533413523, 0.260228538022164 }, + { 0.140263068698347, 0.003332731847933, 1.256323839764616 } + }, + { + { 2.174221102665550, 0.550846463313377, 0.279203836454105 }, + { 0.963674962388503, 2.222317715620410, 1.033020817699738 }, + { 0.733257221615105, 0.559435366776953, 1.508376826849517 } + }, + { + { 3.274163243250245e5, 2.405867301580962e5, 3.034390382218154e5 }, + { 5.886803379935408e5, 4.325844111569120e5, 5.456065763194024e5 }, + { 4.671930521670584e5, 3.433084310645007e5, 4.330101744194682e5 } + } + }; + + double correctColumnNorms[5] = + { + 1.009005583407142, + 1.025788228214852, + 1.765512451893010, + 3.871153286669158, + 1.383289714485624e+06 + }; + + for ( unsigned int i = 0; i < 5; ++i ) + { + assignMat( solver.Q_, testMats[i] ); + expQ = solver.computeMatrixExponential(); + assert( doubleEq( matColNorm( expQ ), correctColumnNorms[i] ) ); + + //Comparing termwise just to be doubly sure. + for( unsigned int j = 0; j < 3; ++j ) + { + for( unsigned int k = 0; k < 3; ++k ) + assert( doubleEq( (*expQ)[j][k], correctExps[i][j][k] ) ); + } + + delete expQ; + } + + ///////////////// + //Testing state space interpolation. + //////////////// + const Cinfo* rateTableCinfo = MarkovRateTable::initCinfo(); + const Cinfo* interpol2dCinfo = Interpol2D::initCinfo(); + const Cinfo* vectorTableCinfo = VectorTable::initCinfo(); + const Cinfo* markovSolverCinfo = MarkovSolver::initCinfo(); + + vector< DimInfo > single; + + Id rateTable2dId = Id::nextId(); + Id rateTable1dId = Id::nextId(); + Id int2dId = Id::nextId(); + Id vecTableId = Id::nextId(); + Id solver2dId = Id::nextId(); + Id solver1dId = Id::nextId(); + + ObjId eRateTable2d = Element( rateTable2dId, rateTableCinfo, + "rateTable2d", single, 1, true ); + Element *eRateTable1d = new Element( rateTable1dId, rateTableCinfo, + "rateTable1d", single, 1, true ); + Element *eInt2d = new Element( int2dId, interpol2dCinfo, "int2d", single, 1 ); + Element *eVecTable = new Element( vecTableId, vectorTableCinfo, "vecTable", + single, 1, true ); + Element *eSolver2d = new Element( solver2dId, markovSolverCinfo, + "solver2d", single, 1, true ); + Element *eSolver1d = new Element( solver1dId, markovSolverCinfo, + "solver1d", single, 1, true ); + + Eref rateTable2dEref( eRateTable2d, 0 ); + Eref rateTable1dEref( eRateTable1d, 0 ); + Eref int2dEref( eInt2d, 0 ); + Eref vecTableEref( eVecTable, 0 ); + Eref solver2dEref( eSolver2d, 0 ); + Eref solver1dEref( eSolver1d, 0 ); + + vector< double > table1d; + vector< vector< double > > table2d; + double v, conc; + + MarkovRateTable* rateTable2d = reinterpret_cast< MarkovRateTable* > + ( rateTable2dEref.data() ); + MarkovRateTable* rateTable1d = reinterpret_cast< MarkovRateTable* > + ( rateTable1dEref.data() ); + VectorTable* vecTable = reinterpret_cast< VectorTable* > + ( vecTableEref.data() ); + Interpol2D* int2d = reinterpret_cast< Interpol2D* > + ( int2dEref.data() ); + MarkovSolver* markovSolver2d = reinterpret_cast< MarkovSolver* > + ( solver2dEref.data() ); + MarkovSolver* markovSolver1d = reinterpret_cast< MarkovSolver* > + ( solver1dEref.data() ); + + rateTable2d->init( 3 ); + rateTable1d->init( 3 ); + + vecTable->setMin( -0.10 ); + vecTable->setMax( 0.10 ); + vecTable->setDiv( 200 ); + + v = vecTable->getMin(); + for ( unsigned int i = 0; i < 201; ++i ) + { + table1d.push_back( 1e3 * exp( 9 * v - 0.45 ) ); + v += 0.001; + } + + vecTable->setTable( table1d ); + + rateTable2d->setVtChildTable( 1, 3, vecTableId, 0 ); + rateTable1d->setVtChildTable( 1, 3, vecTableId, 0 ); + + int2d->setXmin( -0.10 ); + int2d->setXmax( 0.10 ); + int2d->setYmin( 0 ); + int2d->setYmax( 50e-6 ); + + v = int2d->getXmin(); + table2d.resize( 201 ); + for ( unsigned int i = 0; i < 201; ++i ) + { + conc = int2d->getYmin(); + for ( unsigned int j = 0; j < 50; ++j ) + { + table2d[i].push_back( 1e3 * conc * exp( -45 * v + 0.65 ) ); + conc += 1e-6; + } + v += 0.001; + } + + int2d->setTableVector( table2d ); + + rateTable2d->setInt2dChildTable( 2, 3, int2dId ); + + rateTable2d->setConstantRate( 3, 1, 0.652 ); + rateTable2d->setConstantRate( 2, 1, 1.541 ); + + rateTable1d->setConstantRate( 3, 1, 0.652 ); + rateTable1d->setConstantRate( 2, 1, 1.541 ); + + markovSolver2d->init( rateTable2dId, 0.1 ); + markovSolver1d->init( rateTable1dId, 0.1 ); + + markovSolver2d->setVm( 0.0533 ); + markovSolver2d->setLigandConc( 3.41e-6 ); + + markovSolver1d->setVm( 0.0533 ); + + Vector initState; + + initState.push_back( 0.2 ); + initState.push_back( 0.4 ); + initState.push_back( 0.4 ); + + markovSolver2d->setInitialState( initState ); + markovSolver2d->computeState(); + + Vector interpState = markovSolver2d->getState(); + + rateTable2dId.destroy(); + rateTable1dId.destroy(); + int2dId.destroy(); + vecTableId.destroy(); + solver2dId.destroy(); + solver1dId.destroy(); + + cout << "." << flush; } #endif // #if 0 diff --git a/biophysics/testBiophysics.cpp b/biophysics/testBiophysics.cpp index e46b301257..0dde4834ba 100644 --- a/biophysics/testBiophysics.cpp +++ b/biophysics/testBiophysics.cpp @@ -105,8 +105,8 @@ void testIntFireNetwork( unsigned int runsteps = 5 ) // by multiple threads if the above Set call is not complete. vector< double > origVm( size, 0.0 ); - - // NOTE: From + + // NOTE: From moose::mtseed( 5489UL ); for ( unsigned int i = 0; i < size; ++i ) origVm[i] = moose::mtrand() * Vmax; @@ -214,7 +214,7 @@ void testIntFireNetwork( unsigned int runsteps = 5 ) #if 0 cout << endl; - cout << std::setprecision(12) << retVm100 << endl; + cout << std::setprecision(12) << retVm100 << endl; cout << std::setprecision(12) << retVm101 << endl; cout << std::setprecision(12) << retVm102 << endl; cout << std::setprecision(11) << retVm99 << endl; @@ -746,26 +746,26 @@ void testHHChannel() //Sample current obtained from channel in Chapter 20, Sakmann & Neher, Pg. 603. //The current is sampled at intervals of 10 usec. static double sampleCurrent[] = {0.0, // This is to handle the value sent by ChanBase on reinit - 0.0000000e+00, 3.0005743e-26, 1.2004594e-25, 2.7015505e-25, 4.8036751e-25, 7.5071776e-25, - 1.0812402e-24, 1.4719693e-24, 1.9229394e-24, 2.4341850e-24, 3.0057404e-24, 3.6376401e-24, - 4.3299183e-24, 5.0826095e-24, 5.8957481e-24, 6.7693684e-24, 7.7035046e-24, 8.6981913e-24, - 9.7534627e-24, 1.0869353e-23, 1.2045897e-23, 1.3283128e-23, 1.4581082e-23, 1.5939791e-23, - 1.7359292e-23, 1.8839616e-23, 2.0380801e-23, 2.1982878e-23, 2.3645883e-23, 2.5369850e-23, - 2.7154813e-23, 2.9000806e-23, 3.0907863e-23, 3.2876020e-23, 3.4905309e-23, 3.6995766e-23, - 3.9147423e-23, 4.1360317e-23, 4.3634480e-23, 4.5969946e-23, 4.8366751e-23, 5.0824928e-23, - 5.3344511e-23, 5.5925535e-23, 5.8568033e-23, 6.1272040e-23, 6.4037589e-23, 6.6864716e-23, - 6.9753453e-23, 7.2703835e-23, 7.5715897e-23, 7.8789672e-23, 8.1925194e-23, 8.5122497e-23, - 8.8381616e-23, 9.1702584e-23, 9.5085435e-23, 9.8530204e-23, 1.0203692e-22, 1.0560563e-22, - 1.0923636e-22, 1.1292913e-22, 1.1668400e-22, 1.2050099e-22, 1.2438013e-22, 1.2832146e-22, - 1.3232502e-22, 1.3639083e-22, 1.4051894e-22, 1.4470937e-22, 1.4896215e-22, 1.5327733e-22, - 1.5765494e-22, 1.6209501e-22, 1.6659757e-22, 1.7116267e-22, 1.7579032e-22, 1.8048057e-22, - 1.8523345e-22, 1.9004900e-22, 1.9492724e-22, 1.9986821e-22, 2.0487195e-22, 2.0993849e-22, - 2.1506786e-22, 2.2026010e-22, 2.2551524e-22, 2.3083331e-22, 2.3621436e-22, 2.4165840e-22, - 2.4716548e-22, 2.5273563e-22, 2.5836888e-22, 2.6406527e-22, 2.6982483e-22, 2.7564760e-22, - 2.8153360e-22, 2.8748287e-22, 2.9349545e-22, 2.9957137e-22, 3.0571067e-22 - }; - -void testMarkovGslSolver() + 0.0000000e+00, 3.0005743e-26, 1.2004594e-25, 2.7015505e-25, 4.8036751e-25, 7.5071776e-25, + 1.0812402e-24, 1.4719693e-24, 1.9229394e-24, 2.4341850e-24, 3.0057404e-24, 3.6376401e-24, + 4.3299183e-24, 5.0826095e-24, 5.8957481e-24, 6.7693684e-24, 7.7035046e-24, 8.6981913e-24, + 9.7534627e-24, 1.0869353e-23, 1.2045897e-23, 1.3283128e-23, 1.4581082e-23, 1.5939791e-23, + 1.7359292e-23, 1.8839616e-23, 2.0380801e-23, 2.1982878e-23, 2.3645883e-23, 2.5369850e-23, + 2.7154813e-23, 2.9000806e-23, 3.0907863e-23, 3.2876020e-23, 3.4905309e-23, 3.6995766e-23, + 3.9147423e-23, 4.1360317e-23, 4.3634480e-23, 4.5969946e-23, 4.8366751e-23, 5.0824928e-23, + 5.3344511e-23, 5.5925535e-23, 5.8568033e-23, 6.1272040e-23, 6.4037589e-23, 6.6864716e-23, + 6.9753453e-23, 7.2703835e-23, 7.5715897e-23, 7.8789672e-23, 8.1925194e-23, 8.5122497e-23, + 8.8381616e-23, 9.1702584e-23, 9.5085435e-23, 9.8530204e-23, 1.0203692e-22, 1.0560563e-22, + 1.0923636e-22, 1.1292913e-22, 1.1668400e-22, 1.2050099e-22, 1.2438013e-22, 1.2832146e-22, + 1.3232502e-22, 1.3639083e-22, 1.4051894e-22, 1.4470937e-22, 1.4896215e-22, 1.5327733e-22, + 1.5765494e-22, 1.6209501e-22, 1.6659757e-22, 1.7116267e-22, 1.7579032e-22, 1.8048057e-22, + 1.8523345e-22, 1.9004900e-22, 1.9492724e-22, 1.9986821e-22, 2.0487195e-22, 2.0993849e-22, + 2.1506786e-22, 2.2026010e-22, 2.2551524e-22, 2.3083331e-22, 2.3621436e-22, 2.4165840e-22, + 2.4716548e-22, 2.5273563e-22, 2.5836888e-22, 2.6406527e-22, 2.6982483e-22, 2.7564760e-22, + 2.8153360e-22, 2.8748287e-22, 2.9349545e-22, 2.9957137e-22, 3.0571067e-22 + }; + +void testMarkovOdeSolver() { Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() ); unsigned size = 1; @@ -774,7 +774,7 @@ void testMarkovGslSolver() Id comptId = shell->doCreate( "Compartment", nid, "compt", size ); Id rateTabId = shell->doCreate( "MarkovRateTable", comptId, "rateTab", size ); Id mChanId = shell->doCreate( "MarkovChannel", comptId, "mChan", size ); - Id gslSolverId = shell->doCreate( "MarkovGslSolver", comptId, "gslSolver", size ); + Id gslSolverId = shell->doCreate( "MarkovOdeSolver", comptId, "gslSolver", size ); Id tabId = shell->doCreate( "Table", nid, "tab", size ); @@ -902,9 +902,10 @@ void testMarkovGslSolver() { if (!doubleEq( sampleCurrent[i] * 1e25, vec[i] * 1e25 )) { - cout << "testMarkovGslSolver: sample=" << sampleCurrent[i]*1e25 << " calculated=" << vec[i]*1e25 << endl; + cout << "testMarkovOdeSolver: sample=" << sampleCurrent[i]*1e25 << " calculated=" << vec[i]*1e25 << endl; } ASSERT_DOUBLE_EQ("", sampleCurrent[i] * 1e25, vec[i] * 1e25 ); + // cout << sampleCurrent[i] << " " << vec[i] << endl; } //Currents involved here are incredibly small. Scaling them up is necessary //for the doubleEq function to do its job. @@ -914,16 +915,16 @@ void testMarkovGslSolver() } //////////////// -//The testMarkovGslSolver() function includes the MarkovChannel object, but +//The testMarkovOdeSolver() function includes the MarkovChannel object, but //is a rather trivial case, in that the rates are all constant. -//This test simultaneously tests the MarkovChannel, MarkovGslSolver, +//This test simultaneously tests the MarkovChannel, MarkovOdeSolver, //MarkovSolverBase and MarkovSolver classes. //This test involves simulating the 4-state NMDA channel model specified //in the following paper : //"Voltage Dependence of NMDA-Activated Macroscopic Conductances Predicted //by Single-Channel Kinetics", Craig E. Jahr and Charles F. Stevens, The Journal //of Neuroscience, 1990, 10(9), pp. 3178-3182. -//It is expected that the MarkovGslSolver and the MarkovSolver objects will +//It is expected that the MarkovOdeSolver and the MarkovSolver objects will //give the same answer. // //Note that this is different from the NMDAChan test which involves synapses. @@ -944,12 +945,12 @@ void testMarkovChannel() Id exptlRateTableId = shell->doCreate( "MarkovRateTable", exptlComptId, "exptlRateTable", size ); - Id mChanGslId = shell->doCreate( "MarkovChannel", gslComptId, - "mChanGsl", size ); + Id mChanOdeId = shell->doCreate( "MarkovChannel", gslComptId, + "mChanOde", size ); Id mChanExptlId = shell->doCreate( "MarkovChannel", exptlComptId, "mChanExptl", size ); - Id gslSolverId = shell->doCreate( "MarkovGslSolver", gslComptId, + Id gslSolverId = shell->doCreate( "MarkovOdeSolver", gslComptId, "gslSolver", size ); Id exptlSolverId = shell->doCreate( "MarkovSolver", exptlComptId, "exptlSolver", size ); @@ -972,7 +973,7 @@ void testMarkovChannel() //via its ChanBase base class, sends back the conductance and current through //it. ObjId mid = shell->doAddMsg( "Single", ObjId( gslComptId ), "channel", - ObjId( mChanGslId ), "channel" ); + ObjId( mChanOdeId ), "channel" ); assert( !mid.bad() ); mid = shell->doAddMsg( "Single", ObjId( exptlComptId ), "channel", @@ -980,7 +981,7 @@ void testMarkovChannel() assert( !mid.bad() ); //////// - //Connecting up the MarkovGslSolver. + //Connecting up the MarkovOdeSolver. /////// //Connecting Compartment and MarkovRateTable. @@ -996,20 +997,20 @@ void testMarkovChannel() ObjId( gslRateTableId ), "channel" ); assert( !mid.bad() ); - //Connecting the MarkovRateTable with the MarkovGslSolver object. + //Connecting the MarkovRateTable with the MarkovOdeSolver object. //As mentioned earlier, the MarkovRateTable object sends out information - //about Q to the MarkovGslSolver. The MarkovGslSolver then churns out + //about Q to the MarkovOdeSolver. The MarkovOdeSolver then churns out //the state of the system for the next time step. mid = shell->doAddMsg("Single", ObjId( gslRateTableId ), "instratesOut", ObjId( gslSolverId ), "handleQ" ); - //Connecting MarkovGslSolver with MarkovChannel. - //The MarkovGslSolver object, upon computing the state of the channel, + //Connecting MarkovOdeSolver with MarkovChannel. + //The MarkovOdeSolver object, upon computing the state of the channel, //sends this information to the MarkovChannel object. The MarkovChannel //object will compute the expected conductance of the channel and send //this information to the compartment. mid = shell->doAddMsg( "Single", ObjId( gslSolverId ), "stateOut", - ObjId( mChanGslId ), "handleState" ); + ObjId( mChanOdeId ), "handleState" ); assert( !mid.bad() ); ////////// @@ -1035,7 +1036,7 @@ void testMarkovChannel() //Get the current values from the GSL solver based channel. mid = shell->doAddMsg( "Single", ObjId( gslTableId ), "requestOut", - ObjId( mChanGslId ), "getIk" ); + ObjId( mChanOdeId ), "getIk" ); assert( !mid.bad() ); //Get the current values from the matrix exponential solver based channel. @@ -1069,10 +1070,10 @@ void testMarkovChannel() ///////////////// //Number of states and open states. - Field< unsigned int >::set( mChanGslId, "numStates", 4 ); + Field< unsigned int >::set( mChanOdeId, "numStates", 4 ); Field< unsigned int >::set( mChanExptlId, "numStates", 4 ); - Field< unsigned int >::set( mChanGslId, "numOpenStates", 1 ); + Field< unsigned int >::set( mChanOdeId, "numOpenStates", 1 ); Field< unsigned int >::set( mChanExptlId, "numOpenStates", 1 ); vector< string > stateLabels; @@ -1085,7 +1086,7 @@ void testMarkovChannel() stateLabels.push_back( "B2" ); //State 3. stateLabels.push_back( "C" ); //State 4. - Field< vector< string > >::set( mChanGslId, "labels", stateLabels ); + Field< vector< string > >::set( mChanOdeId, "labels", stateLabels ); Field< vector< string > >::set( mChanExptlId, "labels", stateLabels ); //Setting up conductance value for single open state. Value chosen @@ -1094,7 +1095,7 @@ void testMarkovChannel() gBar.push_back( 5.431553e-9 ); - Field< vector< double > >::set( mChanGslId, "gbar", gBar ); + Field< vector< double > >::set( mChanOdeId, "gbar", gBar ); Field< vector< double > >::set( mChanExptlId, "gbar", gBar ); //Initial state of the system. This is really an arbitrary choice. @@ -1105,7 +1106,7 @@ void testMarkovChannel() initState.push_back( 0.80 ); initState.push_back( 0.00 ); - Field< vector< double > >::set( mChanGslId, "initialState", initState ); + Field< vector< double > >::set( mChanOdeId, "initialState", initState ); Field< vector< double > >::set( mChanExptlId, "initialState", initState ); //This initializes the GSL solver object. @@ -1254,7 +1255,7 @@ void testMarkovChannel() "process", 2 ); shell->doUseClock( "/n/gslCompt/gslSolver,/n/exptlCompt/exptlSolver", "process", 3 ); - shell->doUseClock( "/n/gslCompt/mChanGsl,/n/gslTable","process", 4 ); + shell->doUseClock( "/n/gslCompt/mChanOde,/n/gslTable","process", 4 ); shell->doUseClock( "/n/exptlCompt/mChanExptl,/n/exptlTable", "process", 5 ); shell->doReinit(); @@ -1687,13 +1688,15 @@ void testBiophysicsProcess() testIntFireNetwork(); testCompartmentProcess(); // testMarkovGslSolver(); - // testMarkovChannel(); + testMarkovChannel(); #if 0 testHHChannel(); #endif } #else // ifdef DO_UNIT_TESTS + +// Dummy functions to avoid link error. void testBiophysics() { ; diff --git a/builtins/Table.h b/builtins/Table.h index 2d59f67dfd..22b55a3f54 100644 --- a/builtins/Table.h +++ b/builtins/Table.h @@ -10,6 +10,10 @@ #ifndef _TABLE_H #define _TABLE_H + +#if USE_BOOST_FILESYSTEM +#include +#endif /* ----- USE_BOOST_FILESYSTEM ----- */ #include /** diff --git a/devel/issues/issue_69.py b/devel/issues/issue_69.py deleted file mode 100644 index 93388421a8..0000000000 --- a/devel/issues/issue_69.py +++ /dev/null @@ -1,15 +0,0 @@ -# -*- coding: utf-8 -*- -# Github issue #69, Also BhallaLab/moose-gui#3. - -import moose -print( '[INFO] Importing moose from %s' % moose.__file__ ) -print( '[INFO] Version : %s' % moose.__version__ ) -moose.loadModel('../data/acc27.g','/acc27_1',"gsl") -compts = moose.wildcardFind('/acc27_1/##[ISA=ChemCompt]') -for compt in compts: - if moose.exists(compt.path+'/stoich'): - st = moose.element(compt.path+'/stoich') - #print " stoich ksolve ",st.ksolve.path - if moose.exists((st.ksolve).path): - moose.delete(st.ksolve) - moose.delete(st) diff --git a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gees.hpp b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gees.hpp index 472ef3675b..2c2531c930 100644 --- a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gees.hpp +++ b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gees.hpp @@ -140,9 +140,7 @@ namespace boost { namespace numeric { namespace bindings { >::value)); #endif -#ifndef NDEBUG typedef typename MatrA::value_type value_type ; -#endif int const n = traits::matrix_size1 (a); assert (n == traits::matrix_size2 (a)); @@ -187,9 +185,7 @@ namespace boost { namespace numeric { namespace bindings { >::value)); #endif -#ifndef NDEBUG typedef typename MatrA::value_type value_type ; -#endif int const n = traits::matrix_size1 (a); assert (n == traits::matrix_size2 (a)); @@ -269,9 +265,7 @@ namespace boost { namespace numeric { namespace bindings { inline int operator() (char jobvs, MatrA& a, EigVal& w, SchVec& vs, optimal_workspace ) const { typedef typename MatrA::value_type value_type ; -#ifndef NDEBUG typedef typename traits::type_traits< value_type >::real_type real_type ; -#endif int n = traits::matrix_size1( a ); @@ -284,10 +278,7 @@ namespace boost { namespace numeric { namespace bindings { inline int operator() (char jobvs, MatrA& a, EigVal& w, SchVec& vs, minimal_workspace ) const { typedef typename MatrA::value_type value_type ; - -#ifndef NDEBUG typedef typename traits::type_traits< value_type >::real_type real_type ; -#endif int n = traits::matrix_size1( a ); diff --git a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesdd.hpp b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesdd.hpp index b7d39f0fe3..f0ab760457 100644 --- a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesdd.hpp +++ b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesdd.hpp @@ -428,14 +428,9 @@ namespace boost { namespace numeric { namespace bindings { || (jobz == 'A' && traits::leading_dimension (vt) >= n) || (jobz == 'S' && traits::leading_dimension (vt) >= minmn)); -#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS - typedef typename traits::matrix_traits::value_type val_t; -#else - typedef typename MatrA::value_type val_t; -#endif - assert (traits::vector_size (w) - >= detail::gesdd_min_work (val_t(), jobz, m, n)); - assert (traits::vector_size (iw) >= detail::gesdd_iwork (m, n)); + //assert (traits::vector_size (w) + // >= detail::gesdd_min_work (val_t(), jobz, m, n)); + //assert (traits::vector_size (iw) >= detail::gesdd_iwork (m, n)); int info; detail::gesdd (jobz, m, n, @@ -507,16 +502,13 @@ namespace boost { namespace numeric { namespace bindings { && traits::leading_dimension (vt) >= n) || (jobz == 'A' && traits::leading_dimension (vt) >= n) || (jobz == 'S' && traits::leading_dimension (vt) >= minmn)); -#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS - typedef typename traits::matrix_traits::value_type val_t; -#else - typedef typename MatrA::value_type val_t; -#endif +#if 0 assert (traits::vector_size (w) >= detail::gesdd_min_work (val_t(), jobz, m, n)); assert (traits::vector_size (rw) >= detail::gesdd_rwork (val_t(), jobz, m, n)); assert (traits::vector_size (iw) >= detail::gesdd_iwork (m, n)); +#endif int info; detail::gesdd (jobz, m, n, diff --git a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesvd.hpp b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesvd.hpp index eba01080e5..3533b6918d 100644 --- a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesvd.hpp +++ b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/gesvd.hpp @@ -295,9 +295,7 @@ namespace boost { namespace numeric { namespace bindings { || (jobvt == 'A' && traits::leading_dimension (vt) >= n) || (jobvt == 'S' && traits::leading_dimension (vt) >= minmn)); #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS -#ifndef NDEBUG typedef typename traits::matrix_traits::value_type val_t; -#endif #else typedef typename MatrA::value_type val_t; #endif @@ -363,13 +361,9 @@ namespace boost { namespace numeric { namespace bindings { || (jobvt == 'O' && traits::leading_dimension (vt) >= 1) || (jobvt == 'A' && traits::leading_dimension (vt) >= n) || (jobvt == 'S' && traits::leading_dimension (vt) >= minmn)); -#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS - typedef typename traits::matrix_traits::value_type val_t; -#else - typedef typename MatrA::value_type val_t; -#endif - assert (traits::vector_size(w) >= detail::gesvd_min_work(val_t(),m,n)); - assert (traits::vector_size(rw) >= detail::gesvd_rwork(val_t(),m,n)); + + // assert (traits::vector_size(w) >= detail::gesvd_min_work(val_t(),m,n)); + // assert (traits::vector_size(rw) >= detail::gesvd_rwork(val_t(),m,n)); int info; detail::gesvd (jobu, jobvt, m, n, diff --git a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/ormqr.hpp b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/ormqr.hpp index efbec08045..f7be063c2e 100644 --- a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/ormqr.hpp +++ b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/ormqr.hpp @@ -214,7 +214,7 @@ namespace boost { namespace numeric { namespace bindings { template inline int ormqr (char side, char trans, const A& a, const Tau& tau, C& c, detail::workspace1 workspace ) { - typedef typename A::value_type value_type ; + // typedef typename A::value_type value_type ; return detail::ormqr( side, trans, a, tau, c, workspace.w_ ); } diff --git a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/syev.hpp b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/syev.hpp index fdd5253bcb..89e75bee45 100644 --- a/external/boost-numeric-bindings/boost/numeric/bindings/lapack/syev.hpp +++ b/external/boost-numeric-bindings/boost/numeric/bindings/lapack/syev.hpp @@ -137,9 +137,7 @@ namespace boost { namespace numeric { namespace bindings { inline int syev (char jobz, char uplo, A& a, W& w, detail::workspace1 workspace ) { -#ifndef NDEBUG typedef typename A::value_type value_type ; -#endif return detail::syev(jobz, uplo, a, w, workspace.w_); } // syev() @@ -198,9 +196,7 @@ namespace boost { namespace numeric { namespace bindings { template inline int syev (char jobz, A& a, W& w, detail::workspace1 workspace ) { -#ifndef NDEBUG typedef typename A::value_type value_type ; -#endif char uplo = traits::matrix_uplo_tag( a ) ; #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK typedef typename traits::matrix_traits::matrix_structure matrix_structure ; diff --git a/external/libsoda/.travis.yml b/external/libsoda/.travis.yml new file mode 100644 index 0000000000..df47800d9e --- /dev/null +++ b/external/libsoda/.travis.yml @@ -0,0 +1,25 @@ +language: cpp +branches: + only: + - master + - devel + +compiler: + - clang + - gcc + +matrix: + # works on Precise and Trusty + - os: linux + addons: + apt: + sources: + - ubuntu-toolchain-r-test + packages: + - g++-5 + env: + - MATRIX_EVAL="CC=gcc-5 && CXX=g++-5" + +script: + - cmake . && make && ctest --output-on-failure + diff --git a/external/libsoda/CMakeLists.txt b/external/libsoda/CMakeLists.txt new file mode 100644 index 0000000000..960fabdb4c --- /dev/null +++ b/external/libsoda/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 2.8) + +add_library(lsoda ${CMAKE_CURRENT_SOURCE_DIR}/LSODA.cpp) + +find_package(Threads REQUIRED) + +add_executable(test_lsoda test_LSODA.cpp) +target_link_libraries(test_lsoda lsoda) + +add_executable(benchmark_lsoda benchmark_LSODA.cpp) +target_link_libraries(benchmark_lsoda lsoda) +target_link_libraries(benchmark_lsoda ${CMAKE_THREAD_LIBS_INIT}) + +enable_testing() +add_test(NAME test_lsoda COMMAND $) +add_test(NAME test_benchmark COMMAND $) diff --git a/external/libsoda/LICENSE b/external/libsoda/LICENSE new file mode 100644 index 0000000000..2a5304883f --- /dev/null +++ b/external/libsoda/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2018 Dilawar Singh + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/external/libsoda/LSODA.cpp b/external/libsoda/LSODA.cpp new file mode 100644 index 0000000000..6e932a557d --- /dev/null +++ b/external/libsoda/LSODA.cpp @@ -0,0 +1,2268 @@ +/* + * HISTORY: + * This is a CPP version of the LSODA library for integration into MOOSE + somulator. + * The original was aquired from + * http://www.ccl.net/cca/software/SOURCES/C/kinetics2/index.shtml and modified + by + * Heng Li . Heng merged several C files into one and added a + * simpler interface. [Available + here](http://lh3lh3.users.sourceforge.net/download/lsoda.c) + + * The original source code came with no license or copyright + * information. Heng Li released his modification under the MIT/X11 license. I + * maintain the same license. I have removed quite a lot of text/comments from + * this library. Please refer to the standard documentation. + * + * Contact: Dilawar Singh +*/ + +#include "LSODA.h" +#include +#include +#include +#include +#include +#include +#include + +#include "helper.h" + +using namespace std; + +#define ETA 2.2204460492503131e-16 + +LSODA::LSODA() +{ + // Initialize arrays. + mord = {{12, 5}}; + sm1 = {{ + 0., 0.5, 0.575, 0.55, 0.45, 0.35, 0.25, 0.2, 0.15, 0.1, 0.075, 0.05, + 0.025 + } + }; + el = {{0}}; + cm1 = {{0}}; + cm2 = {{0}}; +} + +LSODA::~LSODA() {} + +bool LSODA::abs_compare(double a, double b) +{ + return (std::abs(a) < std::abs(b)); +} + +/* Purpose : Find largest component of double vector dx */ +size_t LSODA::idamax1( const vector& dx, const size_t n, const size_t offset=0) +{ + + size_t v = 0, vmax = 0; + size_t idmax = 1; + for( size_t i = 1; i <= n; i++) + { + v = abs(dx[i+offset]); + if( v > vmax ) + { + vmax = v; + idmax = i; + } + } + return idmax; + + // Following has failed with seg-fault. Probably issue with STL. + // return std::max_element( dx.begin()+1+offset, dx.begin()+1+n, LSODA::abs_compare) - dx.begin() - offset; +} + +/* Purpose : scalar vector multiplication + dx = da * dx +*/ +void LSODA::dscal1(const double da, vector &dx, const size_t n, + const size_t offset = 0) +{ + std::transform( dx.begin()+1+offset, dx.end(), dx.begin()+1+offset, + [&da](double x) -> double { return da*x; } ); +} + +/* Purpose : Inner product dx . dy */ +double LSODA::ddot1(const vector &a, const vector &b, + const size_t n, const size_t offsetA = 0, + const size_t offsetB = 0) +{ + double sum = 0.0; + for (size_t i = 1; i <= n; i++) + sum += a[i+offsetA] * b[i+offsetB]; + return sum; +} + +void LSODA::daxpy1(const double da, const vector &dx, + vector &dy, const size_t n, const size_t offsetX = 0, + const size_t offsetY = 0) +{ + + for (size_t i = 1; i <= n; i++) + dy[i + offsetY] = da * dx[i + offsetX] + dy[i + offsetY]; +} + +// See BLAS documentation. The first argument has been changed to vector. +void LSODA::dgesl(const vector> &a, const size_t n, + vector &ipvt, vector &b, const size_t job) +{ + size_t k, j; + double t; + + /* + Job = 0, solve a * x = b. + */ + if (job == 0) + { + /* + First solve L * y = b. + */ + for (k = 1; k <= n; k++) + { + t = ddot1(a[k], b, k - 1); + b[k] = (b[k] - t) / a[k][k]; + } + /* + Now solve U * x = y. + */ + for (k = n - 1; k >= 1; k--) + { + b[k] = b[k] + ddot1(a[k], b, n - k, k, k); + j = ipvt[k]; + if (j != k) + { + t = b[j]; + b[j] = b[k]; + b[k] = t; + } + } + return; + } + /* + Job = nonzero, solve Transpose(a) * x = b. + + First solve Transpose(U) * y = b. + */ + for (k = 1; k <= n - 1; k++) + { + j = ipvt[k]; + t = b[j]; + if (j != k) + { + b[j] = b[k]; + b[k] = t; + } + daxpy1(t, a[k], b, n - k, k, k); + } + /* + Now solve Transpose(L) * x = y. + */ + for (k = n; k >= 1; k--) + { + b[k] = b[k] / a[k][k]; + t = -b[k]; + daxpy1(t, a[k], b, k - 1); + } +} + +// See BLAS documentation. All double* has been changed to std::vector . +void LSODA::dgefa(vector> &a, const size_t n, vector &ipvt, + size_t *const info) +{ + size_t j = 0, k = 0, i = 0; + double t = 0.0; + + /* Gaussian elimination with partial pivoting. */ + + *info = 0; + for (k = 1; k <= n - 1; k++) + { + /* + Find j = pivot index. Note that a[k]+k-1 is the address of + the 0-th element of the row vector whose 1st element is a[k][k]. + */ + j = idamax1(a[k], n - k + 1, k - 1) + k - 1; + ipvt[k] = j; + /* + Zero pivot implies this row already triangularized. + */ + if (a[k][j] == 0.) + { + *info = k; + continue; + } + /* + Interchange if necessary. + */ + if (j != k) + { + t = a[k][j]; + a[k][j] = a[k][k]; + a[k][k] = t; + } + /* + Compute multipliers. + */ + t = -1. / a[k][k]; + dscal1(t, a[k], n - k, k); + + /* + Column elimination with row indexing. + */ + for (i = k + 1; i <= n; i++) + { + t = a[i][j]; + if (j != k) + { + a[i][j] = a[i][k]; + a[i][k] = t; + } + daxpy1(t, a[k], a[i], n - k, k, k); + } + } /* end k-loop */ + + ipvt[n] = n; + if (a[n][n] == 0.) + *info = n; +} + +/* Terminate lsoda due to illegal input. */ +void LSODA::terminate(int *istate) +{ + if (illin == 5) + cerr << "[lsoda] repeated occurrence of illegal input. run aborted.. " + "apparent infinite loop." + << endl; + else + { + illin++; + *istate = -3; + } +} + +/* Terminate lsoda due to various error conditions. */ +void LSODA::terminate2(vector &y, double *t) +{ + for (size_t i = 1; i <= n; i++) + y[i] = yh_[1][i]; + *t = tn_; + illin = 0; + return; +} + +/* + The following block handles all successful returns from lsoda. + If itask != 1, y is loaded from yh_ and t is set accordingly. + *Istate is set to 2, the illegal input counter is zeroed, and the + optional outputs are loaded into the work arrays before returning. +*/ + +void LSODA::successreturn(vector &y, double *t, int itask, int ihit, + double tcrit, int *istate) +{ + for (size_t i = 1; i <= n; i++) + y[i] = yh_[1][i]; + *t = tn_; + if (itask == 4 || itask == 5) + if (ihit) + *t = tcrit; + *istate = 2; + illin = 0; +} + +/* +c references.. +c 1. alan c. hindmarsh, odepack, a systematized collection of ode +c solvers, in scientific computing, r. s. stepleman et al. (eds.), +c north-holland, amsterdam, 1983, pp. 55-64. +c 2. linda r. petzold, automatic selection of methods for solving +c stiff and nonstiff systems of ordinary differential equations, +c siam j. sci. stat. comput. 4 (1983), pp. 136-148. +c----------------------------------------------------------------------- +*/ +void LSODA::lsoda(LSODA_ODE_SYSTEM_TYPE f, const size_t neq, vector &y, + double *t, double tout, int itask, int *istate, int iopt, + int jt, array &iworks, array &rworks, + void *_data) +{ + assert(tout > *t); + + int mxstp0 = 500, mxhnl0 = 10; + + int iflag = 0, lenyh = 0, ihit = 0; + + double atoli = 0, ayi = 0, big = 0, h0 = 0, hmax = 0, hmx = 0, rh = 0, + rtoli = 0, tcrit = 0, tdist = 0, tnext = 0, tol = 0, tolsf = 0, tp = 0, + size = 0, sum = 0, w0 = 0; + + /* + Block a. + This code block is executed on every call. + It tests *istate and itask for legality and branches appropriately. + If *istate > 1 but the flag init shows that initialization has not + yet been done, an error return occurs. + If *istate = 1 and tout = t, return immediately. + */ + + if (*istate < 1 || *istate > 3) + { + // fprintf(stderr, "[lsoda] illegal istate = %d\n", *istate); + cerr << "[lsoda] illegal istate = " << *istate << endl; + terminate(istate); + return; + } + if (itask < 1 || itask > 5) + { + fprintf(stderr, "[lsoda] illegal itask = %d\n", itask); + terminate(istate); + return; + } + if (init == 0 && (*istate == 2 || *istate == 3)) + { + fprintf(stderr, "[lsoda] istate > 1 but lsoda not initialized\n"); + terminate(istate); + return; + } + + /* + Block b. + The next code block is executed for the initial call ( *istate = 1 ), + or for a continuation call with parameter changes ( *istate = 3 ). + It contains checking of all inputs and various initializations. + + First check legality of the non-optional inputs neq, itol, iopt, + jt, ml, and mu. + */ + + if (*istate == 1 || *istate == 3) + { + ntrep = 0; + if (neq <= 0) + { + cerr << "[lsoda] neq = " << neq << " is less than 1." << endl; + terminate(istate); + return; + } + if (*istate == 3 && neq > n) + { + cerr << "[lsoda] istate = 3 and neq increased" << endl; + terminate(istate); + return; + } + n = neq; + if (itol_ < 1 || itol_ > 4) + { + cerr << "[lsoda] itol = " << itol_ << " illegal" << endl; + terminate(istate); + return; + } + if (iopt < 0 || iopt > 1) + { + cerr << "[lsoda] iopt = " << iopt << " illegal" << endl; + terminate(istate); + return; + } + if (jt == 3 || jt < 1 || jt > 5) + { + cerr << "[lsoda] jt = " << jt << " illegal" << endl; + terminate(istate); + return; + } + jtyp = jt; + if (jt > 2) + { + ml = iworks[0]; + mu = iworks[1]; + if (ml >= n) + { + cerr << "[lsoda] ml = " << ml << " not between 1 and neq" << endl; + terminate(istate); + return; + } + if (mu >= n) + { + cerr << "[lsoda] mu = " << mu << " not between 1 and neq" << endl; + terminate(istate); + return; + } + } + + /* Next process and check the optional inpus. */ + /* Default options. */ + if (iopt == 0) + { + ixpr = 0; + mxstep = mxstp0; + mxhnil = mxhnl0; + hmxi = 0.; + hmin = 0.; + if (*istate == 1) + { + h0 = 0.; + mxordn = mord[0]; + mxords = mord[1]; + } + } + /* end if ( iopt == 0 ) */ + /* Optional inputs. */ + else /* if ( iopt = 1 ) */ + { + ixpr = iworks[2]; + if (ixpr > 1) + { + cerr << "[lsoda] ixpr = " << ixpr << " is illegal" << endl; + terminate(istate); + return; + } + + mxstep = iworks[3]; + if (mxstep == 0) + mxstep = mxstp0; + mxhnil = iworks[4]; + + if (*istate == 1) + { + h0 = rworks[1]; + mxordn = iworks[5]; + + if (mxordn == 0) + mxordn = 100; + + mxordn = min(mxordn, mord[0]); + mxords = iworks[6]; + + // if mxords is not given use 100. + if (mxords == 0) + mxords = 100; + + mxords = min(mxords, mord[1]); + + if ((tout - *t) * h0 < 0.) + { + cerr << "[lsoda] tout = " << tout << " behind t = " << *t + << ". integration direction is given by " << h0 << endl; + terminate(istate); + return; + } + } /* end if ( *istate == 1 ) */ + hmax = rworks[2]; + if (hmax < 0.) + { + cerr << "[lsoda] hmax < 0." << endl; + terminate(istate); + return; + } + hmxi = 0.; + if (hmax > 0) + hmxi = 1. / hmax; + + hmin = rworks[3]; + if (hmin < 0.) + { + cerr << "[lsoda] hmin < 0." << endl; + terminate(istate); + return; + } + } /* end else */ /* end iopt = 1 */ + } /* end if ( *istate == 1 || *istate == 3 ) */ + /* + If *istate = 1, meth_ is initialized to 1. + + Also allocate memory for yh_, wm_, ewt, savf, acor, ipvt. + */ + if (*istate == 1) + { + /* + If memory were not freed, *istate = 3 need not reallocate memory. + Hence this section is not executed by *istate = 3. + */ + sqrteta = sqrt(ETA); + meth_ = 1; + + nyh = n; + lenyh = 1 + max(mxordn, mxords); + + yh_.resize(lenyh + 1, std::vector(nyh + 1, 0.0)); + wm_.resize(nyh + 1, std::vector(nyh + 1, 0.0)); + ewt.resize(1 + nyh, 0); + savf.resize(1 + nyh, 0); + acor.resize(nyh + 1, 0.0); + ipvt.resize(nyh + 1, 0.0); + } + /* + Check rtol and atol for legality. + */ + if (*istate == 1 || *istate == 3) + { + rtoli = rtol_[1]; + atoli = atol_[1]; + for (size_t i = 1; i <= n; i++) + { + if (itol_ >= 3) + rtoli = rtol_[i]; + if (itol_ == 2 || itol_ == 4) + atoli = atol_[i]; + if (rtoli < 0.) + { + fprintf(stderr, "[lsoda] rtol = %g is less than 0.\n", rtoli); + terminate(istate); + return; + } + if (atoli < 0.) + { + fprintf(stderr, "[lsoda] atol = %g is less than 0.\n", atoli); + terminate(istate); + return; + } + } /* end for */ + } /* end if ( *istate == 1 || *istate == 3 ) */ + + /* If *istate = 3, set flag to signal parameter changes to stoda. */ + if (*istate == 3) + { + jstart = -1; + } + /* + Block c. + The next block is for the initial call only ( *istate = 1 ). + It contains all remaining initializations, the initial call to f, + and the calculation of the initial step size. + The error weights in ewt are inverted after being loaded. + */ + if (*istate == 1) + { + tn_ = *t; + tsw = *t; + maxord = mxordn; + if (itask == 4 || itask == 5) + { + tcrit = rworks[0]; + if ((tcrit - tout) * (tout - *t) < 0.) + { + fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tout\n"); + terminate(istate); + return; + } + if (h0 != 0. && (*t + h0 - tcrit) * h0 > 0.) + h0 = tcrit - *t; + } + + jstart = 0; + nhnil = 0; + nst = 0; + nje = 0; + nslast = 0; + hu = 0.; + nqu = 0; + mused = 0; + miter = 0; + ccmax = 0.3; + maxcor = 3; + msbp = 20; + mxncf = 10; + + /* Initial call to f. */ + assert(yh_.size() == lenyh + 1); + assert(yh_[0].size() == nyh + 1); + + (*f)(*t, &y[1], &yh_[2][1], _data); + nfe = 1; + + /* Load the initial value vector in yh_. */ + for (size_t i = 1; i <= n; i++) + yh_[1][i] = y[i]; + + /* Load and invert the ewt array. ( h_ is temporarily set to 1. ) */ + nq = 1; + h_ = 1.; + ewset(y); + for (size_t i = 1; i <= n; i++) + { + if (ewt[i] <= 0.) + { + cerr << "[lsoda] ewt[" << i << "] = " << ewt[i] << " <= 0.\n" << endl; + terminate2(y, t); + return; + } + ewt[i] = 1. / ewt[i]; + } + + /* + The coding below computes the step size, h0, to be attempted on the + first step, unless the user has supplied a value for this. + First check that tout - *t differs significantly from zero. + A scalar tolerance quantity tol is computed, as max(rtol[i]) + if this is positive, or max(atol[i]/fabs(y[i])) otherwise, adjusted + so as to be between 100*ETA and 0.001. + Then the computed value h0 is given by + + h0^(-2) = 1. / ( tol * w0^2 ) + tol * ( norm(f) )^2 + + where w0 = max( fabs(*t), fabs(tout) ), + f = the initial value of the vector f(t,y), and + norm() = the weighted vector norm used throughout, given by + the vmnorm function routine, and weighted by the + tolerances initially loaded into the ewt array. + + The sign of h0 is inferred from the initial values of tout and *t. + fabs(h0) is made < fabs(tout-*t) in any case. + */ + if (h0 == 0.) + { + tdist = fabs(tout - *t); + w0 = max(fabs(*t), fabs(tout)); + if (tdist < 2. * ETA * w0) + { + fprintf(stderr, "[lsoda] tout too close to t to start integration\n "); + terminate(istate); + return; + } + tol = rtol_[1]; + if (itol_ > 2) + { + for (size_t i = 2; i <= n; i++) + tol = max(tol, rtol_[i]); + } + if (tol <= 0.) + { + atoli = atol_[1]; + for (size_t i = 1; i <= n; i++) + { + if (itol_ == 2 || itol_ == 4) + atoli = atol_[i]; + ayi = fabs(y[i]); + if (ayi != 0.) + tol = max(tol, atoli / ayi); + } + } + tol = max(tol, 100. * ETA); + tol = min(tol, 0.001); + sum = vmnorm(n, yh_[2], ewt); + sum = 1. / (tol * w0 * w0) + tol * sum * sum; + h0 = 1. / sqrt(sum); + h0 = min(h0, tdist); + h0 = h0 * ((tout - *t >= 0.) ? 1. : -1.); + } /* end if ( h0 == 0. ) */ + /* + Adjust h0 if necessary to meet hmax bound. + */ + rh = fabs(h0) * hmxi; + if (rh > 1.) + h0 /= rh; + + /* + Load h_ with h0 and scale yh_[2] by h0. + */ + h_ = h0; + for (size_t i = 1; i <= n; i++) + yh_[2][i] *= h0; + } /* if ( *istate == 1 ) */ + /* + Block d. + The next code block is for continuation calls only ( *istate = 2 or 3 ) + and is to check stop conditions before taking a step. + */ + if (*istate == 2 || *istate == 3) + { + nslast = nst; + switch (itask) + { + case 1: + if ((tn_ - tout) * h_ >= 0.) + { + intdy(tout, 0, y, &iflag); + if (iflag != 0) + { + fprintf(stderr, "[lsoda] trouble from intdy, itask = %d, tout = %g\n", + itask, tout); + terminate(istate); + return; + } + *t = tout; + *istate = 2; + illin = 0; + return; + } + break; + case 2: + break; + case 3: + tp = tn_ - hu * (1. + 100. * ETA); + if ((tp - tout) * h_ > 0.) + { + fprintf(stderr, "[lsoda] itask = %d and tout behind tcur - hu\n", + itask); + terminate(istate); + return; + } + if ((tn_ - tout) * h_ < 0.) + break; + successreturn(y, t, itask, ihit, tcrit, istate); + return; + case 4: + tcrit = rworks[0]; + if ((tn_ - tcrit) * h_ > 0.) + { + fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tcur\n"); + terminate(istate); + return; + } + if ((tcrit - tout) * h_ < 0.) + { + fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tout\n"); + terminate(istate); + return; + } + if ((tn_ - tout) * h_ >= 0.) + { + intdy(tout, 0, y, &iflag); + if (iflag != 0) + { + fprintf(stderr, "[lsoda] trouble from intdy, itask = %d, tout = %g\n", + itask, tout); + terminate(istate); + return; + } + *t = tout; + *istate = 2; + illin = 0; + return; + } + case 5: + if (itask == 5) + { + tcrit = rworks[0]; + if ((tn_ - tcrit) * h_ > 0.) + { + fprintf(stderr, "[lsoda] itask = 4 or 5 and tcrit behind tcur\n"); + terminate(istate); + return; + } + } + hmx = fabs(tn_) + fabs(h_); + ihit = fabs(tn_ - tcrit) <= (100. * ETA * hmx); + if (ihit) + { + *t = tcrit; + successreturn(y, t, itask, ihit, tcrit, istate); + return; + } + tnext = tn_ + h_ * (1. + 4. * ETA); + if ((tnext - tcrit) * h_ <= 0.) + break; + h_ = (tcrit - tn_) * (1. - 4. * ETA); + if (*istate == 2) + jstart = -2; + break; + } /* end switch */ + } /* end if ( *istate == 2 || *istate == 3 ) */ + /* + Block e. + The next block is normally executed for all calls and contains + the call to the one-step core integrator stoda. + + This is a looping point for the integration steps. + + First check for too many steps being taken, update ewt ( if not at + start of problem). Check for too much accuracy being requested, and + check for h_ below the roundoff level in *t. + */ + while (1) + { + if (*istate != 1 || nst != 0) + { + if ((nst - nslast) >= mxstep) + { + cerr << "[lsoda] " << mxstep << " steps taken before reaching tout" + << endl; + *istate = -1; + terminate2(y, t); + return; + } + + ewset(yh_[1]); + for (size_t i = 1; i <= n; i++) + { + if (ewt[i] <= 0.) + { + cerr << "[lsoda] ewt[" << i << "] = " << ewt[i] << " <= 0." << endl; + *istate = -6; + terminate2(y, t); + return; + } + ewt[i] = 1. / ewt[i]; + } + } + tolsf = ETA * vmnorm(n, yh_[1], ewt); + if (tolsf > 0.01) + { + tolsf = tolsf * 200.; + if (nst == 0) + { + fprintf(stderr, "lsoda -- at start of problem, too much accuracy\n"); + fprintf(stderr, " requested for precision of machine,\n"); + fprintf(stderr, " suggested scaling factor = %g\n", tolsf); + terminate(istate); + return; + } + fprintf(stderr, "lsoda -- at t = %g, too much accuracy requested\n", *t); + fprintf(stderr, " for precision of machine, suggested\n"); + fprintf(stderr, " scaling factor = %g\n", tolsf); + *istate = -2; + terminate2(y, t); + return; + } + + if ((tn_ + h_) == tn_) + { + nhnil++; + if (nhnil <= mxhnil) + { + fprintf(stderr, "lsoda -- warning..internal t = %g and h_ = %g are\n", + tn_, h_); + fprintf( + stderr, + " such that in the machine, t + h_ = t on the next step\n"); + fprintf(stderr, " solver will continue anyway.\n"); + if (nhnil == mxhnil) + { + cerr << "lsoda -- above warning has been issued " << nhnil + << " times, " << endl + << " it will not be issued again for this problem" << endl; + } + } + } + + /* Call stoda */ + stoda(neq, y, f, _data); + if (kflag == 0) + { + /* + Block f. + The following block handles the case of a successful return from the + core integrator ( kflag = 0 ). + If a method switch was just made, record tsw, reset maxord, + set jstart to -1 to signal stoda to complete the switch, + and do extra printing of data if ixpr = 1. + Then, in any case, check for stop conditions. + */ + init = 1; + if (meth_ != mused) + { + tsw = tn_; + maxord = mxordn; + if (meth_ == 2) + maxord = mxords; + jstart = -1; + if (ixpr) + { + if (meth_ == 2) + cerr << "[lsoda] a switch to the stiff method has occurred " + << endl; + if (meth_ == 1) + cerr << "[lsoda] a switch to the nonstiff method has occurred" + << endl; + } + } /* end if ( meth_ != mused ) */ + /* + itask = 1. + If tout has been reached, interpolate. + */ + if (1 == itask) + { + if ((tn_ - tout) * h_ < 0.) + continue; + + intdy(tout, 0, y, &iflag); + *t = tout; + *istate = 2; + illin = 0; + return; + } + /* + itask = 2. + */ + if (itask == 2) + { + successreturn(y, t, itask, ihit, tcrit, istate); + return; + } + /* + itask = 3. + Jump to exit if tout was reached. + */ + if (itask == 3) + { + if ((tn_ - tout) * h_ >= 0.) + { + successreturn(y, t, itask, ihit, tcrit, istate); + return; + } + continue; + } + /* + itask = 4. + See if tout or tcrit was reached. Adjust h_ if necessary. + */ + if (itask == 4) + { + if ((tn_ - tout) * h_ >= 0.) + { + intdy(tout, 0, y, &iflag); + *t = tout; + *istate = 2; + illin = 0; + return; + } + else + { + hmx = fabs(tn_) + fabs(h_); + ihit = fabs(tn_ - tcrit) <= (100. * ETA * hmx); + if (ihit) + { + successreturn(y, t, itask, ihit, tcrit, istate); + return; + } + tnext = tn_ + h_ * (1. + 4. * ETA); + if ((tnext - tcrit) * h_ <= 0.) + continue; + h_ = (tcrit - tn_) * (1. - 4. * ETA); + jstart = -2; + continue; + } + } /* end if ( itask == 4 ) */ + /* + itask = 5. + See if tcrit was reached and jump to exit. + */ + if (itask == 5) + { + hmx = fabs(tn_) + fabs(h_); + ihit = fabs(tn_ - tcrit) <= (100. * ETA * hmx); + successreturn(y, t, itask, ihit, tcrit, istate); + return; + } + } /* end if ( kflag == 0 ) */ + /* + kflag = -1, error test failed repeatedly or with fabs(h_) = hmin. + kflag = -2, convergence failed repeatedly or with fabs(h_) = hmin. + */ + if (kflag == -1 || kflag == -2) + { + fprintf(stderr, "lsoda -- at t = %g and step size h_ = %g, the\n", tn_, + h_); + if (kflag == -1) + { + fprintf(stderr, " error test failed repeatedly or\n"); + fprintf(stderr, " with fabs(h_) = hmin\n"); + *istate = -4; + } + if (kflag == -2) + { + fprintf(stderr, + " corrector convergence failed repeatedly or\n"); + fprintf(stderr, " with fabs(h_) = hmin\n"); + *istate = -5; + } + big = 0.; + imxer = 1; + for (size_t i = 1; i <= n; i++) + { + size = fabs(acor[i]) * ewt[i]; + if (big < size) + { + big = size; + imxer = i; + } + } + terminate2(y, t); + return; + } /* end if ( kflag == -1 || kflag == -2 ) */ + } /* end while */ +} /* end lsoda */ + +void LSODA::stoda(const size_t neq, vector &y, LSODA_ODE_SYSTEM_TYPE f, + void *_data) +{ + assert(neq + 1 == y.size()); + + size_t corflag = 0, orderflag = 0; + size_t i = 0, i1 = 0, j = 0, m = 0, ncf = 0; + double del = 0.0, delp = 0.0, dsm = 0.0, dup = 0.0, exup = 0.0, r = 0.0, + rh = 0.0, rhup = 0.0, told = 0.0; + double pdh = 0.0, pnorm = 0.0; + + /* + stoda performs one step of the integration of an initial value + problem for a system of ordinary differential equations. + Note.. stoda is independent of the value of the iteration method + indicator miter, when this is != 0, and hence is independent + of the type of chord method used, or the Jacobian structure. + Communication with stoda is done with the following variables: + + jstart = an integer used for input only, with the following + values and meanings: + + 0 perform the first step, + > 0 take a new step continuing from the last, + -1 take the next step with a new value of h_, + n, meth_, miter, and/or matrix parameters. + -2 take the next step with a new value of h_, + but with other inputs unchanged. + + kflag = a completion code with the following meanings: + + 0 the step was successful, + -1 the requested error could not be achieved, + -2 corrector convergence could not be achieved, + -3 fatal error in prja or solsy. + + miter = corrector iteration method: + + 0 functional iteration, + >0 a chord method corresponding to jacobian type jt. + + */ + kflag = 0; + told = tn_; + ncf = 0; + ierpj = 0; + iersl = 0; + jcur = 0; + delp = 0.; + + /* + On the first call, the order is set to 1, and other variables are + initialized. rmax is the maximum ratio by which h_ can be increased + in a single step. It is initially 1.e4 to compensate for the small + initial h_, but then is normally equal to 10. If a filure occurs + (in corrector convergence or error test), rmax is set at 2 for + the next increase. + cfode is called to get the needed coefficients for both methods. + */ + if (jstart == 0) + { + lmax = maxord + 1; + nq = 1; + l = 2; + ialth = 2; + rmax = 10000.; + rc = 0.; + el0 = 1.; + crate = 0.7; + hold = h_; + nslp = 0; + ipup = miter; + /* + Initialize switching parameters. meth_ = 1 is assumed initially. + */ + icount = 20; + irflag = 0; + pdest = 0.; + pdlast = 0.; + ratio = 5.; + cfode(2); + for (i = 1; i <= 5; i++) + cm2[i] = tesco[i][2] * elco[i][i + 1]; + cfode(1); + for (i = 1; i <= 12; i++) + cm1[i] = tesco[i][2] * elco[i][i + 1]; + resetcoeff(); + } /* end if ( jstart == 0 ) */ + /* + The following block handles preliminaries needed when jstart = -1. + ipup is set to miter to force a matrix update. + If an order increase is about to be considered ( ialth = 1 ), + ialth is reset to 2 to postpone consideration one more step. + If the caller has changed meth_, cfode is called to reset + the coefficients of the method. + If h_ is to be changed, yh_ must be rescaled. + If h_ or meth_ is being changed, ialth is reset to l = nq + 1 + to prevent further changes in h_ for that many steps. + */ + if (jstart == -1) + { + ipup = miter; + lmax = maxord + 1; + if (ialth == 1) + ialth = 2; + if (meth_ != mused) + { + cfode(meth_); + ialth = l; + resetcoeff(); + } + if (h_ != hold) + { + rh = h_ / hold; + h_ = hold; + scaleh(&rh, &pdh); + } + } /* if ( jstart == -1 ) */ + if (jstart == -2) + { + if (h_ != hold) + { + rh = h_ / hold; + h_ = hold; + scaleh(&rh, &pdh); + } + } /* if ( jstart == -2 ) */ + + /* + Prediction. + This section computes the predicted values by effectively + multiplying the yh_ array by the pascal triangle matrix. + rc is the ratio of new to old values of the coefficient h_ * el[1]. + When rc differs from 1 by more than ccmax, ipup is set to miter + to force pjac to be called, if a jacobian is involved. + In any case, prja is called at least every msbp steps. + */ + while (1) + { + while (1) + { + if (fabs(rc - 1.) > ccmax) + ipup = miter; + if (nst >= nslp + msbp) + ipup = miter; + tn_ += h_; + for (size_t j = nq; j >= 1; j--) + for (size_t i1 = j; i1 <= nq; i1++) + for (i = 1; i <= n; i++) + yh_[i1][i] += yh_[i1 + 1][i]; + + pnorm = vmnorm(n, yh_[1], ewt); + correction(neq, y, f, &corflag, pnorm, &del, &delp, &told, &ncf, &rh, &m, + _data); + if (corflag == 0) + break; + if (corflag == 1) + { + rh = max(rh, hmin / fabs(h_)); + scaleh(&rh, &pdh); + continue; + } + if (corflag == 2) + { + kflag = -2; + hold = h_; + jstart = 1; + return; + } + } /* end inner while ( corrector loop ) */ + + /* + The corrector has converged. jcur is set to 0 + to signal that the Jacobian involved may need updating later. + The local error test is done now. + */ + jcur = 0; + if (m == 0) + dsm = del / tesco[nq][2]; + if (m > 0) + dsm = vmnorm(n, acor, ewt) / tesco[nq][2]; + + if (dsm <= 1.) + { + /* + After a successful step, update the yh_ array. + Decrease icount by 1, and if it is -1, consider switching methods. + If a method switch is made, reset various parameters, + rescale the yh_ array, and exit. If there is no switch, + consider changing h_ if ialth = 1. Otherwise decrease ialth by 1. + If ialth is then 1 and nq < maxord, then acor is saved for + use in a possible order increase on the next step. + If a change in h_ is considered, an increase or decrease in order + by one is considered also. A change in h_ is made only if it is by + a factor of at least 1.1. If not, ialth is set to 3 to prevent + testing for that many steps. + */ + kflag = 0; + nst++; + hu = h_; + nqu = nq; + mused = meth_; + for (size_t j = 1; j <= l; j++) + { + r = el[j]; + for (i = 1; i <= n; i++) + yh_[j][i] += r * acor[i]; + } + icount--; + if (icount < 0) + { + methodswitch(dsm, pnorm, &pdh, &rh); + if (meth_ != mused) + { + rh = max(rh, hmin / fabs(h_)); + scaleh(&rh, &pdh); + rmax = 10.; + endstoda(); + break; + } + } + /* + No method switch is being made. Do the usual step/order selection. + */ + ialth--; + if (ialth == 0) + { + rhup = 0.; + if (l != lmax) + { + for (i = 1; i <= n; i++) + savf[i] = acor[i] - yh_[lmax][i]; + dup = vmnorm(n, savf, ewt) / tesco[nq][3]; + exup = 1. / (double)(l + 1); + rhup = 1. / (1.4 * pow(dup, exup) + 0.0000014); + } + + orderswitch(&rhup, dsm, &pdh, &rh, &orderflag); + + /* + No change in h_ or nq. + */ + if (orderflag == 0) + { + endstoda(); + break; + } + /* + h_ is changed, but not nq. + */ + if (orderflag == 1) + { + rh = max(rh, hmin / fabs(h_)); + scaleh(&rh, &pdh); + rmax = 10.; + endstoda(); + break; + } + /* + both nq and h_ are changed. + */ + if (orderflag == 2) + { + resetcoeff(); + rh = max(rh, hmin / fabs(h_)); + scaleh(&rh, &pdh); + rmax = 10.; + endstoda(); + break; + } + } /* end if ( ialth == 0 ) */ + if (ialth > 1 || l == lmax) + { + endstoda(); + break; + } + + for (size_t i = 1; i <= n; i++) + yh_[lmax][i] = acor[i]; + + endstoda(); + break; + } + /* end if ( dsm <= 1. ) */ + /* + The error test failed. kflag keeps track of multiple failures. + Restore tn_ and the yh_ array to their previous values, and prepare + to try the step again. Compute the optimum step size for this or + one lower. After 2 or more failures, h_ is forced to decrease + by a factor of 0.2 or less. + */ + else + { + kflag--; + tn_ = told; + for (j = nq; j >= 1; j--) + { + for (i1 = j; i1 <= nq; i1++) + for (i = 1; i <= n; i++) + yh_[i1][i] -= yh_[i1 + 1][i]; + } + rmax = 2.; + if (fabs(h_) <= hmin * 1.00001) + { + kflag = -1; + hold = h_; + jstart = 1; + break; + } + if (kflag > -3) + { + rhup = 0.; + orderswitch(&rhup, dsm, &pdh, &rh, &orderflag); + if (orderflag == 1 || orderflag == 0) + { + if (orderflag == 0) + rh = min(rh, 0.2); + rh = max(rh, hmin / fabs(h_)); + scaleh(&rh, &pdh); + } + if (orderflag == 2) + { + resetcoeff(); + rh = max(rh, hmin / fabs(h_)); + scaleh(&rh, &pdh); + } + continue; + } + /* if ( kflag > -3 ) */ + /* + Control reaches this section if 3 or more failures have occurred. + If 10 failures have occurred, exit with kflag = -1. + It is assumed that the derivatives that have accumulated in the + yh_ array have errors of the wrong order. Hence the first + derivative is recomputed, and the order is set to 1. Then + h_ is reduced by a factor of 10, and the step is retried, + until it succeeds or h_ reaches hmin. + */ + else + { + if (kflag == -10) + { + kflag = -1; + hold = h_; + jstart = 1; + break; + } + else + { + rh = 0.1; + rh = max(hmin / fabs(h_), rh); + h_ *= rh; + for (i = 1; i <= n; i++) + y[i] = yh_[1][i]; + (*f)(tn_, &y[1], &savf[1], _data); + nfe++; + for (i = 1; i <= n; i++) + yh_[2][i] = h_ * savf[i]; + ipup = miter; + ialth = 5; + if (nq == 1) + continue; + nq = 1; + l = 2; + resetcoeff(); + continue; + } + } /* end else -- kflag <= -3 */ + } /* end error failure handling */ + } /* end outer while */ + +} /* end stoda */ + +void LSODA::ewset(const vector &ycur) +{ + switch (itol_) + { + case 1: + for (size_t i = 1; i <= n; i++) + ewt[i] = rtol_[1] * fabs(ycur[i]) + atol_[1]; + break; + case 2: + for (size_t i = 1; i <= n; i++) + ewt[i] = rtol_[1] * fabs(ycur[i]) + atol_[i]; + break; + case 3: + for (size_t i = 1; i <= n; i++) + ewt[i] = rtol_[i] * fabs(ycur[i]) + atol_[1]; + break; + case 4: + for (size_t i = 1; i <= n; i++) + ewt[i] = rtol_[i] * fabs(ycur[i]) + atol_[i]; + break; + } + +} /* end ewset */ + +/* + Intdy computes interpolated values of the k-th derivative of the + dependent variable vector y, and stores it in dky. This routine + is called within the package with k = 0 and *t = tout, but may + also be called by the user for any k up to the current order. + ( See detailed instructions in the usage documentation. ) + + The computed values in dky are gotten by interpolation using the + Nordsieck history array yh_. This array corresponds uniquely to a + vector-valued polynomial of degree nqcur or less, and dky is set + to the k-th derivative of this polynomial at t. + The formula for dky is + + q + dky[i] = sum c[k][j] * ( t - tn_ )^(j-k) * h_^(-j) * yh_[j+1][i] + j=k + + where c[k][j] = j*(j-1)*...*(j-k+1), q = nqcur, tn_ = tcur, h_ = hcur. + The quantities nq = nqcur, l = nq+1, n = neq, tn_, and h_ are declared + static globally. The above sum is done in reverse order. + *iflag is returned negative if either k or t is out of bounds. +*/ +void LSODA::intdy(double t, int k, vector &dky, int *iflag) +{ + int ic, jp1 = 0; + double c, r, s, tp; + + *iflag = 0; + if (k < 0 || k > (int)nq) + { + fprintf(stderr, "[intdy] k = %d illegal\n", k); + *iflag = -1; + return; + } + tp = tn_ - hu - 100. * ETA * (tn_ + hu); + if ((t - tp) * (t - tn_) > 0.) + { + fprintf(stderr, + "intdy -- t = %g illegal. t not in interval tcur - hu to tcur\n", + t); + *iflag = -2; + return; + } + s = (t - tn_) / h_; + ic = 1; + for (size_t jj = l - k; jj <= nq; jj++) + ic *= jj; + c = (double)ic; + for (size_t i = 1; i <= n; i++) + dky[i] = c * yh_[l][i]; + + for (int j = nq - 1; j >= k; j--) + { + jp1 = j + 1; + ic = 1; + for (int jj = jp1 - k; jj <= j; jj++) + ic *= jj; + c = (double)ic; + + for (size_t i = 1; i <= n; i++) + dky[i] = c * yh_[jp1][i] + s * dky[i]; + } + if (k == 0) + return; + r = pow(h_, (double)(-k)); + + for (size_t i = 1; i <= n; i++) + dky[i] *= r; + +} /* end intdy */ + +void LSODA::cfode(int meth_) +{ + int i, nq, nqm1, nqp1; + double agamq, fnq, fnqm1, pc[13], pint, ragq, rqfac, rq1fac, tsign, xpin; + /* + cfode is called by the integrator routine to set coefficients + needed there. The coefficients for the current method, as + given by the value of meth_, are set for all orders and saved. + The maximum order assumed here is 12 if meth_ = 1 and 5 if meth_ = 2. + ( A smaller value of the maximum order is also allowed. ) + cfode is called once at the beginning of the problem, and + is not called again unless and until meth_ is changed. + + The elco array contains the basic method coefficients. + The coefficients el[i], 1 < i < nq+1, for the method of + order nq are stored in elco[nq][i]. They are given by a generating + polynomial, i.e., + + l(x) = el[1] + el[2]*x + ... + el[nq+1]*x^nq. + + For the implicit Adams method, l(x) is given by + + dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0. + + For the bdf methods, l(x) is given by + + l(x) = (x+1)*(x+2)*...*(x+nq)/k, + + where k = factorial(nq)*(1+1/2+...+1/nq). + + The tesco array contains test constants used for the + local error test and the selection of step size and/or order. + At order nq, tesco[nq][k] is used for the selection of step + size at order nq-1 if k = 1, at order nq if k = 2, and at order + nq+1 if k = 3. + */ + if (meth_ == 1) + { + elco[1][1] = 1.; + elco[1][2] = 1.; + tesco[1][1] = 0.; + tesco[1][2] = 2.; + tesco[2][1] = 1.; + tesco[12][3] = 0.; + pc[1] = 1.; + rqfac = 1.; + for (nq = 2; nq <= 12; nq++) + { + /* + The pc array will contain the coefficients of the polynomial + + p(x) = (x+1)*(x+2)*...*(x+nq-1). + + Initially, p(x) = 1. + */ + rq1fac = rqfac; + rqfac = rqfac / (double)nq; + nqm1 = nq - 1; + fnqm1 = (double)nqm1; + nqp1 = nq + 1; + /* + Form coefficients of p(x)*(x+nq-1). + */ + pc[nq] = 0.; + for (i = nq; i >= 2; i--) + pc[i] = pc[i - 1] + fnqm1 * pc[i]; + pc[1] = fnqm1 * pc[1]; + /* + Compute integral, -1 to 0, of p(x) and x*p(x). + */ + pint = pc[1]; + xpin = pc[1] / 2.; + tsign = 1.; + for (i = 2; i <= nq; i++) + { + tsign = -tsign; + pint += tsign * pc[i] / (double)i; + xpin += tsign * pc[i] / (double)(i + 1); + } + /* + Store coefficients in elco and tesco. + */ + elco[nq][1] = pint * rq1fac; + elco[nq][2] = 1.; + for (i = 2; i <= nq; i++) + elco[nq][i + 1] = rq1fac * pc[i] / (double)i; + agamq = rqfac * xpin; + ragq = 1. / agamq; + tesco[nq][2] = ragq; + if (nq < 12) + tesco[nqp1][1] = ragq * rqfac / (double)nqp1; + tesco[nqm1][3] = ragq; + } /* end for */ + return; + } /* end if ( meth_ == 1 ) */ + + /* meth_ = 2. */ + pc[1] = 1.; + rq1fac = 1.; + + /* + The pc array will contain the coefficients of the polynomial + p(x) = (x+1)*(x+2)*...*(x+nq). + Initially, p(x) = 1. + */ + for (nq = 1; nq <= 5; nq++) + { + fnq = (double)nq; + nqp1 = nq + 1; + /* + Form coefficients of p(x)*(x+nq). + */ + pc[nqp1] = 0.; + for (i = nq + 1; i >= 2; i--) + pc[i] = pc[i - 1] + fnq * pc[i]; + pc[1] *= fnq; + /* + Store coefficients in elco and tesco. + */ + for (i = 1; i <= nqp1; i++) + elco[nq][i] = pc[i] / pc[2]; + elco[nq][2] = 1.; + tesco[nq][1] = rq1fac; + tesco[nq][2] = ((double)nqp1) / elco[nq][1]; + tesco[nq][3] = ((double)(nq + 2)) / elco[nq][1]; + rq1fac /= fnq; + } + return; + +} /* end cfode */ + +void LSODA::scaleh(double *rh, double *pdh) +{ + double r; + /* + If h_ is being changed, the h_ ratio rh is checked against rmax, hmin, + and hmxi, and the yh_ array is rescaled. ialth is set to l = nq + 1 + to prevent a change of h_ for that many steps, unless forced by a + convergence or error test failure. + */ + *rh = min(*rh, rmax); + *rh = *rh / max(1., fabs(h_) * hmxi * *rh); + /* + If meth_ = 1, also restrict the new step size by the stability region. + If this reduces h_, set irflag to 1 so that if there are roundoff + problems later, we can assume that is the cause of the trouble. + */ + if (meth_ == 1) + { + irflag = 0; + *pdh = max(fabs(h_) * pdlast, 0.000001); + if ((*rh * *pdh * 1.00001) >= sm1[nq]) + { + *rh = sm1[nq] / *pdh; + irflag = 1; + } + } + r = 1.; + for (size_t j = 2; j <= l; j++) + { + r *= *rh; + for (size_t i = 1; i <= n; i++) + yh_[j][i] *= r; + } + h_ *= *rh; + rc *= *rh; + ialth = l; + +} /* end scaleh */ + +void LSODA::prja(const size_t neq, vector &y, LSODA_ODE_SYSTEM_TYPE f, + void *_data) +{ + size_t i = 0, ier = 0, j = 0; + double fac = 0.0, hl0 = 0.0, r = 0.0, r0 = 0.0, yj = 0.0; + /* + prja is called by stoda to compute and process the matrix + P = I - h_ * el[1] * J, where J is an approximation to the Jacobian. + Here J is computed by finite differencing. + J, scaled by -h_ * el[1], is stored in wm_. Then the norm of J ( the + matrix norm consistent with the weighted max-norm on vectors given + by vmnorm ) is computed, and J is overwritten by P. P is then + subjected to LU decomposition in preparation for later solution + of linear systems with p as coefficient matrix. This is done + by dgefa if miter = 2, and by dgbfa if miter = 5. + */ + nje++; + ierpj = 0; + jcur = 1; + hl0 = h_ * el0; + /* + If miter = 2, make n calls to f to approximate J. + */ + if (miter != 2) + { + fprintf(stderr, "[prja] miter != 2\n"); + return; + } + if (miter == 2) + { + fac = vmnorm(n, savf, ewt); + r0 = 1000. * fabs(h_) * ETA * ((double)n) * fac; + if (r0 == 0.) + r0 = 1.; + for (j = 1; j <= n; j++) + { + yj = y[j]; + r = max(sqrteta * fabs(yj), r0 / ewt[j]); + y[j] += r; + fac = -hl0 / r; + (*f)(tn_, &y[1], &acor[1], _data); + for (i = 1; i <= n; i++) + wm_[i][j] = (acor[i] - savf[i]) * fac; + y[j] = yj; + } + nfe += n; + /* + Compute norm of Jacobian. + */ + pdnorm = fnorm(n, wm_, ewt) / fabs(hl0); + /* + Add identity matrix. + */ + for (i = 1; i <= n; i++) + wm_[i][i] += 1.; + /* + Do LU decomposition on P. + */ + dgefa(wm_, n, ipvt, &ier); + if (ier != 0) + ierpj = 1; + return; + } +} /* end prja */ + +/* + This function routine computes the weighted max-norm + of the vector of length n contained in the array v, with weights + contained in the array w of length n. + + vmnorm = max( i = 1, ..., n ) fabs( v[i] ) * w[i]. +*/ +double LSODA::vmnorm(const size_t n, const vector &v, + const vector &w) +{ + double vm = 0.; + for (size_t i = 1; i <= n; i++) + vm = max(vm, fabs(v[i]) * w[i]); + return vm; +} + +double LSODA::fnorm(int n, const vector> &a, + const vector &w) + +/* + This subroutine computes the norm of a full n by n matrix, + stored in the array a, that is consistent with the weighted max-norm + on vectors, with weights stored in the array w. + + fnorm = max(i=1,...,n) ( w[i] * sum(j=1,...,n) fabs( a[i][j] ) / w[j] ) +*/ + +{ + double an = 0, sum = 0; + + for (size_t i = 1; i <= (size_t)n; i++) + { + sum = 0.; + for (size_t j = 1; j <= (size_t)n; j++) + sum += fabs(a[i][j]) / w[j]; + an = max(an, sum * w[i]); + } + return an; +} + +/* + *corflag = 0 : corrector converged, +1 : step size to be reduced, redo prediction, +2 : corrector cannot converge, failure flag. +*/ +void LSODA::correction(const size_t neq, vector &y, + LSODA_ODE_SYSTEM_TYPE f, size_t *corflag, double pnorm, + double *del, double *delp, double *told, size_t *ncf, + double *rh, size_t *m, void *_data) +{ + double rm = 0.0, rate = 0.0, dcon = 0.0; + + /* + Up to maxcor corrector iterations are taken. A convergence test is + made on the r.m.s. norm of each correction, weighted by the error + weight vector ewt. The sum of the corrections is accumulated in the + vector acor[i]. The yh_ array is not altered in the corrector loop. + */ + + *m = 0; + *corflag = 0; + *del = 0.; + + for (size_t i = 1; i <= n; i++) + y[i] = yh_[1][i]; + + (*f)(tn_, &y[1], &savf[1], _data); + + nfe++; + /* + If indicated, the matrix P = I - h_ * el[1] * J is reevaluated and + preprocessed before starting the corrector iteration. ipup is set + to 0 as an indicator that this has been done. + */ + while (1) + { + if (*m == 0) + { + if (ipup > 0) + { + prja(neq, y, f, _data); + ipup = 0; + rc = 1.; + nslp = nst; + crate = 0.7; + if (ierpj != 0) + { + corfailure(told, rh, ncf, corflag); + return; + } + } + for (size_t i = 1; i <= n; i++) + acor[i] = 0.; + } /* end if ( *m == 0 ) */ + if (miter == 0) + { + /* + In case of functional iteration, update y directly from + the result of the last function evaluation. + */ + for (size_t i = 1; i <= n; i++) + { + savf[i] = h_ * savf[i] - yh_[2][i]; + y[i] = savf[i] - acor[i]; + } + *del = vmnorm(n, y, ewt); + for (size_t i = 1; i <= n; i++) + { + y[i] = yh_[1][i] + el[1] * savf[i]; + acor[i] = savf[i]; + } + } + /* end functional iteration */ + /* + In the case of the chord method, compute the corrector error, + and solve the linear system with that as right-hand side and + P as coefficient matrix. + */ + else + { + for (size_t i = 1; i <= n; i++) + y[i] = h_ * savf[i] - (yh_[2][i] + acor[i]); + + solsy(y); + *del = vmnorm(n, y, ewt); + + for (size_t i = 1; i <= n; i++) + { + acor[i] += y[i]; + y[i] = yh_[1][i] + el[1] * acor[i]; + } + } /* end chord method */ + /* + Test for convergence. If *m > 0, an estimate of the convergence + rate constant is stored in crate, and this is used in the test. + + We first check for a change of iterates that is the size of + roundoff error. If this occurs, the iteration has converged, and a + new rate estimate is not formed. + In all other cases, force at least two iterations to estimate a + local Lipschitz constant estimate for Adams method. + On convergence, form pdest = local maximum Lipschitz constant + estimate. pdlast is the most recent nonzero estimate. + */ + if (*del <= 100. * pnorm * ETA) + break; + if (*m != 0 || meth_ != 1) + { + if (*m != 0) + { + rm = 1024.0; + if (*del <= (1024. * *delp)) + rm = *del / *delp; + rate = max(rate, rm); + crate = max(0.2 * crate, rm); + } + dcon = *del * min(1., 1.5 * crate) / (tesco[nq][2] * conit); + if (dcon <= 1.) + { + pdest = max(pdest, rate / fabs(h_ * el[1])); + if (pdest != 0.) + pdlast = pdest; + break; + } + } + /* + The corrector iteration failed to converge. + If miter != 0 and the Jacobian is out of date, prja is called for + the next try. Otherwise the yh_ array is retracted to its values + before prediction, and h_ is reduced, if possible. If h_ cannot be + reduced or mxncf failures have occured, exit with corflag = 2. + */ + (*m)++; + if (*m == maxcor || (*m >= 2 && *del > 2. * *delp)) + { + if (miter == 0 || jcur == 1) + { + corfailure(told, rh, ncf, corflag); + return; + } + ipup = miter; + /* + Restart corrector if Jacobian is recomputed. + */ + *m = 0; + rate = 0.; + *del = 0.; + for (size_t i = 1; i <= n; i++) + y[i] = yh_[1][i]; + + (*f)(tn_, &y[1], &savf[1], _data); + + nfe++; + } + /* + Iterate corrector. + */ + else + { + *delp = *del; + (*f)(tn_, &y[1], &savf[1], _data); + nfe++; + } + } /* end while */ +} /* end correction */ + +void LSODA::corfailure(double *told, double *rh, size_t *ncf, size_t *corflag) +{ + ncf++; + rmax = 2.; + tn_ = *told; + for (size_t j = nq; j >= 1; j--) + for (size_t i1 = j; i1 <= nq; i1++) + for (size_t i = 1; i <= n; i++) + yh_[i1][i] -= yh_[i1 + 1][i]; + + if (fabs(h_) <= hmin * 1.00001 || *ncf == mxncf) + { + *corflag = 2; + return; + } + *corflag = 1; + *rh = 0.25; + ipup = miter; +} + +/* + This routine manages the solution of the linear system arising from + a chord iteration. It is called if miter != 0. + If miter is 2, it calls dgesl to accomplish this. + If miter is 5, it calls dgbsl. + + y = the right-hand side vector on input, and the solution vector + on output. +*/ +void LSODA::solsy(vector &y) +{ + iersl = 0; + if (miter != 2) + { + printf("solsy -- miter != 2\n"); + return; + } + if (miter == 2) + dgesl(wm_, n, ipvt, y, 0); + return; +} + +void LSODA::methodswitch(double dsm, double pnorm, double *pdh, double *rh) +{ + int lm1, lm1p1, lm2, lm2p1, nqm1, nqm2; + double rh1, rh2, rh1it, exm2, dm2, exm1, dm1, alpha, exsm; + + /* + We are current using an Adams method. Consider switching to bdf. + If the current order is greater than 5, assume the problem is + not stiff, and skip this section. + If the Lipschitz constant and error estimate are not polluted + by roundoff, perform the usual test. + Otherwise, switch to the bdf methods if the last step was + restricted to insure stability ( irflag = 1 ), and stay with Adams + method if not. When switching to bdf with polluted error estimates, + in the absence of other information, double the step size. + + When the estimates are ok, we make the usual test by computing + the step size we could have (ideally) used on this step, + with the current (Adams) method, and also that for the bdf. + If nq > mxords, we consider changing to order mxords on switching. + Compare the two step sizes to decide whether to switch. + The step size advantage must be at least ratio = 5 to switch. + */ + if (meth_ == 1) + { + if (nq > 5) + return; + if (dsm <= (100. * pnorm * ETA) || pdest == 0.) + { + if (irflag == 0) + return; + rh2 = 2.; + nqm2 = min(nq, mxords); + } + else + { + exsm = 1. / (double)l; + rh1 = 1. / (1.2 * pow(dsm, exsm) + 0.0000012); + rh1it = 2. * rh1; + *pdh = pdlast * fabs(h_); + if ((*pdh * rh1) > 0.00001) + rh1it = sm1[nq] / *pdh; + rh1 = min(rh1, rh1it); + if (nq > mxords) + { + nqm2 = mxords; + lm2 = mxords + 1; + exm2 = 1. / (double)lm2; + lm2p1 = lm2 + 1; + dm2 = vmnorm(n, yh_[lm2p1], ewt) / cm2[mxords]; + rh2 = 1. / (1.2 * pow(dm2, exm2) + 0.0000012); + } + else + { + dm2 = dsm * (cm1[nq] / cm2[nq]); + rh2 = 1. / (1.2 * pow(dm2, exsm) + 0.0000012); + nqm2 = nq; + } + if (rh2 < ratio * rh1) + return; + } + /* + The method switch test passed. Reset relevant quantities for bdf. + */ + *rh = rh2; + icount = 20; + meth_ = 2; + miter = jtyp; + pdlast = 0.; + nq = nqm2; + l = nq + 1; + return; + } /* end if ( meth_ == 1 ) */ + + /* + We are currently using a bdf method, considering switching to Adams. + Compute the step size we could have (ideally) used on this step, + with the current (bdf) method, and also that for the Adams. + If nq > mxordn, we consider changing to order mxordn on switching. + Compare the two step sizes to decide whether to switch. + The step size advantage must be at least 5/ratio = 1 to switch. + If the step size for Adams would be so small as to cause + roundoff pollution, we stay with bdf. + */ + exsm = 1. / (double)l; + if (mxordn < nq) + { + nqm1 = mxordn; + lm1 = mxordn + 1; + exm1 = 1. / (double)lm1; + lm1p1 = lm1 + 1; + dm1 = vmnorm(n, yh_[lm1p1], ewt) / cm1[mxordn]; + rh1 = 1. / (1.2 * pow(dm1, exm1) + 0.0000012); + } + else + { + dm1 = dsm * (cm2[nq] / cm1[nq]); + rh1 = 1. / (1.2 * pow(dm1, exsm) + 0.0000012); + nqm1 = nq; + exm1 = exsm; + } + rh1it = 2. * rh1; + *pdh = pdnorm * fabs(h_); + if ((*pdh * rh1) > 0.00001) + rh1it = sm1[nqm1] / *pdh; + rh1 = min(rh1, rh1it); + rh2 = 1. / (1.2 * pow(dsm, exsm) + 0.0000012); + if ((rh1 * ratio) < (5. * rh2)) + return; + alpha = max(0.001, rh1); + dm1 *= pow(alpha, exm1); + if (dm1 <= 1000. * ETA * pnorm) + return; + /* + The switch test passed. Reset relevant quantities for Adams. + */ + *rh = rh1; + icount = 20; + meth_ = 1; + miter = 0; + pdlast = 0.; + nq = nqm1; + l = nq + 1; +} /* end methodswitch */ + +/* + This routine returns from stoda to lsoda. Hence freevectors() is + not executed. +*/ +void LSODA::endstoda() +{ + double r = 1. / tesco[nqu][2]; + for (size_t i = 1; i <= n; i++) + acor[i] *= r; + hold = h_; + jstart = 1; +} + +/* + Regardless of the success or failure of the step, factors + rhdn, rhsm, and rhup are computed, by which h_ could be multiplied + at order nq - 1, order nq, or order nq + 1, respectively. + In the case of a failure, rhup = 0. to avoid an order increase. + The largest of these is determined and the new order chosen + accordingly. If the order is to be increased, we compute one + additional scaled derivative. + + orderflag = 0 : no change in h_ or nq, + 1 : change in h_ but not nq, + 2 : change in both h_ and nq. +*/ +void LSODA::orderswitch(double *rhup, double dsm, double *pdh, double *rh, + size_t *orderflag) +{ + size_t newq = 0; + double exsm, rhdn, rhsm, ddn, exdn, r; + + *orderflag = 0; + + exsm = 1. / (double)l; + rhsm = 1. / (1.2 * pow(dsm, exsm) + 0.0000012); + + rhdn = 0.; + if (nq != 1) + { + ddn = vmnorm(n, yh_[l], ewt) / tesco[nq][1]; + exdn = 1. / (double)nq; + rhdn = 1. / (1.3 * pow(ddn, exdn) + 0.0000013); + } + /* + If meth_ = 1, limit rh accordinfg to the stability region also. + */ + if (meth_ == 1) + { + *pdh = max(fabs(h_) * pdlast, 0.000001); + if (l < lmax) + *rhup = min(*rhup, sm1[l] / *pdh); + rhsm = min(rhsm, sm1[nq] / *pdh); + if (nq > 1) + rhdn = min(rhdn, sm1[nq - 1] / *pdh); + pdest = 0.; + } + if (rhsm >= *rhup) + { + if (rhsm >= rhdn) + { + newq = nq; + *rh = rhsm; + } + else + { + newq = nq - 1; + *rh = rhdn; + if (kflag < 0 && *rh > 1.) + *rh = 1.; + } + } + else + { + if (*rhup <= rhdn) + { + newq = nq - 1; + *rh = rhdn; + if (kflag < 0 && *rh > 1.) + *rh = 1.; + } + else + { + *rh = *rhup; + if (*rh >= 1.1) + { + r = el[l] / (double)l; + nq = l; + l = nq + 1; + for (size_t i = 1; i <= n; i++) + yh_[l][i] = acor[i] * r; + + *orderflag = 2; + return; + } + else + { + ialth = 3; + return; + } + } + } + /* + If meth_ = 1 and h_ is restricted by stability, bypass 10 percent test. + */ + if (1 == meth_) + { + if ((*rh * *pdh * 1.00001) < sm1[newq]) + if (kflag == 0 && *rh < 1.1) + { + ialth = 3; + return; + } + } + else + { + if (kflag == 0 && *rh < 1.1) + { + ialth = 3; + return; + } + } + if (kflag <= -2) + *rh = min(*rh, 0.2); + /* + If there is a change of order, reset nq, l, and the coefficients. + In any case h_ is reset according to rh and the yh_ array is rescaled. + Then exit or redo the step. + */ + if (newq == nq) + { + *orderflag = 1; + return; + } + nq = newq; + l = nq + 1; + *orderflag = 2; + +} /* end orderswitch */ + +void LSODA::resetcoeff() +/* + The el vector and related constants are reset + whenever the order nq is changed, or at the start of the problem. +*/ +{ + array ep1; + + ep1 = elco[nq]; + for (size_t i = 1; i <= l; i++) + el[i] = ep1[i]; + rc = rc * el[1] / el0; + el0 = el[1]; + conit = 0.5 / (double)(nq + 2); +} + +void LSODA::_freevectors(void) +{ + // Does nothing. USE c++ memory mechanism here. +} + +/* --------------------------------------------------------------------------*/ +/** + * @Synopsis Simpler interface. + * + * @Param f System + * @Param neq, size of system. + * @Param y, init values of size neq + * @Param yout, results vector for size neq+1, ignore yout[0] + * @Param t, start time. + * @Param tout, stop time. + * @Param _data + * @Param rtol, relative tolerance. + * @Param atol, absolute tolerance. + */ +/* ----------------------------------------------------------------------------*/ +void LSODA::lsoda_update(LSODA_ODE_SYSTEM_TYPE f, const size_t neq, + vector &y, vector &yout, double *t, + const double tout, int *istate, void * _data, + double rtol, double atol + ) +{ + array iworks = {{0}}; + array rworks = {{0.0}}; + + int itask, iopt, jt; + + // cout << "Debug : rtol " << rtol << ". atol " << atol << endl; + + itask = 1; + iopt = 0; + jt = 2; + + yout.resize(neq + 1); + + // Set the tolerance. We should do it only once. + rtol_.resize(neq + 1, rtol); + atol_.resize(neq + 1, atol); + rtol_[0] = 0; + atol_[0] = 0; + + // Fill-in values. + for (size_t i = 1; i <= neq; i++) + yout[i] = y[i - 1]; + + lsoda(f, neq, yout, t, tout, itask, istate, iopt, jt, iworks, rworks, _data); +} diff --git a/external/libsoda/LSODA.h b/external/libsoda/LSODA.h new file mode 100644 index 0000000000..f8b4e13960 --- /dev/null +++ b/external/libsoda/LSODA.h @@ -0,0 +1,158 @@ +/*** + * Created: 2018-08-14 + + * Author: Dilawar Singh + * Organization: NCBS Bangalore + * License: MIT License + */ + +#ifndef LSODE_H +#define LSODE_H + +#include +#include +#include +#include + +using namespace std; + +/* --------------------------------------------------------------------------*/ +/** + * @Synopsis Type definition of LSODA ode system. See the file test_LSODA.cpp + * for an example. + * + * @Param time, double + * @Param y, array of double. + * @Param dydt, array of double + * @Param data, void* + * + * @Returns void + */ +/* ----------------------------------------------------------------------------*/ +typedef void (*LSODA_ODE_SYSTEM_TYPE)(double t, double *y, double *dydt, void *); + +class LSODA +{ + +public: + LSODA(); + ~LSODA(); + + size_t idamax1(const vector &dx, const size_t n, const size_t offset); + + void dscal1(const double da, vector &dx, const size_t n, + const size_t offset); + + double ddot1(const vector &a, const vector &b, const size_t n, + const size_t offsetA, const size_t offsetB); + + void daxpy1(const double da, const vector &dx, vector &dy, + const size_t n, const size_t offsetX, + const size_t offsetY); + + void dgesl(const vector> &a, const size_t n, vector &ipvt, + vector &b, const size_t job); + + void dgefa(vector> &a, const size_t n, vector &ipvt, + size_t *const info); + + void prja(const size_t neq, vector &y, LSODA_ODE_SYSTEM_TYPE f, + void *_data); + + void lsoda(LSODA_ODE_SYSTEM_TYPE f, const size_t neq, vector &y, + double *t, double tout, int itask, int *istate, int iopt, int jt, + array &iworks, array &rworks, void *_data); + + void correction(const size_t neq, vector &y, LSODA_ODE_SYSTEM_TYPE f, + size_t *corflag, double pnorm, double *del, double *delp, + double *told, size_t *ncf, double *rh, size_t *m, + void *_data); + + void stoda(const size_t neq, vector &y, LSODA_ODE_SYSTEM_TYPE f, + void *_data); + + // We call this function in VoxelPools:: + void lsoda_update(LSODA_ODE_SYSTEM_TYPE f, const size_t neq, + vector &y, std::vector &yout, double *t, + const double tout, int *istate, void *const _data, + double rtol = 1e-3, double atol = 1e-6 // Tolerance + ); + + void terminate(int *istate); + void terminate2(vector &y, double *t); + void successreturn(vector &y, double *t, int itask, int ihit, + double tcrit, int *istate); + void _freevectors(void); + void ewset(const vector &ycur); + void resetcoeff(void); + void solsy(vector &y); + void endstoda(void); + void orderswitch(double *rhup, double dsm, double *pdh, double *rh, + size_t *orderflag); + void intdy(double t, int k, vector &dky, int *iflag); + void corfailure(double *told, double *rh, size_t *ncf, size_t *corflag); + void methodswitch(double dsm, double pnorm, double *pdh, double *rh); + void cfode(int meth_); + void scaleh(double *rh, double *pdh); + double fnorm(int n, const vector> &a, const vector &w); + double vmnorm(const size_t n, const vector &v, + const vector &w); + + static bool abs_compare(double a, double b); + +private: + size_t ml, mu, imxer; + double sqrteta; + + // NOTE: initialize in default constructor. Older compiler e.g. 4.8.4 would + // produce error if these are initialized here. With newer compiler, + // initialization can be done here. + array mord; + array sm1; + + array el; // = {0}; + array cm1; // = {0}; + array cm2; // = {0}; + + array, 13> elco; + array, 13> tesco; + + size_t illin, init, ierpj, iersl, jcur, l, miter, maxord, maxcor, msbp, mxncf; + + int kflag, jstart; + + size_t ixpr = 0, jtyp, mused, mxordn, mxords = 12; + size_t meth_; + + size_t n, nq, nst, nfe, nje, nqu; + size_t mxstep, mxhnil; + size_t nslast, nhnil, ntrep, nyh; + + double ccmax, el0, h_ = .0; + double hmin, hmxi, hu, rc, tn_ = 0.0; + double tsw, pdnorm; + double conit, crate, hold, rmax; + + size_t ialth, ipup, lmax; + size_t nslp; + double pdest, pdlast, ratio; + int icount, irflag; + + vector ewt; + vector savf; + vector acor; + vector> yh_; + vector> wm_; + + vector ipvt; + +private: + int itol_ = 2; + std::vector rtol_; + std::vector atol_; + +public: + void *param = nullptr; +}; + +#endif /* end of include guard: LSODE_H */ diff --git a/external/libsoda/README.md b/external/libsoda/README.md new file mode 100644 index 0000000000..cd2020e6a7 --- /dev/null +++ b/external/libsoda/README.md @@ -0,0 +1,13 @@ +[![Build Status](https://travis-ci.org/dilawar/libsoda.svg?branch=master)](https://travis-ci.org/dilawar/libsoda) + +## libsoda++ + +This is `c++11` version of LSODA library. See the `./test_LSODA.cpp` file for example. + +## CREDIT + +This work is based on http://lh3lh3.users.sourceforge.net/download/lsoda.c . + +Related projects: + +1. https://github.com/sdwfrost/liblsoda diff --git a/external/libsoda/benchmark_LSODA.cpp b/external/libsoda/benchmark_LSODA.cpp new file mode 100644 index 0000000000..82692b63f3 --- /dev/null +++ b/external/libsoda/benchmark_LSODA.cpp @@ -0,0 +1,248 @@ +/*** + * Description: benchamark ODE system. + * + * Created: 2018-08-14 + + * Author: Dilawar Singh + * Organization: NCBS Bangalore + * License: MIT License + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "LSODA.h" +#include "helper.h" + +using namespace std; +using namespace std::chrono; + +// Describe the system. +static void fex(double t, double *y, double *ydot, void *data) +{ + ydot[0] = 1.0E4 * y[1] * y[2] - .04E0 * y[0]; + // Don't swap ydot[1] and ydot[2]. The order will change and test will fail. + ydot[2] = 3.0E7 * y[1] * y[1]; + ydot[1] = -1.0 * (ydot[0] + ydot[2]); +} + +static void system_scipy( double t, double* y, double* ydot, void* data) +{ + double mu = 1E4; + ydot[0] = y[1]; + ydot[1] = mu * (1- y[0]*y[0])*y[1] - y[0]; +} + +// This system is described here +// https://github.com/sdwfrost/liblsoda/issues/10 +static void system_github_issue_10(double t, double* y, double* ydot, void* data) +{ + ydot[0]= 9*y[0] + 24*y[1] + 5*cos(t)-(1/3)*sin(t); + ydot[1]= -24*y[0] -51*y[1] -95*cos(t) + (1/3)*sin(t); +} + +double test_github_system( void ) +{ + // cout << "Running test given https://github.com/sdwfrost/liblsoda/issues/10" << endl; + double t = 0e0, tout = 0.5; + + vector y = {4.0/3.0, 2.0/3.0}; + int istate = 1; + + LSODA lsoda; + + chrono::steady_clock::time_point begin = chrono::steady_clock::now(); + vector yout; + vector res; + for (size_t i = 0; i < 10; i++) + { + lsoda.lsoda_update( system_github_issue_10, 2, y, yout, &t, tout, &istate, nullptr ); + res.push_back( yout[1] ); + res.push_back( yout[2] ); + tout += 0.5; + + y[0] = yout[1]; + y[1] = yout[2]; + } + + chrono::steady_clock::time_point end = chrono::steady_clock::now(); + + if( !areEqual(-11.9400786, res[0])) cout << 'x'; + if( !areEqual( 3.8608262, res[1])) cout << 'x'; + double dt = duration_cast(end-begin).count(); + return dt; +} + +double test_scipy_sys( void ) +{ + // cout << "Running test scipy sys" << endl; + double t, tout; + t = 0e0; tout = 10; + + vector y = {10,0}; + int istate = 1; + + LSODA lsoda; + + vector yout; + chrono::steady_clock::time_point begin = chrono::steady_clock::now(); + lsoda.lsoda_update( system_scipy, 2, y, yout, &t, tout, &istate, nullptr ); + chrono::steady_clock::time_point end = chrono::steady_clock::now(); + + // printf(" at t= %12.4e y= %14.6e %14.6e\n", t, yout[1], yout[2]); + + if( ! areEqual(9.999899e+00, yout[1])) cout << 'x'; + if( !areEqual(-1.010111e-05, yout[2])) cout << 'x'; + + if (istate <= 0) + { + cerr << "error istate = " << istate << endl; + exit(0); + } + + double dt = duration_cast(end-begin).count(); + return dt; +} + +double test_fex(void) +{ + // cout << "Running test fex." << endl; + int neq = 3; + double t, tout; + t = 0e0; + tout = 0.4e0; + + vector y = { 1e0, 0e0, 0.0 }; + int istate = 1; + + LSODA lsoda; + setprecision( 12 ); + + chrono::steady_clock::time_point begin = chrono::steady_clock::now(); + + vector res; + vector yout; + for (size_t iout = 1; iout <= 12; iout++) + { + lsoda.lsoda_update( fex, neq, y, yout, &t, tout, &istate, nullptr, 1e-4, 1e-8 ); + // cerr << " at t " << t << " y= " << yout[1] << ' ' << yout[2] << ' ' << yout[3] << endl; + // Update the y for next iteration. + y[0] = yout[1]; + y[1] = yout[2]; + y[2] = yout[3]; + + res.push_back( y[0] ); + res.push_back( y[1] ); + res.push_back( y[2] ); + + if (istate <= 0) + { + cerr << "error istate = " << istate << endl; + exit(0); + } + tout = tout * 10.0E0; + } + + chrono::steady_clock::time_point end= chrono::steady_clock::now(); + double dt = chrono::duration_cast(end - begin).count(); + // cout << "|| Time taken (us)= " << dt << endl; + + vector expected = { 0.985172 ,3.3864e-05 , 0.0147939 + , 0.905514 ,2.24042e-05, 0.0944634 + , 0.715803 ,9.18446e-06, 0.284188 + , 0.450479 ,3.22234e-06, 0.549517 + , 0.183171 ,8.94046e-07, 0.816828 + , 0.0389738 ,1.62135e-07, 0.961026 + , 0.00493686 ,1.98442e-08, 0.995063 + , 0.00051665 ,2.06765e-09, 0.999483 + , 5.20075e-05 ,2.08041e-10, 0.999948 + , 5.20168e-06 ,2.08068e-11, 0.999995 + , 5.19547e-07 ,2.07819e-12, 0.999999 + }; + + // Assert here. + for (size_t i = 0; i < expected.size(); i++) + { + double err = abs( expected[i] - res[i] ); + if( err > 1e-6 ) + { + cerr << "FAILED: Expected " << expected[i] << ". Got " << res[i] << endl; + assert( false ); + } + } + + return dt; +} + +int run_serial( ) +{ + vector t1s, t2s, t3s; + size_t N = 10000; + cout << "Running " << N << " iterations. " << endl; + + // Launch all three tests in parallel. + chrono::steady_clock::time_point begin = chrono::steady_clock::now(); + for (size_t i = 0; i < N; i++) + { + double t1 = test_scipy_sys(); + double t2 = test_fex(); + double t3 = test_github_system( ); + t1s.push_back( t1 ); + t2s.push_back( t2 ); + t3s.push_back( t3 ); + } + + cout << "Avg time for 1 " << std::accumulate(t1s.begin(), t1s.end(), 0)/N + << " us per loop."<< endl; + cout << "Avg time for 2 " << std::accumulate(t2s.begin(), t2s.end(), 0)/N + << " us per loop."<< endl; + cout << "Avg time for 3 " << std::accumulate(t3s.begin(), t3s.end(), 0)/N + << " us per loop."<< endl; + + chrono::steady_clock::time_point end = chrono::steady_clock::now(); + double dt = chrono::duration_cast(end - begin).count(); + cout << "Total time taken " << dt/N << " us. (per loop)" << endl; + + + return 0; +} + +int run_multithreaded( ) +{ + vector t1s, t2s, t3s; + size_t N = 10000; + cout << "Running " << N << " iterations. " << endl; + + // Launch all three tests in parallel. + chrono::steady_clock::time_point begin = chrono::steady_clock::now(); + for (size_t i = 0; i < N; i++) + { + std::thread th1(test_scipy_sys); + std::thread th2(test_fex); + std::thread th3(test_github_system); + th1.join(); + th2.join(); + th3.join(); + } + + chrono::steady_clock::time_point end = chrono::steady_clock::now(); + double dt = chrono::duration_cast(end - begin).count(); + cout << "Total time taken " << dt / N << " us (per loop)." << endl; + return 0; +} + +int main(int argc, const char *argv[]) +{ + cout <<"|| running serial version" << endl; + run_serial(); + + cout <<"|| running multithreaded version" << endl; + run_multithreaded(); + return 0; +} diff --git a/external/libsoda/helper.h b/external/libsoda/helper.h new file mode 100644 index 0000000000..cef7f08ba0 --- /dev/null +++ b/external/libsoda/helper.h @@ -0,0 +1,24 @@ +#ifndef HELPER_H +#define HELPER_H + +template bool areEqual(T a, T b) { + return std::fabs(a - b) < 1e-6; +} + +template +void print_vec(const vector &v, const string prefix = "") { + cout << prefix << " SIZE=" << v.size() << " : "; + for (auto x : v) + cout << x << ','; + cout << endl; +} + +template +void print_arr( const double* a, const size_t n, const string prefix="" ) { + cout << prefix << " SIZE=" << n << " : "; + for (size_t i = 0; i < n; i++) + cout << a[i] << ','; + cout << endl; +} + +#endif /* end of include guard: HELPER_H */ diff --git a/external/libsoda/run_benchmark.sh b/external/libsoda/run_benchmark.sh new file mode 100755 index 0000000000..04320cc393 --- /dev/null +++ b/external/libsoda/run_benchmark.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env bash +set -e +set -x +valgrind --tool=callgrind --callgrind-out-file=callgrind.out ./test_lsoda +gprof2dot -f callgrind ./callgrind.out > callgrind.prof.dot +dot -Tpng ./callgrind.prof.dot > lsoda_prof.png diff --git a/external/libsoda/test_LSODA.cpp b/external/libsoda/test_LSODA.cpp new file mode 100644 index 0000000000..e00dd5a3a2 --- /dev/null +++ b/external/libsoda/test_LSODA.cpp @@ -0,0 +1,199 @@ +/*** + * Description: test libsoda. + * + * Created: 2018-08-14 + + * Author: Dilawar Singh + * Organization: NCBS Bangalore + * License: MIT License + */ + +#include +#include +#include +#include +#include +#include +#include + +#include "LSODA.h" +#include "helper.h" + +using namespace std; + +// Describe the system. +static void fex(double t, double *y, double *ydot, void *data) +{ + ydot[0] = 1.0E4 * y[1] * y[2] - .04E0 * y[0]; + // Don't swap ydot[1] and ydot[2]. The order will change and test will fail. + ydot[2] = 3.0E7 * y[1] * y[1]; + ydot[1] = -1.0 * (ydot[0] + ydot[2]); +} + +static void system_scipy( double t, double* y, double* ydot, void* data) +{ + double mu = 1E4; + ydot[0] = y[1]; + ydot[1] = mu * (1- y[0]*y[0]) * y[1] - y[0]; +} + +// This system is described here +// https://github.com/sdwfrost/liblsoda/issues/10 +static void system_github_issue_10(double t, double* y, double* ydot, void* data) +{ + ydot[0]= 9*y[0] + 24*y[1] + 5*cos(t)-(1/3)*sin(t); + ydot[1]= -24*y[0] -51*y[1] -95*cos(t) + (1/3)*sin(t); +} + +int test_github_system( void ) +{ + // cout << "Running test given https://github.com/sdwfrost/liblsoda/issues/10" << endl; + double t = 0e0, tout = 0.5; + + vector y = {4.0/3.0, 2.0/3.0}; + int istate = 1; + + LSODA lsoda; + + vector yout; + vector res; + + for (size_t i = 0; i < 10; i++) + { + lsoda.lsoda_update( system_github_issue_10, 2, y, yout, &t, tout, &istate, nullptr, 1e-6, 1e-6 ); + res.push_back( yout[1] ); + res.push_back( yout[2] ); + tout += 0.5; + + y[0] = yout[1]; + y[1] = yout[2]; + + cout << t << ' ' << setprecision(8) << y[0] << ' ' << y[1] << endl; + } + + cout << res[0] << ' ' << res[1] << endl; + assert( areEqual(-11.94045, res[0]) ); + assert( areEqual( 3.8610102, res[1]) ); + + if (istate <= 0) + { + cerr << "error istate = " << istate << endl; + throw runtime_error( "Failed to compute the solution." ); + } + + return 0; +} + +int test_scipy_sys( void ) +{ + cout << "Running test scipy sys" << endl; + double t, tout; + t = 0e0; tout = 10; + + vector y = {10,0}; + int istate = 1; + + LSODA lsoda; + + // Create vector to store results. NOTE THAT yout[0] will be ignored. + vector yout; + lsoda.lsoda_update( system_scipy, 2, y, yout, &t, tout, &istate, nullptr ); + + areEqual( 9.999899e+00, yout[1]); + areEqual( -1.010111e-05, yout[2] ); + + if (istate <= 0) + { + cerr << "error istate = " << istate << endl; + exit(0); + } + return 0; +} + +int test_fex(void) +{ + // cout << "Running test fex." << endl; + int neq = 3; + double t, tout; + t = 0e0; + tout = 0.4e0; + vector y = { 1e0, 0e0, 0.0 }; + int istate = 1; + + LSODA lsoda; + setprecision( 12 ); + + vector res; + + vector yout; + for (size_t iout = 1; iout <= 12; iout++) + { + lsoda.lsoda_update( fex, neq, y, yout, &t, tout, &istate, nullptr, 1e-4, 1e-8 ); + // cerr << " at t " << t << " y= " << yout[1] << ' ' << yout[2] << ' ' << yout[3] << endl; + // Update the y for next iteration. + y[0] = yout[1]; + y[1] = yout[2]; + y[2] = yout[3]; + + res.push_back( y[0] ); + res.push_back( y[1] ); + res.push_back( y[2] ); + + if (istate <= 0) + { + cerr << "error istate = " << istate << endl; + exit(0); + } + tout = tout * 10.0E0; + } + + + vector expected = { 0.985172 ,3.3864e-05 , 0.0147939 + , 0.905514 ,2.24042e-05, 0.0944634 + , 0.715803 ,9.18446e-06, 0.284188 + , 0.450479 ,3.22234e-06, 0.549517 + , 0.183171 ,8.94046e-07, 0.816828 + , 0.0389738 ,1.62135e-07, 0.961026 + , 0.00493686 ,1.98442e-08, 0.995063 + , 0.00051665 ,2.06765e-09, 0.999483 + , 5.20075e-05 ,2.08041e-10, 0.999948 + , 5.20168e-06 ,2.08068e-11, 0.999995 + , 5.19547e-07 ,2.07819e-12, 0.999999 + }; + + // Assert here. + for (size_t i = 0; i < expected.size(); i++) + { + double err = abs( expected[i] - res[i] ); + if( err > 1e-6 ) + { + cerr << "FAILED: Expected " << expected[i] << ". Got " << res[i] << endl; + assert( false ); + } + } + + return 0; +} + +int main(int argc, const char *argv[]) +{ + chrono::steady_clock::time_point begin = chrono::steady_clock::now(); + test_scipy_sys(); + chrono::steady_clock::time_point end= chrono::steady_clock::now(); + cout << "|| Time taken (us)= " << chrono::duration_cast(end - begin).count() < #include "VoxelPools.h" diff --git a/ksolve/CMakeLists.txt b/ksolve/CMakeLists.txt index f5d8cf840a..43aa6e84cb 100644 --- a/ksolve/CMakeLists.txt +++ b/ksolve/CMakeLists.txt @@ -4,7 +4,7 @@ include( ${CMAKE_CURRENT_SOURCE_DIR}/../CheckCXXCompiler.cmake) if(WITH_BOOST) find_package(Boost 1.53 REQUIRED COMPONENTS thread) - add_definitions( -DUSE_BOOST_ASYNC ) + add_definitions(-DUSE_BOOST ) set(WITH_BOOST_ODE ON) include_directories( ${Boost_INCLUDE_DIRS} ) # This is still not part of official bindings. @@ -13,12 +13,15 @@ endif() # if boost ode is being used, don't use GSL. if(WITH_BOOST_ODE) + # This is still not part of official bindings. + set(WITH_GSL OFF) + include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../external/boost-numeric-bindings) add_definitions(-DUSE_BOOST_ODE -UUSE_GSL) include_directories(${Boost_INCLUDE_DIRS}) include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../external/boost-numeric-bindings) elseif(WITH_GSL) include_directories( ${GSL_INCLUDE_DIRS} ) - add_definitions(-DUSE_GSL -UUSE_BOOST_ODE) + add_definitions(-DUSE_GSL -UUSE_BOOST_ODE -UUSE_BOOST) endif(WITH_BOOST_ODE) set(KSOLVE_SRCS diff --git a/ksolve/FuncRateTerm.h b/ksolve/FuncRateTerm.h index bfb7db56d7..720026fe19 100644 --- a/ksolve/FuncRateTerm.h +++ b/ksolve/FuncRateTerm.h @@ -32,7 +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. - assert(std::isfinite(v)); + assert(! std::isnan(v)); return v; } @@ -68,7 +68,7 @@ class FuncRate: public ExternReac RateTerm* copyWithVolScaling(double vol, double sub, double prd ) const { - assert(std::isfinite(vol)); + assert(! std::isnan(vol)); double ratio = sub * pow( NA * vol, funcVolPower_ ); FuncRate* ret = new FuncRate( k_ / ratio, func_->getTarget() ); ret->funcVolPower_ = funcVolPower_; diff --git a/ksolve/KinSparseMatrix.cpp b/ksolve/KinSparseMatrix.cpp index 2b0c9d3701..ac70328d71 100644 --- a/ksolve/KinSparseMatrix.cpp +++ b/ksolve/KinSparseMatrix.cpp @@ -40,7 +40,7 @@ double KinSparseMatrix::computeRowRate( for ( const int* i = entry; i != end; ++i ) { ret += *i * v[ *colIndex++ ]; - assert(std::isfinite(ret)); + assert(! std::isnan(ret)); } //assert ( !std::isnan(ret) ); return ret; diff --git a/ksolve/Ksolve.cpp b/ksolve/Ksolve.cpp index 5106460061..e000adc367 100644 --- a/ksolve/Ksolve.cpp +++ b/ksolve/Ksolve.cpp @@ -59,7 +59,8 @@ const Cinfo* Ksolve::initCinfo() "rk4: The Runge-Kutta 4th order fixed dt method" "rk2: The Runge-Kutta 2,3 embedded fixed dt method" "rkck: The Runge-Kutta Cash-Karp (4,5) method" - "rk8: The Runge-Kutta Prince-Dormand (8,9) method" , + "rk8: The Runge-Kutta Prince-Dormand (8,9) method" + "lsoda: LSODA method", &Ksolve::setMethod, &Ksolve::getMethod ); @@ -110,7 +111,7 @@ const Cinfo* Ksolve::initCinfo() static ValueFinfo< Ksolve, unsigned int > numThreads ( "numThreads", - "Number of threads to use (applicable in deterministic case)", + "Number of threads to use", &Ksolve::setNumThreads, &Ksolve::getNumThreads ); @@ -257,20 +258,33 @@ string Ksolve::getMethod() const void Ksolve::setMethod( string method ) { + std::transform(method.begin(), method.end(), method.begin(), ::tolower); + // If user is trying to set ksolve method after ksolve has been initialized, + // show a warning. + if( isBuilt_ ) + { + moose::showWarn( + "You are trying to set Ksolve::method after moose::Stoich has been " + " initialized. This will be ignored. Please do before ksolve is assigned to " + " moose::Stoich." + ); + return; + } + #if USE_GSL if ( method == "rk5" || method == "gsl" ) { method_ = "rk5"; } else if ( method == "rk4" || method == "rk2" || - method == "rk8" || method == "rkck" ) + method == "rk8" || method == "rkck" || method == "lsoda" ) { method_ = method; } else { cout << "Warning: Ksolve::setMethod: '" << method << - "' not known, using rk5\n"; + "' is not known, using default rk5\n"; method_ = "rk5"; } #elif USE_BOOST_ODE diff --git a/ksolve/Ksolve.h b/ksolve/Ksolve.h index 293f1e90b8..3c1eba96f2 100644 --- a/ksolve/Ksolve.h +++ b/ksolve/Ksolve.h @@ -37,6 +37,10 @@ class Ksolve: public ZombiePoolInterface double getEpsRel() const; void setEpsRel( double val ); + // To make API consistent with GssaVoxelPools + double getRelativeAccuracy( ) const; + double getAbsoluteAccuracy( ) const; + /// Assigns Stoich object to Ksolve. Id getStoich() const; void setStoich( Id stoich ); /// Inherited from ZombiePoolInterface. diff --git a/ksolve/VoxelPools.cpp b/ksolve/VoxelPools.cpp index 8c7c670e37..f28199d54b 100644 --- a/ksolve/VoxelPools.cpp +++ b/ksolve/VoxelPools.cpp @@ -1,13 +1,15 @@ -/* -* This program is part of 'MOOSE', the -* Messaging Object Oriented Simulation Environment. -* Copyright (C) 2003-2010 Upinder S. Bhalla. and NCBS -* It is made available under the terms of the -* GNU Lesser General Public License version 2.1 -* See the file COPYING.LIB for the full notice. -*/ +/*** + * Description: Class VoxelPool. + * Authors: Upinder Bhalla , + * Dilawar Singh + * Organization: NCBS Bangalore + * License: GNU GPL3 + */ + +#include #include "../basecode/header.h" +#include "../basecode/SparseMatrix.h" #ifdef USE_GSL #include @@ -23,17 +25,16 @@ using namespace boost::numeric; #include "VoxelPools.h" #include "RateTerm.h" #include "FuncTerm.h" -#include "../basecode/SparseMatrix.h" #include "KinSparseMatrix.h" -#include "../mesh/VoxelJunction.h" #include "XferInfo.h" #include "ZombiePoolInterface.h" +#include "Ksolve.h" #include "Stoich.h" ////////////////////////////////////////////////////////////// // Class definitions -VoxelPools::VoxelPools() +VoxelPools::VoxelPools() : pLSODA(nullptr) { #ifdef USE_GSL driver_ = 0; @@ -62,6 +63,14 @@ void VoxelPools::reinit( double dt ) gsl_odeiv2_driver_reset( driver_ ); gsl_odeiv2_driver_reset_hstart( driver_, dt / 10.0 ); #endif + + // If method is LDODA create lSODA object and save the address of this as + // param (void*). + if( getMethod() == "lsoda" ) + { + pLSODA.reset(new LSODA()); + pLSODA->param = (void *) this; + } } void VoxelPools::setStoich( Stoich* s, const OdeSystem* ode ) @@ -82,105 +91,204 @@ void VoxelPools::setStoich( Stoich* s, const OdeSystem* ode ) gsl_odeiv2_driver_free( driver_ ); driver_ = gsl_odeiv2_driver_alloc_y_new( &sys_, ode->gslStep - , ode->initStepSize, ode->epsAbs, ode->epsRel - ); + , ode->initStepSize, ode->epsAbs, ode->epsRel); } #endif VoxelPoolsBase::reinit(); } +const string VoxelPools::getMethod( ) +{ + Ksolve* k = reinterpret_cast( stoichPtr_->getKsolve().eref().data() ); + return k->getMethod( ); +} + void VoxelPools::advance( const ProcInfo* p ) { double t = p->currTime - p->dt; + Ksolve* k = reinterpret_cast( stoichPtr_->getKsolve().eref().data() ); -#ifdef USE_GSL - int status = gsl_odeiv2_driver_apply( driver_, &t, p->currTime, varS()); - if ( status != GSL_SUCCESS ) + if( getMethod() == "lsoda" ) { - cout << "Error: VoxelPools::advance: GSL integration error at time " - << t << "\n"; - cout << "Error info: " << status << ", " << - gsl_strerror( status ) << endl; - if ( status == GSL_EMAXITER ) - cout << "Max number of steps exceeded\n"; - else if ( status == GSL_ENOPROG ) - cout << "Timestep has gotten too small\n"; - else if ( status == GSL_EBADFUNC ) - cout << "Internal error\n"; - assert( 0 ); - } + vector yout(size()+1); + pLSODA->lsoda_update( &VoxelPools::lsodaSys, size() + , Svec(), yout , &t + , p->currTime, &lsodaState, this + ); -#elif USE_BOOST_ODE + // Now update the y from yout. This is different thant normal GSL or + // BOOST based approach. + for (size_t i = 0; i < size(); i++) + varS()[i] = yout[i+1]; - /*----------------------------------------------------------------------------- - We need to call updateFuncs here (unlike in GSL solver) because there - is no way we can update const vector_type_& y in evalRatesUsingBoost - function. In gsl implmentation one could do it, because const_cast can - take away the constantness of double*. This probably makes the call bit - cleaner. - *-----------------------------------------------------------------------------*/ - stoichPtr_->updateFuncs( &Svec()[0], p->currTime ); - - /*----------------------------------------------------------------------------- - * Using integrate function works with with default stepper type. - * - * NOTICE to developer: - * If you are planning your own custom typdedef of stepper_type_ (see - * file BoostSystem.h), the you may run into troble. Have a look at this - * http://boostw.boost.org/doc/libs/1_56_0/boost/numeric/odeint/integrate/integrate.hpp - *----------------------------------------------------------------------------- - */ - - /** - * @brief Default step size for fixed size iterator. - */ - - // Variout stepper times are listed here: - // https://www.boost.org/doc/libs/1_68_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.explicit_steppers - - // Describe system to be used in boost solver calls. - auto sys = [this](const vector_type_& dy, vector_type_& dydt, const double t) - { - VoxelPools::evalRates(this, dy, dydt); - }; - - // This is usually the default method. It works well in practice. Tested - // with steady-state solver. Closest to GSL rk5 . - if( method_ == "rk5" || method_ == "gsl" || method_ == "boost" ) - odeint::integrate_adaptive( - make_dense_output( epsAbs_, epsRel_, odeint::runge_kutta_dopri5() ) - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt - ); - else if( method_ == "rk5a" || method_ == "adaptive" ) - odeint::integrate_adaptive( odeint::make_controlled( epsAbs_, epsRel_ ) - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt ); - else if( method_ == "rk2" ) - odeint::integrate_const( rk_midpoint_stepper_type_() - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt); - else if( method_ == "rk4" ) - odeint::integrate_const( rk4_stepper_type_() - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt ); - else if ("rk54" == method_ ) - odeint::integrate_const( rk_karp_stepper_type_() - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt); - else if ("rkck" == method_ ) - odeint::integrate_adaptive( odeint::make_controlled( epsAbs_, epsRel_ ) - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt); - else if( method_ == "rk8" ) - odeint::integrate_const( rk_felhberg_stepper_type_() - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt); - else if( method_ == "rk8a" ) - odeint::integrate_adaptive( rk_felhberg_stepper_type_() - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt); + if( lsodaState == 0 ) + { + cerr << "Error: VoxelPools::advance: LSODA integration error at time " + << t << "\n"; + assert(0); + } + } else { - cerr << "Ksolve: Unknow method " << method_ << ", using default!" << endl; - odeint::integrate_const( - make_dense_output( epsAbs_, epsRel_, odeint::runge_kutta_dopri5() ) - , sys, Svec(), p->currTime - p->dt, p->currTime, p->dt - ); + +#ifdef USE_GSL + int status = gsl_odeiv2_driver_apply( driver_, &t, p->currTime, varS()); + if ( status != GSL_SUCCESS ) + { + cerr << "Error: VoxelPools::advance: GSL integration error at time " + << t << "\n"; + cerr << "Error info: " << status << ", " << + gsl_strerror( status ) << endl; + if ( status == GSL_EMAXITER ) + cerr << "Max number of steps exceeded\n"; + else if ( status == GSL_ENOPROG ) + cerr << "Timestep has gotten too small\n"; + else if ( status == GSL_EBADFUNC ) + cerr << "Internal error\n"; + assert( 0 ); + } + +#elif USE_BOOST_ODE + // NOTE: Make sure to assing vp to BoostSys vp. In next call, it will be used by + // updateRates func. Unlike gsl call, we can't pass extra void* to gslFunc. + VoxelPools* vp = reinterpret_cast< VoxelPools* >( this ); + + /*----------------------------------------------------------------------------- +NOTE: 04/21/2016 11:31:42 AM + +We need to call updateFuncs here (unlike in GSL solver) because there +is no way we can update const vector_type_& y in evalRatesUsingBoost +function. In gsl implmentation one could do it, because const_cast can +take away the constantness of double*. This probably makes the call bit +cleaner. + *-----------------------------------------------------------------------------*/ + stoichPtr_->updateFuncs( &Svec()[0], p->currTime ); + + /*----------------------------------------------------------------------------- + * Using integrate function works with with default stepper type. + * + * NOTICE to developer: + * If you are planning your own custom typdedef of stepper_type_ (see + * file BoostSystem.h), the you may run into troble. Have a look at this + * http://boostw.boost.org/doc/libs/1_56_0/boost/numeric/odeint/integrate/integrate.hpp + * + * To make numerical results comparable with "gsl" solvers, + * by default, we pick adaptive step-size control. Use a suffix 'c' to + * force a constant step size. The numerical accuracy of constant step + * size solver is low. + * More details can be found here: + * https://www.boost.org/doc/libs/1_72_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html + */ + + const double fixedDt = 0.1; + + if( method_ == "rk2" ) + odeint::integrate_const( rk_midpoint_stepper_type_() + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt ) + ); + else if( method_ == "rk4c" ) + odeint::integrate_const( rk4_stepper_type_() + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt ) + ); + else if( method_ == "rk5c") + odeint::integrate_const( rk_karp_stepper_type_() + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt ) + ); + else if( method_ == "rk5ck" ) + odeint::integrate_adaptive( + odeint::make_controlled( epsAbs_, epsRel_ ) + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt + , p->currTime + , p->dt + ); + else if ("rk54c" == method_ ) + odeint::integrate_const( rk_karp_stepper_type_() + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt ) + ); + else if ("rk54" == method_ ) + odeint::integrate_adaptive( + odeint::make_controlled( epsAbs_, epsRel_ ) + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt + , p->currTime + , p->dt + ); + else if ("rk5c" == method_ ) + odeint::integrate_const( rk_dopri_stepper_type_() + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt + , p->currTime + , std::min( p->dt, fixedDt ) + ); + else if ("rk5" == method_ || "gsl" == method_) + odeint::integrate_adaptive( + odeint::make_controlled( epsAbs_, epsRel_ ) + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt + , p->currTime + , p->dt + ); + else if( method_ == "rk8c" ) + odeint::integrate_const( rk_felhberg_stepper_type_() + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt, p->currTime, std::min( p->dt, fixedDt ) + ); + else if( method_ == "rk8" ) + odeint::integrate_adaptive( + odeint::make_controlled( epsAbs_, epsRel_ ) + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt + , p->currTime + , p->dt + ); + else + odeint::integrate_adaptive( + odeint::make_controlled( epsAbs_, epsRel_ ) + , [this](const vector_type_& dy, vector_type_& dydt, const double t) { + VoxelPools::evalRates(this, dy, dydt ); + } + , Svec() + , p->currTime - p->dt + , p->currTime + , p->dt + ); +#endif // USE_GSL } -#endif if ( !stoichPtr_->getAllowNegative() ) // clean out negatives { @@ -188,7 +296,7 @@ void VoxelPools::advance( const ProcInfo* p ) double* vs = varS(); for ( unsigned int i = 0; i < nv; ++i ) { - if ( signbit(vs[i]) ) + if ( std::signbit(vs[i]) ) vs[i] = 0.0; } } @@ -202,18 +310,14 @@ void VoxelPools::setInitDt( double dt ) } #ifdef USE_GSL -// static func. This goes into the Gsl solver. -int VoxelPools::gslFunc(double t, const double* y, double *dydt, void* params) +// static func. This is the function that goes into the Gsl solver. +int VoxelPools::gslFunc( double t, const double* y, double *dydt, void* params ) { VoxelPools* vp = reinterpret_cast< VoxelPools* >( params ); double* q = const_cast< double* >( y ); // Assign the func portion. vp->stoichPtr_->updateFuncs( q, t ); vp->updateRates( y, dydt ); -#ifdef USE_GSL return GSL_SUCCESS; -#else - return 0; -#endif } #elif USE_BOOST_ODE // NOT GSL @@ -223,7 +327,23 @@ void VoxelPools::evalRates( VoxelPools* vp, const vector_type_& y, vector_type_ vp->updateRates( &y[0], &dydt[0] ); } -#endif +#endif // USE_BOOST_ODE + +/* --------------------------------------------------------------------------*/ +/** + * @Synopsis Function to pass LSODA::lsoda_update function. Since it is + a * static function, we have to make sure void* param holds the value of + pointer * to VoxelPools. * * @Param t * @Param y * @Param dydt * @Param + params Address of VoxelPools as void*. */ +/* +----------------------------------------------------------------------------*/ +void VoxelPools::lsodaSys( double t, double* y, double* dydt, void* param) +{ + VoxelPools* vp = reinterpret_cast< VoxelPools* >( param ); + // Fill in the values. + vp->stoichPtr_->updateFuncs( y, t ); + vp->updateRates( y, dydt ); +} /////////////////////////////////////////////////////////////////////// // Here are the internal reaction rate calculation functions @@ -276,13 +396,10 @@ void VoxelPools::updateRates( const double* s, double* yprime ) const const KinSparseMatrix& N = stoichPtr_->getStoichiometryMatrix(); vector< double > v( N.nColumns(), 0.0 ); vector< double >::iterator j = v.begin(); - // totVar should include proxyPools only if this voxel uses them unsigned int totVar = stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools(); - // totVar should include proxyPools if this voxel does not use them unsigned int totInvar = stoichPtr_->getNumBufPools(); - assert( N.nColumns() == 0 || N.nRows() == stoichPtr_->getNumAllPools() ); assert( N.nColumns() == rates_.size() ); @@ -291,12 +408,11 @@ void VoxelPools::updateRates( const double* s, double* yprime ) const for (unsigned int i = 0; i < totVar; ++i) { auto rate = N.computeRowRate( i, v ); - assert(std::isfinite(rate)); + assert(! std::isnan(rate)); *yprime++ = rate; } for (unsigned int i = 0; i < totInvar ; ++i) *yprime++ = 0.0; - } /** @@ -304,7 +420,7 @@ void VoxelPools::updateRates( const double* s, double* yprime ) const * This is a utility function for programs like SteadyState that need * to analyze velocity. */ -void VoxelPools::updateReacVelocities(const double* s, vector& v) const +void VoxelPools::updateReacVelocities(const double* s, vector< double >& v) const { const KinSparseMatrix& N = stoichPtr_->getStoichiometryMatrix(); assert( N.nColumns() == rates_.size() ); @@ -318,7 +434,7 @@ void VoxelPools::updateReacVelocities(const double* s, vector& v) const for ( i = rates_.begin(); i != rates_.end(); i++) { *j++ = (**i)( s ); - assert(std::isfinite(*(j-1))); + assert(! std::isnan(*(j-1))); } } diff --git a/ksolve/VoxelPools.h b/ksolve/VoxelPools.h index 04c5bb8c7e..4c79db7f4e 100644 --- a/ksolve/VoxelPools.h +++ b/ksolve/VoxelPools.h @@ -12,12 +12,12 @@ #include "OdeSystem.h" #include "VoxelPoolsBase.h" +#include "external/libsoda/LSODA.h" #ifdef USE_BOOST_ODE #include "BoostSys.h" #endif - class Stoich; class ProcInfo; @@ -45,18 +45,23 @@ class VoxelPools: public VoxelPoolsBase const Stoich* getStoich( ); + const string getMethod( ); + /// Do the numerical integration. Advance the simulation. void advance( const ProcInfo* p ); /// Set initial timestep to use by the solver. void setInitDt( double dt ); -#ifdef USE_GSL /* ----- not USE_BOOST_ODE ----- */ +#ifdef USE_GSL /* ----- not USE_BOOST ----- */ static int gslFunc( double t, const double* y, double *dydt, void* params); #elif USE_BOOST_ODE static void evalRates( VoxelPools* vp, const vector_type_& y, vector_type_& dydt ); #endif /* ----- not USE_BOOST_ODE ----- */ + // System of LSODA. + static void lsodaSys( double t, double* y, double* dydt, void* params); + ////////////////////////////////////////////////////////////////// // Rate manipulation and calculation functions ////////////////////////////////////////////////////////////////// @@ -87,13 +92,17 @@ class VoxelPools: public VoxelPoolsBase * This is a utility function for programs like SteadyState that * need to analyze velocity. */ - void updateReacVelocities( - const double* s, vector< double >& v ) const; + void updateReacVelocities( const double* s, vector< double >& v ) const; /// Used for debugging. void print() const; + private: + std::shared_ptr pLSODA; + LSODA_ODE_SYSTEM_TYPE lsodaSystem; + int lsodaState = 1; + #ifdef USE_GSL gsl_odeiv2_driver* driver_; gsl_odeiv2_system sys_; diff --git a/ksolve/helper.h b/ksolve/helper.h new file mode 100644 index 0000000000..4535e002db --- /dev/null +++ b/ksolve/helper.h @@ -0,0 +1,18 @@ +#ifndef HELPER_H +#define HELPER_H + +template +bool areEqual(T a, T b) +{ + return std::fabs(a - b) < 1e-6; +} + +template +void print_vec( const vector& v, const string prefix="" ) +{ + cout << prefix << " SIZE=" << v.size() << " : "; + for( auto x : v ) cout << x << ','; + cout << endl; +} + +#endif /* end of include guard: HELPER_H */ diff --git a/ksolve/testKsolve.cpp b/ksolve/testKsolve.cpp index c3788bdec9..167ff30800 100644 --- a/ksolve/testKsolve.cpp +++ b/ksolve/testKsolve.cpp @@ -8,6 +8,7 @@ **********************************************************************/ #include "../basecode/header.h" #include "../shell/Shell.h" + #include "RateTerm.h" #include "FuncTerm.h" #include "../basecode/SparseMatrix.h" @@ -18,6 +19,7 @@ #include "XferInfo.h" #include "ZombiePoolInterface.h" #include "Stoich.h" +#include "../mesh/VoxelJunction.h" #include "../builtins/MooseParser.h" #include "../utility/testing_macros.hpp" @@ -263,6 +265,39 @@ void testRunKsolve() cout << "." << flush; } +void testRunKsolveWithLSODA() +{ + double simDt = 0.1; + // double plotDt = 0.1; + Shell* s = reinterpret_cast< Shell* >( Id().eref().data() ); + Id kin = makeReacTest(); + + Id ksolve = s->doCreate( "Ksolve", kin, "ksolve", 1 ); + Field< string >::set( ksolve, "method", "lsoda" ); + + Id stoich = s->doCreate( "Stoich", ksolve, "stoich", 1 ); + Field< Id >::set( stoich, "compartment", kin ); + Field< Id >::set( stoich, "ksolve", ksolve ); + Field< string >::set( stoich, "path", "/kinetics/##" ); + s->doUseClock( "/kinetics/ksolve", "process", 4 ); + s->doSetClock( 4, simDt ); + + s->doReinit(); + s->doStart( 20.0 ); + Id plots( "/kinetics/plots" ); + for ( unsigned int i = 0; i < 7; ++i ) + { + stringstream ss; + ss << "plot." << i; + SetGet2< string, string >::set( ObjId( plots, i ) + , "xplot", "tsr2lsoda.plot", ss.str() + ); + } + s->doDelete( kin ); + cout << "." << flush; +} + + void testRunGsolve() { double simDt = 0.1; diff --git a/pymoose/moosemodule.cpp b/pymoose/moosemodule.cpp index ca25ec3f5f..3084d47234 100644 --- a/pymoose/moosemodule.cpp +++ b/pymoose/moosemodule.cpp @@ -2543,7 +2543,6 @@ int defineLookupFinfos(const Cinfo * cinfo) return 1; } - int defineClass(PyObject * module_dict, const Cinfo * cinfo) { const string& className = cinfo->name(); diff --git a/pymoose/pymooseinit.cpp b/pymoose/pymooseinit.cpp index add0fee13b..00674e2fa0 100644 --- a/pymoose/pymooseinit.cpp +++ b/pymoose/pymooseinit.cpp @@ -59,7 +59,7 @@ extern void testBiophysicsProcess(); extern unsigned int initMsgManagers(); extern void destroyMsgManagers(); // void regressionTests(); -#endif +#endif // DO_UNIT_TESTS extern void speedTestMultiNodeIntFireNetwork( unsigned int size, unsigned int runsteps ); diff --git a/python/moose/chemUtil/add_Delete_ChemicalSolver.py b/python/moose/chemUtil/add_Delete_ChemicalSolver.py index 9b11a5c148..ac69367568 100644 --- a/python/moose/chemUtil/add_Delete_ChemicalSolver.py +++ b/python/moose/chemUtil/add_Delete_ChemicalSolver.py @@ -113,9 +113,9 @@ def setCompartmentSolver(modelRoot, solver): for compt in compts: if solver != 'ee': - if (solver == 'gsl'): + if solver.lower() in [ 'gsl', 'runge kutta', 'lsoda' ]: ksolve = moose.Ksolve(compt.path + '/ksolve') - if (solver == 'gssa') : + elif solver.lower() in ['gssa', 'gillespie']: ksolve = moose.Gsolve(compt.path + '/gsolve') if (len(compts) > 1): diff --git a/python/moose/graph_utils.py b/python/moose/graph_utils.py deleted file mode 100644 index bec44a74a5..0000000000 --- a/python/moose/graph_utils.py +++ /dev/null @@ -1,85 +0,0 @@ -# -*- coding: utf-8 -*- -"""graph_utils.py: Graph related utilties. It does not require networkx library. -It writes files to be used with graphviz. - -Last modified: Sat Jan 18, 2014 05:01PM - -""" - -__author__ = "Dilawar Singh" -__copyright__ = "Copyright 2013, NCBS Bangalore" -__credits__ = ["NCBS Bangalore", "Bhalla Lab"] -__license__ = "GPL" -__version__ = "1.0.0" -__maintainer__ = "Dilawar Singh" -__email__ = "dilawars@ncbs.res.in" -__status__ = "Development" - -import sys -from . import _moose -import inspect -from . import print_utils as debug -import re - -pathPat = re.compile(r'.+?\[\d+\]$') - -def getMoosePaths(pat, isRoot=True): - ''' Return a list of paths for a given pattern. ''' - if type(pat) != str: - pat = pat.path - assert type(pat) == str - moose_paths = [x.path for x in _moose.wildcardFind(pat)] - return moose_paths - -def writeGraphviz(filename=None, pat='/##[TYPE=Compartment]'): - '''This is a generic function. It takes the the pattern, search for paths - and write a graphviz file. - ''' - - def fix(path): - '''Fix a given path so it can be written to a graphviz file''' - # If no [0] is at end of the path then append it. - global pathPat - if not pathPat.match(path): - path = path + '[0]' - return path - - - pathList = getMoosePaths(pat) - compList = _moose.wildcardFind(pat) - if not compList: - debug.dump("WARN" - , "No compartment found" - , frame = inspect.currentframe() - ) - - dot = [] - dot.append("digraph G {") - dot.append("\tconcentrate=true;") - for c in compList: - if c.neighbors['raxial']: - for n in c.neighbors['raxial']: - lhs = fix(c.path) - rhs = fix(n.path) - dot.append('\t"{}" -> "{}";'.format(lhs, rhs)) - elif c.neighbors['axial']: - for n in c.neighbors['axial']: - lhs = fix(c.path) - rhs = fix(n.path) - dot.append('\t"{}" -> "{}" [dir=back];'.format(lhs, rhs)) - else: - p = fix(c.path) - dot.append('\t"{}"'.format(p)) - - dot.append('}') - dot = '\n'.join(dot) - if not filename: - print(dot) - else: - with open(filename, 'w') as graphviz: - debug.dump("INFO" - , "Writing compartment topology to file {}".format(filename) - ) - graphviz.write(dot) - return True - diff --git a/python/moose/graphutil.py b/python/moose/graphutil.py deleted file mode 100644 index ed9e9060ed..0000000000 --- a/python/moose/graphutil.py +++ /dev/null @@ -1,193 +0,0 @@ -# -*- coding: utf-8 -*- -# graphutil.py --- -# -# Filename: graphutil.py -# Description: -# Author: Subhasis Ray -# Maintainer: -# Created: Sun Mar 18 13:42:28 2012 (+0530) -# Version: -# Last-Updated: Wed Mar 28 18:27:26 2012 (+0530) -# By: subha -# Update #: 103 -# URL: -# Keywords: -# Compatibility: -# -# - -# Commentary: -# -# -# -# - -# Change log: -# -# - -# Code: - -from collections import defaultdict -import re -import networkx as nx -import numpy as np -from . import moose - -def moosegraph(element, ies=['childMsg'], ied=['parentMsg'], iv=[], keep_solitary=False): - """Create a graph out of all objects under teh element - specified. Ignore edges going outside the subtree rooted at - element. - - If element is a string, it can be a wildcard path in MOOSE - - ies is a list of sourcefields to be ignored when constructing edges. - - ied is a list of destination fields to be ignored when - constructing edges. - - iv is a list of paths, possibly wildcards, to be ignored when - building vertex set. - - keep_solitary -- if True solitary nodes are not discarded - - Note that this is a temporary solution. We rely on neighbours - field to create the graph. Ideally we should use the messages - between objects. A bug in Messages does not allow us to access the - fields on both sides properly in the general case. We shall switch - to messgaes once that bug is fixed.""" - - path = None - if isinstance(element, str): - path = element - elif isinstance(element, moose.Id): - path = element.getPath() - elif isinstance(element, moose.ObjId): - path = element.getField('path') - elif isinstance(element, moose.Neutral): - path = element.path - else: - raise TypeError('Require Id, ObjId, Neutral or string object for element') - iv_re = [re.compile(pattern) for pattern in iv] - ies_re = [re.compile(pattern) for pattern in ies] - ied_re = [re.compile(pattern) for pattern in ied] - if '#' in path: - all_v = moose.wildcardFind(path) - else: - if path.endswith('/'): - path = '%s##' % (path) - else: - path = '%s/##' % (path) - all_v = moose.wildcardFind(path) - valid_v = [] - # Collect all not-to-be-ignored vertices - for vv in all_v: - include = True - for iv in iv_re: - if iv.search(vv.path): - include = False - break - if include: - valid_v.append(vv) - graph = nx.DiGraph() - graph.add_nodes_from([vv.id_ for vv in valid_v]) - for vv in graph.nodes(): - for fname in vv[0].getFieldNames('srcFinfo'): - matches = [True for regex in ies_re if regex.search(fname)] - if matches: - continue - nlist = vv[0].getNeighbors(fname) - for nn in nlist: - if nn in graph.nodes(): - try: - src = graph.edge[vv][nn]['src'] - src.add(fname) - except KeyError: - graph.add_edge(vv, nn, src=set([fname])) - - for fname in vv[0].getFieldNames('destFinfo'): - matches = [True for regex in ied_re if regex.search(fname)] - if matches: - continue - nlist = vv[0].getNeighbors(fname) - for nn in nlist: - if nn in graph.nodes(): - try: - dest = graph.edge[vv][nn]['dest'] - dest.add(fname) - except KeyError: - graph.add_edge(nn, vv, dest=set([fname])) - - for fname in vv[0].getFieldNames('sharedFinfo'): - nlist = vv[0].getNeighbors(fname) - for nn in nlist: - if nn in graph.nodes(): - try: - src = graph.edge[vv][nn]['src'] - src.add(fname) - except KeyError: - graph.add_edge(vv, nn, src=set([fname])) - try: - dest = graph.edge[nn][vv]['dest'] - dest.add(fname) - except KeyError: - graph.add_edge(nn, vv, dest=set([fname])) - if not keep_solitary: - solitary=[ n for n,d in graph.degree_iter() if d==0 ] - graph.remove_nodes_from(solitary) - return graph - -def draw_moosegraph(graph, pos=None, label_edges=True): - if pos is None: - pos = nx.spring_layout(graph) - edge_labels = {} - if label_edges: - for ee in graph.edges(): - try: - src = ','.join(graph.edge[ee[0]][ee[1]]['src']) - except KeyError: - src = '?' - try: - dest = ','.join(graph.edge[ee[0]][ee[1]]['dest']) - except KeyError: - dest = '?' - label = '%s-%s' % (src, dest) - edge_labels[ee] = label - nx.draw(graph, pos) - nx.draw_networkx_edge_labels(graph, pos, edge_labels) - -import unittest -import pylab - -class TestMooseGraph(unittest.TestCase): - def setUp(self): - cell = moose.Neutral('/n') - comp1 = moose.Compartment('/n/c1') - chan = moose.HHChannel('/n/c1/ch') - moose.connect(comp1, 'channel', chan, 'channel') - comp2 = moose.Compartment('/n/c2') - comp2.connect('raxial', comp1, 'axial') - self.graph = moosegraph('/n') - - def test_edgelabels(self): - for ee in self.graph.edges(): - keys = self.graph.edge[ee[0]][ee[1]].keys() - self.assertTrue('src' in keys) - self.assertTrue('dest' in keys) - - def test_vertexid(self): - for vv in self.graph.nodes(): - self.assertTrue(isinstance(vv, moose.Id)) - - def test_plot(self): - draw_moosegraph(self.graph) - pylab.show() - - -if '__main__' == __name__: - unittest.main() - - - -# -# graphutil.py ends here diff --git a/python/moose/moose_constants.py b/python/moose/moose_constants.py index 804bc7ec21..848213b310 100644 --- a/python/moose/moose_constants.py +++ b/python/moose/moose_constants.py @@ -21,13 +21,9 @@ """ __author__ = "Dilawar Singh" -__copyright__ = "Copyright 2013, Dilawar Singh, NCBS Bangalore" -__credits__ = ["NCBS Bangalore"] __license__ = "GNU GPL" -__version__ = "1.0.0" __maintainer__ = "Dilawar Singh" __email__ = "dilawars@ncbs.res.in" -__status__ = "Development" ## for Ca Pool # FARADAY = 96154.0 @@ -66,4 +62,3 @@ LOOKUPCLOCK = 6 STIMCLOCK = 7 PLOTCLOCK = 8 - diff --git a/python/moose/moose_test.py b/python/moose/moose_test.py index 59d6dfea40..3238bafe5c 100644 --- a/python/moose/moose_test.py +++ b/python/moose/moose_test.py @@ -62,6 +62,10 @@ def run(self, timeout, **kwargs): try: timer.start() output, errors = self.process.communicate() + try: + errors = errors.decode('utf8') + except Exception: + pass logging.warn('%s %s' % (output, errors)) finally: timer.cancel() diff --git a/python/moose/nettwork_utils.py b/python/moose/nettwork_utils.py new file mode 100644 index 0000000000..a87219f2ba --- /dev/null +++ b/python/moose/nettwork_utils.py @@ -0,0 +1,129 @@ +# -*- coding: utf-8 -*- +from __future__ import print_function, division, absolute_import + +# Some network analysis utilities. + +__author__ = "Dilawar Singh" +__copyright__ = "Copyright 2018-19, NCBS Bangalore" +__maintainer__ = "Dilawar Singh" +__email__ = "dilawars@ncbs.res.in" + +import sys +import hashlib +import moose._moose as _moose +import re + +import logging +logger_ = logging.getLogger('moose.utils.graph') + +pathPat_ = re.compile(r'.+?\[\d+\]$') + +def morphologyToGraphviz(filename=None, pat='/##[TYPE=Compartment]'): + '''Write Electrical network to a dot graph. + + Params: + + :filename: Default None. Write graphviz file to this path. If None, write to + stdout. + :pat: Compartment path. By default, search for all moose.Compartment. + ''' + + def fix(path): + '''Fix a given path so it can be written to a graphviz file''' + # If no [0] is at end of the path then append it. + global pathPat_ + if not pathPat_.match(path): + path = path + '[0]' + return path + + compList = _moose.wildcardFind(pat) + if not compList: + logger_.warn("No compartment found") + + dot = [] + dot.append("digraph G {") + dot.append("\tconcentrate=true;") + for c in compList: + if c.neighbors['raxial']: + for n in c.neighbors['raxial']: + lhs = fix(c.path) + rhs = fix(n.path) + dot.append('\t"{}" -> "{}";'.format(lhs, rhs)) + elif c.neighbors['axial']: + for n in c.neighbors['axial']: + lhs = fix(c.path) + rhs = fix(n.path) + dot.append('\t"{}" -> "{}" [dir=back];'.format(lhs, rhs)) + else: + p = fix(c.path) + dot.append('\t"{}"'.format(p)) + + dot.append('}') + dot = '\n'.join(dot) + if not filename: + print(dot) + else: + with open(filename, 'w') as graphviz: + logger_.info("Writing compartment topology to file {}".format(filename)) + graphviz.write(dot) + return True + + +def chemicalReactionNetworkToGraphviz(compt, path=None): + """Write chemical reaction network to a graphviz file. + + :param compt: Given compartment. + :param filepath: Save to this filepath. If None, write to stdout. + """ + dot = _crn(compt) + if path is None: + print(dot, file=sys.stdout) + return + with open(path, 'w') as f: + f.write(dot) + +# aliases +crnToDot = chemicalReactionNetworkToGraphviz +crnToGraphviz = chemicalReactionNetworkToGraphviz + +def _fixLabel(name): + name = name.replace('*', 'star') + name = name.replace('.', '_') + return "\"{}\"".format(name) + +def _addNode(n, nodes, dot): + node = hashlib.sha224(n.path.encode()).hexdigest() + nodeType = 'pool' + if isinstance(n, _moose.Reac) or isinstance(n, _moose.ZombieReac): + node = 'r'+node + nodeType = 'reac' + else: + node = 'p'+node + if node in nodes: + return node + + nLabel = n.name + if nodeType == 'reac': + nLabel = "kf=%g kb=%g"%(n.Kf, n.Kb) + dot.append('\t%s [label="%s", kf=%g, kb=%g, shape=rect]' % ( + node, nLabel, n.Kf, n.Kb)) + else: + dot.append('\t%s [label="%s", concInit=%g]' % ( + node, nLabel, n.concInit)) + return node + +def _crn(compt): + nodes = {} + reacs = _moose.wildcardFind(compt.path+'/##[TYPE=Reac]') + reacs += _moose.wildcardFind(compt.path+'/##[TYPE=ZombieReac]') + dot = ['digraph %s {\n\t overlap=false' % compt.name ] + for r in reacs: + rNode = _addNode(r, nodes, dot) + for s in r.neighbors['sub']: + sNode = _addNode(s, nodes, dot) + dot.append('\t%s -> %s' % (sNode, rNode)) + for p in r.neighbors['prd']: + pNode = _addNode(p, nodes, dot) + dot.append('\t%s -> %s' % (rNode, pNode)) + dot.append('}') + return '\n'.join(dot) diff --git a/python/moose/plot_utils.py b/python/moose/plot_utils.py index eeeec611e9..97d3ee0fe8 100644 --- a/python/moose/plot_utils.py +++ b/python/moose/plot_utils.py @@ -125,14 +125,23 @@ def plotTable(table, **kwargs): plt.legend(loc='best') def plotTables(tables, outfile=None, **kwargs): - """Plot a list of tables onto one figure only. + """Plot a dict of tables. + + Each axis will be labeled with dict keys. The dict values must be + moose.Table/moose.Table2. + + :param outfile: Default `None`. If given, plot will be saved to this filepath. + :param grid: A tuple with (cols, rows), default is (len(tables)//2+1, 2) + :param figsize: Size of figure (W, H) in inches. Default (10, 1.5*len(tables)) """ assert type(tables) == dict, "Expected a dict of moose.Table" - plt.figure(figsize=(10, 1.5*len(tables))) + plt.figure(figsize=kwargs.get('figsize', (10, 1.5*len(tables)))) subplot = kwargs.get('subplot', True) + gridSize = kwargs.get('grid', (len(tables)//2+1, 2)) for i, tname in enumerate(tables): if subplot: - plt.subplot(len(tables), 1, i+1) + assert gridSize[0] <= 9 + plt.subplot(100*gridSize[0]+10*gridSize[1]+(i+1)) yvec = tables[tname].vector xvec = np.linspace(0, moose.Clock('/clock').currentTime, len(yvec)) plt.plot(xvec, yvec, label=tname) @@ -149,9 +158,7 @@ def plotTables(tables, outfile=None, **kwargs): try: plt.savefig(outfile, transparent=True) except Exception as e: - pu.dump("WARN" - , "Failed to save figure, plotting onto a window" - ) + pu.dump("WARN", "Failed to save figure. Errror %s"%e) plt.show() else: plt.show() diff --git a/python/moose/utils.py b/python/moose/utils.py index 5facab9f0f..4ed1402e7d 100644 --- a/python/moose/utils.py +++ b/python/moose/utils.py @@ -17,16 +17,21 @@ import re import logging -logger = logging.getLogger('moose') +logger_ = logging.getLogger('moose.utils') from moose.moose_constants import * from moose.print_utils import * +try: + from moose.network_utils import * +except Exception as e: + logger_.warn("Netowrk utilities are not loaded due to %s" % e) + # Print and Plot utilities. try: from moose.plot_utils import * except Exception as e: - info( "Plot utilities are not loaded due to '%s'" % e ) + logger_.warn( "Plot utilities are not loaded due to '%s'" % e ) def create_table_path(model, graph, element, field): @@ -523,12 +528,15 @@ def assignDefaultTicks(modelRoot='/model', dataRoot='/data', solver='hsolve'): if len(tab.neighbors['input']) == 0: moose.useClock(9, tab.path, 'process') -def stepRun(simtime, steptime, verbose=True): +def stepRun(simtime, steptime, verbose=True, logger=None): """Run the simulation in steps of `steptime` for `simtime`.""" + global logger_ + if logger is None: + logger = logger_ clock = moose.element('/clock') if verbose: msg = 'Starting simulation for %g' % (simtime) - logger.info(msg) + logger_.info(msg) ts = datetime.now() while clock.currentTime < simtime - steptime: ts1 = datetime.now() @@ -537,20 +545,20 @@ def stepRun(simtime, steptime, verbose=True): td = te - ts1 if verbose: msg = 'Simulated till %g. Left: %g. %g of simulation took: %g s' % (clock.currentTime, simtime - clock.currentTime, steptime, td.days * 86400 + td.seconds + 1e-6 * td.microseconds) - logger.info(msg) + logger_.info(msg) remaining = simtime - clock.currentTime if remaining > 0: if verbose: msg = 'Running the remaining %g.' % (remaining) - logger.info(msg) + logger_.info(msg) moose.start(remaining) te = datetime.now() td = te - ts dt = min([t for t in moose.element('/clock').dts if t > 0.0]) if verbose: msg = 'Finished simulation of %g with minimum dt=%g in %g s' % (simtime, dt, td.days * 86400 + td.seconds + 1e-6 * td.microseconds) - logger.info(msg) + logger_.info(msg) diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py index d9b06a04ec..ee6105d700 100644 --- a/python/rdesigneur/rdesigneur.py +++ b/python/rdesigneur/rdesigneur.py @@ -112,6 +112,7 @@ def __init__(self, stimList = [], plotList = [], # elecpath, geom_expr, object, field, title ['wave' [min max]] moogList = [], + ode_method = "gsl", # gsl, lsoda, gssa, gillespie params = None ): """ Constructor of the rdesigner. This just sets up internal fields @@ -121,6 +122,7 @@ def __init__(self, self.modelPath = modelPath self.turnOffElec = turnOffElec self.useGssa = useGssa + self.ode_method = ode_method self.combineSegments = combineSegments self.stealCellFromLibrary = stealCellFromLibrary self.verbose = verbose @@ -1348,6 +1350,7 @@ def _configureSolvers( self ) : raise BuildError( "configureSolvers: no chem meshes defined." ) fixXreacs.fixXreacs( self.chemid.path ) dmksolve = moose.Ksolve( self.dendCompt.path + '/ksolve' ) + dmksolve.method = self.ode_method dmdsolve = moose.Dsolve( self.dendCompt.path + '/dsolve' ) dmstoich = moose.Stoich( self.dendCompt.path + '/stoich' ) dmstoich.compartment = self.dendCompt @@ -1361,6 +1364,7 @@ def _configureSolvers( self ) : smksolve = moose.Gsolve( self.spineCompt.path + '/ksolve' ) else: smksolve = moose.Ksolve( self.spineCompt.path + '/ksolve' ) + smksolve.method = self.ode_method smdsolve = moose.Dsolve( self.spineCompt.path + '/dsolve' ) smstoich = moose.Stoich( self.spineCompt.path + '/stoich' ) smstoich.compartment = self.spineCompt @@ -1372,6 +1376,7 @@ def _configureSolvers( self ) : pmksolve = moose.Gsolve( self.psdCompt.path + '/ksolve' ) else: pmksolve = moose.Ksolve( self.psdCompt.path + '/ksolve' ) + pmksolve.method = self.ode_method pmdsolve = moose.Dsolve( self.psdCompt.path + '/dsolve' ) pmstoich = moose.Stoich( self.psdCompt.path + '/stoich' ) pmstoich.compartment = self.psdCompt diff --git a/python/rdesigneur/rmoogli.py b/python/rdesigneur/rmoogli.py index f673dad7ca..c34060f060 100644 --- a/python/rdesigneur/rmoogli.py +++ b/python/rdesigneur/rmoogli.py @@ -1,4 +1,4 @@ -from __future__ import absolute_import, print_function +from __future__ import absolute_import, print_function, division # rmoogli.py: rdesigneur Moogli interface # This is a fallback version designed to work with moogul but using almost # the same interface as far as rdesigneur is concerned. diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp index 944a5bc2f9..e905d705c7 100644 --- a/scheduling/Clock.cpp +++ b/scheduling/Clock.cpp @@ -700,13 +700,10 @@ void Clock::buildTicks( const Eref& e ) void Clock::handleStart( const Eref& e, double runtime, bool notify ) { notify_ = notify; - if ( stride_ == 0 || stride_ == ~0U ) stride_ = 1; unsigned long n = round( runtime / ( stride_ * dt_ ) ); - handleStep( e, n ); - } void Clock::handleStep( const Eref& e, unsigned long numSteps ) @@ -911,7 +908,7 @@ void Clock::buildDefaultTick() defaultTick_["HHChannel2D"] = 2; defaultTick_["Leakage"] = 2; defaultTick_["MarkovChannel"] = 2; - defaultTick_["MarkovGslSolver"] = 2; + defaultTick_["MarkovOdeSolver"] = 2; defaultTick_["MarkovRateTable"] = 2; defaultTick_["MarkovSolver"] = 2; defaultTick_["MarkovSolverBase"] = 2; diff --git a/setup.py b/setup.py index 8c3e46be45..fb2d67c6dd 100644 --- a/setup.py +++ b/setup.py @@ -83,7 +83,6 @@ def run(self): self.spawn(["ctest", "--output-on-failure", '-j%d'%numCores_]) os.chdir(sdir_) - class build_ext(_build_ext): user_options = [ ('with-boost', None, 'Use Boost Libraries (OFF)') diff --git a/tests/data/acc27.g b/tests/data/acc27.g deleted file mode 100644 index c2105ce77f..0000000000 --- a/tests/data/acc27.g +++ /dev/null @@ -1,325 +0,0 @@ -//genesis -// kkit Version 8 flat dumpfile - -// Saved on Tue Jan 22 17:41:13 2002 - -include kkit {argv 1} - -FASTDT = 5e-05 -SIMDT = 0.005 -CONTROLDT = 10 -PLOTDT = 5 -MAXTIME = 6000 -TRANSIENT_TIME = 2 -VARIABLE_DT_FLAG = 1 -DEFAULT_VOL = 1.6667e-21 -VERSION = 8.0 -setfield /file/modpath value /home2/bhalla/scripts/modules -kparms - -//genesis - -initdump -version 3 -ignoreorphans 1 -simobjdump table input output alloced step_mode stepsize x y z -simobjdump xtree path script namemode sizescale -simobjdump xcoredraw xmin xmax ymin ymax -simobjdump xtext editable -simobjdump xgraph xmin xmax ymin ymax overlay -simobjdump xplot pixflags script fg ysquish do_slope wy -simobjdump group xtree_fg_req xtree_textfg_req plotfield expanded movealone \ - link savename file version md5sum mod_save_flag x y z -simobjdump kpool CoTotal CoInit Co n nInit nTotal nMin vol slave_enable notes \ - xtree_fg_req xtree_textfg_req x y z -simobjdump kreac kf kb notes xtree_fg_req xtree_textfg_req x y z -simobjdump kenz CoComplexInit CoComplex nComplexInit nComplex vol k1 k2 k3 \ - keepconc usecomplex notes xtree_fg_req xtree_textfg_req link x y z -simobjdump stim level1 width1 delay1 level2 width2 delay2 baselevel trig_time \ - trig_mode notes xtree_fg_req xtree_textfg_req is_running x y z -simobjdump xtab input output alloced step_mode stepsize notes editfunc \ - xtree_fg_req xtree_textfg_req baselevel last_x last_y is_running x y z -simobjdump kchan perm gmax Vm is_active use_nernst notes xtree_fg_req \ - xtree_textfg_req x y z -simobjdump proto x y z -simobjdump linkinfo xtree_fg_req xtree_textfg_req uplink downlink x y z -simobjdump uplink xtree_fg_req xtree_textfg_req x y z -simobjdump downlink xtree_fg_req xtree_textfg_req x y z -simobjdump mirror notes xtree_fg_req x y z -simundump group /kinetics/MAPK 0 yellow black x 0 0 "" MAPK \ - /home2/bhalla/scripts/modules/MAPK_0.g 0 0 0 1 10 0 -simundump kpool /kinetics/MAPK/MAPK 0 1 0.3 0.3 0.3 0.3 1 0 1 0 "" 35 yellow \ - -8 -7 0 -simundump kpool /kinetics/MAPK/MKKK 0 0.1 0.1 0.1 0.1 0.1 0.1 0 1 0 "" 16 \ - yellow -8 5 0 -simundump kpool /kinetics/MAPK/MKK 0 0.3 0.3 0.3 0.3 0.3 0.3 0 1 0 "" 60 \ - yellow -8 -1 0 -simundump kpool /kinetics/MAPK/int1 0 1 0.001 0.001 0.001 0.001 1 0 1 0 "" 30 \ - yellow -4 4 0 -simundump kenz /kinetics/MAPK/int1/2 0 0 0 0 0 0.001 156.25 1 0.25 0 1 "" red \ - 30 "" -4 5 0 -simundump kpool /kinetics/MAPK/MKKK-P 0 0.05 0 0 0 0 0.05 0 1 0 "" 51 yellow \ - 0 5 0 -simundump kenz /kinetics/MAPK/MKKK-P/3 0 0 0 0 0 0.001 8.3333 0.1 0.025 0 1 \ - "" red 51 "" -4 2 0 -simundump kenz /kinetics/MAPK/MKKK-P/4 0 0 0 0 0 0.001 8.3333 0.1 0.025 0 1 \ - "" red 51 "" 4 2 0 -simundump kpool /kinetics/MAPK/int3 0 1 0.001 0.001 0.001 0.001 1 0 1 0 "" \ - blue yellow -4 -2 0 -simundump kenz /kinetics/MAPK/int3/6 0 0 0 0 0 0.001 250 3 0.75 0 1 "" red \ - blue "" -4 -1 0 -simundump kpool /kinetics/MAPK/int5 0 1 0.001 0.001 0.001 0.001 1 0 1 0 "" 1 \ - yellow -4 -8 0 -simundump kenz /kinetics/MAPK/int5/10 0 0 0 0 0 0.001 166.67 2 0.5 0 1 "" red \ - 1 "" -4 -7 0 -simundump kpool /kinetics/MAPK/MKK-P 0 0.1 0 0 0 0 0.1 0 1 0 "" 5 yellow 0 -1 \ - 0 -simundump kpool /kinetics/MAPK/MAPK-P 0 0.1 0 0 0 0 0.1 0 1 0 "" 55 yellow 0 \ - -7 0 -simundump kpool /kinetics/MAPK/int2 0 1 0.001 0.001 0.001 0.001 1 0 1 0 "" 2 \ - yellow 4 -2 0 -simundump kenz /kinetics/MAPK/int2/5 0 0 0 0 0 0.001 250 3 0.75 0 1 "" red 2 \ - "" 4 -1 0 -simundump kpool /kinetics/MAPK/int4 0 1 0.001 0.001 0.001 0.001 1 0 1 0 "" 17 \ - yellow 4 -8 0 -simundump kenz /kinetics/MAPK/int4/9 0 0 0 0 0 0.001 166.67 2 0.5 0 1 "" red \ - 17 "" 4 -7 0 -simundump kpool /kinetics/MAPK/Ras-MKKKK 0 0.1 0.001 0.001 0.001 0.001 0.1 0 \ - 1 0 "" 47 yellow 6 8 0 -simundump kenz /kinetics/MAPK/Ras-MKKKK/1 0 0 0 0 0 0.001 1250 10 2.5 0 1 "" \ - red 47 "" -4 8 0 -simundump kpool /kinetics/MAPK/inactiveRas-MKKK 0 0 0 0 0 0 0 0 1 0 "" 30 \ - yellow 11 8 0 -simundump kreac /kinetics/MAPK/Neg_feedback 0 1 0.009 "" white yellow 11 2 0 -simundump kpool /kinetics/MAPK/MKK-PP 0 0.1 0 0 0 0 0.1 0 1 0 "" 60 yellow 8 \ - -1 0 -simundump kenz /kinetics/MAPK/MKK-PP/7 0 0 0 0 0 0.001 8.3333 0.1 0.025 0 1 \ - "" red 60 "" -4 -4 0 -simundump kenz /kinetics/MAPK/MKK-PP/8 0 0 0 0 0 0.001 8.3333 0.1 0.025 0 1 \ - "" red 60 "" 4 -4 0 -simundump kpool /kinetics/MAPK/MAPK-PP 0 0.1 0 0 0 0 0.1 0 1 0 "" 46 yellow 8 \ - -7 0 -simundump xgraph /graphs/conc1 0 0 6000 0 0.3 0 -simundump xgraph /graphs/conc2 0 0 6000 4.5157e-05 0.3 0 -simundump xplot /graphs/conc1/MAPK-PP.Co 3 524288 \ - "delete_plot.w ; edit_plot.D " 46 0 0 1 -simundump xplot /graphs/conc1/MAPK.Co 3 524288 \ - "delete_plot.w ; edit_plot.D " 35 0 0 1 -simundump xplot /graphs/conc2/Ras-MKKKK.Co 3 524288 \ - "delete_plot.w ; edit_plot.D " 47 0 0 1 -simundump xplot /graphs/conc2/MAPK.Co 3 524288 \ - "delete_plot.w ; edit_plot.D " 35 0 0 1 -simundump xplot /graphs/conc2/MKKK.Co 3 524288 \ - "delete_plot.w ; edit_plot.D " 16 0 0 1 -simundump xplot /graphs/conc2/MKK.Co 3 524288 \ - "delete_plot.w ; edit_plot.D " 60 0 0 1 -simundump xgraph /moregraphs/conc3 0 0 6000 0 1 0 -simundump xgraph /moregraphs/conc4 0 0 6000 0 1 0 -simundump xcoredraw /edit/draw 0 -6.5262 17.834 -8.7578 11.542 -simundump xtree /edit/draw/tree 0 \ - /kinetics/#[],/kinetics/#[]/#[],/kinetics/#[]/#[]/#[][TYPE!=proto],/kinetics/#[]/#[]/#[][TYPE!=linkinfo]/##[] \ - "edit_elm.D ; drag_from_edit.w " auto 0.6 -simundump xtext /file/notes 0 1 -xtextload /file/notes \ -"22 Jan 2002" \ -" " \ -" This model is based on Kholodenko, B.N." \ -" Eur. J. Biochem. 267, 1583-1588(2000)" \ -"" -addmsg /kinetics/MAPK/MKK-PP/7 /kinetics/MAPK/MAPK REAC sA B -addmsg /kinetics/MAPK/int5/10 /kinetics/MAPK/MAPK MM_PRD pA -addmsg /kinetics/MAPK/Ras-MKKKK/1 /kinetics/MAPK/MKKK REAC sA B -addmsg /kinetics/MAPK/int1/2 /kinetics/MAPK/MKKK MM_PRD pA -addmsg /kinetics/MAPK/MKKK-P/3 /kinetics/MAPK/MKK REAC sA B -addmsg /kinetics/MAPK/int3/6 /kinetics/MAPK/MKK MM_PRD pA -addmsg /kinetics/MAPK/int1 /kinetics/MAPK/int1/2 ENZYME n -addmsg /kinetics/MAPK/MKKK-P /kinetics/MAPK/int1/2 SUBSTRATE n -addmsg /kinetics/MAPK/Ras-MKKKK/1 /kinetics/MAPK/MKKK-P MM_PRD pA -addmsg /kinetics/MAPK/int1/2 /kinetics/MAPK/MKKK-P REAC sA B -addmsg /kinetics/MAPK/MKKK-P /kinetics/MAPK/MKKK-P/3 ENZYME n -addmsg /kinetics/MAPK/MKK /kinetics/MAPK/MKKK-P/3 SUBSTRATE n -addmsg /kinetics/MAPK/MKKK-P /kinetics/MAPK/MKKK-P/4 ENZYME n -addmsg /kinetics/MAPK/MKK-P /kinetics/MAPK/MKKK-P/4 SUBSTRATE n -addmsg /kinetics/MAPK/int3 /kinetics/MAPK/int3/6 ENZYME n -addmsg /kinetics/MAPK/MKK-P /kinetics/MAPK/int3/6 SUBSTRATE n -addmsg /kinetics/MAPK/int5 /kinetics/MAPK/int5/10 ENZYME n -addmsg /kinetics/MAPK/MAPK-P /kinetics/MAPK/int5/10 SUBSTRATE n -addmsg /kinetics/MAPK/MKKK-P/4 /kinetics/MAPK/MKK-P REAC sA B -addmsg /kinetics/MAPK/MKKK-P/3 /kinetics/MAPK/MKK-P MM_PRD pA -addmsg /kinetics/MAPK/int3/6 /kinetics/MAPK/MKK-P REAC sA B -addmsg /kinetics/MAPK/int2/5 /kinetics/MAPK/MKK-P MM_PRD pA -addmsg /kinetics/MAPK/MKK-PP/8 /kinetics/MAPK/MAPK-P REAC sA B -addmsg /kinetics/MAPK/MKK-PP/7 /kinetics/MAPK/MAPK-P MM_PRD pA -addmsg /kinetics/MAPK/int5/10 /kinetics/MAPK/MAPK-P REAC sA B -addmsg /kinetics/MAPK/int4/9 /kinetics/MAPK/MAPK-P MM_PRD pA -addmsg /kinetics/MAPK/int2 /kinetics/MAPK/int2/5 ENZYME n -addmsg /kinetics/MAPK/MKK-PP /kinetics/MAPK/int2/5 SUBSTRATE n -addmsg /kinetics/MAPK/int4 /kinetics/MAPK/int4/9 ENZYME n -addmsg /kinetics/MAPK/MAPK-PP /kinetics/MAPK/int4/9 SUBSTRATE n -addmsg /kinetics/MAPK/Neg_feedback /kinetics/MAPK/Ras-MKKKK REAC A B -addmsg /kinetics/MAPK/Ras-MKKKK /kinetics/MAPK/Ras-MKKKK/1 ENZYME n -addmsg /kinetics/MAPK/MKKK /kinetics/MAPK/Ras-MKKKK/1 SUBSTRATE n -addmsg /kinetics/MAPK/Neg_feedback /kinetics/MAPK/inactiveRas-MKKK REAC B A -addmsg /kinetics/MAPK/MAPK-PP /kinetics/MAPK/Neg_feedback SUBSTRATE n -addmsg /kinetics/MAPK/Ras-MKKKK /kinetics/MAPK/Neg_feedback SUBSTRATE n -addmsg /kinetics/MAPK/inactiveRas-MKKK /kinetics/MAPK/Neg_feedback PRODUCT n -addmsg /kinetics/MAPK/MKKK-P/4 /kinetics/MAPK/MKK-PP MM_PRD pA -addmsg /kinetics/MAPK/int2/5 /kinetics/MAPK/MKK-PP REAC sA B -addmsg /kinetics/MAPK/MKK-PP /kinetics/MAPK/MKK-PP/7 ENZYME n -addmsg /kinetics/MAPK/MAPK /kinetics/MAPK/MKK-PP/7 SUBSTRATE n -addmsg /kinetics/MAPK/MKK-PP /kinetics/MAPK/MKK-PP/8 ENZYME n -addmsg /kinetics/MAPK/MAPK-P /kinetics/MAPK/MKK-PP/8 SUBSTRATE n -addmsg /kinetics/MAPK/MKK-PP/8 /kinetics/MAPK/MAPK-PP MM_PRD pA -addmsg /kinetics/MAPK/int4/9 /kinetics/MAPK/MAPK-PP REAC sA B -addmsg /kinetics/MAPK/Neg_feedback /kinetics/MAPK/MAPK-PP REAC A B -addmsg /kinetics/MAPK/MAPK-PP /graphs/conc1/MAPK-PP.Co PLOT Co *MAPK-PP.Co *46 -addmsg /kinetics/MAPK/MAPK /graphs/conc1/MAPK.Co PLOT Co *MAPK.Co *35 -addmsg /kinetics/MAPK/Ras-MKKKK /graphs/conc2/Ras-MKKKK.Co PLOT Co *Ras-MKKKK.Co *47 -addmsg /kinetics/MAPK/MAPK /graphs/conc2/MAPK.Co PLOT Co *MAPK.Co *35 -addmsg /kinetics/MAPK/MKKK /graphs/conc2/MKKK.Co PLOT Co *MKKK.Co *16 -addmsg /kinetics/MAPK/MKK /graphs/conc2/MKK.Co PLOT Co *MKK.Co *60 -enddump -// End of dump - -call /kinetics/MAPK/notes LOAD \ -"This is the oscillatory MAPK model from Kholodenko 2000" \ -"Eur J. Biochem 267:1583-1588" \ -"The original model is formulated in terms of idealized" \ -"Michaelis-Menten enzymes and the enzyme-substrate complex" \ -"concentrations are therefore assumed negligible. The" \ -"current implementation of the model uses explicit enzyme" \ -"reactions involving substrates and is therefore an" \ -"approximation to the Kholodenko model. The approximation is" \ -"greatly improved if the enzyme is flagged as Available" \ -"which is an option in Kinetikit. This flag means that the" \ -"enzyme protein concentration is not reduced even when it" \ -"is involved in a complex. However, the substrate protein" \ -"continues to participate in enzyme-substrate complexes" \ -"and its concentration is therefore affected. Overall," \ -"this model works almost the same as the Kholodenko model" \ -"but the peak MAPK-PP amplitudes are a little reduced and" \ -"the period of oscillations is about 10% longer." \ -"If the enzymes are not flagged as Available then the" \ -"oscillations commence only when the Km for enzyme 1" \ -"is set to 0.1 uM." -call /kinetics/MAPK/MAPK/notes LOAD \ -"The total concn. of MAPK is 300nM " \ -"from" \ -"Kholodenko, 2000." -call /kinetics/MAPK/MKKK/notes LOAD \ -"The total concn. of MKKK is 100nM " \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/MKK/notes LOAD \ -"The total concn. of MKK is 300nM " \ -"from" \ -"Kholodenko,2000" -call /kinetics/MAPK/int1/notes LOAD \ -"This is the intermediate enzyme which catalyses the " \ -"dephosphorylation of MKKK-P to MKKK. The concentration" \ -"is set to 1 nM based on" \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/int1/2/notes LOAD \ -"Km is 8nM and Vmax is 0.25nM.s-1 " \ -"from" \ -"Kholodenko, 2000." -call /kinetics/MAPK/MKKK-P/notes LOAD \ -"This is the phosphorylated form of MKKK which converts MKK" \ -"to MKK-P and then to MKK-PP" \ -"from" \ -"Kholodenko, 2000." -call /kinetics/MAPK/MKKK-P/3/notes LOAD \ -"Km is 15 nM and Vmax is 0.025s-1" \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/MKKK-P/4/notes LOAD \ -"Km is 15nM and Vmax is 0.025s-1" \ -"from " \ -"Kholodenko, 2000." -call /kinetics/MAPK/int3/notes LOAD \ -"This intermediate enzyme catalyses the dephosphorylation of" \ -"MKK-P to MKK. The concentration is 1nM" \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/int3/6/notes LOAD \ -"The Km is 15nM and the Vmax is 0.75nM.s-1" \ -"from" \ -"Kholodenko 2000." -call /kinetics/MAPK/int5/notes LOAD \ -"This catalyses the conversion of MAPK-P to MAPK. The " \ -"concenration is 1nM." \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/int5/10/notes LOAD \ -"The Km is 15nM and Vmax is 0.5nM.s-1" \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/MKK-P/notes LOAD \ -"This is the single phoshorylated form of MKK." \ -"from" \ -"Kholodenko, 2000." -call /kinetics/MAPK/MAPK-P/notes LOAD \ -"This is the single phopshorylated form of MAPK" \ -"from" \ -"Kholodenko, 2000." -call /kinetics/MAPK/int2/notes LOAD \ -"This intermediate enzyme which catalyses the dephosphorylation of" \ -"MKK-PP to MKK-P. The concentration is 1nM." \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/int2/5/notes LOAD \ -"The Km is 15nM and Vmax is 0.75nM.s-1 " \ -"from" \ -"Kholodenko, 2000" \ -"" -call /kinetics/MAPK/int4/notes LOAD \ -"This intermediate enzyme catalyses the dephosphorylation of" \ -"MAPK-PP to MAPK-P. The concentration is 1nM." \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/int4/9/notes LOAD \ -"The Km is 15nM and Vmax is 0.5nM.s-1 " \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/Ras-MKKKK/notes LOAD \ -"The concn. of Ras-MKKKK* is set to 1 nM implicitly" \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/Ras-MKKKK/1/notes LOAD \ -"The Km is 10nM and Vmax is 2.5nM sec^-1. We assume that" \ -"there is 1 nM of the Ras-MKKKK." \ -"From Kholodenko, 2000." \ -"" \ -"If the enzymes are not flagged as Available, then this" \ -"Km should be set to 0.1 to obtain oscillations." -call /kinetics/MAPK/inactiveRas-MKKK/notes LOAD \ -"This is the inactive form of Ras-MKKK. Based on the" \ -"reaction scheme from Kholodenko 2000, this is equivalent" \ -"to a binding of the MAPK-PP to the Ras. The amount of" \ -"Ras in the model is small enough that negligible amounts" \ -"of MAPK are involved in this reaction. So it is a fair" \ -"approximation to the negative feedback mechanism from" \ -"Kholodenko, 2000." -call /kinetics/MAPK/Neg_feedback/notes LOAD \ -"From Kholodenko, 2000 Eur J Biochem 267" \ -"the Kd is 9 nM. We use a rather fast Kf of 1/sec/uM" \ -"so that equilibrium is maintained." \ -"" -call /kinetics/MAPK/MKK-PP/notes LOAD \ -"This is the double phosphorylated and active form of MKK" \ -"from" \ -"Kholodenko, 2000" -call /kinetics/MAPK/MKK-PP/7/notes LOAD \ -"The Km is 15nM which is 0.015uM Vmax is 0.025s-1" \ -"from" \ -"Kholodenko, 2000." \ -"" -call /kinetics/MAPK/MKK-PP/8/notes LOAD \ -"The Km is 15nM which is 0.015uM and Vmax is 0.025s-1" \ -"from" \ -"Kholodenko, 2000" \ -"" -call /kinetics/MAPK/MAPK-PP/notes LOAD \ -"This is the double phosphorylated and active form of MAPK." \ -"from" \ -"Kholodenko, 2000." -complete_loading diff --git a/tests/py_moose/test_ksolve.py b/tests/py_moose/test_ksolve.py index 4edabc933e..92f4ba2ff9 100644 --- a/tests/py_moose/test_ksolve.py +++ b/tests/py_moose/test_ksolve.py @@ -37,7 +37,7 @@ moose.reinit() print( '[INFO] Using method = %s' % ksolve.method ) t1 = time.time() -moose.start( 10 ) +moose.start( 100 ) print('[INFO] Time taken %s' % (time.time() - t1 )) expected = [ 7.77859 , 2.77858 , 2.27541 , 6.12141 , 7.77858 , 2.77858 , 2.27541 , 6.12141 , 7.77858 , 2.77858 , 2.27541 , 6.12141 , 7.77858 diff --git a/tests/py_moose/test_moose_attribs.py b/tests/py_moose/test_moose_attribs.py index 6bbc51bc85..636f77f4d6 100644 --- a/tests/py_moose/test_moose_attribs.py +++ b/tests/py_moose/test_moose_attribs.py @@ -24,7 +24,7 @@ 'HHGate2D', 'HSolve', 'INFINITE', 'IntFire', 'IntFireBase', 'Interpol', 'Interpol2D', 'IzhIF', 'IzhikevichNrn', 'Ksolve', 'LIF', 'Leakage', 'LookupField', 'MMPump', 'MMenz', 'MarkovChannel' , - # 'MarkovGslSolver', # This is GSL specific. + 'MarkovOdeSolver', 'MarkovRateTable', 'MarkovSolver', 'MarkovSolverBase', 'MeshEntry', 'MgBlock', 'Msg', 'Mstring', 'NMDAChan', 'Nernst', 'NeuroMesh', 'Neuron', 'Neutral', 'OneToAllMsg', 'OneToOneDataIndexMsg', @@ -59,8 +59,8 @@ def main(): global attribs for at in attribs: + print( "\tTesting for attrib %s" % at ) assert hasattr( moose, at ), 'Attrib %s not found' % at - print( getattr(moose, at ) ) if __name__ == '__main__': main() diff --git a/tests/python/test_cylinder_diffusion_gsolve+dsolve.py b/tests/python/test_cylinder_diffusion_gsolve+dsolve.py new file mode 100644 index 0000000000..db3dd77344 --- /dev/null +++ b/tests/python/test_cylinder_diffusion_gsolve+dsolve.py @@ -0,0 +1,146 @@ +######################################################################### +## This program is part of 'MOOSE', the +## Messaging Object Oriented Simulation Environment. +## Copyright (C) 2014 Upinder S. Bhalla. and NCBS +## It is made available under the terms of the +## GNU Lesser General Public License version 2.1 +## See the file COPYING.LIB for the full notice. +######################################################################### + + +import sys +import math +import numpy as np +import moose +print('[INFO] Using moose from %s' % moose.__file__ ) +import os + +moose.seed( 1 ) + +concA = 0.005 # millimolar +def makeModel(): + # create container for model + r0 = 2e-6 # m + r1 = 1e-6 # m + num = 100 + diffLength = 1e-6 # m + L = num * diffLength # m + diffConst = 10e-12 + #motorRate = 1e-6 + #diffConst = 0 + motorRate = 0 + + model = moose.Neutral( 'model' ) + compartment = moose.CylMesh( '/model/compartment' ) + compartment.r0 = r0 + compartment.r1 = r1 + compartment.x0 = 0 + compartment.x1 = L + compartment.diffLength = diffLength + assert( compartment.numDiffCompts == num ) + + # create molecules and reactions + a = moose.Pool( '/model/compartment/a' ) + b = moose.Pool( '/model/compartment/b' ) + c = moose.Pool( '/model/compartment/c' ) + d = moose.Pool( '/model/compartment/d' ) + r1 = moose.Reac( '/model/compartment/r1' ) + moose.connect( r1, 'sub', b, 'reac' ) + moose.connect( r1, 'sub', d, 'reac' ) + moose.connect( r1, 'prd', c, 'reac' ) + r1.Kf = 1000.0 # 1/(mM.sec) + r1.Kb = 1 # 1/sec + + # Assign parameters + a.diffConst = diffConst + b.diffConst = diffConst / 2.0 + b.motorConst = motorRate + c.diffConst = diffConst + d.diffConst = diffConst + + # Make solvers + ksolve = moose.Gsolve( '/model/compartment/ksolve' ) + dsolve = moose.Dsolve( '/model/compartment/dsolve' ) + stoich = moose.Stoich( '/model/compartment/stoich' ) + stoich.compartment = compartment + stoich.ksolve = ksolve + stoich.dsolve = dsolve + stoich.path = "/model/compartment/##" + assert( dsolve.numPools == 4 ) + a.vec.concInit = concA + b.vec.concInit = concA / 5.0 + c.vec.concInit = concA + d.vec.concInit = concA / 5.0 + for i in range( num ): + d.vec[i].concInit = concA * 2 * i / num + +def main(): + """ + This example illustrates how to set up a diffusion/transport model with + a simple reaction-diffusion system in a tapering cylinder: + + | Molecule **a** diffuses with diffConst of 10e-12 m^2/s. + | Molecule **b** diffuses with diffConst of 5e-12 m^2/s. + | Molecule **b** also undergoes motor transport with a rate of 10e-6 m/s + | Thus it 'piles up' at the end of the cylinder. + | Molecule **c** does not move: diffConst = 0.0 + | Molecule **d** does not move: diffConst = 10.0e-12 but it is buffered. + | Because it is buffered, it is treated as non-diffusing. + + All molecules other than **d** start out only in the leftmost (first) + voxel, with a concentration of 1 mM. **d** is present throughout + at 0.2 mM, except in the last voxel, where it is at 1.0 mM. + + The cylinder has a starting radius of 2 microns, and end radius of + 1 micron. So when the molecule undergoing motor transport gets to the + narrower end, its concentration goes up. + + There is a little reaction in all compartments: ``b + d <===> c`` + + As there is a high concentration of **d** in the last compartment, + when the molecule **b** reaches the end of the cylinder, the reaction + produces lots of **c**. + + Note that molecule **a** does not participate in this reaction. + + The concentrations of all molecules are displayed in an animation. + """ + runtime = 20.0 + diffdt = 0.005 + plotdt = 0.1 + makeModel() + # Set up clocks. The dsolver to know before assigning stoich + moose.setClock( 10, diffdt ) # 10 is the standard clock for Dsolve. + moose.setClock( 16, plotdt ) # 16 is the standard clock for Ksolve. + + a = moose.element( '/model/compartment/a' ) + b = moose.element( '/model/compartment/b' ) + c = moose.element( '/model/compartment/c' ) + d = moose.element( '/model/compartment/d' ) + + moose.reinit() + atot = sum( a.vec.n ) + btot = sum( b.vec.n ) + ctot = sum( c.vec.n ) + dtot = sum( d.vec.n ) + for t in np.arange( 0, runtime, plotdt ): + moose.start( plotdt ) + + atot2 = sum( a.vec.n ) + btot2 = sum( b.vec.n ) + ctot2 = sum( c.vec.n ) + dtot2 = sum( d.vec.n ) + + msg = 'Ratio of initial to final total numbers of ' + got = np.array((atot2/atot, btot2/btot, ctot2/ctot, dtot2/dtot)) + msg += 'a=%f b=%f, c=%f, d=%f' % (tuple(got)) + print(msg) + print('Initial to final (b+c)=%f' % (float(btot2 + ctot2) / (btot + ctot ))) + expected = np.array((1.00003087, 1.39036644, 0.92191184, 1.11427514)) + error = got - expected + rerror = np.abs( error ) / expected + assert np.allclose(got, expected, atol=1e-3), "Got %s, expected %s" % (got, expected) + +# Run the 'main' if this script is executed standalone. +if __name__ == '__main__': + main() diff --git a/tests/python/test_gsolve_parallel.py b/tests/python/test_gsolve_parallel.py new file mode 100644 index 0000000000..e7887ce853 --- /dev/null +++ b/tests/python/test_gsolve_parallel.py @@ -0,0 +1,83 @@ +import numpy as np +import moose +print( '[INFO] Using moose from %s' % moose.__file__ ) +import time +moose.seed( 10 ) + +def main( nT ): + """ + This example implements a reaction-diffusion like system which is + bistable and propagates losslessly. It is based on the NEURON example + rxdrun.py, but incorporates more compartments and runs for a longer time. + The system is implemented in a function rather than as a proper system + of chemical reactions. Please see rxdReacDiffusion.py for a variant that + uses a reaction plus a function object to control its rates. + """ + print( 'Using %d threads' % nT ) + + dt = 0.1 + + # define the geometry + compt = moose.CylMesh( '/cylinder' ) + compt.r0 = compt.r1 = 100e-9 + compt.x1 = 200e-9 + compt.diffLength = 0.2e-9 + assert( compt.numDiffCompts == compt.x1/compt.diffLength ) + + #define the molecule. Its geometry is defined by its parent volume, cylinder + c = moose.Pool( '/cylinder/pool' ) + c.diffConst = 1e-13 # define diffusion constant + + + # Here we set up a function calculation + func = moose.Function( '/cylinder/pool/func' ) + func.expr = "(-x0 * (30e-9 - x0) * (100e-9 - x0))*0.0001" + func.x.num = 1 #specify number of input variables. + + #Connect the molecules to the func + moose.connect( c, 'nOut', func.x[0], 'input' ) + #Connect the function to the pool + moose.connect( func, 'valueOut', c, 'increment' ) + + #Set up solvers + ksolve = moose.Gsolve( '/cylinder/Gsolve' ) + try: + ksolve.numThreads = nT + except Exception as e: + print( 'OLD MOOSE. Does not support multithreading' ) + dsolve = moose.Dsolve( '/cylinder/dsolve' ) + stoich = moose.Stoich( '/cylinder/stoich' ) + stoich.compartment = compt + stoich.ksolve = ksolve + stoich.dsolve = dsolve + stoich.path = '/cylinder/##' + + #initialize + x = np.arange( 0, compt.x1, compt.diffLength ) + c.vec.nInit = [ 1000.0 for q in x ] + + # Run and plot it. + moose.reinit() + updateDt = 50 + runtime = updateDt * 10 + t1 = time.time() + res = [] + for t in range( 0, runtime-1, updateDt ): + moose.start( updateDt ) + y = c.vec.n + res.append( (np.mean(y), np.std(y)) ) + + expected = [(9.0, 0.0), (6.0, 0.0), (5.0, 0.0), (3.0, 0.0), (2.0, 0.0), + (2.0, 0.0), (2.0, 0.0), (1.0, 0.0), (1.0, 0.0), (1.0, 0.0)] + print(("Time = ", time.time() - t1)) + print( res ) + assert res == expected + +if __name__ == '__main__': + import sys + import multiprocessing + nT = int(multiprocessing.cpu_count() / 2) + if len(sys.argv) > 1: + nT = int(sys.argv[1]) + main( nT ) + diff --git a/tests/python/test_ksolve_parallel.py b/tests/python/test_ksolve_parallel.py new file mode 100644 index 0000000000..72da77d781 --- /dev/null +++ b/tests/python/test_ksolve_parallel.py @@ -0,0 +1,105 @@ +import os +import numpy as np +import moose +print( 'Using moose from %s' % moose.__file__ ) +import time +os.environ['MOOSE_SHOW_PROFILING_INFO'] = '1' + +def main( nthreads = 1 ): + """ + This example implements a reaction-diffusion like system which is + bistable and propagates losslessly. It is based on the NEURON example + rxdrun.py, but incorporates more compartments and runs for a longer time. + The system is implemented as a hybrid of a reaction and a function which + sets its rates. Please see rxdFuncDiffusion.py for a variant that uses + just a function object to set up the system. + """ + + dt = 0.1 + + # define the geometry + compt = moose.CylMesh( '/cylinder' ) + compt.r0 = compt.r1 = 1 + # compt.diffLength = 1e-6 + compt.x1 = 1e-2 + assert compt.numDiffCompts == int(compt.x1/compt.diffLength), ( + compt.numDiffCompts, compt.x1 / compt.diffLength ) + print( 'No of compartment %d' % compt.numDiffCompts ) + + #define the molecule. Its geometry is defined by its parent volume, cylinder + c = moose.Pool( '/cylinder/pool' ) + c.diffConst = 1 # define diffusion constant + # There is an implicit reaction substrate/product. MOOSE makes it explicit. + buf = moose.BufPool( '/cylinder/buf' ) + buf.nInit = 1 + + # The reaction is something entirely peculiar, not a chemical thing. + reaction = moose.Reac( '/cylinder/reac' ) + reaction.Kb = 0 + + # so here we set up a function calculation to do the same thing. + func = moose.Function( '/cylinder/reac/func' ) + func.expr = "(1 - x0) * (0.3 - x0)" + func.x.num = 1 #specify number of input variables. + + #Connect the reaction to the pools + moose.connect( reaction, 'sub', c, 'reac' ) + moose.connect( reaction, 'prd', buf, 'reac' ) + + #Connect the function to the reaction + moose.connect( func, 'valueOut', reaction, 'setNumKf' ) + + #Connect the molecules to the func + moose.connect( c, 'nOut', func.x[0], 'input' ) + + #Set up solvers + ksolve = moose.Ksolve( '/cylinder/ksolve' ) + ksolve.numThreads = nthreads + dsolve = moose.Dsolve( '/cylinder/dsolve' ) + stoich = moose.Stoich( '/cylinder/stoich' ) + stoich.compartment = compt + stoich.ksolve = ksolve + stoich.dsolve = dsolve + stoich.path = '/cylinder/##' + for i in range( 10, 18 ): + moose.setClock( i, dt ) + + #initialize + x = np.arange( 0, compt.x1, compt.diffLength ) + c.vec.nInit = [ (q < 0.2 * compt.x1) for q in x ] + + expected = [ (0.2, 0.40000000000000013) + , (2.6704795776286974e-07, 1.2678976830753021e-17) + , (8.167639617309419e-14, 3.8777269301457245e-24) + , (2.498062905267963e-20, 1.1860363878961374e-30) + , (7.64029581501609e-27, 3.6273808003690943e-37) + ] + + # Run and plot it. + moose.reinit() + updateDt = 50 + runtime = updateDt * 4 + # plt = pylab.plot( x, c.vec.n, label='t = 0 ') + yvec = c.vec.n + u1, m1 = np.mean( yvec ), np.std( yvec ) + print( u1, m1 ) + assert np.isclose( (u1, m1), expected[0], atol=1e-5).all(), expected[0] + t1 = time.time() + for i, t in enumerate(range( 0, runtime-1, updateDt)): + moose.start( updateDt ) + # plt = pylab.plot( x, c.vec.n, label='t = '+str(t + updateDt) ) + yvec = c.vec.n + u1, m1 = np.mean( yvec ), np.std( yvec ) + print( u1, m1 ) + np.isclose( (u1,m1), expected[i+1], atol=1e-5 ).all(), expected[i+1] + return time.time() - t1 + + +if __name__ == '__main__': + import multiprocessing + import sys + nT = int(multiprocessing.cpu_count()/2) + if len(sys.argv) > 1: + nT = int(sys.argv[1]) + t2 = main( nT ) + print( 'With threads=%d, Time=%g' % (nT, t2)) diff --git a/tests/python/test_lsoda.py b/tests/python/test_lsoda.py new file mode 100644 index 0000000000..094af5da50 --- /dev/null +++ b/tests/python/test_lsoda.py @@ -0,0 +1,53 @@ +# -*- coding: utf-8 -*- +import sys +import numpy as np +import time +import moose +print('Using moose from %s' % moose.__file__ ) + +compt = moose.CubeMesh( '/compt' ) +compt.volume = 1e-20 + +pools = [] +for r in range( 10 ): + a1 = moose.Pool( '/compt/a1%s' % r ) + a1.concInit = 10 + a2 = moose.Pool( '/compt/a2%s' % r ) + a2.concInit = 5 + b1 = moose.Pool( '/compt/b1%s' % r ) + b1.concInit = 0.054 + b2 = moose.Pool( '/compt/b2%s' % r ) + b2.concInit = 3.9 + r = moose.Reac( '/compt/reac%s'% r ) + moose.connect( r, 'sub', a1, 'reac' ) + moose.connect( r, 'sub', a2, 'reac' ) + moose.connect( r, 'prd', b1, 'reac' ) + moose.connect( r, 'prd', b2, 'reac' ) + r.Kf = 2.9 + r.Kb = 4.5 + pools += [ a1, a2, b1, b2 ] + +ksolve = moose.Ksolve( '/compt/ksolve' ) +ksolve.method = 'lsoda' +stoich = moose.Stoich( '/compt/stoich' ) +stoich.compartment = compt +stoich.ksolve = ksolve +stoich.path = '/compt/##' +moose.reinit() +print( '[INFO] Using method = %s' % ksolve.method ) +t1 = time.time() +moose.start( 100 ) +print('[INFO] Time taken %s' % (time.time() - t1 )) +expected = [ 7.77859 , 2.77858 , 2.27541 , 6.12141 , 7.77858 , 2.77858 + , 2.27541 , 6.12141 , 7.77858 , 2.77858 , 2.27541 , 6.12141 , 7.77858 + , 2.77858 , 2.27541 , 6.12141 , 7.77858 , 2.77858 , 2.27541 , 6.12141 + , 7.77858 , 2.77858 , 2.27541 , 6.12141 , 7.77858 , 2.77858 , 2.27541 + , 6.12141 , 7.77858 , 2.77858 , 2.27541 , 6.12141 , 7.77858 , 2.77858 + , 2.27541 , 6.12141 , 7.77858 , 2.77858 , 2.27541 , 6.12141 + ] +concs = [ p.conc for p in pools ] +if(not np.isclose( concs, expected ).all() ): + print( " Expected %s" % expected ) + print( " Got %s" % concs ) + quit(1) +print( 'Test passed' ) diff --git a/utility/boost_ode.h b/utility/boost_ode.h new file mode 100644 index 0000000000..28818230f7 --- /dev/null +++ b/utility/boost_ode.h @@ -0,0 +1,39 @@ +/* + * ===================================================================================== + * + * Filename: boost_ode.h + * + * Description: Utility function for Boost ODEINT library. + * + * Version: 1.0 + * Created: Saturday 24 November 2018 04:54:32 IST + * Revision: none + * Compiler: gcc + * + * Author: Dilawar Singh (), dilawars@ncbs.res.in + * Organization: NCBS Bangalore + * + * ===================================================================================== + */ +#ifndef BOOST_ODE_H +#define BOOST_ODE_H + +#include +#include +#include + +typedef double value_type_; +typedef std::vector vector_type_; + +typedef boost::numeric::odeint::runge_kutta4< vector_type_ > rk4_stepper_type_; +typedef boost::numeric::odeint::runge_kutta_dopri5< vector_type_ > rk_dopri_stepper_type_; +typedef boost::numeric::odeint::modified_midpoint< vector_type_ > rk_midpoint_stepper_type_; + +/*----------------------------------------------------------------------------- + * This stepper type found to be most suitable for adaptive solver. The gsl + * implementation has runge_kutta_fehlberg78 solver. + *-----------------------------------------------------------------------------*/ +typedef boost::numeric::odeint::runge_kutta_cash_karp54< vector_type_ > rk_karp_stepper_type_; +typedef boost::numeric::odeint::runge_kutta_fehlberg78< vector_type_ > rk_felhberg_stepper_type_; + +#endif /* end of include guard: BOOST_ODE_H */ diff --git a/utility/print_function.hpp b/utility/print_function.hpp index 25d9fef04c..761009616e 100644 --- a/utility/print_function.hpp +++ b/utility/print_function.hpp @@ -214,8 +214,8 @@ namespace moose { /* Be safe than sorry */ if(!reset) ss << T_RESET; + cout << ss.str() << endl; - cout.flush( ); } /* @@ -307,5 +307,20 @@ namespace moose { logF.close(); } + template + void print_array( T* a, size_t n, const string prefix = "" ) + { + stringstream ss; + ss << prefix; + for (size_t i = 0; i < n; i++) + { + ss << a[i] << ' '; + if( (i+1) % 20 == 0 ) + ss << endl; + } + ss << endl; + cerr << ss.str(); + } + } #endif /* ----- #ifndef print_function_INC ----- */ diff --git a/utility/testing_macros.hpp b/utility/testing_macros.hpp index 21fd8a1254..c3c51a79e9 100644 --- a/utility/testing_macros.hpp +++ b/utility/testing_macros.hpp @@ -132,7 +132,7 @@ static ostringstream assertStream; #define ASSERT_EQ(a, b, token) \ if( ! doubleEq((a), (b)) ) { \ assertStream.str(""); \ - assertStream.precision( 9 ); \ + assertStream.precision( 12 ); \ LOCATION(assertStream) \ assertStream << "Expected " << a << ", received " << b << endl; \ assertStream << token << endl; \ @@ -143,7 +143,8 @@ static ostringstream assertStream; if(! doubleEq(a, b) ) { \ assertStream.str(""); \ LOCATION(assertStream); \ - assertStream << "Expected " << std::fixed << b << ", received " << a << endl; \ + assertStream.precision( 12 ); \ + assertStream << "Expected " << a << ", received " << b << endl; \ assertStream << token; \ moose::__dump__(assertStream.str(), moose::failed); \ throw std::runtime_error( "float equality test failed" ); \