diff --git a/moose-core/.travis/travis_prepare_osx.sh b/moose-core/.travis/travis_prepare_osx.sh index 8c67277ad623b7f9e83e21d63a21f369a9179110..61fdea5b0a4d2f23b16d19ee1d8b1348b69142e5 100755 --- a/moose-core/.travis/travis_prepare_osx.sh +++ b/moose-core/.travis/travis_prepare_osx.sh @@ -23,8 +23,8 @@ rvm get head brew update #brew outdated cmake || brew install cmake brew install gsl -brew outdated hdf5 || brew install homebrew/science/hdf5 -brew outdated libsbml || brew install homebrew/science/libsbml +brew install hdf5 +brew install homebrew/science/libsbml #brew outdated python || brew install python #brew outdated numpy || brew install homebrew/python/numpy #brew unlink numpy && brew link numpy || echo "Failed to link numpy" diff --git a/moose-core/CMakeLists.txt b/moose-core/CMakeLists.txt index 7c4b3226309fc1b138eb866a5f3f8e9c160d156c..068c20b2316d6bc32641955441f1b9d82dad7de0 100644 --- a/moose-core/CMakeLists.txt +++ b/moose-core/CMakeLists.txt @@ -151,11 +151,11 @@ set(STATIC_LIBRARIES "" ) set(SYSTEM_SHARED_LIBS ${LibXML2_LIBRARIES}) if(WITH_GSL) - find_package(GSL 1.16) + find_package(GSL 1.16 REQUIRED) if(NOT GSL_FOUND) message(FATAL_ERROR "=====================================================================\n" - " FATAL gsl(>1.16) not found.\n\n" + " FATAL gsl(>=1.16) not found.\n\n" " MOOSE requires Gnu Scientific Library (GSL) 1.16 or higher. \n" " Please install the `dev` or `devel` package e.g. \n" " $ sudo apt-get install libgsl0-dev \n" @@ -451,8 +451,35 @@ endif() ############################ CTEST ###################################### +include( CTest ) +# If CTEST_OUTPUT_ON_FAILURE environment variable is set, the output is printed +# onto the console if a test fails. +set(ENV{CTEST_OUTPUT_ON_FAILURE} ON) + +# Run this test in debug mode. In Release mode, this does not do anything. +set(MOOSE_BIN_LOCATION $<TARGET_FILE:moose.bin>) +message(STATUS "Executable moose.bin will be at ${MOOSE_BIN_LOCATION}" ) +if(DEBUG) + add_test(NAME moose.bin-raw-run COMMAND moose.bin -u -q) +endif(DEBUG) + +## PyMOOSE tests. +set(PYMOOSE_TEST_DIRECTORY ${CMAKE_SOURCE_DIR}/tests/python) + +# Collect all python script in tests folder and run them. +file(GLOB PY_TEST_SCRIPTS "${PYMOOSE_TEST_DIRECTORY}/test_*.py" ) -include(${CMAKE_CURRENT_SOURCE_DIR}/MooseTests.cmake) +# Run python tests. +foreach( _test_script ${PY_TEST_SCRIPTS} ) + get_filename_component( _test_name ${_test_script} NAME_WE) + add_test( NAME ${_test_name} + COMMAND ${PYTHON_EXECUTABLE} ${_test_script} + WORKING_DIRECTORY ${PYMOOSE_TEST_DIRECTORY} + ) + set_tests_properties( ${_test_name} + PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" + ) +endforeach( ) ############################ CPACK ###################################### include(${CMAKE_CURRENT_SOURCE_DIR}/cmake_moose_cpack.cmake) diff --git a/moose-core/CheckCXXCompiler.cmake b/moose-core/CheckCXXCompiler.cmake index 15a4d730b98679a10bf42e9c5f8df04d6de3fd58..92fbfdfa3c57c433baf213930c59c61300b088a7 100644 --- a/moose-core/CheckCXXCompiler.cmake +++ b/moose-core/CheckCXXCompiler.cmake @@ -34,17 +34,12 @@ if(COMPILER_SUPPORTS_CXX11) if(APPLE) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++" ) endif(APPLE) -elseif(COMPILER_SUPPORTS_CXX0X) - message(STATUS "Your compiler supports c++0x features. Enabling it") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x") - add_definitions( -DENABLE_CXX11 ) +else(COMPILER_SUPPORTS_CXX11) add_definitions( -DBOOST_NO_CXX11_SCOPED_ENUMS -DBOOST_NO_SCOPED_ENUMS ) - if(APPLE) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++" ) - endif(APPLE) -else() - add_definitions( -DBOOST_NO_CXX11_SCOPED_ENUMS -DBOOST_NO_SCOPED_ENUMS ) - message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support.") -endif() + message(FATAL_ERROR "\ + The compiler ${CMAKE_CXX_COMPILER} is too old. \ + Please use a compiler which has c++11 support. + ") +endif(COMPILER_SUPPORTS_CXX11) diff --git a/moose-core/MooseTests.cmake b/moose-core/MooseTests.cmake deleted file mode 100644 index 6026946ded7d8a9ff912baf6c56a7e43f8d46f28..0000000000000000000000000000000000000000 --- a/moose-core/MooseTests.cmake +++ /dev/null @@ -1,154 +0,0 @@ -ENABLE_TESTING() - -FIND_PACKAGE(PythonInterp REQUIRED) - -# If CTEST_OUTPUT_ON_FAILURE environment variable is set, the output is printed -# onto the console if a test fails. -SET(ENV{CTEST_OUTPUT_ON_FAILURE} ON) - -ADD_TEST(NAME moose.bin-raw-run - COMMAND moose.bin -u -q - ) - -## PyMOOSE tests. - -SET(PYMOOSE_TEST_DIRECTORY ${CMAKE_SOURCE_DIR}/tests/python) - -# This test does not work with python-2.6 because of unittest changed API. -#ADD_TEST(NAME pymoose-test-pymoose -# COMMAND ${PYTHON_EXECUTABLE} test_pymoose.py -# WORKING_DIRECTORY ${PYMOOSE_TEST_DIRECTORY} -# ) - -IF(WITH_MPI) - SET(TEST_COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 - ${MPIEXEC_PREFLAGS} ${PYTHON_EXECUTABLE} ${MPIEXEC_POSTFLAGS} - ) -else(WITH_MPI) - SET(TEST_COMMAND ${PYTHON_EXECUTABLE}) -endif(WITH_MPI) - -ADD_TEST(NAME pymoose-test-synchan - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_synchan.py - ) -set_tests_properties(pymoose-test-synchan - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -ADD_TEST(NAME pymoose-test-function - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_function.py - ) -set_tests_properties(pymoose-test-function - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -ADD_TEST(NAME pymoose-test-vec - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_vec.py - ) -set_tests_properties(pymoose-test-vec PROPERTIES ENVIRONMENT - "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -ADD_TEST(NAME pymoose-pyrun - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_pyrun.py - ) -set_tests_properties(pymoose-pyrun - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# Do not run this test after packaging. -ADD_TEST(NAME pymoose-neuroml-reader-test - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_neuroml.py - ) -set_tests_properties(pymoose-neuroml-reader-test - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -ADD_TEST(NAME pymoose-nsdf-sanity-test - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_nsdf.py - ) -set_tests_properties(pymoose-nsdf-sanity-test - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# Test Ksolve -ADD_TEST(NAME pymoose-ksolve-test - COMMAND ${TEST_COMMAND} ${PROJECT_SOURCE_DIR}/tests/python/test_ksolve.py - ) -set_tests_properties(pymoose-ksolve-test - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -## Test basic SBML support. Only for python2. -if( PYTHON_VERSION_STRING VERSION_LESS "3.0.0" ) - ADD_TEST(NAME pymoose-test-basic-sbml-support - COMMAND ${TEST_COMMAND} - ${PROJECT_SOURCE_DIR}/tests/python/test_sbml.py - ) - set_tests_properties(pymoose-test-basic-sbml-support - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) -endif( ) - -ADD_TEST(NAME pymoose-test-rng - COMMAND ${PROJECT_SOURCE_DIR}/tests/python/test_random_gen.sh - ${PYTHON_EXECUTABLE} ${PROJECT_BINARY_DIR}/python - ) -set_tests_properties(pymoose-test-rng PROPERTIES - ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python PYTHON=${PYTHON_EXECUTABLE}" - ) - -# Test Streamer class -ADD_TEST( NAME pymoose-test-streamer - COMMAND ${TEST_COMMAND} - ${PROJECT_SOURCE_DIR}/tests/python/test_streamer.py - ) -set_tests_properties(pymoose-test-streamer - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# Test streaming supports in tables. -ADD_TEST( NAME pymoose-test-streaming_in_tables - COMMAND ${TEST_COMMAND} - ${PROJECT_SOURCE_DIR}/tests/python/test_table_streaming_support.py - ) -set_tests_properties(pymoose-test-streaming_in_tables - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# Test kkit support. -ADD_TEST( NAME pymoose-test-kkit - COMMAND ${TEST_COMMAND} - ${PROJECT_SOURCE_DIR}/tests/python/test_kkit.py - ) -set_tests_properties(pymoose-test-kkit - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# Test Calcium hsolve support. -ADD_TEST( NAME pymoose-test-calcium-hsolve - COMMAND ${TEST_COMMAND} - ${PROJECT_SOURCE_DIR}/tests/python/test_hsolve_externalCalcium.py - ) -set_tests_properties(pymoose-test-calcium-hsolve - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# Test accessing existing path without moose.element. -ADD_TEST( NAME pymoose-test-path - COMMAND ${TEST_COMMAND} - ${PROJECT_SOURCE_DIR}/tests/python/test_paths.py - ) -set_tests_properties(pymoose-test-path - PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" - ) - -# # NOTE: These tests are not enabled yet. They take lot of time. Not all scripts -# # are fixed yet. -# # Test moose-examples with very short timeout. -# ADD_TEST( NAME pymoose-test-moose-examples -# COMMAND ${TEST_COMMAND} -c "import moose; moose.test( timeout = 10 );" -# ) -# set_tests_properties(pymoose-test-moose-examples -# PROPERTIES ENVIRONMENT "PYTHONPATH=${PROJECT_BINARY_DIR}/python" -# ) diff --git a/moose-core/biophysics/Neuron.cpp b/moose-core/biophysics/Neuron.cpp index 4ca30e107cdb0b6b090453c41e29aa9c43bff606..4ae85145ffa04acc1521e7bfa8c9319fa809f44a 100644 --- a/moose-core/biophysics/Neuron.cpp +++ b/moose-core/biophysics/Neuron.cpp @@ -420,6 +420,16 @@ const Cinfo* Neuron::initCinfo() &Neuron::getParentCompartmentOfSpine ); + static ReadOnlyLookupElementValueFinfo< Neuron, vector< ObjId >, vector< ObjId > > + spineIdsFromCompartmentIds( + "spineIdsFromCompartmentIds", + "Vector of ObjIds of spine entries (FieldElements on this Neuron, " + "used for scaling) that map to the the specified " + "electrical compartments. If a bad compartment Id is given, the" + "corresponding spine entry is the root Id.", + &Neuron::getSpineIdsFromCompartmentIds + ); + ///////////////////////////////////////////////////////////////////// // DestFinfos ///////////////////////////////////////////////////////////////////// @@ -498,6 +508,7 @@ const Cinfo* Neuron::initCinfo() &spinesFromExpression, // ReadOnlyLookupValueFinfo &spinesOnCompartment, // ReadOnlyLookupValueFinfo &parentCompartmentOfSpine, // ReadOnlyLookupValueFinfo + &spineIdsFromCompartmentIds, // ReadOnlyLookupValueFinfo &buildSegmentTree, // DestFinfo &setSpineAndPsdMesh, // DestFinfo &setSpineAndPsdDsolve, // DestFinfo @@ -1058,6 +1069,37 @@ ObjId Neuron::getParentCompartmentOfSpine( return ObjId(); } +vector< ObjId > Neuron::getSpineIdsFromCompartmentIds( + const Eref& e, vector< ObjId > compt ) const +{ + vector< ObjId > ret; + map< Id, unsigned int > lookupSpine; + Id spineBase = Id( e.id().value() + 1 ); + for ( unsigned int i = 0; i < spines_.size(); ++i ) { + for ( vector< Id >::const_iterator j = spines_[i].begin(); j != spines_[i].end(); ++j ) { + lookupSpine[ *j ] = i; + } + } + // cout << "################## " << lookupSpine.size() << endl; + for ( map< Id, unsigned int >::const_iterator k = lookupSpine.begin(); k != lookupSpine.end(); ++k ) { + // cout << "spine[" << k->second << "] has " << k->first.element()->getName() << endl; + // cout << "spine[" << k->second << "] has " << k->first << endl; + + } + for ( vector< ObjId >::const_iterator j = compt.begin(); j != compt.end(); ++j ) + { + // cout << "compt: " << *j << " " << j->element()->getName() << endl; + map< Id, unsigned int >::const_iterator k = lookupSpine.find( j->id ); + if ( k != lookupSpine.end() ) { + ret.push_back( ObjId( spineBase, e.dataIndex(), k->second ) ); + // cout << "spine[" << k->second << "] has " << j->element()->getName() << endl; + } else { + ret.push_back( ObjId() ); + } + } + return ret; +} + void Neuron::buildElist( const Eref& e, const vector< string >& line, vector< ObjId >& elist, @@ -1974,6 +2016,9 @@ const vector< Id >& Neuron::spineIds( unsigned int index ) const void Neuron::scaleBufAndRates( unsigned int spineNum, double lenScale, double diaScale ) const { + double volScale = lenScale * diaScale * diaScale; + if ( doubleEq( volScale, 1.0 ) ) + return; if ( spineStoich_.size() == 0 ) // Perhaps no chem stuff in model, but user could have forgotten // to assign psd and spine meshes. @@ -1999,7 +2044,6 @@ void Neuron::scaleBufAndRates( unsigned int spineNum, // The chem system for the spine may not have been defined. return; } - double volScale = lenScale * diaScale * diaScale; SetGet2< unsigned int, double >::set( ss, "scaleBufsAndRates", spineToMeshOrdering_[spineNum], volScale ); volScale = diaScale * diaScale; @@ -2024,7 +2068,9 @@ void Neuron::scaleHeadDiffusion( unsigned int spineNum, double len, double dia) const { double vol = len * dia * dia * PI * 0.25; - double diffScale = dia * dia * 0.25 * PI / len; + // Note that the diffusion scale for the PSD uses half the length + // of the head compartment. I'm explicitly putting this in below. + double diffScale = dia * dia * 0.25 * PI / (len/2.0); unsigned int meshIndex = spineToMeshOrdering_[ spineNum ]; Id headCompt = Field< Id >::get( headDsolve_, "compartment" ); LookupField< unsigned int, double >::set( headCompt, "oneVoxelVolume", diff --git a/moose-core/biophysics/Neuron.h b/moose-core/biophysics/Neuron.h index 3cef4dab7edf5a39bfee37e17b41e678c0d1fa82..777c5448600b2f56209dc016c951a359881e4da6 100644 --- a/moose-core/biophysics/Neuron.h +++ b/moose-core/biophysics/Neuron.h @@ -50,6 +50,8 @@ class Neuron const Eref& e, ObjId compt ) const; ObjId getParentCompartmentOfSpine( const Eref& e, ObjId compt ) const; + vector< ObjId > getSpineIdsFromCompartmentIds( + const Eref& e, vector< ObjId > compt ) const; void setChannelDistribution( const Eref& e, vector< string > v ); vector< string > getChannelDistribution( const Eref& e ) const; void setPassiveDistribution( const Eref& e, vector< string > v ); diff --git a/moose-core/biophysics/Spine.cpp b/moose-core/biophysics/Spine.cpp index 6909b54555d332e3bf3b72aa47a401871428af71..03ca6eacab84e0f3a1b7d48137c301ead12b80c3 100644 --- a/moose-core/biophysics/Spine.cpp +++ b/moose-core/biophysics/Spine.cpp @@ -59,7 +59,7 @@ const Cinfo* Spine::initCinfo() "scaled by the user, both the diameter and the length of the " "spine head scale by the cube root of the ratio to the " "previous volume. The diameter of the PSD is pegged to the " - "diameter fo the spine head. \n" + "diameter of the spine head. \n" "This is useful to scale total # of molecules in the head. ", &Spine::setHeadVolume, &Spine::getHeadVolume diff --git a/moose-core/builtins/Table.cpp b/moose-core/builtins/Table.cpp index e07e4172ad6e40f8f7c41c30b6a7ef4a1eaed551..c82c6a8d145ff1c16bcd9d928334263bd1dca54d 100644 --- a/moose-core/builtins/Table.cpp +++ b/moose-core/builtins/Table.cpp @@ -59,6 +59,15 @@ const Cinfo* Table::initCinfo() , &Table::getUseStreamer ); + static ValueFinfo< Table, bool > useSpikeMode( + "useSpikeMode" + , "When set to true, look for spikes in a time-series." + " Normally used for monitoring Vm for action potentials." + " Could be used for any thresholded event. Default is False." + , &Table::setUseSpikeMode + , &Table::getUseSpikeMode + ); + static ValueFinfo< Table, string > outfile( "outfile" , "Set the name of file to which data is written to. If set, " @@ -131,6 +140,7 @@ const Cinfo* Table::initCinfo() &columnName, // Value &outfile, // Value &useStreamer, // Value + &useSpikeMode, // Value handleInput(), // DestFinfo &spike, // DestFinfo requestOut(), // SrcFinfo @@ -200,7 +210,13 @@ const Cinfo* Table::initCinfo() static const Cinfo* tableCinfo = Table::initCinfo(); -Table::Table() : threshold_( 0.0 ) , lastTime_( 0.0 ) , input_( 0.0 ), dt_( 0.0 ) +Table::Table() : + threshold_( 0.0 ) , + lastTime_( 0.0 ) , + input_( 0.0 ), + fired_(false), + useSpikeMode_(false), + dt_( 0.0 ) { // Initialize the directory to which each table should stream. rootdir_ = "_tables"; @@ -238,7 +254,13 @@ void Table::process( const Eref& e, ProcPtr p ) // Copy incoming data to ret and insert into vector. vector< double > ret; requestOut()->send( e, &ret ); - vec().insert( vec().end(), ret.begin(), ret.end() ); + if (useSpikeMode_) { + for ( vector< double >::const_iterator + i = ret.begin(); i != ret.end(); ++i ) + spike( *i ); + } else { + vec().insert( vec().end(), ret.begin(), ret.end() ); + } /* If we are streaming to a file, let's write to a file. And clean the * vector. @@ -268,6 +290,7 @@ void Table::reinit( const Eref& e, ProcPtr p ) unsigned int numTick = e.element()->getTick(); Clock* clk = reinterpret_cast<Clock*>(Id(1).eref().data()); dt_ = clk->getTickDt( numTick ); + fired_ = false; /** Create the default filepath for this table. */ if( useStreamer_ ) @@ -290,7 +313,13 @@ void Table::reinit( const Eref& e, ProcPtr p ) lastTime_ = 0; vector< double > ret; requestOut()->send( e, &ret ); - vec().insert( vec().end(), ret.begin(), ret.end() ); + if (useSpikeMode_) { + for ( vector< double >::const_iterator + i = ret.begin(); i != ret.end(); ++i ) + spike( *i ); + } else { + vec().insert( vec().end(), ret.begin(), ret.end() ); + } if( useStreamer_ ) { @@ -313,8 +342,15 @@ void Table::input( double v ) void Table::spike( double v ) { - if ( v > threshold_ ) - vec().push_back( lastTime_ ); + if ( fired_ ) { // Wait for it to go below threshold + if ( v < threshold_ ) + fired_ = false; + } else { + if ( v > threshold_ ) { // wait for it to go above threshold. + fired_ = true; + vec().push_back( lastTime_ ); + } + } } ////////////////////////////////////////////////////////////// @@ -371,6 +407,18 @@ bool Table::getUseStreamer( void ) const return useStreamer_; } +/* Enable/disable spike mode. */ +void Table::setUseSpikeMode( bool useSpikeMode ) +{ + useSpikeMode_ = useSpikeMode; +} + +bool Table::getUseSpikeMode( void ) const +{ + return useSpikeMode_; +} + + /* set/get outfile_ */ void Table::setOutfile( string outpath ) { diff --git a/moose-core/builtins/Table.h b/moose-core/builtins/Table.h index 5a21b2d31c13d9d26ef53a2d9d9c0e6366691aab..9fb78c8b69f5416c999d3383708c58f127e50be8 100644 --- a/moose-core/builtins/Table.h +++ b/moose-core/builtins/Table.h @@ -43,6 +43,9 @@ public: void setUseStreamer ( bool status ); bool getUseStreamer ( void ) const; + void setUseSpikeMode ( bool status ); + bool getUseSpikeMode ( void ) const; + void setOutfile ( string outfilepath ); string getOutfile ( void ) const; @@ -75,6 +78,8 @@ private: double threshold_; double lastTime_; double input_; + bool fired_; + bool useSpikeMode_; /** * @brief Keep the data, each entry is preceeded by time value. diff --git a/moose-core/cmake_modules/FindGSL.cmake b/moose-core/cmake_modules/FindGSL.cmake index 93c2686f78e6fbbdde540e86241d0995934c8813..b45aca3a8c776e4f64e1a28cba727c3b3e0df80b 100644 --- a/moose-core/cmake_modules/FindGSL.cmake +++ b/moose-core/cmake_modules/FindGSL.cmake @@ -27,7 +27,7 @@ # Set this envrionment variable to search in this path first. SET(GSL_ROOT_DIR $ENV{GSL_ROOT_DIR}) if(GSL_ROOT_DIR) - message( STATUS "Debug: GSL_ROOT_DIR value ${GSL_ROOT_DIR}") + message( STATUS "Debug: GSL_ROOT_DIR = ${GSL_ROOT_DIR}") endif(GSL_ROOT_DIR) IF(WIN32) @@ -54,7 +54,7 @@ IF(WIN32) ENDIF (GSL_LIB AND GSLCBLAS_LIB) ELSE(WIN32) - + # UNIX IF(GSL_USE_STATIC_LIBRARIES) SET(GSL_LIB_NAMES libgsl.a) SET(GSL_CBLAS_LIB_NAMES libgslcblas.a) @@ -63,31 +63,40 @@ ELSE(WIN32) SET(GSL_CBLAS_LIB_NAMES gslcblas) ENDIF( ) - FIND_LIBRARY(GSL_LIB - NAMES ${GSL_LIB_NAMES} - PATHS - ${GSL_ROOT_DIR}/lib ${GSL_ROOT_DIR}/lib64 - /opt/lib /opt/lib64 - ) + if(GSL_ROOT_DIR) + FIND_LIBRARY(GSL_LIB + NAMES ${GSL_LIB_NAMES} + PATHS ${GSL_ROOT_DIR}/lib NO_DEFAULT_PATH + ) - FIND_LIBRARY(GSLCBLAS_LIB - NAMES ${GSL_CBLAS_LIB_NAMES} - PATHS - ${GSL_ROOT_DIR}/lib ${GSL_ROOT_DIR}/lib64 - /opt/lib /opt/lib64 - ) - IF (GSL_LIB AND GSLCBLAS_LIB) - SET (GSL_LIBRARIES ${GSL_LIB} ${GSLCBLAS_LIB}) - ENDIF (GSL_LIB AND GSLCBLAS_LIB) + FIND_LIBRARY(GSLCBLAS_LIB + NAMES ${GSL_CBLAS_LIB_NAMES} + PATHS ${GSL_ROOT_DIR}/lib NO_DEFAULT_PATH + ) + IF (GSL_LIB AND GSLCBLAS_LIB) + SET (GSL_LIBRARIES ${GSL_LIB} ${GSLCBLAS_LIB}) + ENDIF (GSL_LIB AND GSLCBLAS_LIB) - FIND_PATH(GSL_INCLUDE_DIRS NAMES gsl/gsl_blas.h - PATHS ${GSL_ROOT_DIR}/include /opt/include - ) + FIND_PATH(GSL_INCLUDE_DIRS NAMES gsl/gsl_blas.h + PATHS ${GSL_ROOT_DIR}/include NO_DEFAULT_PATH + ) + else(GSL_ROOT_DIR) + FIND_LIBRARY(GSL_LIB NAMES ${GSL_LIB_NAMES} ) + FIND_LIBRARY(GSLCBLAS_LIB NAMES ${GSL_CBLAS_LIB_NAMES}) + + IF (GSL_LIB AND GSLCBLAS_LIB) + SET (GSL_LIBRARIES ${GSL_LIB} ${GSLCBLAS_LIB}) + ENDIF (GSL_LIB AND GSLCBLAS_LIB) + FIND_PATH(GSL_INCLUDE_DIRS NAMES gsl/gsl_blas.h + PATHS ${GSL_ROOT_DIR}/include /opt/include + ) + endif( ) ENDIF(WIN32) # FIND version +# message(STATUS "Searching in ${GSL_INCLUDE_DIRS}") if(GSL_INCLUDE_DIRS) file(READ "${GSL_INCLUDE_DIRS}/gsl/gsl_version.h" GSL_VERSION_TEXT) string(REGEX REPLACE ".*define[ ]+GSL_MAJOR_VERSION[ ]*([0-9]+).*" "\\1" diff --git a/moose-core/kinetics/BufPool.cpp b/moose-core/kinetics/BufPool.cpp index 230be5d6cebbb52abe1df35e25b93fcac4faee25..54e183780a6f76fcb1614878e70e02981b28ae25 100644 --- a/moose-core/kinetics/BufPool.cpp +++ b/moose-core/kinetics/BufPool.cpp @@ -78,7 +78,6 @@ void BufPool::vSetConcInit( const Eref& e, double conc ) vSetConc( e, conc ); } - ////////////////////////////////////////////////////////////// // MsgDest Definitions ////////////////////////////////////////////////////////////// @@ -97,3 +96,4 @@ void BufPool::vReinit( const Eref& e, ProcPtr p ) ////////////////////////////////////////////////////////////// // Field Definitions ////////////////////////////////////////////////////////////// + diff --git a/moose-core/kinetics/Pool.cpp b/moose-core/kinetics/Pool.cpp index 86508cf4429f043e47483bbd419c88f2917d496b..8b0eb62940d7b361d6e3b9db3236a70e310aedb3 100644 --- a/moose-core/kinetics/Pool.cpp +++ b/moose-core/kinetics/Pool.cpp @@ -20,54 +20,21 @@ const Cinfo* Pool::initCinfo() ////////////////////////////////////////////////////////////// // Field Definitions: All inherited from PoolBase. ////////////////////////////////////////////////////////////// - static ElementValueFinfo< Pool, bool > isBuffered( - "isBuffered", - "Flag: True if Pool is buffered. This field changes the " - "type of the Pool object to BufPool, or vice versa. " - "None of the messages are affected. " - "This object class flip can only be done in the non-zombified " - "form of the Pool.", - &Pool::setIsBuffered, - &Pool::getIsBuffered - ); - ////////////////////////////////////////////////////////////// - // MsgDest Definitions: All but increment and decrement inherited + // MsgDest Definitions: All inherited ////////////////////////////////////////////////////////////// - static DestFinfo increment( "increment", - "Increments mol numbers by specified amount. Can be +ve or -ve", - new OpFunc1< Pool, double >( &Pool::increment ) - ); - - static DestFinfo decrement( "decrement", - "Decrements mol numbers by specified amount. Can be +ve or -ve", - new OpFunc1< Pool, double >( &Pool::decrement ) - ); - - static DestFinfo nIn( "nIn", - "Set the number of molecules by specified amount", - new OpFunc1< Pool, double >( &Pool::nIn ) - ); - ////////////////////////////////////////////////////////////// // SrcFinfo Definitions: All inherited. ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// // SharedMsg Definitions: All inherited. ////////////////////////////////////////////////////////////// - static Finfo* poolFinfos[] = { - &isBuffered, // ElementValueFinfo - &increment, // DestFinfo - &decrement, // DestFinfo - &nIn, // DestFinfo - }; - static Dinfo< Pool > dinfo; static Cinfo poolCinfo ( "Pool", PoolBase::initCinfo(), - poolFinfos, - sizeof( poolFinfos ) / sizeof ( Finfo* ), + 0, + 0, &dinfo ); @@ -98,10 +65,10 @@ Pool::~Pool() * It uses a low-level replaceCinfo call to just change the * identity of the Cinfo used, leaving everything else as is. */ -void Pool::setIsBuffered( const Eref& e, bool v ) +void Pool::vSetIsBuffered( const Eref& e, bool v ) { static const Cinfo* bufPoolCinfo = Cinfo::find( "BufPool" ); - if (getIsBuffered( e ) == v) + if (vGetIsBuffered( e ) == v) return; if (v) { e.element()->replaceCinfo( bufPoolCinfo ); @@ -110,8 +77,11 @@ void Pool::setIsBuffered( const Eref& e, bool v ) } } -bool Pool::getIsBuffered( const Eref& e ) const +bool Pool::vGetIsBuffered( const Eref& e ) const { + /// We need this explicit check because when the moose class is + /// flipped, the internal C++ class isn't. + /// Inherited by BufPool. return e.element()->cinfo()->name() == "BufPool"; } @@ -125,14 +95,6 @@ void Pool::vProcess( const Eref& e, ProcPtr p ) { // double A = e.sumBuf( aSlot ); // double B = e.sumBuf( bSlot ); - /* - if ( n_ < 0 ) - cout << "nugh" << e.objId().path() << endl; - if ( B_ < 0 ) - cout << "bugh" << e.objId().path() << endl; - if ( p->dt < 0 ) - cout << "tugh" << e.objId().path() << endl; - */ if ( n_ > EPSILON && B_ > EPSILON ) { double C = exp( -B_ * p->dt / n_ ); @@ -162,7 +124,7 @@ void Pool::vReac( double A, double B ) B_ += B; } -void Pool::increment( double val ) +void Pool::vIncrement( double val ) { if ( val > 0 ) A_ += val; @@ -170,7 +132,7 @@ void Pool::increment( double val ) B_ -= val; } -void Pool::decrement( double val ) +void Pool::vDecrement( double val ) { if ( val < 0 ) A_ -= val; @@ -178,37 +140,13 @@ void Pool::decrement( double val ) B_ += val; } -void Pool::nIn( double val) +void Pool::vnIn( double val) { n_ = val; B_ = 0; A_ = 0; } -/* -void Pool::vRemesh( const Eref& e, - double oldvol, - unsigned int numTotalEntries, unsigned int startEntry, - const vector< unsigned int >& localIndices, - const vector< double >& vols ) -{ - if ( e.dataIndex() != 0 ) - return; - Neutral* n = reinterpret_cast< Neutral* >( e.data() ); - assert( vols.size() > 0 ); - double concInit = nInit_ / ( NA * oldvol ); - if ( vols.size() != e.element()->dataHandler()->localEntries() ) - n->setLastDimension( e, q, vols.size() ); - // Note that at this point the Pool pointer may be invalid! - // But we need to update the concs anyway. - assert( e.element()->dataHandler()->localEntries() == vols.size() ); - Pool* pooldata = reinterpret_cast< Pool* >( e.data() ); - for ( unsigned int i = 0; i < vols.size(); ++i ) { - pooldata[i].nInit_ = pooldata[i].n_ = concInit * vols[i] * NA; - } -} -*/ - void Pool::vHandleMolWt( const Eref& e, double v ) { ; // Here I should update DiffConst too. diff --git a/moose-core/kinetics/Pool.h b/moose-core/kinetics/Pool.h index d55696aae486670b864fe1226693317ddf05470a..bc4e88770a781fc67b59b2a6a6f442f4e32f46dd 100644 --- a/moose-core/kinetics/Pool.h +++ b/moose-core/kinetics/Pool.h @@ -57,8 +57,8 @@ class Pool: public PoolBase /** * Functions to examine and change class between Pool and BufPool. */ - void setIsBuffered( const Eref& e, bool v ); - bool getIsBuffered( const Eref& e ) const; + void vSetIsBuffered( const Eref& e, bool v ); + bool vGetIsBuffered( const Eref& e) const; ////////////////////////////////////////////////////////////////// // Dest funcs. These too override virtual funcs in the Pool base @@ -69,13 +69,13 @@ class Pool: public PoolBase void vProcess( const Eref& e, ProcPtr p ); void vReinit( const Eref& e, ProcPtr p ); void vReac( double A, double B ); + void vIncrement( double val ); + void vDecrement( double val ); + void vnIn( double val ); ////////////////////////////////////////////////////////////////// // Novel Dest funcs not present in Pool base class. ////////////////////////////////////////////////////////////////// - void increment( double val ); - void decrement( double val ); - void nIn( double val ); ////////////////////////////////////////////////////////////////// static const Cinfo* initCinfo(); diff --git a/moose-core/kinetics/PoolBase.cpp b/moose-core/kinetics/PoolBase.cpp index b67e70cf7001372b166554d67311daadbfe7475e..f1d7e0f44ba4984981c3a1b2f48e8fd6d0678843 100644 --- a/moose-core/kinetics/PoolBase.cpp +++ b/moose-core/kinetics/PoolBase.cpp @@ -84,6 +84,18 @@ const Cinfo* PoolBase::initCinfo() &PoolBase::getSpecies ); + static ElementValueFinfo< PoolBase, bool > isBuffered( + "isBuffered", + "Flag: True if Pool is buffered. " + "In the case of Pool and BufPool the field can be assigned, to " + "change the type of the Pool object to BufPool, or vice versa. " + "None of the messages are affected. " + "This object class flip can only be done in the non-zombified " + "form of the Pool/BufPool. In Zombies it is read-only.", + &PoolBase::setIsBuffered, + &PoolBase::getIsBuffered + ); + ////////////////////////////////////////////////////////////// // MsgDest Definitions ////////////////////////////////////////////////////////////// @@ -104,6 +116,24 @@ const Cinfo* PoolBase::initCinfo() "Should only be used in SharedMsg with species.", new EpFunc1< PoolBase, double >( &PoolBase::handleMolWt ) ); + ////////////////////////////////////////////////////////////// + // MsgDest Definitions: These three are used for non-reaction + // calculations involving algebraically defined rate terms. + ////////////////////////////////////////////////////////////// + static DestFinfo increment( "increment", + "Increments mol numbers by specified amount. Can be +ve or -ve", + new OpFunc1< PoolBase, double >( &PoolBase::increment ) + ); + + static DestFinfo decrement( "decrement", + "Decrements mol numbers by specified amount. Can be +ve or -ve", + new OpFunc1< PoolBase, double >( &PoolBase::decrement ) + ); + + static DestFinfo nIn( "nIn", + "Assigns the number of molecules in Pool to specified value", + new OpFunc1< PoolBase, double >( &PoolBase::nIn ) + ); ////////////////////////////////////////////////////////////// // SrcFinfo Definitions @@ -155,6 +185,10 @@ const Cinfo* PoolBase::initCinfo() &concInit, // Value &volume, // Readonly Value &speciesId, // Value + &isBuffered, // Value + &increment, // DestFinfo + &decrement, // DestFinfo + &nIn, // DestFinfo &reac, // SharedFinfo &proc, // SharedFinfo &species, // SharedFinfo @@ -208,6 +242,21 @@ void PoolBase::reinit( const Eref& e, ProcPtr p ) vReinit( e, p ); } +void PoolBase::increment( double val ) +{ + vIncrement(val); +} + +void PoolBase::decrement( double val ) +{ + vDecrement( val ); +} + +void PoolBase::nIn( double val) +{ + vnIn(val); +} + void PoolBase::reac( double A, double B ) { vReac( A, B ); @@ -234,6 +283,15 @@ void PoolBase::vReac( double A, double B ) void PoolBase::vHandleMolWt( const Eref& e, double v ) {;} +void PoolBase::vIncrement( double val ) +{;} + +void PoolBase::vDecrement( double val ) +{;} + +void PoolBase::vnIn( double val) +{;} + ////////////////////////////////////////////////////////////// // Field Definitions ////////////////////////////////////////////////////////////// @@ -328,6 +386,19 @@ unsigned int PoolBase::getSpecies( const Eref& e ) const return vGetSpecies( e ); } +/** + * setIsBuffered is active only for Pool and BufPool. Otherwise ignored. + */ +void PoolBase::setIsBuffered( const Eref& e, bool v ) +{ + vSetIsBuffered( e, v ); +} + +bool PoolBase::getIsBuffered( const Eref& e ) const +{ + return vGetIsBuffered( e ); +} + ////////////////////////////////////////////////////////////// // Virtual Field Definitions ////////////////////////////////////////////////////////////// @@ -341,6 +412,10 @@ double PoolBase::vGetMotorConst(const Eref& e ) const return 0.0; } +/// Dummy default function for most pool subclasses. +void PoolBase::vSetIsBuffered( const Eref& e, bool v ) +{;} + ////////////////////////////////////////////////////////////// // Zombie conversion routine: Converts Pool subclasses. There // will typically be a target specific follow-up function, for example, diff --git a/moose-core/kinetics/PoolBase.h b/moose-core/kinetics/PoolBase.h index 3821569ab7cfa84050472e23034199115a69da08..b5aaa1b8d545e363f1c413786ba7c73b5ada9400 100644 --- a/moose-core/kinetics/PoolBase.h +++ b/moose-core/kinetics/PoolBase.h @@ -68,6 +68,11 @@ class PoolBase void setSpecies( const Eref& e, SpeciesId v ); SpeciesId getSpecies( const Eref& e ) const; + /** + * Functions to examine and change class between Pool and BufPool. + */ + void setIsBuffered( const Eref& e, bool v ); + bool getIsBuffered( const Eref& e ) const; ////////////////////////////////////////////////////////////////// // Here are the inner virtual funcs for fields. @@ -91,6 +96,9 @@ class PoolBase virtual void vSetVolume( const Eref& e, double v ) = 0; virtual void vSetSpecies( const Eref& e, SpeciesId v ) = 0; virtual SpeciesId vGetSpecies( const Eref& e ) const = 0; + /// I put in a default empty function for vSetIsBuffered. + virtual void vSetIsBuffered( const Eref& e, bool v ); + virtual bool vGetIsBuffered( const Eref& e) const = 0; /** * Assign whatever info is needed by the zombie based on the * solver Element. Encapsulates some unpleasant field extraction, @@ -126,6 +134,9 @@ class PoolBase void reinit( const Eref& e, ProcPtr p ); void reac( double A, double B ); void handleMolWt( const Eref& e, double v ); + void increment( double val ); + void decrement( double val ); + void nIn( double val ); ////////////////////////////////////////////////////////////////// // Virtual Dest funcs. Most of these have a generic do-nothing @@ -136,6 +147,9 @@ class PoolBase virtual void vReinit( const Eref& e, ProcPtr p ); virtual void vReac( double A, double B ); virtual void vHandleMolWt( const Eref& e, double v); + virtual void vIncrement( double val ); + virtual void vDecrement( double val ); + virtual void vnIn( double val ); ////////////////////////////////////////////////////////////////// static const Cinfo* initCinfo(); diff --git a/moose-core/ksolve/GssaVoxelPools.cpp b/moose-core/ksolve/GssaVoxelPools.cpp index 0b09a418712f143bba241b6c7dc780abed4779f6..438652d8af1f748a9226d4befabbd04a0de234ce 100644 --- a/moose-core/ksolve/GssaVoxelPools.cpp +++ b/moose-core/ksolve/GssaVoxelPools.cpp @@ -70,7 +70,7 @@ void GssaVoxelPools::updateDependentMathExpn( // The lower commented block has all funcs updated every time step, // but this too // doesn't update the cascading reacs. So this is now handled by the - // useClockedUpdate flag, and we use the upper block here instead. + // useClockedUpdate flag, and we use the lower block here instead. /* const vector< unsigned int >& deps = g->dependentMathExpn[ rindex ]; for( vector< unsigned int >::const_iterator @@ -78,11 +78,15 @@ void GssaVoxelPools::updateDependentMathExpn( g->stoich->funcs( *i )->evalPool( varS(), time ); } */ + /* unsigned int numFuncs = g->stoich->getNumFuncs(); for( unsigned int i = 0; i < numFuncs; ++i ) { g->stoich->funcs( i )->evalPool( varS(), time ); } + */ + // This function is equivalent to the loop above. + g->stoich->updateFuncs( varS(), time ); } void GssaVoxelPools::updateDependentRates( @@ -135,6 +139,7 @@ void GssaVoxelPools::setNumReac( unsigned int n ) */ bool GssaVoxelPools::refreshAtot( const GssaSystem* g ) { + g->stoich->updateFuncs( varS(), t_ ); updateReacVelocities( g, S(), v_ ); atot_ = 0; for ( vector< double >::const_iterator @@ -156,7 +161,6 @@ bool GssaVoxelPools::refreshAtot( const GssaSystem* g ) */ void GssaVoxelPools::recalcTime( const GssaSystem* g, double currTime ) { - updateDependentMathExpn( g, 0, currTime ); refreshAtot( g ); assert( t_ > currTime ); t_ = currTime; @@ -174,6 +178,8 @@ void GssaVoxelPools::advance( const ProcInfo* p, const GssaSystem* g ) if ( atot_ <= 0.0 ) // reac system is stuck, will not advance. { t_ = nextt; + g->stoich->updateFuncs( varS(), t_ ); + // updateDependentMathExpn( g, 0, t_ ); return; } unsigned int rindex = pickReac(); @@ -185,6 +191,8 @@ void GssaVoxelPools::advance( const ProcInfo* p, const GssaSystem* g ) if ( !refreshAtot( g ) ) // Stuck state. { t_ = nextt; + g->stoich->updateFuncs( varS(), t_ ); + // updateDependentMathExpn( g, 0, t_ ); return; } // We had a roundoff error, fixed it, but now need to be sure @@ -210,8 +218,8 @@ void GssaVoxelPools::advance( const ProcInfo* p, const GssaSystem* g ) r = rng_.uniform(); } t_ -= ( 1.0 / atot_ ) * log( r ); - // g->stoich->updateFuncs( varS(), t_ ); // Handled next line. - updateDependentMathExpn( g, rindex, t_ ); + g->stoich->updateFuncs( varS(), t_ ); + // updateDependentMathExpn( g, rindex, t_ ); updateDependentRates( g->dependency[ rindex ], g->stoich ); } } @@ -371,7 +379,7 @@ void GssaVoxelPools::xferIn( XferInfo& xf, // cout << x << " i = " << *i << *j << " m = " << *m << endl; double dx = *i++ - *j++; double base = floor( dx ); - if ( rng_.uniform() > dx - base ) + if ( rng_.uniform() >= (dx - base) ) x += base; else x += base + 1.0; @@ -425,7 +433,7 @@ void GssaVoxelPools::xferInOnlyProxies( if ( *k >= stoichPtr_->getNumVarPools() && *k < proxyEndIndex ) { double base = floor( *i ); - if ( rng_.uniform() > *i - base ) + if ( rng_.uniform() >= (*i - base) ) varSinit()[*k] = (varS()[*k] += base ); else varSinit()[*k] = (varS()[*k] += base + 1.0 ); diff --git a/moose-core/ksolve/KinSparseMatrix.cpp b/moose-core/ksolve/KinSparseMatrix.cpp index 59783b29269fa71c24fdabdd397af606da799b79..1b1c1d31e23e1ec85726fe71f1701db83930ca63 100644 --- a/moose-core/ksolve/KinSparseMatrix.cpp +++ b/moose-core/ksolve/KinSparseMatrix.cpp @@ -164,21 +164,6 @@ void makeVecUnique( vector< unsigned int >& v ) v.resize( pos - v.begin() ); } -vector< int > KinSparseMatrix::matrixEntry() const -{ - return N_; -} - -vector< unsigned int > KinSparseMatrix::colIndex() const -{ - return colIndex_; -} - -vector< unsigned int > KinSparseMatrix::rowStart() const -{ - return rowStart_; -} - #ifdef DO_UNIT_TESTS #include "header.h" diff --git a/moose-core/ksolve/KinSparseMatrix.h b/moose-core/ksolve/KinSparseMatrix.h index 4cc48aecd31e37a6b06b0a7e690a4557336fce99..52132dcc0ac9597691c20f4aa34865dbe89d52d2 100644 --- a/moose-core/ksolve/KinSparseMatrix.h +++ b/moose-core/ksolve/KinSparseMatrix.h @@ -61,10 +61,6 @@ class KinSparseMatrix: public SparseMatrix< int > */ void truncateRow( unsigned int maxColumnIndex ); - /// Here we expose the sparse matrix for MOOSE use. - vector< int > matrixEntry() const; - vector< unsigned int > colIndex() const; - vector< unsigned int > rowStart() const; private: /** * End colIndex for rows (molecules in the transposed matrix) diff --git a/moose-core/ksolve/Stoich.cpp b/moose-core/ksolve/Stoich.cpp index 425c19ddbd0b0a07b71a39b3b5843c6a316927b7..198f1978f1acce36a84d62b2f8c24de970341ea0 100644 --- a/moose-core/ksolve/Stoich.cpp +++ b/moose-core/ksolve/Stoich.cpp @@ -631,6 +631,12 @@ const FuncTerm* Stoich::funcs( unsigned int i ) const return funcs_[i]; } +bool Stoich::isFuncTarget( unsigned int poolIndex ) const +{ + assert( poolIndex < funcTarget_.size() ); + return ( funcTarget_[poolIndex] != ~0 ); +} + vector< int > Stoich::getMatrixEntry() const { return N_.matrixEntry(); @@ -968,6 +974,10 @@ void Stoich::resizeArrays() species_.resize( totNumPools, 0 ); + funcTarget_.clear(); + // Only the pools controlled by a func (targets) have positive indices. + funcTarget_.resize( totNumPools, ~0 ); + unsigned int totNumRates = ( reacVec_.size() + offSolverReacVec_.size() ) * (1+useOneWay_) + ( enzVec_.size() + offSolverEnzVec_.size() ) * (2 + useOneWay_ ) + @@ -1118,10 +1128,13 @@ void Stoich::installAndUnschedFunc( Id func, Id pool, double volScale ) string expr = Field< string >::get( func, "expr" ); ft->setExpr( expr ); // Tie the output of the FuncTerm to the pool it controls. - ft->setTarget( convertIdToPoolIndex( pool ) ); + unsigned int targetIndex = convertIdToPoolIndex( pool ); + ft->setTarget( targetIndex ); ft->setVolScale( volScale ); unsigned int funcIndex = convertIdToFuncIndex( func ); assert( funcIndex != ~0U ); + // funcTarget_ vector tracks which pools are controlled by which func. + funcTarget_[targetIndex] = funcIndex; funcs_[ funcIndex ] = ft; } diff --git a/moose-core/ksolve/Stoich.h b/moose-core/ksolve/Stoich.h index 2951e3655b2f50accd2c25bed97b040867fdef27..6dbe51de34bea6d86df19a2c32005ed2fd58b978 100644 --- a/moose-core/ksolve/Stoich.h +++ b/moose-core/ksolve/Stoich.h @@ -126,6 +126,8 @@ class Stoich unsigned int getNumFuncs() const; const FuncTerm* funcs( unsigned int i ) const; + /// Returns true if the specified pool is controlled by a func + bool isFuncTarget( unsigned int poolIndex ) const; vector< int > getMatrixEntry() const; vector< unsigned int > getColIndex() const; @@ -620,6 +622,13 @@ class Stoich */ vector< Id > poolFuncVec_; + /** + * vector tracks which pool is controlled by which func. + * Unused entries are flagged by ~0. + * funcTarget_[ poolIndex ] == funcIndex + */ + vector< unsigned int > funcTarget_; + /** * Vector of funcs controlling pool increment, that is dN/dt * This is handled as a rateTerm. diff --git a/moose-core/ksolve/VoxelPoolsBase.cpp b/moose-core/ksolve/VoxelPoolsBase.cpp index b9e4e3065bc03ae6288fee36dc2a39b47ed11864..22831c66d816d09784d4dc0cabfc960bf68031e8 100644 --- a/moose-core/ksolve/VoxelPoolsBase.cpp +++ b/moose-core/ksolve/VoxelPoolsBase.cpp @@ -118,8 +118,11 @@ void VoxelPoolsBase::scaleVolsBufsRates( unsigned int start = stoichPtr_->getNumVarPools(); unsigned int end = start + stoichPtr_->getNumBufPools(); assert( end == Sinit_.size() ); - for ( unsigned int i = start; i < end; ++i ) - S_[i] = Sinit_[i]; + for ( unsigned int i = start; i < end; ++i ) { + // Must not reassign pools that are controlled by functions. + if ( !stoichPtr->isFuncTarget(i) ) + S_[i] = Sinit_[i]; + } // Scale rates. Start by clearing out old rates if any for ( unsigned int i = 0; i < rates_.size(); ++i ) diff --git a/moose-core/ksolve/ZombieBufPool.cpp b/moose-core/ksolve/ZombieBufPool.cpp index a00a6499c7fb5a5351a48ba6270c8ed21078287c..826655299061397f79ccd60177a83c4310831960 100644 --- a/moose-core/ksolve/ZombieBufPool.cpp +++ b/moose-core/ksolve/ZombieBufPool.cpp @@ -72,3 +72,7 @@ void ZombieBufPool::vSetConcInit( const Eref& e, double conc ) vSetConc( e, conc ); } +bool ZombieBufPool::vGetIsBuffered( const Eref& e ) const +{ + return true; +} diff --git a/moose-core/ksolve/ZombieBufPool.h b/moose-core/ksolve/ZombieBufPool.h index 8fda4df18c82e4d82656f2ab5695766bbbab142f..a00344e3cabce63b970c6583dc5a6eb8681f8126 100644 --- a/moose-core/ksolve/ZombieBufPool.h +++ b/moose-core/ksolve/ZombieBufPool.h @@ -21,6 +21,7 @@ class ZombieBufPool: public ZombiePool void vSetNinit( const Eref& e, double v ); void vSetConc( const Eref& e, double v ); void vSetConcInit( const Eref& e, double v ); + bool vGetIsBuffered( const Eref& e ) const; static const Cinfo* initCinfo(); private: diff --git a/moose-core/ksolve/ZombiePool.cpp b/moose-core/ksolve/ZombiePool.cpp index b87beb6e61e8416e8c1159dbf611cb5d12932176..6d612c8332bdbc6e1bae085238874ed2f5b81456 100644 --- a/moose-core/ksolve/ZombiePool.cpp +++ b/moose-core/ksolve/ZombiePool.cpp @@ -179,6 +179,11 @@ double ZombiePool::vGetVolume( const Eref& e ) const return lookupVolumeFromMesh( e ); } +bool ZombiePool::vGetIsBuffered( const Eref& e ) const +{ + return false; +} + ////////////////////////////////////////////////////////////// // Zombie conversion functions. ////////////////////////////////////////////////////////////// diff --git a/moose-core/ksolve/ZombiePool.h b/moose-core/ksolve/ZombiePool.h index d0cefe8fd467fe4733c319fca0fc47af1608165f..147823fcca7aa4019699b42f58a60ea3a46ad35b 100644 --- a/moose-core/ksolve/ZombiePool.h +++ b/moose-core/ksolve/ZombiePool.h @@ -46,6 +46,7 @@ class ZombiePool: public PoolBase void vSetMotorConst( const Eref& e, double v ); double vGetMotorConst( const Eref& e ) const; + bool vGetIsBuffered( const Eref& e ) const; ////////////////////////////////////////////////////////////////// // Dest funcs ////////////////////////////////////////////////////////////////// diff --git a/moose-core/python/moose/OrderedDict.py b/moose-core/python/moose/OrderedDict.py new file mode 100644 index 0000000000000000000000000000000000000000..cf7e171642da148f9ac792e0dd03a11009a47919 --- /dev/null +++ b/moose-core/python/moose/OrderedDict.py @@ -0,0 +1,261 @@ +# -*- coding: utf-8 -*- +# Backport of OrderedDict() class that runs on Python 2.4, 2.5, 2.6, 2.7 and pypy. +# Passes Python2.7's test suite and incorporates all the latest updates. + +import sys + +if sys.version_info > (3, 0 ): + from _thread import get_ident as _get_ident +else: + from thread import get_ident as _get_ident + +try: + from _abcoll import KeysView, ValuesView, ItemsView +except ImportError: + pass + + +class OrderedDict(dict): + 'Dictionary that remembers insertion order' + # An inherited dict maps keys to values. + # The inherited dict provides __getitem__, __len__, __contains__, and get. + # The remaining methods are order-aware. + # Big-O running times for all methods are the same as for regular dictionaries. + + # The internal self.__map dictionary maps keys to links in a doubly linked list. + # The circular doubly linked list starts and ends with a sentinel element. + # The sentinel element never gets deleted (this simplifies the algorithm). + # Each link is stored as a list of length three: [PREV, NEXT, KEY]. + + def __init__(self, *args, **kwds): + '''Initialize an ordered dictionary. Signature is the same as for + regular dictionaries, but keyword arguments are not recommended + because their insertion order is arbitrary. + + ''' + if len(args) > 1: + raise TypeError('expected at most 1 arguments, got %d' % len(args)) + try: + self.__root + except AttributeError: + self.__root = root = [] # sentinel node + root[:] = [root, root, None] + self.__map = {} + self.__update(*args, **kwds) + + def __setitem__(self, key, value, dict_setitem=dict.__setitem__): + 'od.__setitem__(i, y) <==> od[i]=y' + # Setting a new item creates a new link which goes at the end of the linked + # list, and the inherited dictionary is updated with the new key/value pair. + if key not in self: + root = self.__root + last = root[0] + last[1] = root[0] = self.__map[key] = [last, root, key] + dict_setitem(self, key, value) + + def __delitem__(self, key, dict_delitem=dict.__delitem__): + 'od.__delitem__(y) <==> del od[y]' + # Deleting an existing item uses self.__map to find the link which is + # then removed by updating the links in the predecessor and successor nodes. + dict_delitem(self, key) + link_prev, link_next, key = self.__map.pop(key) + link_prev[1] = link_next + link_next[0] = link_prev + + def __iter__(self): + 'od.__iter__() <==> iter(od)' + root = self.__root + curr = root[1] + while curr is not root: + yield curr[2] + curr = curr[1] + + def __reversed__(self): + 'od.__reversed__() <==> reversed(od)' + root = self.__root + curr = root[0] + while curr is not root: + yield curr[2] + curr = curr[0] + + def clear(self): + 'od.clear() -> None. Remove all items from od.' + try: + for node in self.__map.itervalues(): + del node[:] + root = self.__root + root[:] = [root, root, None] + self.__map.clear() + except AttributeError: + pass + dict.clear(self) + + def popitem(self, last=True): + '''od.popitem() -> (k, v), return and remove a (key, value) pair. + Pairs are returned in LIFO order if last is true or FIFO order if false. + + ''' + if not self: + raise KeyError('dictionary is empty') + root = self.__root + if last: + link = root[0] + link_prev = link[0] + link_prev[1] = root + root[0] = link_prev + else: + link = root[1] + link_next = link[1] + root[1] = link_next + link_next[0] = root + key = link[2] + del self.__map[key] + value = dict.pop(self, key) + return key, value + + # -- the following methods do not depend on the internal structure -- + + def keys(self): + 'od.keys() -> list of keys in od' + return list(self) + + def values(self): + 'od.values() -> list of values in od' + return [self[key] for key in self] + + def items(self): + 'od.items() -> list of (key, value) pairs in od' + return [(key, self[key]) for key in self] + + def iterkeys(self): + 'od.iterkeys() -> an iterator over the keys in od' + return iter(self) + + def itervalues(self): + 'od.itervalues -> an iterator over the values in od' + for k in self: + yield self[k] + + def iteritems(self): + 'od.iteritems -> an iterator over the (key, value) items in od' + for k in self: + yield (k, self[k]) + + def update(*args, **kwds): + '''od.update(E, **F) -> None. Update od from dict/iterable E and F. + + If E is a dict instance, does: for k in E: od[k] = E[k] + If E has a .keys() method, does: for k in E.keys(): od[k] = E[k] + Or if E is an iterable of items, does: for k, v in E: od[k] = v + In either case, this is followed by: for k, v in F.items(): od[k] = v + + ''' + if len(args) > 2: + raise TypeError('update() takes at most 2 positional ' + 'arguments (%d given)' % (len(args),)) + elif not args: + raise TypeError('update() takes at least 1 argument (0 given)') + self = args[0] + # Make progressively weaker assumptions about "other" + other = () + if len(args) == 2: + other = args[1] + if isinstance(other, dict): + for key in other: + self[key] = other[key] + elif hasattr(other, 'keys'): + for key in other.keys(): + self[key] = other[key] + else: + for key, value in other: + self[key] = value + for key, value in kwds.items(): + self[key] = value + + __update = update # let subclasses override update without breaking __init__ + + __marker = object() + + def pop(self, key, default=__marker): + '''od.pop(k[,d]) -> v, remove specified key and return the corresponding value. + If key is not found, d is returned if given, otherwise KeyError is raised. + + ''' + if key in self: + result = self[key] + del self[key] + return result + if default is self.__marker: + raise KeyError(key) + return default + + def setdefault(self, key, default=None): + 'od.setdefault(k[,d]) -> od.get(k,d), also set od[k]=d if k not in od' + if key in self: + return self[key] + self[key] = default + return default + + def __repr__(self, _repr_running={}): + 'od.__repr__() <==> repr(od)' + call_key = id(self), _get_ident() + if call_key in _repr_running: + return '...' + _repr_running[call_key] = 1 + try: + if not self: + return '%s()' % (self.__class__.__name__,) + return '%s(%r)' % (self.__class__.__name__, self.items()) + finally: + del _repr_running[call_key] + + def __reduce__(self): + 'Return state information for pickling' + items = [[k, self[k]] for k in self] + inst_dict = vars(self).copy() + for k in vars(OrderedDict()): + inst_dict.pop(k, None) + if inst_dict: + return (self.__class__, (items,), inst_dict) + return self.__class__, (items,) + + def copy(self): + 'od.copy() -> a shallow copy of od' + return self.__class__(self) + + @classmethod + def fromkeys(cls, iterable, value=None): + '''OD.fromkeys(S[, v]) -> New ordered dictionary with keys from S + and values equal to v (which defaults to None). + + ''' + d = cls() + for key in iterable: + d[key] = value + return d + + def __eq__(self, other): + '''od.__eq__(y) <==> od==y. Comparison to another OD is order-sensitive + while comparison to a regular mapping is order-insensitive. + + ''' + if isinstance(other, OrderedDict): + return len(self)==len(other) and self.items() == other.items() + return dict.__eq__(self, other) + + def __ne__(self, other): + return not self == other + + # -- the following methods are only used in Python 2.7 -- + + def viewkeys(self): + "od.viewkeys() -> a set-like object providing a view on od's keys" + return KeysView(self) + + def viewvalues(self): + "od.viewvalues() -> an object providing a view on od's values" + return ValuesView(self) + + def viewitems(self): + "od.viewitems() -> a set-like object providing a view on od's items" + return ItemsView(self) diff --git a/moose-core/python/moose/SBML/readSBML.py b/moose-core/python/moose/SBML/readSBML.py index aa3216526126dc6c1a63ff97f4e4a778dfef178e..c5b00c528641e20764fec236712e8ea70aa167ff 100644 --- a/moose-core/python/moose/SBML/readSBML.py +++ b/moose-core/python/moose/SBML/readSBML.py @@ -1,52 +1,71 @@ -# -*- coding: utf-8 -*- - ''' ******************************************************************* -* File: readSBML.py -* Description: -* Author: HarshaRani -* E-mail: hrani@ncbs.res.in -* This program is part of 'MOOSE', the -* Messaging Object Oriented Simulation Environment, -* also known as GENESIS 3 base code. -* copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS -* -* Created : Thu May 12 10:19:00 2016(+0530) -* Version -* Last-Updated: Wed Jan 11 2017 -* By: + * File: readSBML.py + * Description: + * Author: HarshaRani + * E-mail: hrani@ncbs.res.in + ********************************************************************/ + +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment, +** also known as GENESIS 3 base code. +** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS +Created : Thu May 13 10:19:00 2016(+0530) +Version +Last-Updated: Wed Oct 4 14:50:00 2017(+0530) + + By:HarshaRani **********************************************************************/ - -TODO in -- Compartment - --Need to add group - --Need to deal with compartment outside -- Molecule - -- Need to add group - -- mathML only AssisgmentRule is taken partly I have checked addition and multiplication, - --, need to do for other calculation. - -- In Assisgment rule one of the variable is a function, in moose since assignment is done using function, - function can't get input from another function (model 000740 in l3v1) -- Loading Model from SBML - --Tested 1-30 testcase example model provided by l3v1 and l2v4 std. - ---These are the models that worked (sbml testcase)1-6,10,14-15,17-21,23-25,34,35,58 -- Need to check - ----what to do when boundarycondition is true i.e., - differential equation derived from the reaction definitions - should not be calculated for the species(7-9,11-13,16) - ----kineticsLaw, Math fun has fraction,ceiling,reminder,power 28etc. - ----Events to be added 26 - ----initial Assisgment for compartment 27 - ----when stoichiometry is rational number 22 - ---- For Michaelis Menten kinetics km is not defined which is most of the case need to calculate - +2017 +Oct 4 : - loadpath is cleaned up +Sep 13: - After EnzymaticReaction's k2 is set, explicity ratio is set to 4 to make sure it balance. + - If units are defined in the rate law for the reaction then check is made and if not in milli mole the base unit + then converted to milli unit + - Moose doesn't allow Michaelis-Menten Enz to have more than one substrates/product +Sep 12: - Now mooseReadSBML return model and errorFlag + - check's are made if model is valid if its not errorFlag is set + - check if model has atleast one compartment, if not errorFlag is set + - errorFlag is set for Rules (for now piecewise is set which is not read user are warned) + - rateLaw are also calculated depending on units and number of substrates/product + +Sep 8 : -functionDefinitions is read, + - if Kf and Kb unit are not defined then checked if substance units is defined and depending on this unit Kf and Kb is calculated + -kf and kb units is not defined and if substance units is also not defined validator fails +Aug 9 : a check made to for textColor while adding to Annotator +Aug 8 : removed "findCompartment" function to chemConnectUtil and imported the function from the same file + + TODO in + -Compartment + --Need to deal with compartment outside + -Molecule + -- mathML only AssisgmentRule is taken partly I have checked addition and multiplication, + -- concentration as piecewise (like tertiary operation with time ) + -- need to do for other calculation. + -- In Assisgment rule one of the variable is a function, in moose since assignment is done using function, + function can't get input from another function (model 000740 in l3v1) + -Loading Model from SBML + --Tested 1-30 testcase example model provided by l3v1 and l2v4 std. + ---These are the models that worked (sbml testcase)1-6,10,14-15,17-21,23-25,34,35,58 + ---Need to check + ----what to do when boundarycondition is true i.e., + differential equation derived from the reaction definitions + should not be calculated for the species(7-9,11-13,16) + ----kineticsLaw, Math fun has fraction,ceiling,reminder,power 28etc. + ----Events to be added 26 + ----initial Assisgment for compartment 27 + ----when stoichiometry is rational number 22 + ---- For Michaelis Menten kinetics km is not defined which is most of the case need to calculate ''' + import sys -import os.path import collections import moose -from .validation import validateModel +from moose.chemUtil.chemConnectUtil import * +from moose.SBML.validation import validateModel +import re +import os foundLibSBML_ = False try: @@ -60,15 +79,19 @@ def mooseReadSBML(filepath, loadpath, solver="ee"): if not foundLibSBML_: print('No python-libsbml found.' '\nThis module can be installed by following command in terminal:' - '\n\t easy_install python-libsbl' + '\n\t easy_install python-libsbml' + '\n\t apt-get install python-libsbml' ) - return None + return moose.element('/') if not os.path.isfile(filepath): print('%s is not found ' % filepath) - return None + return moose.element('/') + with open(filepath, "r") as filep: + loadpath = loadpath[loadpath.find('/')+1:] + loaderror = None filep = open(filepath, "r") document = libsbml.readSBML(filepath) tobecontinue = False @@ -111,13 +134,14 @@ def mooseReadSBML(filepath, loadpath, solver="ee"): print("\n") if (model.getNumCompartments() == 0): - return moose.element('/') + return moose.element('/'), "Atleast one compartment is needed" else: + loadpath ='/'+loadpath baseId = moose.Neutral(loadpath) + basePath = baseId # All the model will be created under model as # a thumbrule - basePath = moose.Neutral( - baseId.path + '/model') + basePath = moose.Neutral(baseId.path) # Map Compartment's SBML id as key and value is # list of[ Moose ID and SpatialDimensions ] global comptSbmlidMooseIdMap @@ -125,41 +149,103 @@ def mooseReadSBML(filepath, loadpath, solver="ee"): warning = " " global msg msg = " " + msgRule = "" + msgReac = "" + noRE = "" + groupInfo = {} + funcDef = {} + modelAnnotaInfo = {} comptSbmlidMooseIdMap = {} - #print(("modelPath:" + basePath.path)) globparameterIdValue = {} - modelAnnotaInfo = {} + mapParameter(model, globparameterIdValue) errorFlag = createCompartment( basePath, model, comptSbmlidMooseIdMap) + + groupInfo = checkGroup(basePath,model) + funcDef = checkFuncDef(model) if errorFlag: specInfoMap = {} errorFlag,warning = createSpecies( - basePath, model, comptSbmlidMooseIdMap, specInfoMap, modelAnnotaInfo) - #print(errorFlag,warning) + basePath, model, comptSbmlidMooseIdMap, specInfoMap, modelAnnotaInfo,groupInfo) if errorFlag: - errorFlag = createRules( - model, specInfoMap, globparameterIdValue) + msgRule = createRules( + model, specInfoMap, globparameterIdValue) if errorFlag: - errorFlag, msg = createReaction( - model, specInfoMap, modelAnnotaInfo, globparameterIdValue) + errorFlag, msgReac = createReaction( + model, specInfoMap, modelAnnotaInfo, globparameterIdValue,funcDef,groupInfo) + if len(moose.wildcardFind(moose.element(loadpath).path+"/##[ISA=ReacBase],/##[ISA=EnzBase]")) == 0: + errorFlag = False + noRE = ("Atleast one reaction should be present ") getModelAnnotation( model, baseId, basePath) if not errorFlag: - print(msg) # Any time in the middle if SBML does not read then I # delete everything from model level This is important # as while reading in GUI the model will show up untill # built which is not correct print "Deleted rest of the # model" moose.delete(basePath) - basePath = moose.Shell('/') - return basePath + loadpath = moose.Shell('/') + #return basePath, "" + loaderror = msgRule+msgReac+noRE + if loaderror != "": + loaderror = loaderror +" to display in the widget" + return moose.element(loadpath), loaderror else: print("Validation failed while reading the model.") - return moose.element('/') + return moose.element('/'), "This document is not valid SBML" + +def checkFuncDef(model): + funcDef = {} + for findex in range(0,model.getNumFunctionDefinitions()): + bvar = [] + funMathML = "" + foundbvar = False + foundfuncMathML = False + f = model.getFunctionDefinition(findex) + fmath = f.getMath() + for i in range(0,fmath.getNumBvars()): + bvar.append(fmath.getChild(i).getName()) + foundbvar = True + funcMathML = fmath.getRightChild() + if fmath.getRightChild(): + foundfuncMathML = True + if foundbvar and foundfuncMathML: + funcDef[f.getName()] = {'bvar':bvar, "MathML": fmath.getRightChild()} + return funcDef +def checkGroup(basePath,model): + groupInfo = {} + modelAnnotaInfo = {} + if model.getPlugin("groups") != None: + mplugin = model.getPlugin("groups") + modelgn = mplugin.getNumGroups() + for gindex in range(0, mplugin.getNumGroups()): + p = mplugin.getGroup(gindex) + groupAnnoInfo = {} + groupAnnoInfo = getObjAnnotation(p, modelAnnotaInfo) + if groupAnnoInfo != {}: + if moose.exists(basePath.path+'/'+groupAnnoInfo["Compartment"]): + if not moose.exists(basePath.path+'/'+groupAnnoInfo["Compartment"]+'/'+p.getId()): + moosegrp = moose.Neutral(basePath.path+'/'+groupAnnoInfo["Compartment"]+'/'+p.getId()) + else: + moosegrp = moose.element(basePath.path+'/'+groupAnnoInfo["Compartment"]+'/'+p.getId()) + moosegrpinfo = moose.Annotator(moosegrp.path+'/info') + moosegrpinfo.color = groupAnnoInfo["bgColor"] + else: + print ("Compartment not found") + if p.getKind() == 2: + if p.getId() not in groupInfo: + #groupInfo[p.getId()] + for gmemIndex in range(0,p.getNumMembers()): + mem = p.getMember(gmemIndex) + if p.getId() in groupInfo: + groupInfo[p.getId()].append(mem.getIdRef()) + else: + groupInfo[p.getId()] =[mem.getIdRef()] + return groupInfo def setupEnzymaticReaction(enz, groupName, enzName, specInfoMap, modelAnnotaInfo): enzPool = (modelAnnotaInfo[groupName]["enzyme"]) @@ -168,7 +254,6 @@ def setupEnzymaticReaction(enz, groupName, enzName, cplx = (modelAnnotaInfo[groupName]["complex"]) cplx = str(idBeginWith(cplx)) complx = moose.element(specInfoMap[cplx]["Mpath"].path) - enzyme_ = moose.Enz(enzParent.path + '/' + enzName) moose.move(complx, enzyme_) moose.connect(enzyme_, "cplx", complx, "reac") @@ -203,7 +288,10 @@ def addSubPrd(reac, reName, type, reactSBMLIdMooseId, specInfoMap): for rt in range(0, reac.getNumReactants()): rct = reac.getReactant(rt) sp = rct.getSpecies() - rctMapIter[sp] = rct.getStoichiometry() + if rct.isSetStoichiometry(): + rctMapIter[sp] = rct.getStoichiometry() + else: + rctMapIter[sp] = 1 if rct.getStoichiometry() > 1: pass # print " stoich ",reac.name,rct.getStoichiometry() @@ -223,7 +311,11 @@ def addSubPrd(reac, reName, type, reactSBMLIdMooseId, specInfoMap): for rt in range(0, reac.getNumProducts()): rct = reac.getProduct(rt) sp = rct.getSpecies() - rctMapIter[sp] = rct.getStoichiometry() + if rct.isSetStoichiometry(): + rctMapIter[sp] = rct.getStoichiometry() + else: + rctMapIter[sp] = 1 + if rct.getStoichiometry() > 1: pass # print " stoich prd",reac.name,rct.getStoichiometry() @@ -284,7 +376,7 @@ def getModelAnnotation(obj, baseId, basepath): for plots in plotlist: plots = plots.replace(" ", "") plotorg = plots - if moose.exists(basepath.path + plotorg): + if( moose.exists(basepath.path + plotorg) and isinstance(moose.element(basepath.path+plotorg),moose.PoolBase)) : plotSId = moose.element( basepath.path + plotorg) # plotorg = convertSpecialChar(plotorg) @@ -299,8 +391,7 @@ def getModelAnnotation(obj, baseId, basepath): if not fullPath in tablelistname: tab = moose.Table2(fullPath) tablelistname.append(fullPath) - moose.connect( - tab, "requestOut", plotSId, "getConc") + moose.connect(tab, "requestOut", plotSId, "getConc") def getObjAnnotation(obj, modelAnnotationInfo): @@ -309,11 +400,11 @@ def getObjAnnotation(obj, modelAnnotationInfo): # modelAnnotaInfo= {} annotateMap = {} if (obj.getAnnotation() is not None): + annoNode = obj.getAnnotation() for ch in range(0, annoNode.getNumChildren()): childNode = annoNode.getChild(ch) - if (childNode.getPrefix() == "moose" and (childNode.getName() == - "ModelAnnotation" or childNode.getName() == "EnzymaticReaction")): + if (childNode.getPrefix() == "moose" and (childNode.getName() in["ModelAnnotation","EnzymaticReaction","GroupAnnotation"])): sublist = [] for gch in range(0, childNode.getNumChildren()): grandChildNode = childNode.getChild(gch) @@ -324,7 +415,6 @@ def getObjAnnotation(obj, modelAnnotationInfo): else: print( "Error: expected exactly ONE child of ", nodeName) - if nodeName == "xCord": annotateMap[nodeName] = nodeValue if nodeName == "yCord": @@ -333,11 +423,15 @@ def getObjAnnotation(obj, modelAnnotationInfo): annotateMap[nodeName] = nodeValue if nodeName == "textColor": annotateMap[nodeName] = nodeValue + if nodeName == "Group": + annotateMap[nodeName] = nodeValue + if nodeName == "Compartment": + annotateMap[nodeName] = nodeValue return annotateMap def getEnzAnnotation(obj, modelAnnotaInfo, rev, - globparameterIdValue, specInfoMap): + globparameterIdValue, specInfoMap,funcDef): name = obj.getId() name = name.replace(" ", "_space_") # modelAnnotaInfo= {} @@ -383,8 +477,6 @@ def getEnzAnnotation(obj, modelAnnotaInfo, rev, groupName = list(annotateMap["grpName"])[0] klaw = obj.getKineticLaw() mmsg = "" - errorFlag, mmsg, k1, k2 = getKLaw( - obj, klaw, rev, globparameterIdValue, specInfoMap) if 'substrates' in annotateMap: sublist = list(annotateMap["substrates"]) @@ -394,6 +486,10 @@ def getEnzAnnotation(obj, modelAnnotaInfo, rev, prdlist = list(annotateMap["product"]) else: prdlist = {} + noOfsub = len(sublist) + noOfprd = len(prdlist) + errorFlag, mmsg, k1, k2 = getKLaw( + obj, klaw, noOfsub, noOfprd, rev, globparameterIdValue, funcDef,specInfoMap) if list(annotateMap["stage"])[0] == '1': if groupName in modelAnnotaInfo: @@ -438,7 +534,7 @@ def getEnzAnnotation(obj, modelAnnotaInfo, rev, return(groupName) -def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): +def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue,funcDef,groupInfo): # print " reaction " # Things done for reaction # --Reaction is not created, if substrate and product is missing @@ -451,15 +547,25 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): errorFlag = True reactSBMLIdMooseId = {} msg = "" - rName = "" reaction_ = None for ritem in range(0, model.getNumReactions()): reactionCreated = False groupName = "" + rName = "" + rId = "" reac = model.getReaction(ritem) + group = "" + reacAnnoInfo = {} + reacAnnoInfo = getObjAnnotation(reac, modelAnnotaInfo) + # if "Group" in reacAnnoInfo: + # group = reacAnnoInfo["Group"] + if (reac.isSetId()): rId = reac.getId() + groups = [k for k, v in groupInfo.iteritems() if rId in v] + if groups: + group = groups[0] if (reac.isSetName()): rName = reac.getName() rName = rName.replace(" ", "_space_") @@ -474,39 +580,45 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): "\"") if (reac.getAnnotation() is not None): groupName = getEnzAnnotation( - reac, modelAnnotaInfo, rev, globparameterIdValue, specInfoMap) + reac, modelAnnotaInfo, rev, globparameterIdValue, specInfoMap,funcDef) if (groupName != "" and list( modelAnnotaInfo[groupName]["stage"])[0] == 3): reaction_, reactionCreated = setupEnzymaticReaction( reac, groupName, rName, specInfoMap, modelAnnotaInfo) reaction_.k3 = modelAnnotaInfo[groupName]['k3'] - reaction_.k2 = modelAnnotaInfo[groupName]['k2'] reaction_.concK1 = modelAnnotaInfo[groupName]['k1'] + reaction_.k2 = modelAnnotaInfo[groupName]['k2'] + reaction_.ratio = 4 if reactionCreated: if (reac.isSetNotes): pullnotes(reac, reaction_) reacAnnoInfo = {} reacAnnoInfo = getObjAnnotation(reac, modelAnnotaInfo) - if reacAnnoInfo: - if not moose.exists(reaction_.path + '/info'): - reacInfo = moose.Annotator(reaction_.path + '/info') - else: - reacInfo = moose.element(reaction_.path + '/info') - for k, v in list(reacAnnoInfo.items()): - if k == 'xCord': - reacInfo.x = float(v) - elif k == 'yCord': - reacInfo.y = float(v) - elif k == 'bgColor': - reacInfo.color = v - else: - reacInfo.textColor = v + + #if reacAnnoInfo.keys() in ['xCord','yCord','bgColor','Color']: + if not moose.exists(reaction_.path + '/info'): + reacInfo = moose.Annotator(reaction_.path + '/info') + else: + reacInfo = moose.element(reaction_.path + '/info') + for k, v in list(reacAnnoInfo.items()): + if k == 'xCord': + reacInfo.x = float(v) + elif k == 'yCord': + reacInfo.y = float(v) + elif k == 'bgColor': + reacInfo.color = v + elif k == 'Color': + reacInfo.textColor = v elif(groupName == ""): + numRcts = reac.getNumReactants() numPdts = reac.getNumProducts() nummodifiers = reac.getNumModifiers() + # if (nummodifiers > 0 and (numRcts > 1 or numPdts >1)): + # print("Warning: %s" %(rName)," : Enzymatic Reaction has more than one Substrate or Product which is not allowed in moose, we will be skiping creating this reaction in MOOSE") + # reactionCreated = False if not (numRcts and numPdts): print("Warning: %s" %(rName)," : Substrate or Product is missing, we will be skiping creating this reaction in MOOSE") @@ -527,6 +639,11 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): sp = react.getSpecies() sp = str(idBeginWith(sp)) speCompt = specInfoMap[sp]["comptId"].path + if group: + if moose.exists(speCompt+'/'+group): + speCompt = speCompt+'/'+group + else: + speCompt = (moose.Neutral(speCompt+'/'+group)).path reaction_ = moose.Reac(speCompt + '/' + rName) reactionCreated = True reactSBMLIdMooseId[rName] = { @@ -548,20 +665,21 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): pullnotes(reac, reaction_) reacAnnoInfo = {} reacAnnoInfo = getObjAnnotation(reac, modelAnnotaInfo) - if reacAnnoInfo: - if not moose.exists(reaction_.path + '/info'): - reacInfo = moose.Annotator(reaction_.path + '/info') - else: - reacInfo = moose.element(reaction_.path + '/info') - for k, v in list(reacAnnoInfo.items()): - if k == 'xCord': - reacInfo.x = float(v) - elif k == 'yCord': - reacInfo.y = float(v) - elif k == 'bgColor': - reacInfo.color = v - else: - reacInfo.textColor = v + + #if reacAnnoInfo.keys() in ['xCord','yCord','bgColor','Color']: + if not moose.exists(reaction_.path + '/info'): + reacInfo = moose.Annotator(reaction_.path + '/info') + else: + reacInfo = moose.element(reaction_.path + '/info') + for k, v in list(reacAnnoInfo.items()): + if k == 'xCord': + reacInfo.x = float(v) + elif k == 'yCord': + reacInfo.y = float(v) + elif k == 'bgColor': + reacInfo.color = v + elif k == 'Color': + reacInfo.textColor = v addSubPrd(reac, rName, "sub", reactSBMLIdMooseId, specInfoMap) addSubPrd(reac, rName, "prd", reactSBMLIdMooseId, specInfoMap) @@ -569,7 +687,7 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): klaw = reac.getKineticLaw() mmsg = "" errorFlag, mmsg, kfvalue, kbvalue = getKLaw( - model, klaw, rev, globparameterIdValue, specInfoMap) + model, klaw, reac.num_reactants,reac.num_products, rev, globparameterIdValue,funcDef,specInfoMap) if not errorFlag: msg = "Error while importing reaction \"" + \ rName + "\"\n Error in kinetics law " @@ -577,9 +695,6 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): msg = msg + mmsg return(errorFlag, msg) else: - # print " reactSBMLIdMooseId - # ",reactSBMLIdMooseId[rName]["nSub"], " prd - # ",reactSBMLIdMooseId[rName]["nPrd"] if reaction_.className == "Reac": subn = reactSBMLIdMooseId[rName]["nSub"] prdn = reactSBMLIdMooseId[rName]["nPrd"] @@ -591,7 +706,7 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue): return (errorFlag, msg) -def getKLaw(model, klaw, rev, globparameterIdValue, specMapList): +def getKLaw(model, klaw,noOfsub, noOfprd,rev, globparameterIdValue, funcDef, specMapList): parmValueMap = {} amt_Conc = "amount" value = 0.0 @@ -610,8 +725,8 @@ def getKLaw(model, klaw, rev, globparameterIdValue, specMapList): kbparm = "" kfvalue = 0 kbvalue = 0 - kfp = None - kbp = None + kfp = "" + kbp = "" mssgstr = "" for i in ruleMemlist: if i in parmValueMap or i in globparameterIdValue: @@ -632,24 +747,127 @@ def getKLaw(model, klaw, rev, globparameterIdValue, specMapList): kbvalue = globparameterIdValue[i] kbp = model.getParameter(kbparm) index += 1 - elif not (i in specMapList or i in comptSbmlidMooseIdMap): mssgstr = "\"" + i + "\" is not defined " - return (False, mssgstr) + return (False, mssgstr, 0,0) + if kfp != "": - # print " unit set for rate law kfp ",kfparm, " ",kfp.isSetUnits() + lvalue =1.0 if kfp.isSetUnits(): kfud = kfp.getDerivedUnitDefinition() - # print " kfud ",kfud + lvalue = transformUnits( 1,kfud ,"substance",True) + else: + unitscale = 1 + if (noOfsub >1): + #converting units to milli M for Moose + #depending on the order of reaction,millifactor will calculated + unitscale = unitsforRates(model) + unitscale = unitscale*1000 + lvalue = pow(unitscale,1-noOfsub) + kfvalue = kfvalue*lvalue; + if kbp != "": - pass - # print " unit set for rate law kbp ",kbparm, " ",kbp.isSetUnits() - + lvalue = 1.0; + if kbp.isSetUnits(): + kbud = kbp.getDerivedUnitDefinition() + else: + unitscale = 1 + if (noOfprd >1): + unitscale = unitsforRates(model) + unitscale = unitscale*1000 + lvalue = pow(unitscale,1-noOfprd) + kbvalue = kbvalue*lvalue; return (True, mssgstr, kfvalue, kbvalue) - +def transformUnits( mvalue, ud, type, hasonlySubUnit ): + lvalue = mvalue; + if (type == "compartment"): + if(ud.getNumUnits() == 0): + unitsDefined = False + else: + for ut in range(0, ud.getNumUnits()): + unit = ud.getUnit(ut) + if ( unit.isLitre() ): + exponent = unit.getExponent() + multiplier = unit.getMultiplier() + scale = unit.getScale() + offset = unit.getOffset() + lvalue *= pow( multiplier * pow(10.0,scale), exponent ) + offset + # Need to check if spatial dimension is less than 3 then, + # then volume conversion e-3 to convert cubicmeter shd not be done. + #lvalue *= pow(1e-3,exponent) + + elif(type == "substance"): + exponent = 1.0 + if(ud.getNumUnits() == 0): + unitsDefined = False + else: + for ut in range(0, ud.getNumUnits()): + unit = ud.getUnit(ut) + if ( unit.isMole() ): + exponent = unit.getExponent() + multiplier = unit.getMultiplier() + scale = unit.getScale() + offset = unit.getOffset() + #scale+3 is to convert to milli moles for moose + lvalue *= pow( multiplier * pow(10.0,scale+3), exponent ) + offset + + elif (unit.isItem()): + exponent = unit.getExponent() + multiplier = unit.getMultiplier() + scale = unit.getScale() + offset = unit.getOffset() + lvalue *= pow( multiplier * pow(10.0,scale), exponent ) + offset + else: + + exponent = unit.getExponent() + multiplier = unit.getMultiplier() + scale = unit.getScale() + offset = unit.getOffset() + lvalue *= pow( multiplier * pow(10.0,scale), exponent ) + offset + return lvalue + +def unitsforRates(model): + lvalue =1; + if model.getNumUnitDefinitions(): + for n in range(0,model.getNumUnitDefinitions()): + ud = model.getUnitDefinition(n) + for ut in range(0,ud.getNumUnits()): + unit = ud.getUnit(ut) + if (ud.getId() == "substance"): + if ( unit.isMole() ): + exponent = unit.getExponent(); + multiplier = unit.getMultiplier(); + scale = unit.getScale(); + offset = unit.getOffset(); + lvalue *= pow( multiplier * pow(10.0,scale), exponent ) + offset; + return lvalue + else: + return lvalue def getMembers(node, ruleMemlist): - if node.getType() == libsbml.AST_PLUS: + msg = "" + found = True + if node.getType() == libsbml.AST_LAMBDA: + #In lambda get Bvar values and getRighChild which will be kineticLaw + if node.getNumBvars() == 0: + print ("0") + return False + + for i in range (0,node.getNumBvars()): + ruleMemlist.append(node.getChild(i).getName()) + #funcD[funcName] = {"bvar" : bvar, "MathML":node.getRightChild()} + elif node.getType() == libsbml.AST_FUNCTION: + #funcName = node.getName() + #funcValue = [] + #functionfound = False + for i in range(0,node.getNumChildren()): + #functionfound = True + #print " $$ ",node.getChild(i).getName() + #funcValue.append(node.getChild(i).getName()) + getMembers(node.getChild(i),ruleMemlist) + #funcKL[node.getName()] = funcValue + + elif node.getType() == libsbml.AST_PLUS: if node.getNumChildren() == 0: print ("0") return False @@ -662,6 +880,7 @@ def getMembers(node, ruleMemlist): pass elif node.getType() == libsbml.AST_NAME: # This will be the ci term" + ruleMemlist.append(node.getName()) elif node.getType() == libsbml.AST_MINUS: if node.getNumChildren() == 0: @@ -692,81 +911,139 @@ def getMembers(node, ruleMemlist): # Multiplication getMembers(node.getChild(i), ruleMemlist) + elif node.getType() == libsbml.AST_LAMBDA: + if node.getNumChildren() == 0: + print ("0") + return False + getMembers(node.getChild(0), ruleMemlist) + for i in range(1, node.getNumChildren()): + getMembers(node.getChild(i), ruleMemlist) + elif node.getType() == libsbml.AST_FUNCTION_POWER: - pass + msg = msg + "\n moose is yet to handle \""+node.getName() + "\" operator" + found = False + elif node.getType() == libsbml.AST_FUNCTION_PIECEWISE: + #print " piecewise ", libsbml.formulaToL3String(node) + msg = msg + "\n moose is yet to handle \""+node.getName() + "\" operator" + found = False + ''' + if node.getNumChildren() == 0: + print("0") + return False + else: + lchild = node.getLeftChild() + getMembers(lchild, ruleMemlist) + rchild = node.getRightChild() + getMembers(rchild, ruleMemlist) + ''' else: - - print(" this case need to be handled", node.getType()) + #print(" this case need to be handled", node.getName()) + msg = msg + "\n moose is yet to handle \""+node.getName() + "\" operator" + found = False # if len(ruleMemlist) > 2: # print "Sorry! for now MOOSE cannot handle more than 2 parameters" # return True - + return found, msg def createRules(model, specInfoMap, globparameterIdValue): + #This section where assigment, Algebraic, rate rules are converted + #For now assignment rules are written and that too just summation of pools for now. + msg = "" + found = True for r in range(0, model.getNumRules()): rule = model.getRule(r) + comptvolume = [] + if (rule.isAssignment()): rule_variable = rule.getVariable() - rule_variable = parentSp = str(idBeginWith(rule_variable)) - poolList = specInfoMap[rule_variable]["Mpath"].path - poolsCompt = findCompartment(moose.element(poolList)) - if not isinstance(moose.element(poolsCompt), moose.ChemCompt): - return -2 - else: - if poolsCompt.name not in comptvolume: - comptvolume.append(poolsCompt.name) - - funcId = moose.Function(poolList + '/func') - - objclassname = moose.element(poolList).className - if objclassname == "BufPool" or objclassname == "ZombieBufPool": - moose.connect(funcId, 'valueOut', poolList, 'setN') - elif objclassname == "Pool" or objclassname == "ZombiePool": - # moose.connect( funcId, 'valueOut', poolList ,'increament' ) - moose.connect(funcId, 'valueOut', poolList, 'setN') - elif objclassname == "Reac" or objclassname == "ZombieReac": - moose.connect(funcId, 'valueOut', poolList, 'setNumkf') - - ruleMath = rule.getMath() - ruleMemlist = [] - speFunXterm = {} - getMembers(ruleMath, ruleMemlist) - - for i in ruleMemlist: - - if (i in specInfoMap): - i = str(idBeginWith(i)) - specMapList = specInfoMap[i]["Mpath"] - poolsCompt = findCompartment(moose.element(specMapList)) - if not isinstance(moose.element( - poolsCompt), moose.ChemCompt): - return -2 - else: - if poolsCompt.name not in comptvolume: - comptvolume.append(poolsCompt.name) - numVars = funcId.numVars - x = funcId.path + '/x[' + str(numVars) + ']' - speFunXterm[i] = 'x' + str(numVars) - moose.connect(specMapList, 'nOut', x, 'input') - funcId.numVars = numVars + 1 - - elif not(i in globparameterIdValue): - print("check the variable type ", i) - exp = rule.getFormula() - for mem in ruleMemlist: - if (mem in specInfoMap): - exp1 = exp.replace(mem, str(speFunXterm[mem])) - exp = exp1 - elif(mem in globparameterIdValue): - exp1 = exp.replace(mem, str(globparameterIdValue[mem])) - exp = exp1 + + if specInfoMap.has_key(rule_variable): + #In assignment rule only if pool exist, then that is conveted to moose as + # this can be used as summation of pool's, P1+P2+P3 etc + rule_variable = parentSp = str(idBeginWith(rule_variable)) + poolList = specInfoMap[rule_variable]["Mpath"].path + poolsCompt = findCompartment(moose.element(poolList)) + #If pool comes without a compartment which is not allowed moose + #then returning with -2 + if not isinstance(moose.element(poolsCompt), moose.ChemCompt): + return -2 else: - print("Math expression need to be checked") - exp = exp.replace(" ", "") - funcId.expr = exp.strip(" \t\n\r") - # return True + if poolsCompt.name not in comptvolume: + comptvolume.append(poolsCompt.name) + + ruleMath = rule.getMath() + #libsbml.writeMathMLToString(ruleMath) + ruleMemlist = [] + speFunXterm = {} + found, msg = getMembers(ruleMath, ruleMemlist) + + if found: + allPools = True + for i in ruleMemlist: + if specInfoMap.has_key(i): + pass + else: + allPools = False + if allPools: + #only if addition then summation works, only then I create a function in moose + # which is need to get the summation's output to a pool + funcId = moose.Function(poolList + '/func') + objclassname = moose.element(poolList).className + if objclassname == "BufPool" or objclassname == "ZombieBufPool": + moose.connect(funcId, 'valueOut', poolList, 'setN') + elif objclassname == "Pool" or objclassname == "ZombiePool": + # moose.connect( funcId, 'valueOut', poolList ,'increament' ) + moose.connect(funcId, 'valueOut', poolList, 'setN') + elif objclassname == "Reac" or objclassname == "ZombieReac": + moose.connect(funcId, 'valueOut', poolList, 'setNumkf') + for i in ruleMemlist: + if (i in specInfoMap): + i = str(idBeginWith(i)) + specMapList = specInfoMap[i]["Mpath"] + poolsCompt = findCompartment(moose.element(specMapList)) + if not isinstance(moose.element( + poolsCompt), moose.ChemCompt): + return -2 + else: + if poolsCompt.name not in comptvolume: + comptvolume.append(poolsCompt.name) + numVars = funcId.numVars + x = funcId.path + '/x[' + str(numVars) + ']' + speFunXterm[i] = 'x' + str(numVars) + moose.connect(specMapList, 'nOut', x, 'input') + funcId.numVars = numVars + 1 + + elif not(i in globparameterIdValue): + msg = msg + "check the variable name in mathML, this object neither pool or a constant \"" + i+"\" in assignmentRule " +rule.getVariable() + + exp = rule.getFormula() + exprOK = True + #print " specFunXTerm ",speFunXterm + for mem in ruleMemlist: + if (mem in specInfoMap): + #exp1 = exp.replace(mem, str(speFunXterm[mem])) + #print " mem ",mem, "$ ", speFunXterm[mem], "$$ ",exp + exp1 = re.sub(r'\b%s\b'% (mem), speFunXterm[mem], exp) + exp = exp1 + elif(mem in globparameterIdValue): + #exp1 = exp.replace(mem, str(globparameterIdValue[mem])) + exp1 = re.sub(r'\b%s\b'% (mem), globparameterIdValue[mem], exp) + exp = exp1 + else: + msg = msg +"Math expression need to be checked", exp + exprOK = False + if exprOK: + exp = exp.replace(" ", "") + funcId.expr = exp.strip(" \t\n\r") + else: + msg = msg +"\nAssisgment Rule has parameter as variable, currently moose doesn't have this capability so ignoring."\ + + rule.getVariable() + " is not converted to moose." + found = False + + + elif(rule.isRate()): print( "Warning : For now this \"", @@ -781,7 +1058,7 @@ def createRules(model, specInfoMap, globparameterIdValue): if len(comptvolume) > 1: warning = "\nFunction ", moose.element( poolList).name, " has input from different compartment which is depricated in moose and running this model cause moose to crash" - return True + return msg def pullnotes(sbmlId, mooseId): @@ -798,7 +1075,7 @@ def pullnotes(sbmlId, mooseId): def createSpecies(basePath, model, comptSbmlidMooseIdMap, - specInfoMap, modelAnnotaInfo): + specInfoMap, modelAnnotaInfo,groupInfo): # ToDo: # - Need to add group name if exist in pool # - Notes @@ -808,8 +1085,18 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap, else: for sindex in range(0, model.getNumSpecies()): spe = model.getSpecies(sindex) + group = "" + specAnnoInfo = {} + specAnnoInfo = getObjAnnotation(spe, modelAnnotaInfo) + # if "Group" in specAnnoInfo: + # group = specAnnoInfo["Group"] + sName = None sId = spe.getId() + + groups = [k for k, v in groupInfo.iteritems() if sId in v] + if groups: + group = groups[0] if spe.isSetName(): sName = spe.getName() sName = sName.replace(" ", "_space_") @@ -826,31 +1113,33 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap, hasonlySubUnit = spe.getHasOnlySubstanceUnits() # "false": is {unit of amount}/{unit of size} (i.e., concentration or density). # "true": then the value is interpreted as having a unit of amount only. - + if group: + if moose.exists(comptEl+'/'+group): + comptEl = comptEl+'/'+group + else: + comptEl = (moose.Neutral(comptEl+'/'+group)).path if (boundaryCondition): poolId = moose.BufPool(comptEl + '/' + sName) else: poolId = moose.Pool(comptEl + '/' + sName) - if (spe.isSetNotes): pullnotes(spe, poolId) - specAnnoInfo = {} - specAnnoInfo = getObjAnnotation(spe, modelAnnotaInfo) - if specAnnoInfo: - if not moose.exists(poolId.path + '/info'): - poolInfo = moose.Annotator(poolId.path + '/info') - else: - poolInfo = moose.element(poolId.path + '/info') - for k, v in list(specAnnoInfo.items()): - if k == 'xCord': - poolInfo.x = float(v) - elif k == 'yCord': - poolInfo.y = float(v) - elif k == 'bgColor': - poolInfo.color = v - else: - poolInfo.textColor = v + #if specAnnoInfo.keys() in ['xCord','yCord','bgColor','textColor']: + if not moose.exists(poolId.path + '/info'): + poolInfo = moose.Annotator(poolId.path + '/info') + else: + poolInfo = moose.element(poolId.path + '/info') + + for k, v in list(specAnnoInfo.items()): + if k == 'xCord': + poolInfo.x = float(v) + elif k == 'yCord': + poolInfo.y = float(v) + elif k == 'bgColor': + poolInfo.color = v + elif k == 'Color': + poolInfo.textColor = v specInfoMap[sId] = { "Mpath": poolId, "const": constant, @@ -888,8 +1177,8 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap, if not unitset: # print " unit is not set" unitfactor = pow(10, -3) - initvalue = initvalue * unitfactor + poolId.concInit = initvalue else: nr = model.getNumRules() @@ -902,6 +1191,7 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap, if (rule_variable == sId): found = True break + if not (found): print( "Invalid SBML: Either initialConcentration or initialAmount must be set or it should be found in assignmentRule but non happening for ", @@ -950,12 +1240,14 @@ def transformUnit(unitForObject, hasonlySubUnit=False): lvalue = lvalue * pow(6.0221409e23, 1) elif hasonlySubUnit == False: # Pool units in moose is mM - if scale > 0: - lvalue *= pow(multiplier * pow(10.0, - scale - 3), exponent) + offset - elif scale <= 0: - lvalue *= pow(multiplier * pow(10.0, - scale + 3), exponent) + offset + + lvalue = lvalue * pow(multiplier*pow(10.0,scale+3),exponent)+offset + # if scale > 0: + # lvalue *= pow(multiplier * pow(10.0, + # scale - 3), exponent) + offset + # elif scale <= 0: + # lvalue *= pow(multiplier * pow(10.0, + # scale + 3), exponent) + offset unitset = True unittype = "Mole" return (lvalue, unitset, unittype) @@ -989,8 +1281,9 @@ def transformUnit(unitForObject, hasonlySubUnit=False): def createCompartment(basePath, model, comptSbmlidMooseIdMap): # ToDoList : Check what should be done for the spaitialdimension is 2 or # 1, area or length + #print " createCompartment ",model.getNumCompartments() if not(model.getNumCompartments()): - return False, + return False, "Model has no compartment, atleast one compartment should exist to display the widget" else: for c in range(0, model.getNumCompartments()): compt = model.getCompartment(c) @@ -1055,20 +1348,21 @@ def setupMMEnzymeReaction(reac, rName, specInfoMap, reactSBMLIdMooseId, pullnotes(reac, MMEnz) reacAnnoInfo = {} reacAnnoInfo = getObjAnnotation(reac, modelAnnotaInfo) - if reacAnnoInfo: - if not moose.exists(MMEnz.path + '/info'): - reacInfo = moose.Annotator(MMEnz.path + '/info') - else: - reacInfo = moose.element(MMEnz.path + '/info') - for k, v in list(reacAnnoInfo.items()): - if k == 'xCord': - reacInfo.x = float(v) - elif k == 'yCord': - reacInfo.y = float(v) - elif k == 'bgColor': - reacInfo.color = v - else: - reacInfo.textColor = v + + #if reacAnnoInfo.keys() in ['xCord','yCord','bgColor','Color']: + if not moose.exists(MMEnz.path + '/info'): + reacInfo = moose.Annotator(MMEnz.path + '/info') + else: + reacInfo = moose.element(MMEnz.path + '/info') + for k, v in list(reacAnnoInfo.items()): + if k == 'xCord': + reacInfo.x = float(v) + elif k == 'yCord': + reacInfo.y = float(v) + elif k == 'bgColor': + reacInfo.color = v + elif k == 'Color': + reacInfo.textColor = v return(reactionCreated, MMEnz) @@ -1093,22 +1387,12 @@ def idBeginWith(name): def convertSpecialChar(str1): d = {"&": "_and", "<": "_lessthan_", ">": "_greaterthan_", "BEL": "°", "-": "_minus_", "'": "_prime_", "+": "_plus_", "*": "_star_", "/": "_slash_", "(": "_bo_", ")": "_bc_", - "[": "_sbo_", "]": "_sbc_", ".": "_dot_", " ": "_" + "[": "_sbo_", "]": "_sbc_", " ": "_" } for i, j in list(d.items()): str1 = str1.replace(i, j) return str1 - -def mooseIsInstance(element, classNames): - return moose.element(element).__class__.__name__ in classNames - - -def findCompartment(element): - while not mooseIsInstance(element, ["CubeMesh", "CyclMesh"]): - element = element.parent - return element - if __name__ == "__main__": try: sys.argv[1] diff --git a/moose-core/python/moose/SBML/validation.py b/moose-core/python/moose/SBML/validation.py index 8b6569229483931e0a7c57ab05883c704610570f..ee2b2cae90c9a239784e17710def616bdf826bc8 100644 --- a/moose-core/python/moose/SBML/validation.py +++ b/moose-core/python/moose/SBML/validation.py @@ -1,4 +1,24 @@ # -*- coding: utf-8 -*- +''' +******************************************************************* + * File: validation.py + * Description: + * Author: HarshaRani + * E-mail: hrani@ncbs.res.in + ********************************************************************/ + +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment, +** also known as GENESIS 3 base code. +** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS +Created : Thu May 12 10:19:00 2016(+0530) +Version +Last-Updated: Fri Jul 28 15:50:00 2017(+0530) + By: +**********************************************************************/ + +''' foundLibSBML_ = False try: from libsbml import * @@ -93,4 +113,4 @@ def validateModel(sbmlDoc): if __name__ == '__main__': sbmlDoc = libsbml.readSBML('00001-sbml-l3v1.xml') - validateModel(sbmlDoc) + validateModel(sbmlDoc) \ No newline at end of file diff --git a/moose-core/python/moose/SBML/writeSBML.py b/moose-core/python/moose/SBML/writeSBML.py index 82ce209f38c1a08272c2d154980e22ebd0fcf094..99fb3f21e7b678197e8b5d20336d7b6ce78073c9 100644 --- a/moose-core/python/moose/SBML/writeSBML.py +++ b/moose-core/python/moose/SBML/writeSBML.py @@ -13,25 +13,28 @@ ** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS Created : Friday May 27 12:19:00 2016(+0530) Version -Last-Updated: Wed Jan 11 15:20:00 2017(+0530) - By: +Last-Updated: Tue 8 Aug 11:10:00 2017(+0530) + By: HarshaRani **********************************************************************/ /**************************** +2017 +Aug 8 : removed "findCompartment" function to chemConnectUtil and imported the function from the same file + convertSpecialChar for setId and convertSpecialCharshot for setName. + specialChar like /,\,[,],space are not allowed as moose doesn't take +Aug 3 : Added recalculatecoordinates,cleanup in groupName + ''' import sys import re -from collections import Counter - import moose -from .validation import validateModel +from moose.SBML.validation import validateModel from moose.chemUtil.chemConnectUtil import * from moose.chemUtil.graphUtils import * # ToDo: # Table should be written -# Group's should be added # boundary condition for buffer pool having assignment statment constant # shd be false @@ -48,11 +51,12 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): if not foundLibSBML_: print('No python-libsbml found.' '\nThis module can be installed by following command in terminal:' - '\n\t easy_install python-libsbml' + '\n\t easy_install python-libsbml or' + '\n\t apt-get install python-libsbml' ) - return -1, msg + return -2, "Could not save the model in to SBML file. \nThis module can be installed by following command in terminal: \n\t easy_install python-libsbml or \n\t apt-get install python-libsbml",'' - sbmlDoc = SBMLDocument(3, 1) + #sbmlDoc = SBMLDocument(3, 1) filepath, filenameExt = os.path.split(filename) if filenameExt.find('.') != -1: filename = filenameExt[:filenameExt.find('.')] @@ -65,17 +69,36 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): spe_constTrue = [] global nameList_ nameList_ = [] + + xcord,ycord = [],[] + if moose.exists(moose.element(modelpath).path): + mObj = moose.wildcardFind(moose.element(modelpath).path+'/##[ISA=PoolBase]'+','+ + moose.element(modelpath).path+'/##[ISA=ReacBase]'+','+ + moose.element(modelpath).path+'/##[ISA=EnzBase]'+','+ + moose.element(modelpath).path+'/##[ISA=StimulusTable]') + for p in mObj: + if not isinstance(moose.element(p.parent),moose.CplxEnzBase): + if moose.exists(p.path+'/info'): + xcord.append(moose.element(p.path+'/info').x) + ycord.append(moose.element(p.path+'/info').y) + recalculatecoordinates(mObj,xcord,ycord) positionInfoexist = False - xmlns = XMLNamespaces() - xmlns.add("http://www.sbml.org/sbml/level3/version1") - xmlns.add("http://www.moose.ncbs.res.in", "moose") - xmlns.add("http://www.w3.org/1999/xhtml", "xhtml") - sbmlDoc.setNamespaces(xmlns) + + xmlns = SBMLNamespaces(3, 1) + xmlns.addNamespace("http://www.w3.org/1999/xhtml", "xhtml") + xmlns.addNamespace("http://www.moose.ncbs.res.in", "moose") + xmlns.addNamespace("http://www.sbml.org/sbml/level3/version1/groups/version1", "groups") + sbmlDoc = SBMLDocument(xmlns) + sbmlDoc.setPackageRequired("groups",bool(0)) + cremodel_ = sbmlDoc.createModel() cremodel_.setId(filename) - cremodel_.setTimeUnits("second") + cremodel_.setTimeUnits("time") cremodel_.setExtentUnits("substance") cremodel_.setSubstanceUnits("substance") + cremodel_.setVolumeUnits("volume") + cremodel_.setAreaUnits("area") + cremodel_.setLengthUnits("length") neutralNotes = "" specieslist = moose.wildcardFind(modelpath + '/##[ISA=PoolBase]') @@ -91,25 +114,41 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): cremodel_.setNotes(notesString) srcdesConnection = {} - if not bool(sceneitems): - meshEntry,xmin,xmax,ymin,ymax,positionInfoexist,sitem = setupMeshObj(modelpath) - #moose.coordinateUtil.setupItem(modelpath,srcdesConnection) - setupItem(modelpath,srcdesConnection) - if not positionInfoexist: - autoCoordinates(meshEntry,srcdesConnection) writeUnits(cremodel_) modelAnno = writeSimulationAnnotation(modelpath) if modelAnno: cremodel_.setAnnotation(modelAnno) + groupInfo = {} + compartexist, groupInfo = writeCompt(modelpath, cremodel_) - compartexist = writeCompt(modelpath, cremodel_) if compartexist == True: - species = writeSpecies( modelpath,cremodel_,sbmlDoc,sceneitems) + species = writeSpecies( modelpath,cremodel_,sbmlDoc,sceneitems,groupInfo) if species: writeFunc(modelpath, cremodel_) - writeReac(modelpath, cremodel_, sceneitems) - writeEnz(modelpath, cremodel_, sceneitems) + reacGroup = {} + writeReac(modelpath, cremodel_, sceneitems,groupInfo) + writeEnz(modelpath, cremodel_, sceneitems,groupInfo) + if groupInfo: + for key,value in groupInfo.items(): + mplugin = cremodel_.getPlugin("groups") + group = mplugin.createGroup() + name = str(idBeginWith(moose.element(key).name)) + group.setId(name) + group.setKind("collection") + ginfo = moose.element(key.path+'/info') + groupCompartment = findCompartment(key) + if ginfo.color != '': + grpAnno = "<moose:GroupAnnotation>" + grpAnno = grpAnno + "<moose:Compartment>" + groupCompartment.name + "</moose:Compartment>\n" + if ginfo.color: + grpAnno = grpAnno + "<moose:bgColor>" + ginfo.color + "</moose:bgColor>\n" + grpAnno = grpAnno + "</moose:GroupAnnotation>" + group.setAnnotation(grpAnno) + + for values in value: + member = group.createMember() + member.setIdRef(values) consistencyMessages = "" SBMLok = validateModel(sbmlDoc) if (SBMLok): @@ -118,7 +157,8 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): return True, consistencyMessages, writeTofile if (not SBMLok): - cerr << "Errors encountered " << endl + #cerr << "Errors encountered " << endl + consistencyMessages = "Errors encountered" return -1, consistencyMessages else: return False, "Atleast one compartment should exist to write SBML" @@ -127,7 +167,7 @@ def calPrime(x): prime = int((20 * (float(x - cmin) / float(cmax - cmin))) - 10) return prime -def writeEnz(modelpath, cremodel_, sceneitems): +def writeEnz(modelpath, cremodel_, sceneitems,groupInfo): for enz in moose.wildcardFind(modelpath + '/##[ISA=EnzBase]'): enzannoexist = False enzGpnCorCol = " " @@ -135,17 +175,26 @@ def writeEnz(modelpath, cremodel_, sceneitems): enzSubt = () compt = "" notesE = "" + groupName = moose.element("/") + if moose.exists(enz.path + '/info'): + groupName = moose.element("/") Anno = moose.Annotator(enz.path + '/info') notesE = Anno.notes element = moose.element(enz) ele = getGroupinfo(element) - if element.className == "Neutral" or sceneitems or Anno.x or Anno.y: + ele = findGroup_compt(element) + if ele.className == "Neutral" or sceneitems or Anno.x or Anno.y: enzannoexist = True if enzannoexist: enzAnno = "<moose:ModelAnnotation>\n" if ele.className == "Neutral": - enzGpnCorCol = "<moose:Group>" + ele.name + "</moose:Group>\n" + groupName = ele + #enzGpnCorCol = "<moose:Group>" + ele.name + "</moose:Group>\n" + # if ele.name not in groupInfo: + # groupInfo[ele.name]=[setId] + # else: + # groupInfo[ele.name].append(setId) if sceneitems: #Saved from GUI, then scene co-ordinates are passed enzGpnCorCol = enzGpnCorCol + "<moose:xCord>" + \ @@ -165,7 +214,6 @@ def writeEnz(modelpath, cremodel_, sceneitems): if Anno.textColor: enzGpnCorCol = enzGpnCorCol + "<moose:textColor>" + \ Anno.textColor + "</moose:textColor>\n" - if (enz.className == "Enz" or enz.className == "ZombieEnz"): enzyme = cremodel_.createReaction() if notesE != "": @@ -180,15 +228,22 @@ def writeEnz(modelpath, cremodel_, sceneitems): compt = comptVec.name + "_" + \ str(comptVec.getId().value) + "_" + \ str(comptVec.getDataIndex()) + "_" - - enzyme.setId(str(idBeginWith(cleanEnzname + + enzsetId = str(idBeginWith(cleanEnzname + "_" + str(enz.getId().value) + "_" + str(enz.getDataIndex()) + "_" + - "Complex_formation_"))) - enzyme.setName(cleanEnzname) + "Complex_formation_")) + enzyme.setId(enzsetId) + + if groupName != moose.element('/'): + if groupName not in groupInfo: + groupInfo[groupName]=[enzsetId] + else: + groupInfo[groupName].append(enzsetId) + enzyme.setName(str(idBeginWith(convertSpecialCharshot(enz.name)))) + #enzyme.setName(cleanEnzname) enzyme.setFast(False) enzyme.setReversible(True) k1 = enz.concK1 @@ -228,7 +283,7 @@ def writeEnz(modelpath, cremodel_, sceneitems): noofSub, sRateLaw = getSubprd(cremodel_, True, "sub", enzSubt) #rec_order = rec_order + noofSub rec_order = noofSub - rate_law = compt + " * " + rate_law + "*" + sRateLaw + rate_law = compt + " * ( " + rate_law + " * " + sRateLaw enzPrd = enz.neighbors["cplxDest"] if not enzPrd: @@ -238,7 +293,7 @@ def writeEnz(modelpath, cremodel_, sceneitems): for i in range(0, len(nameList_)): enzAnno = enzAnno + "<moose:product>" + \ nameList_[i] + "</moose:product>\n" - rate_law = rate_law + " - " + compt + "* k2" + '*' + sRateLaw + rate_law = rate_law + " - " + " k2 " + ' * ' + sRateLaw +" )" prd_order = noofPrd enzAnno = enzAnno + "<moose:groupName>" + cleanEnzname + "_" + \ @@ -261,14 +316,22 @@ def writeEnz(modelpath, cremodel_, sceneitems): unit = parmUnit(rec_order - 1, cremodel_) printParameters(kl, "k1", k1, unit) enzyme = cremodel_.createReaction() - enzyme.setId(str(idBeginWith(cleanEnzname + + enzsetIdP = str(idBeginWith(cleanEnzname + "_" + str(enz.getId().value) + "_" + str(enz.getDataIndex()) + "_" + - "Product_formation_"))) - enzyme.setName(cleanEnzname) + "Product_formation_")) + enzyme.setId(enzsetIdP) + enzyme.setName(str(idBeginWith(convertSpecialCharshot(enz.name)))) + #enzyme.setName(cleanEnzname) + if groupName != moose.element('/'): + if groupName not in groupInfo: + groupInfo[groupName]=[enzsetIdP] + else: + groupInfo[groupName].append(enzsetIdP) + enzyme.setFast(False) enzyme.setReversible(False) enzAnno2 = "<moose:EnzymaticReaction>" @@ -318,6 +381,7 @@ def writeEnz(modelpath, cremodel_, sceneitems): printParameters(kl, "k3", k3, unit) elif(enz.className == "MMenz" or enz.className == "ZombieMMenz"): + enzSub = enz.neighbors["sub"] enzPrd = enz.neighbors["prd"] if (len(enzSub) != 0 and len(enzPrd) != 0): @@ -335,13 +399,20 @@ def writeEnz(modelpath, cremodel_, sceneitems): notesStringE = "<body xmlns=\"http://www.w3.org/1999/xhtml\">\n \t \t" + \ cleanNotesE + "\n\t </body>" enzyme.setNotes(notesStringE) - enzyme.setId(str(idBeginWith(cleanEnzname + + mmenzsetId = str(idBeginWith(cleanEnzname + "_" + str(enz.getId().value) + "_" + str(enz.getDataIndex()) + - "_"))) - enzyme.setName(cleanEnzname) + "_")) + enzyme.setId(mmenzsetId) + if groupName != moose.element('/'): + if groupName not in groupInfo: + groupInfo[groupName]=[mmenzsetId] + else: + groupInfo[groupName].append(mmenzsetId) + enzyme.setName(str(idBeginWith(convertSpecialCharshot(enz.name)))) + #enzyme.setName(cleanEnzname) enzyme.setFast(False) enzyme.setReversible(True) if enzannoexist: @@ -360,14 +431,15 @@ def writeEnz(modelpath, cremodel_, sceneitems): enzPrd = enz.neighbors["prd"] noofPrd, sRateLawP = getSubprd(cremodel_, False, "prd", enzPrd) kl = enzyme.createKineticLaw() - fRate_law = "kcat *" + sRateLawS + "*" + sRateLawM + \ - "/(" + compt + " * (" + "Km" + "+" + sRateLawS + "))" + fRate_law = compt + " * ( kcat * " + sRateLawS + " * " + sRateLawM + \ + " / ( Km" + " + " + sRateLawS + "))" kl.setFormula(fRate_law) kl.setNotes( "<body xmlns=\"http://www.w3.org/1999/xhtml\">\n\t\t" + fRate_law + "\n \t </body>") - printParameters(kl, "Km", Km, "substance") + KmUnit(cremodel_) + printParameters(kl, "Km", Km, "mmole_per_litre") kcatUnit = parmUnit(0, cremodel_) printParameters(kl, "kcat", kcat, kcatUnit) @@ -378,6 +450,29 @@ def printParameters(kl, k, kvalue, unit): para.setValue(kvalue) para.setUnits(unit) +def KmUnit(cremodel_): + unit_stream = "mmole_per_litre" + lud = cremodel_.getListOfUnitDefinitions() + flag = False + for i in range(0, len(lud)): + ud = lud.get(i) + if (ud.getId() == unit_stream): + flag = True + break + if (not flag): + unitdef = cremodel_.createUnitDefinition() + unitdef.setId(unit_stream) + unit = unitdef.createUnit() + unit.setKind(UNIT_KIND_LITRE) + unit.setExponent(-1) + unit.setMultiplier(1) + unit.setScale(0) + unit = unitdef.createUnit() + unit.setKind(UNIT_KIND_MOLE) + unit.setExponent(1) + unit.setMultiplier(1) + unit.setScale(-3) + return unit_stream def parmUnit(rct_order, cremodel_): order = rct_order @@ -386,7 +481,7 @@ def parmUnit(rct_order, cremodel_): elif order == 1: unit_stream = "litre_per_mmole_per_second" elif order == 2: - unit_stream = "litre_per_mmole_sq_per_second" + unit_stream = "sq_litre_per_mmole_sq_per_second" else: unit_stream = "litre_per_mmole_" + str(rct_order) + "_per_second" @@ -405,7 +500,7 @@ def parmUnit(rct_order, cremodel_): if order != 0: unit = unitdef.createUnit() unit.setKind(UNIT_KIND_LITRE) - unit.setExponent(1) + unit.setExponent(order) unit.setMultiplier(1) unit.setScale(0) unit = unitdef.createUnit() @@ -421,6 +516,8 @@ def parmUnit(rct_order, cremodel_): unit.setScale(0) return unit_stream +def Counter(items): + return dict((i, items.count(i)) for i in items) def getSubprd(cremodel_, mobjEnz, type, neighborslist): if type == "sub": @@ -507,7 +604,7 @@ def listofname(reacSub, mobjEnz): nameList_.append(clean_name) -def writeReac(modelpath, cremodel_, sceneitems): +def writeReac(modelpath, cremodel_, sceneitems,reacGroup): for reac in moose.wildcardFind(modelpath + '/##[ISA=ReacBase]'): reacSub = reac.neighbors["sub"] reacPrd = reac.neighbors["prd"] @@ -517,13 +614,15 @@ def writeReac(modelpath, cremodel_, sceneitems): reacannoexist = False reacGpname = " " cleanReacname = convertSpecialChar(reac.name) - reaction.setId(str(idBeginWith(cleanReacname + + setId = str(idBeginWith(cleanReacname + "_" + str(reac.getId().value) + "_" + str(reac.getDataIndex()) + - "_"))) - reaction.setName(cleanReacname) + "_")) + reaction.setId(setId) + reaction.setName(str(idBeginWith(convertSpecialCharshot(reac.name)))) + #reaction.setName(cleanReacname) #Kf = reac.numKf #Kb = reac.numKb Kf = reac.Kf @@ -548,7 +647,11 @@ def writeReac(modelpath, cremodel_, sceneitems): if reacannoexist: reacAnno = "<moose:ModelAnnotation>\n" if ele.className == "Neutral": - reacAnno = reacAnno + "<moose:Group>" + ele.name + "</moose:Group>\n" + #reacAnno = reacAnno + "<moose:Group>" + ele.name + "</moose:Group>\n" + if ele not in reacGroup: + reacGroup[ele]=[setId] + else: + reacGroup[ele].append(setId) if sceneitems: #Saved from GUI, then scene co-ordinates are passed reacAnno = reacAnno + "<moose:xCord>" + \ @@ -574,7 +677,7 @@ def writeReac(modelpath, cremodel_, sceneitems): kl_s, sRL, pRL, compt = "", "", "", "" if not reacSub and not reacPrd: - print(" Reaction ", reac.name, "missing substrate and product") + print(" Reaction ", reac.name, "missing substrate and product, this is not allowed in SBML which will not be written") else: kfl = reaction.createKineticLaw() if reacSub: @@ -629,8 +732,9 @@ def writeReac(modelpath, cremodel_, sceneitems): kl_s + "\n \t </body>") kfl.setFormula(kl_s) + else: - print(" Reaction ", reac.name, "missing substrate and product") + print(" Reaction ", reac.name, "missing substrate and product, this is not allowed in SBML which will not be written") def writeFunc(modelpath, cremodel_): @@ -667,7 +771,6 @@ def writeFunc(modelpath, cremodel_): rule.setVariable(fName) rule.setFormula(expr) - def convertNotesSpecialChar(str1): d = {"&": "_and", "<": "_lessthan_", ">": "_greaterthan_", "BEL": "°"} for i, j in d.items(): @@ -676,7 +779,6 @@ def convertNotesSpecialChar(str1): str1 = str1.strip(' \t\n\r') return str1 - def getGroupinfo(element): # Note: At this time I am assuming that if group exist (incase of Genesis) # 1. for 'pool' its between compartment and pool, /modelpath/Compartment/Group/pool @@ -685,18 +787,7 @@ def getGroupinfo(element): # if /modelpath/Compartment/Group/Group1/Pool, then I check and get Group1 # And /modelpath is also a NeutralObject,I stop till I find Compartment - while not mooseIsInstance(element, ["Neutral"]) and not mooseIsInstance( - element, ["CubeMesh", "CyclMesh"]): - element = element.parent - return element - - -def mooseIsInstance(element, classNames): - return moose.element(element).__class__.__name__ in classNames - - -def findCompartment(element): - while not mooseIsInstance(element, ["CubeMesh", "CyclMesh"]): + while not mooseIsInstance(element, ["Neutral", "CubeMesh", "CyclMesh"]): element = element.parent return element @@ -707,6 +798,23 @@ def idBeginWith(name): changedName = "_" + name return changedName +def findGroup_compt(melement): + while not (mooseIsInstance(melement, ["Neutral","CubeMesh", "CyclMesh"])): + melement = melement.parent + return melement + +def convertSpecialCharshot(str1): + d = { "BEL" : "°", + "'" : "_prime_", + "\\" : "_slash_", + "/" : "_slash_", + "[" : "_sbo_", + "]" : "_sbc_", + ": " : "_" , + " " : "_" } + for i, j in d.items(): + str1 = str1.replace(i, j) + return str1 def convertSpecialChar(str1): d = {"&": "_and", "<": "_lessthan_", ">": "_greaterthan_", "BEL": "°", "-": "_minus_", "'": "_prime_", @@ -718,7 +826,7 @@ def convertSpecialChar(str1): return str1 -def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems): +def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems,speGroup): # getting all the species for spe in moose.wildcardFind(modelpath + '/##[ISA=PoolBase]'): sName = convertSpecialChar(spe.name) @@ -749,7 +857,8 @@ def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems): sName = convertSpecialChar( enzPool + "_" + enzname + "_" + sName) - s1.setName(sName) + #s1.setName(sName) + s1.setName(str(idBeginWith(convertSpecialCharshot(spe.name)))) # s1.setInitialAmount(spe.nInit) s1.setInitialConcentration(spe.concInit) s1.setCompartment(compt) @@ -798,7 +907,13 @@ def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems): if speciannoexist: speciAnno = "<moose:ModelAnnotation>\n" if ele.className == "Neutral": - speciAnno = speciAnno + "<moose:Group>" + ele.name + "</moose:Group>\n" + #speciAnno = speciAnno + "<moose:Group>" + ele.name + "</moose:Group>\n" + if ele not in speGroup: + speGroup[ele]=[spename] + else: + speGroup[ele].append(spename) + + if sceneitems: #Saved from GUI, then scene co-ordinates are passed speciAnno = speciAnno + "<moose:xCord>" + \ @@ -826,6 +941,7 @@ def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems): def writeCompt(modelpath, cremodel_): # getting all the compartments compts = moose.wildcardFind(modelpath + '/##[ISA=ChemCompt]') + groupInfo = {} for compt in compts: comptName = convertSpecialChar(compt.name) # converting m3 to litre @@ -843,10 +959,20 @@ def writeCompt(modelpath, cremodel_): c1.setSize(size) c1.setSpatialDimensions(ndim) c1.setUnits('volume') + #For each compartment get groups information along + for grp in moose.wildcardFind(compt.path+'/##[TYPE=Neutral]'): + grp_cmpt = findGroup_compt(grp.parent) + try: + value = groupInfo[moose.element(grp)] + except KeyError: + # Grp is not present + groupInfo[moose.element(grp)] = [] + + if compts: - return True + return True,groupInfo else: - return False + return False,groupInfo # write Simulation runtime,simdt,plotdt def writeSimulationAnnotation(modelpath): modelAnno = "" @@ -877,17 +1003,41 @@ def writeSimulationAnnotation(modelpath): graphSpefound = True if graphSpefound: if not plots: - #plots = ori[ori.find(q.name)-1:len(ori)] - plots = '/' + q.name + '/' + name + plots = ori[ori.find(q.name)-1:len(ori)] + #plots = '/' + q.name + '/' + name else: - #plots = plots + "; "+ori[ori.find(q.name)-1:len(ori)] - plots = plots + "; /" + q.name + '/' + name + plots = plots + "; "+ori[ori.find(q.name)-1:len(ori)] + #plots = plots + "; /" + q.name + '/' + name if plots != " ": modelAnno = modelAnno + "<moose:plots> " + plots + "</moose:plots>\n" modelAnno = modelAnno + "</moose:ModelAnnotation>" return modelAnno +def xyPosition(objInfo,xory): + try: + return(float(moose.element(objInfo).getField(xory))) + except ValueError: + return (float(0)) + +def recalculatecoordinates(mObjlist,xcord,ycord): + positionInfoExist = not(len(np.nonzero(xcord)[0]) == 0 \ + and len(np.nonzero(ycord)[0]) == 0) + + if positionInfoExist: + #Here all the object has been taken now recalculate and reassign back x and y co-ordinates + xmin = min(xcord) + xmax = max(xcord) + ymin = min(ycord) + ymax = max(ycord) + for merts in mObjlist: + objInfo = merts.path+'/info' + if moose.exists(objInfo): + Ix = (xyPosition(objInfo,'x')-xmin)/(xmax-xmin) + Iy = (ymin-xyPosition(objInfo,'y'))/(ymax-ymin) + moose.element(objInfo).x = Ix + moose.element(objInfo).y = Iy + def writeUnits(cremodel_): unitVol = cremodel_.createUnitDefinition() unitVol.setId("volume") @@ -905,6 +1055,30 @@ def writeUnits(cremodel_): unit.setExponent(1.0) unit.setScale(-3) + unitLen = cremodel_.createUnitDefinition() + unitLen.setId("length") + unit = unitLen.createUnit() + unit.setKind(UNIT_KIND_METRE) + unit.setMultiplier(1.0) + unit.setExponent(1.0) + unit.setScale(0) + + unitArea = cremodel_.createUnitDefinition() + unitArea.setId("area") + unit = unitArea.createUnit() + unit.setKind(UNIT_KIND_METRE) + unit.setMultiplier(1.0) + unit.setExponent(2.0) + unit.setScale(0) + + unitTime = cremodel_.createUnitDefinition() + unitTime.setId("time") + unit = unitTime.createUnit() + unit.setKind(UNIT_KIND_SECOND) + unit.setExponent(1) + unit.setMultiplier(1) + unit.setScale(0) + if __name__ == "__main__": try: sys.argv[1] diff --git a/moose-core/python/moose/chemMerge/__init__.py b/moose-core/python/moose/chemMerge/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..0116e5dbe066dd86f67a71deecb70994c4a7cd1f --- /dev/null +++ b/moose-core/python/moose/chemMerge/__init__.py @@ -0,0 +1,2 @@ +# -*- coding: utf-8 -*- +from .merge import * diff --git a/moose-core/python/moose/merge/merge.py b/moose-core/python/moose/chemMerge/merge.py similarity index 55% rename from moose-core/python/moose/merge/merge.py rename to moose-core/python/moose/chemMerge/merge.py index cb1ab62d7ecbeb4c1b8b648100fa5dc07d8a4cb0..5b406a25a99ef8428987ac88ab68cb8b2e9a3368 100644 --- a/moose-core/python/moose/merge/merge.py +++ b/moose-core/python/moose/chemMerge/merge.py @@ -12,185 +12,288 @@ #** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS #Created : Friday Dec 16 23:19:00 2016(+0530) #Version -#Last-Updated: Thursday Jan 12 17:30:33 2017(+0530) +#Last-Updated: Wed Oct 14 00:25:33 2017(+0530) # By: Harsha #**********************************************************************/ -# This program is used to merge models -# -- Model B is merged to modelA -#Rules are -#--In this models are mergered at group level (if exists) - - +# This program is used to merge chem models from src to destination +#Rules are : +# -- If Compartment from the src model doesn't exist in destination model, +# then entire compartment and its children are copied over including groups +# -- Models are mergered at group level (if exists) +# (Group is Neutral object in moose, which may represent pathway in network model) +# -- Pool's are copied from source to destination if it doesn't exist, if exist nothing is done +# -- Reaction (Reac), Enzyme (Enz) are copied +# --- if any danglling Reac or Enz exist then that is not copied +# +# --- if Reac Name's is different for a given path (group level) +# then copy the entire Reac along with substrate/product +# --- if same Reac Name and same sub and prd then nothing is copied +# --- if same Reac Name but sub or prd is different then duplicated and copied +# +# --- if Enz Name's is different for a given parent pool path +# then copy the entire Enz along with substrate/product +# --- if same Enz Name and same sub and prd then nothing is copied +# --- if same Enz Name but sub or prd is different then duplicated and copied +# -- Function are copied only if destination pool to which its suppose to connect doesn't exist with function of its own +# +''' +Change log +Oct 14: absolute import with mtypes just for python3 + +Oct 12: clean way of cheking the type of path provided, filepath,moose obj, moose path are taken, + if source is empty then nothing to copy, + if destination was empty list is update with new object + +Oct 11: missing group are copied instead of creating one in new path which also copies Annotator info + earlier if user asked to save the model, it was saving default to kkit format, now user need to run the command to save (if this is run in command) + To check: When Gui is allowed to merge 2 models, need to see what happens + +''' import sys import os #from . import _moose as moose import moose -import mtypes +from moose.chemMerge import mtypes from moose.chemUtil.chemConnectUtil import * from moose.chemUtil.graphUtils import * +#from moose.genesis import mooseWriteKkit -def mergeChemModel(A,B): - """ Merges two model or the path """ +def checkFile_Obj_str(file_Obj_str): + model = moose.element('/') + loaded = False + found = False + if isinstance(file_Obj_str, str): + if os.path.isfile(file_Obj_str) == True: + model,loaded = loadModels(file_Obj_str) + found = True + elif file_Obj_str.find('/') != -1 : + if not isinstance(file_Obj_str,moose.Neutral): + if moose.exists(file_Obj_str): + model = file_Obj_str + loaded = True + found = True + elif isinstance(file_Obj_str, moose.Neutral): + if moose.exists(file_Obj_str.path): + model = file_Obj_str.path + loaded = True + found = True + elif isinstance(file_Obj_str, moose.Neutral): + if moose.exists(file_Obj_str.path): + model = file_Obj_str.path + loaded = True + found = True + if not found: + print ("%s path or filename doesnot exist. " % (file_Obj_str)) + return model,loaded - modelA,loadedA = loadModels(A) - modelB,loadedB = loadModels(B) - if not loadedA or not loadedB: - if not loadedA: - modelB = moose.Shell('/') - if not loadedB: - modelA = moose.Shell('/') - else: - directory, bfname = os.path.split(B) - global grpNotcopiedyet,poolListina +def mergeChemModel(src,des): + """ Merges two model or the path """ + A = src + B = des + sfile = src + dfile = des + loadedA = False + loadedB = False + modelA = moose.element('/') + modelB = moose.element('/') + modelA,loadedA = checkFile_Obj_str(A) + modelB,loadedB = checkFile_Obj_str(B) + + if loadedA and loadedB: + ## yet deleteSolver is called to make sure all the moose object are off from solver + deleteSolver(modelA) + deleteSolver(modelB) + + global poolListina poolListina = {} grpNotcopiedyet = [] dictComptA = dict( [ (i.name,i) for i in moose.wildcardFind(modelA+'/##[ISA=ChemCompt]') ] ) dictComptB = dict( [ (i.name,i) for i in moose.wildcardFind(modelB+'/##[ISA=ChemCompt]') ] ) - poolNotcopiedyet = [] - - for key in list(dictComptB.keys()): - if key not in dictComptA: - # if compartmentname from modelB does not exist in modelA, then copy - copy = moose.copy(dictComptB[key],moose.element(modelA)) + if len(dictComptA): + for key in list(dictComptA.keys()): + if key not in dictComptB: + # if compartmentname from modelB does not exist in modelA, then copy + copy = moose.copy(dictComptA[key],moose.element(modelB)) + dictComptB[key]=moose.element(copy) + else: + #if compartmentname from modelB exist in modelA, + #volume is not same, then change volume of ModelB same as ModelA + if abs(dictComptB[key].volume - dictComptA[key].volume): + #hack for now + while (abs(dictComptB[key].volume - dictComptA[key].volume) != 0.0): + dictComptA[key].volume = float(dictComptB[key].volume) + dictComptB = dict( [ (i.name,i) for i in moose.wildcardFind(modelB+'/##[ISA=ChemCompt]') ] ) + + #Mergering pool + poolMerge(dictComptB[key],dictComptA[key],poolNotcopiedyet) + + if grpNotcopiedyet: + # objA = moose.element(comptA).parent.name + # if not moose.exists(objA+'/'+comptB.name+'/'+bpath.name): + # print bpath + # moose.copy(bpath,moose.element(objA+'/'+comptB.name)) + pass + comptBdict = comptList(modelB) + poolListinb = {} + poolListinb = updatePoolList(comptBdict) + R_Duplicated, R_Notcopiedyet,R_Daggling = [], [], [] + E_Duplicated, E_Notcopiedyet,E_Daggling = [], [], [] + for key in list(dictComptA.keys()): + funcExist, funcNotallowed = [], [] + funcExist,funcNotallowed = functionMerge(dictComptB,dictComptA,key) + + poolListinb = updatePoolList(dictComptB) + R_Duplicated,R_Notcopiedyet,R_Daggling = reacMerge(dictComptB,dictComptA,key,poolListinb) + + poolListinb = updatePoolList(dictComptB) + E_Duplicated,E_Notcopiedyet,E_Daggling = enzymeMerge(dictComptB,dictComptA,key,poolListinb) + ''' + if isinstance(src, str): + if os.path.isfile(src) == True: + spath, sfile = os.path.split(src) + else: + sfile = src else: - #if compartmentname from modelB exist in modelA, - #volume is not same, then change volume of ModelB same as ModelA - if abs(dictComptA[key].volume - dictComptB[key].volume): - #hack for now - while (abs(dictComptA[key].volume - dictComptB[key].volume) != 0.0): - dictComptB[key].volume = float(dictComptA[key].volume) - dictComptA = dict( [ (i.name,i) for i in moose.wildcardFind(modelA+'/##[ISA=ChemCompt]') ] ) - - #Mergering pool - poolMerge(dictComptA[key],dictComptB[key],poolNotcopiedyet) - - if grpNotcopiedyet: - # objA = moose.element(comptA).parent.name - # if not moose.exists(objA+'/'+comptB.name+'/'+bpath.name): - # print bpath - # moose.copy(bpath,moose.element(objA+'/'+comptB.name)) - pass - - comptAdict = comptList(modelA) - poolListina = {} - poolListina = updatePoolList(comptAdict) - funcNotallowed = [] - R_Duplicated, R_Notcopiedyet,R_Daggling = [], [], [] - E_Duplicated, E_Notcopiedyet,E_Daggling = [], [], [] - for key in list(dictComptB.keys()): - funcNotallowed = [] - funcNotallowed = functionMerge(dictComptA,dictComptB,key) - - poolListina = updatePoolList(dictComptA) - R_Duplicated,R_Notcopiedyet,R_Daggling = reacMerge(dictComptA,dictComptB,key,poolListina) - - poolListina = updatePoolList(dictComptA) - E_Duplicated,E_Notcopiedyet,E_Daggling = enzymeMerge(dictComptA,dictComptB,key,poolListina) - - print("\n Model is merged to %s" %modelA) - - if funcNotallowed: - print( "\nPool already connected to a function, this function is not to connect to same pool, since no two function are allowed to connect to same pool:") - for fl in list(funcNotallowed): - print("\t [Pool]: %s [Function]: %s \n" %(str(fl.parent.name), str(fl.path))) - - if R_Duplicated or E_Duplicated: - print ("Reaction / Enzyme are Duplicate" - "\n 1. The once whoes substrate / product names are different for a give reaction name " - "\n 2. its compartment to which it belongs to may be is different" - "\n Models have to decide to keep or delete these reaction/enzyme") - if E_Duplicated: - print("Reaction: ") - for rd in list(R_Duplicated): - print ("%s " %str(rd.name)) - - if E_Duplicated: - print ("Enzyme:") - for ed in list(E_Duplicated): - print ("%s " %str(ed.name)) - - if R_Notcopiedyet or E_Notcopiedyet: - - print ("\nThese Reaction/Enzyme in model are not dagging but while copying the associated substrate or product is missing") - if R_Notcopiedyet: - print("Reaction: ") - for rd in list(R_Notcopiedyet): - print ("%s " %str(rd.name)) - if E_Notcopiedyet: - print ("Enzyme:") - for ed in list(E_Notcopiedyet): - print ("%s " %str(ed.name)) - - if R_Daggling or E_Daggling: - print ("\n Daggling reaction/enzyme are not not allowed in moose, these are not merged") - if R_Daggling: - print("Reaction: ") - for rd in list(R_Daggling): + sfile = src + if isinstance(des, str): + print " A str",des + if os.path.isfile(des) == True: + dpath, dfile = os.path.split(des) + else: + dfile = des + else: + dfile = des + ''' + print("\nThe content of %s (src) model is merged to %s (des)." %(sfile, dfile)) + # Here any error or warning during Merge is written it down + if funcExist: + print( "\nIn model \"%s\" pool already has connection from a function, these function from model \"%s\" is not allowed to connect to same pool,\n since no two function are allowed to connect to same pool:"%(dfile, sfile)) + for fl in list(funcExist): + print("\t [Pool]: %s [Function]: %s \n" %(str(fl.parent.name), str(fl.path))) + if funcNotallowed: + print( "\nThese functions is not to copied, since pool connected to function input are from different compartment:") + for fl in list(funcNotallowed): + print("\t [Pool]: %s [Function]: %s \n" %(str(fl.parent.name), str(fl.path))) + if R_Duplicated or E_Duplicated: + print ("These Reaction / Enzyme are \"Duplicated\" into destination file \"%s\", due to " + "\n 1. If substrate / product name's are different for a give reaction/Enzyme name " + "\n 2. If product belongs to different compartment " + "\n Models have to decide to keep or delete these reaction/enzyme in %s" %(dfile, dfile)) + if E_Duplicated: + print("Reaction: ") + for rd in list(R_Duplicated): print ("%s " %str(rd.name)) - if E_Daggling: - print ("Enzyme:") - for ed in list(E_Daggling): - print ("%s " %str(ed.name)) + if E_Duplicated: + print ("Enzyme:") + for ed in list(E_Duplicated): + print ("%s " %str(ed.name)) + + if R_Notcopiedyet or E_Notcopiedyet: + + print ("\nThese Reaction/Enzyme in model are not dagging but while copying the associated substrate or product is missing") + if R_Notcopiedyet: + print("Reaction: ") + for rd in list(R_Notcopiedyet): + print ("%s " %str(rd.name)) + if E_Notcopiedyet: + print ("Enzyme:") + for ed in list(E_Notcopiedyet): + print ("%s " %str(ed.name)) + + if R_Daggling or E_Daggling: + print ("\n Daggling reaction/enzyme are not allowed in moose, these are not merged to %s from %s" %(dfile, sfile)) + if R_Daggling: + print("Reaction: ") + for rd in list(R_Daggling): + print ("%s " %str(rd.name)) + if E_Daggling: + print ("Enzyme:") + for ed in list(E_Daggling): + print ("%s " %str(ed.name)) + ## Model is saved + print ("\n ") + print ('\nMerged model is available under moose.element(\'%s\')' %(modelB)) + print (' From the python terminal itself \n to save the model in to genesis format use \n >moose.mooseWriteKkit(\'%s\',\'filename.g\')' %(modelB)) + print (' to save into SBML format \n >moose.mooseWriteSBML(\'%s\',\'filename.xml\')' %(modelB)) + return modelB + # savemodel = raw_input("Do you want to save the model? \"YES\" \"NO\" ") + # if savemodel.lower() == 'yes' or savemodel.lower() == 'y': + # mergeto = raw_input("Enter File name ") + # if mergeto and mergeto.strip(): + # filenameto = 'merge.g' + # else: + # if str(mergeto).rfind('.') != -1: + # mergeto = mergeto[:str(mergeto).rfind('.')] + # if str(mergeto).rfind('/'): + # mergeto = mergeto+'merge' + + # filenameto = mergeto+'.g' + + # error,written = moose.mooseWriteKkit(modelB, filenameto) + # if written == False: + # print('Could not save the Model, check the files') + # else: + # if error == "": + # print(" \n The merged model is saved into \'%s\' " %(filenameto)) + # else: + # print('Model is saved but these are not written\n %s' %(error)) + + # else: + # print ('\nMerged model is available under moose.element(\'%s\')' %(modelB)) + # print (' If you are in python terminal you could save \n >moose.mooseWriteKkit(\'%s\',\'filename.g\')' %(modelB)) + # print (' If you are in python terminal you could save \n >moose.mooseWriteSBML(\'%s\',\'filename.g\')' %(modelB)) + #return modelB + else: + print ('\nSource file has no objects to copy(\'%s\')' %(modelA)) + return moose.element('/') def functionMerge(comptA,comptB,key): - funcNotallowed = [] + funcNotallowed, funcExist = [], [] comptApath = moose.element(comptA[key]).path comptBpath = moose.element(comptB[key]).path - funcListina = moose.wildcardFind(comptApath+'/##[ISA=PoolBase]') - funcListinb = moose.wildcardFind(comptBpath+'/##[ISA=Function]') objA = moose.element(comptApath).parent.name objB = moose.element(comptBpath).parent.name - #For function mergering pool name is taken as reference - funcNotcopiedyet = [] - for fb in funcListinb: - if fb.parent.className in ['ZombiePool','Pool','ZombieBufPool','BufPool']: - objA = moose.element(comptApath).parent.name - fbpath = fb.path - #funcpath = fbpath[fbpath.find(findCompartment(fb).name)-1:len(fbpath)] - funcparentB = fb.parent.path - funcpath = fbpath.replace(objB,objA) - funcparentA = funcparentB.replace(objB,objA) - tt = moose.element(funcparentA).neighbors['setN'] - if tt: - funcNotallowed.append(fb) - else: - if len(moose.element(fb.path+'/x').neighbors["input"]): - #inputB = moose.element(fb.path+'/x').neighbors["input"] - inputB = subprdList(moose.element(fb.path+'/x'),"input") - inputB_expr = fb.expr - if moose.exists(funcpath): - #inputA = moose.element(objA+funcpath+'/x').neighbors["input"] - inputA = subprdList(moose.element(funcpath+'/x'),"input") - inputA_expr = moose.element(funcpath).expr - hassameExpr = False - if inputA_expr == inputB_expr: - hassameExpr = True - hassameLen,hassameS,hassameVols = same_len_name_vol(inputA,inputB) - if not all((hassameLen,hassameS,hassameVols,hassameExpr)): - fb.name = fb.name+'_duplicatedF' - createFunction(fb,inputB,objB,objA) + #This will give us all the function which exist in modelB + funcListinb = moose.wildcardFind(comptBpath+'/##[ISA=Function]') + for fb in funcListinb: + #This will give us all the pools that its connected to, for this function + fvalueOut = moose.element(fb).neighbors['valueOut'] + for poolinB in fvalueOut: + poolinBpath = poolinB.path + poolinA = poolinBpath.replace(objB,objA) + connectionexist = [] + if moose.exists(poolinA): + #This is give us if pool which is to be connected already exist any connection + connectionexist = moose.element(poolinA).neighbors['setN']+moose.element(poolinA).neighbors['setConc']+ moose.element(poolinA).neighbors['increment'] + if len(connectionexist) == 0: + #This pool in model A doesnot exist with any function + inputs = moose.element(fb.path+'/x').neighbors['input'] + volumes = [] + for ins in inputs: + volumes.append((findCompartment(moose.element(ins))).volume) + if len(set(volumes)) == 1: + # If all the input connected belongs to one compartment then copy + createFunction(fb,poolinA,objB,objA) else: - #function doesnot exist then copy - if len(inputB): - volinput = [] - for inb in inputB: - volinput.append(findCompartment(moose.element(inb)).volume) - if len(set(volinput)) == 1: - # If all the input connected belongs to one compartment then copy - createFunction(fb,inputB,objB,objA) - else: - # moose doesn't allow function input to come from different compartment - funcNotallowed.append(fb) - return funcNotallowed - -def createFunction(fb,inputB,objB,objA): + # moose doesn't allow function's input to come from different compartment + funcNotallowed.append(fb) + else: + #Pool in model 'A' already exist function " + funcExist.append(fb) + else: + print(" Path in model A doesn't exists %s" %(poolinA)) + + return funcExist,funcNotallowed + +def createFunction(fb,setpool,objB,objA): fapath1 = fb.path.replace(objB,objA) fapath = fapath1.replace('[0]','') + if not moose.exists(fapath): # if fb.parent.className in ['CubeMesh','CyclMesh']: # des = moose.Function('/'+objA+'/'+fb.parent.name+'/'+fb.name) @@ -199,7 +302,13 @@ def createFunction(fb,inputB,objB,objA): # if fb.parent.name == akey.name: # des = moose.Function(akey.path+'/'+fb.name) des = moose.Function(fapath) - moose.connect(des, 'valueOut', moose.element(fapath).parent,'setN' ) + else: + des = moose.element(fapath) + inputB = moose.element(fb.path+'/x').neighbors["input"] + moose.connect(des, 'valueOut', moose.element(setpool),'setN' ) + inputA = [] + inputA = moose.element(fapath+'/x').neighbors["input"] + if not inputA: for src in inputB: pool = ((src.path).replace(objB,objA)).replace('[0]','') numVariables = des.numVars @@ -210,46 +319,47 @@ def createFunction(fb,inputB,objB,objA): des.expr = expr moose.connect( pool, 'nOut', des.x[numVariables], 'input' ) - #if fb.expr != des.expr: - # print "Function ",des, " which is duplicated from modelB, expression is different, this is tricky in moose to know what those constants are connected to " - # print "ModelB ", fb, fb.expr, "\nModelA ",des, des.expr - def comptList(modelpath): comptdict = {} for ca in moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]'): comptdict[ca.name] = ca return comptdict -def loadModels(filename): +def loadModels(filepath): """ load models into moose if file, if moosepath itself it passes back the path and delete solver if exist """ modelpath = '/' loaded = False - if os.path.isfile(filename) : - modelpath = filename[filename.rfind('/'): filename.rfind('.')] - ext = os.path.splitext(filename)[1] - filename = filename.strip() - modeltype = mtypes.getType(filename) - subtype = mtypes.getSubtype(filename, modeltype) + if os.path.isfile(filepath) : + fpath, filename = os.path.split(filepath) + # print " head and tail ",head, " ",tail + # modelpath = filename[filename.rfind('/'): filename.rfind('.')] + # print "modelpath ",modelpath + # ext = os.path.splitext(filename)[1] + # filename = filename.strip() + modelpath = '/'+filename[:filename.rfind('.')] + modeltype = mtypes.getType(filepath) + subtype = mtypes.getSubtype(filepath, modeltype) + if subtype == 'kkit' or modeltype == "cspace": - moose.loadModel(filename,modelpath) + if moose.exists(modelpath): + moose.delete(modelpath) + moose.loadModel(filepath,modelpath) loaded = True elif subtype == 'sbml': - #moose.ReadSBML() + #moose.mooseReadSBML(filename,modelpath) + #loaded = True pass else: print("This file is not supported for mergering") modelpath = moose.Shell('/') - elif moose.exists(filename): - modelpath = filename + + elif moose.exists(filepath): + modelpath = filepath loaded = True - ## default is 'ee' solver while loading the model using moose.loadModel, - ## yet deleteSolver is called just to be assured - if loaded: - deleteSolver(modelpath) return modelpath,loaded @@ -281,7 +391,11 @@ def poolMerge(comptA,comptB,poolNotcopiedyet): bpath.name = bpath.name+"_grp" l = moose.Neutral(grp_cmpt) else: - moose.Neutral(grp_cmpt) + #moose.Neutral(grp_cmpt) + src = bpath + srcpath = (bpath.parent).path + des = srcpath.replace(objB,objA) + moose.copy(bpath,moose.element(des)) apath = moose.element(bpath.path.replace(objB,objA)) @@ -301,15 +415,15 @@ def copy_deleteUnlyingPoolObj(pool,path): # which will automatically copie's the pool copied = False - if pool.parent.className not in ["Enz","ZombieEnz"]: + if pool.parent.className not in ["Enz","ZombieEnz","MMenz","ZombieMMenz"]: poolcopied = moose.copy(pool,path) copied = True # deleting function and enzyme which gets copied if exist under pool # This is done to ensure daggling function / enzyme not copied. funclist = [] for types in ['setConc','setN','increment']: - funclist.extend(moose.element(poolcopied).neighbors[types]) + for fl in funclist: moose.delete(fl) enzlist = moose.element(poolcopied).neighbors['reac'] @@ -365,20 +479,15 @@ def enzymeMerge(comptA,comptB,key,poolListina): allclean = True else: # didn't find sub or prd for this Enzyme - #print ("didn't find sub or prd for this reaction" ) RE_Notcopiedyet.append(eb) else: # -- it is dagging reaction RE_Daggling.append(eb) - #print ("This reaction \""+eb.path+"\" has no substrate/product daggling reaction are not copied") - #war_msg = war_msg+"\nThis reaction \""+eb.path+"\" has no substrate/product daggling reaction are not copied" else: #Same Enzyme name # -- Same substrate and product including same volume then don't copy # -- different substrate/product or if sub/prd's volume is different then DUPLICATE the Enzyme allclean = False - #ea = moose.element('/'+obj+'/'+enzcompartment.name+'/'+enzparent.name+'/'+eb.name) - #ea = moose.element(pA.path+'/'+eb.name) ea = moose.element(eb.path.replace(objB,objA)) eAsubname = subprdList(ea,"sub") eBsubname = subprdList(eb,"sub") @@ -408,7 +517,6 @@ def enzymeMerge(comptA,comptB,key,poolListina): moose.connect(moose.element(enz).parent,"nOut",moose.element(enz),"enzDest") #moose.connect(moose.element(enz),"enz",moose.element(enz).parent,"reac") - #moose.connect( cplxItem, 'reac', enz, 'cplx' ) connectObj(enz,eBsubname,"sub",comptA,war_msg) connectObj(enz,eBprdname,"prd",comptA,war_msg) RE_Duplicated.append(enz) @@ -423,12 +531,9 @@ def enzymeMerge(comptA,comptB,key,poolListina): # -- it may be connected Enzyme cplx if eBsubname and eBprdname: RE_Notcopiedyet.append(eb) - #print ("This Enzyme \""+eb.path+"\" has no substrate/product must be connect to cplx") - #war_msg = war_msg+ "\nThis Enzyme \""+rb.path+"\" has no substrate/product must be connect to cplx" else: RE_Daggling.append(eb) - #print ("This enzyme \""+eb.path+"\" has no substrate/product daggling reaction are not copied") - #war_msg = war_msg+"\nThis reaction \""+eb.path+"\" has no substrate/product daggling reaction are not copied" + return RE_Duplicated,RE_Notcopiedyet,RE_Daggling def reacMerge(comptA,comptB,key,poolListina): @@ -512,17 +617,12 @@ def reacMerge(comptA,comptB,key,poolListina): # -- it may be connected Enzyme cplx if rBsubname and rBprdname: RE_Notcopiedyet.append(rb) - #print ("This reaction \""+rb.path+"\" has no substrate/product must be connect to cplx") - #war_msg = war_msg+ "\nThis reaction \""+rb.path+"\" has no substrate/product must be connect to cplx" else: RE_Daggling.append(rb) - #print ("This reaction \""+rb.path+"\" has no substrate/product daggling reaction are not copied") - #war_msg = war_msg+"\nThis reaction \""+rb.path+"\" has no substrate/product daggling reaction are not copied" return RE_Duplicated,RE_Notcopiedyet,RE_Daggling def subprdList(reac,subprd): - #print "Reac ",reac rtype = moose.element(reac).neighbors[subprd] rname = [] for rs in rtype: @@ -573,8 +673,6 @@ def connectObj(reac,spList,spType,comptA,war_msg): allclean = True else: #It should not come here unless the sub/prd is connected to enzyme cplx pool - #print ("This pool \""+rsp.name+"\" doesnot exists in this "+comptName+" compartment to connect to this reaction \""+reac.name+"\"") - #war_msg = war_msg+ "This pool \""+rsp.name+"\" doesnot exists in this "+comptName+" compartment to connect to this reaction \""+reac.name+"\"" allclean = False return allclean @@ -604,7 +702,40 @@ def mooseIsInstance(element, classNames): if __name__ == "__main__": - - modelA = '/home/harsha/genesis_files/gfile/acc92.g' - modelB = '/home/harsha/genesis_files/gfile/acc50.g' - mergered = mergeChemModel(modelA,modelB) + try: + sys.argv[1] + except IndexError: + print("Source filename or path not given") + exit(0) + else: + src = sys.argv[1] + if not os.path.exists(src): + print("Filename or path does not exist",src) + else: + try: + sys.argv[2] + except IndexError: + print("Destination filename or path not given") + exit(0) + else: + des = sys.argv[2] + if not os.path.exists(src): + print("Filename or path does not exist",des) + exit(0) + else: + mergered = mergeChemModel(src,des) + ''' + try: + sys.argv[3] + except IndexError: + print ("Merge to save not specified") + mergeto = "merge" + else: + mergeto = sys.argv[3] + if str(mergeto).rfind('.') != -1: + mergeto = mergeto[:str(mergeto).rfind('.')] + if str(mergeto).rfind('/'): + mergeto = mergeto+'merge' + + mergered = mergeChemModel(src,des,mergeto) + ''' diff --git a/moose-core/python/moose/merge/mtypes.py b/moose-core/python/moose/chemMerge/mtypes.py similarity index 100% rename from moose-core/python/moose/merge/mtypes.py rename to moose-core/python/moose/chemMerge/mtypes.py diff --git a/moose-core/python/moose/chemUtil/chemConnectUtil.py b/moose-core/python/moose/chemUtil/chemConnectUtil.py index be0c4ba6f36fd8d60aaaa3c2a75864f9601cdd2f..48e06e9482e16c922f9124746ca189af27bbcfc0 100644 --- a/moose-core/python/moose/chemUtil/chemConnectUtil.py +++ b/moose-core/python/moose/chemUtil/chemConnectUtil.py @@ -1,7 +1,22 @@ # -*- coding: utf-8 -*- +""" ChemconnectUtil.py Some of the command function are written """ + +#Created : Friday May 27 12:19:00 2016(+0530) + +__author__ = "Harsha Rani" +__copyright__ = "Copyright 2017, Harsha Rani and NCBS Bangalore" +__credits__ = ["NCBS Bangalore"] +__license__ = "GNU GPL" +__version__ = "1.0.0" +__maintainer__ = "Harsha Rani" +__email__ = "hrani@ncbs.res.in" +__status__ = "Development" +__updated__ = "Aug 8 2017" + +#Aug 8 : added findCompartment function here + import moose import numpy as np -from collections import Counter def xyPosition(objInfo,xory): try: @@ -119,7 +134,7 @@ def setupItem(modelPath,cntDict): uniqItem,countuniqItem = countitems(items,'prd') prdNo = uniqItem if (len(subNo) == 0 or len(prdNo) == 0): - print ("Substrate Product is empty ",path, " ",items) + print ("\nSubstrate Product is empty for "+items.path) for prd in uniqItem: prdlist.append((moose.element(prd),'p',countuniqItem[prd])) @@ -178,5 +193,18 @@ def countitems(mitems,objtype): items = [] items = moose.element(mitems).neighbors[objtype] uniqItems = set(items) - countuniqItems = Counter(items) + #countuniqItems = Counter(items) + countuniqItems = dict((i, items.count(i)) for i in items) return(uniqItems,countuniqItems) + +def findCompartment(element): + if element.path == '/': + return moose.element('/') + elif mooseIsInstance(element, ["CubeMesh", "CyclMesh"]): + return (element) + else: + return findCompartment(moose.element(element.parent)) + + +def mooseIsInstance(element, classNames): + return moose.element(element).__class__.__name__ in classNames \ No newline at end of file diff --git a/moose-core/python/moose/chemUtil/graphUtils.py b/moose-core/python/moose/chemUtil/graphUtils.py index 061b8ea0111ae46ec6b5b22ce80b72352c0f5af4..2ca109f30110f866dd80f5381e7ce3f09e8fea7b 100644 --- a/moose-core/python/moose/chemUtil/graphUtils.py +++ b/moose-core/python/moose/chemUtil/graphUtils.py @@ -1,4 +1,26 @@ # -*- coding: utf-8 -*- +''' +******************************************************************* + * File: chemConnectUtil.py + * Description: + * Author: HarshaRani + * E-mail: hrani@ncbs.res.in + ********************************************************************/ +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment, +** also known as GENESIS 3 base code. +** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS +Created : Friday May 27 12:19:00 2016(+0530) +Version +Last-Updated: Mon Aug 7 17:02:00 2017(+0530) + By: HarshaRani +**********************************************************************/ +/**************************** + +Aug 7: cleaned up space +''' + import moose pygraphvizFound_ = True @@ -47,7 +69,7 @@ def autoCoordinates(meshEntry,srcdesConnection): for items in (items for items in out[0] ): G.add_edge(moose.element(items[0]).path,inn.path) if len(out[1]) == 0: - print(inn.className + ':' + inn.name + "doesn't have output mssg") + print(inn.className + ':' + inn.name + " doesn't have output mssg") else: for items in (items for items in out[1] ): G.add_edge(inn.path,moose.element(items[0]).path) diff --git a/moose-core/python/moose/genesis/writeKkit.py b/moose-core/python/moose/genesis/writeKkit.py index 2e25bc7ed974242abfe42a9d978315a57a39d973..c94ca3daddfc1e0a0215ccb04e84986c2c0fad8a 100644 --- a/moose-core/python/moose/genesis/writeKkit.py +++ b/moose-core/python/moose/genesis/writeKkit.py @@ -9,16 +9,25 @@ __version__ = "1.0.0" __maintainer__ = "Harsha Rani" __email__ = "hrani@ncbs.res.in" __status__ = "Development" +__updated__ = "Aug 8 2017" + +# 2017 +# Aug 8 : All the moose object which doesn't have compartment are not written. Error message are appended import sys import random import re -import matplotlib import moose - from moose.chemUtil.chemConnectUtil import * from moose.chemUtil.graphUtils import * +foundmatplotlib_ = False +try: + import matplotlib + foundmatplotlib_ = True +except Exception as e: + pass + GENESIS_COLOR_SEQUENCE = ((248, 0, 255), (240, 0, 255), (232, 0, 255), (224, 0, 255), (216, 0, 255), (208, 0, 255), (200, 0, 255), (192, 0, 255), (184, 0, 255), (176, 0, 255), (168, 0, 255), (160, 0, 255), (152, 0, 255), (144, 0, 255), (136, 0, 255), (128, 0, 255), (120, 0, 255), (112, 0, 255), (104, 0, 255), (96, 0, 255), (88, 0, 255), (80, 0, 255), @@ -40,78 +49,100 @@ GENESIS_COLOR_SEQUENCE = ((248, 0, 255), (240, 0, 255), (232, 0, 255), (224, 0, # --StimulusTable def mooseWriteKkit( modelpath, filename, sceneitems={}): - if filename.rfind('.') != -1: - filename = filename[:filename.rfind('.')] + global foundmatplotlib_ + if not foundmatplotlib_: + print('No maplotlib found.' + '\nThis module can be installed by following command in terminal:' + '\n\t sudo apt install python-maplotlib', "") + return False else: - filename = filename[:len(filename)] - filename = filename+'.g' - global NA - NA = 6.0221415e23 - global cmin,cmax,xmin,xmax,ymin,ymax - cmin, xmin, ymin = 0, 0, 0 - cmax, xmax, ymax = 1, 1, 1 - - compt = moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]') - maxVol = estimateDefaultVol(compt) - positionInfoExist = True - if compt: - if bool(sceneitems): - cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems) - elif not bool(sceneitems): - srcdesConnection = {} - setupItem(modelpath,srcdesConnection) - meshEntry,xmin,xmax,ymin,ymax,positionInfoExist,sceneitems = setupMeshObj(modelpath) + error = "" + + if filename.rfind('.') != -1: + filename = filename[:filename.rfind('.')] + else: + filename = filename[:len(filename)] + filename = filename+'.g' + global NA + NA = 6.0221415e23 + global cmin,cmax,xmin,xmax,ymin,ymax + cmin, xmin, ymin = 0, 0, 0 + cmax, xmax, ymax = 1, 1, 1 + + compt = moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]') + maxVol = estimateDefaultVol(compt) + positionInfoExist = True + if compt: + if bool(sceneitems): + cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems) + elif not bool(sceneitems): + srcdesConnection = {} + setupItem(modelpath,srcdesConnection) + meshEntry,xmin,xmax,ymin,ymax,positionInfoExist,sceneitems = setupMeshObj(modelpath) + if not positionInfoExist: + #cmin,cmax,sceneitems = autoCoordinates(meshEntry,srcdesConnection) + sceneitems = autoCoordinates(meshEntry,srcdesConnection) + if not positionInfoExist: - #cmin,cmax,sceneitems = autoCoordinates(meshEntry,srcdesConnection) - sceneitems = autoCoordinates(meshEntry,srcdesConnection) - - if not positionInfoExist: - # if position are not from kkit, then zoom factor is applied while - # writing to genesis. Like if position is from pyqtSceneItem or auto-coordinates - cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems) - for k,v in list(sceneitems.items()): - anno = moose.element(k.path+'/info') - x1 = calPrime(v['x']) - y1 = calPrime(v['y']) - sceneitems[k]['x'] = x1 - sceneitems[k]['y'] = y1 - - f = open(filename, 'w') - writeHeader (f,maxVol) - - gtId_vol = writeCompartment(modelpath,compt,f) - writePool(modelpath,f,gtId_vol,sceneitems) - reacList = writeReac(modelpath,f,sceneitems) - enzList = writeEnz(modelpath,f,sceneitems) - writeSumtotal(modelpath,f) - f.write("simundump xgraph /graphs/conc1 0 0 99 0.001 0.999 0\n" - "simundump xgraph /graphs/conc2 0 0 100 0 1 0\n") - tgraphs = moose.wildcardFind(modelpath+'/##[ISA=Table2]') - first, second = " ", " " - if tgraphs: - first,second = writeplot(tgraphs,f) - if first: - f.write(first) - f.write("simundump xgraph /moregraphs/conc3 0 0 100 0 1 0\n" - "simundump xgraph /moregraphs/conc4 0 0 100 0 1 0\n") - if second: - f.write(second) - f.write("simundump xcoredraw /edit/draw 0 -6 4 -2 6\n" - "simundump xtree /edit/draw/tree 0 \\\n" - " /kinetics/#[],/kinetics/#[]/#[],/kinetics/#[]/#[]/#[][TYPE!=proto],/kinetics/#[]/#[]/#[][TYPE!=linkinfo]/##[] \"edit_elm.D <v>; drag_from_edit.w <d> <S> <x> <y> <z>\" auto 0.6\n" - "simundump xtext /file/notes 0 1\n") - storeReacMsg(reacList,f) - storeEnzMsg(enzList,f) - if tgraphs: - storePlotMsgs(tgraphs,f) - writeFooter1(f) - writeNotes(modelpath,f) - writeFooter2(f) - print('Written to file '+filename) - return True - else: - print(("Warning: writeKkit:: No model found on " , modelpath)) - return False + # if position are not from kkit, then zoom factor is applied while + # writing to genesis. Like if position is from pyqtSceneItem or auto-coordinates + cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems) + for k,v in list(sceneitems.items()): + anno = moose.element(k.path+'/info') + #x1 = calPrime(v['x']) + #y1 = calPrime(v['y']) + #sceneitems[k]['x'] = x1 + #sceneitems[k]['y'] = y1 + + f = open(filename, 'w') + writeHeader (f,maxVol) + + gtId_vol = writeCompartment(modelpath,compt,f) + errors = "" + error = writePool(modelpath,f,gtId_vol,sceneitems) + errors = errors+error + + reacList,error= writeReac(modelpath,f,sceneitems) + errors = errors+error + + enzList,error = writeEnz(modelpath,f,sceneitems) + errors = errors+error + + error = writeSumtotal(modelpath,f) + errors = errors+error + + error = writeStimulus(modelpath,f) + errors = errors+error + + #writeSumtotal(modelpath,f) + f.write("simundump xgraph /graphs/conc1 0 0 99 0.001 0.999 0\n" + "simundump xgraph /graphs/conc2 0 0 100 0 1 0\n") + tgraphs = moose.wildcardFind(modelpath+'/##[ISA=Table2]') + first, second = " ", " " + if tgraphs: + first,second = writeplot(tgraphs,f) + if first: + f.write(first) + f.write("simundump xgraph /moregraphs/conc3 0 0 100 0 1 0\n" + "simundump xgraph /moregraphs/conc4 0 0 100 0 1 0\n") + if second: + f.write(second) + f.write("simundump xcoredraw /edit/draw 0 -6 4 -2 6\n" + "simundump xtree /edit/draw/tree 0 \\\n" + " /kinetics/#[],/kinetics/#[]/#[],/kinetics/#[]/#[]/#[][TYPE!=proto],/kinetics/#[]/#[]/#[][TYPE!=linkinfo]/##[] \"edit_elm.D <v>; drag_from_edit.w <d> <S> <x> <y> <z>\" auto 0.6\n" + "simundump xtext /file/notes 0 1\n") + storeReacMsg(reacList,f) + storeEnzMsg(enzList,f) + if tgraphs: + storePlotMsgs(tgraphs,f) + writeFooter1(f) + writeNotes(modelpath,f) + writeFooter2(f) + print('Written to file '+filename) + return errors, True + else: + print("Warning: writeKkit:: No model found on " , modelpath) + return False def findMinMax(sceneitems): cmin = 0.0 @@ -176,75 +207,82 @@ def storeEnzMsg( enzList, f): storeCplxEnzMsgs( enz, f ) def writeEnz( modelpath,f,sceneitems): + error = "" enzList = moose.wildcardFind(modelpath+'/##[ISA=EnzBase]') for enz in enzList: - x = random.randrange(0,10) - y = random.randrange(0,10) - textcolor = "" - color = "" - k1 = 0; - k2 = 0; - k3 = 0; - nInit = 0; - concInit = 0; - n = 0; - conc = 0; - enzParent = enz.parent - if (isinstance(enzParent.className,moose.Pool)) or (isinstance(enzParent.className,moose.ZombiePool)): - print(" raise exception enz doesn't have pool as parent") - return False + if findCompartment(enz) == moose.element('/'): + error = error + " \n "+enz.path+ " doesn't have compartment ignored to write to genesis" else: - vol = enzParent.volume * NA * 1e-3; - isMichaelisMenten = 0; - enzClass = enz.className - if (enzClass == "ZombieMMenz" or enzClass == "MMenz"): - k1 = enz.numKm - k3 = enz.kcat - k2 = 4.0*k3; - k1 = (k2 + k3) / k1; - isMichaelisMenten = 1; - - elif (enzClass == "ZombieEnz" or enzClass == "Enz"): - k1 = enz.k1 - k2 = enz.k2 - k3 = enz.k3 - cplx = enz.neighbors['cplx'][0] - nInit = cplx.nInit[0]; - if sceneitems != None: - # value = sceneitems[enz] - # x = calPrime(value['x']) - # y = calPrime(value['y']) - x = sceneitems[enz]['x'] - y = sceneitems[enz]['y'] - - einfo = enz.path+'/info' - if moose.exists(einfo): - color = moose.Annotator(einfo).getField('color') - color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) - - textcolor = moose.Annotator(einfo).getField('textColor') - textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) - - if color == "" or color == " ": - color = getRandColor() - if textcolor == "" or textcolor == " ": - textcolor = getRandColor() - - f.write("simundump kenz /kinetics/" + trimPath(enz) + " " + str(0)+ " " + - str(concInit) + " " + - str(conc) + " " + - str(nInit) + " " + - str(n) + " " + - str(vol) + " " + - str(k1) + " " + - str(k2) + " " + - str(k3) + " " + - str(0) + " " + - str(isMichaelisMenten) + " " + - "\"\"" + " " + - str(textcolor) + " " + str(color) + " \"\"" + - " " + str(int(x)) + " " + str(int(y)) + " "+str(0)+"\n") - return enzList + x = random.randrange(0,10) + y = random.randrange(0,10) + textcolor = "" + color = "" + k1 = 0; + k2 = 0; + k3 = 0; + nInit = 0; + concInit = 0; + n = 0; + conc = 0; + enzParent = enz.parent + if (isinstance(enzParent.className,moose.Pool)) or (isinstance(enzParent.className,moose.ZombiePool)): + print(" raise exception enz doesn't have pool as parent") + return False + else: + vol = enzParent.volume * NA * 1e-3; + isMichaelisMenten = 0; + enzClass = enz.className + if (enzClass == "ZombieMMenz" or enzClass == "MMenz"): + k1 = enz.numKm + k3 = enz.kcat + k2 = 4.0*k3; + k1 = (k2 + k3) / k1; + isMichaelisMenten = 1; + + elif (enzClass == "ZombieEnz" or enzClass == "Enz"): + k1 = enz.k1 + k2 = enz.k2 + k3 = enz.k3 + if enz.neighbors['cplx']: + cplx = enz.neighbors['cplx'][0] + nInit = cplx.nInit[0]; + else: + cplx = moose.Pool(enz.path+"/cplx") + moose.Annotator(cplx.path+'/info') + moose.connect( enz, 'cplx', cplx, 'reac' ) + nInit = cplx.nInit + + einfo = enz.path+'/info' + if moose.exists(einfo): + x = sceneitems[enz]['x'] + y = sceneitems[enz]['y'] + color = moose.Annotator(einfo).getField('color') + color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) + + textcolor = moose.Annotator(einfo).getField('textColor') + textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) + else: + error = error + "\n x and y co-ordinates are not specified for `" + enz.name+ "` zero will be assigned \n " + if color == "" or color == " ": + color = getRandColor() + if textcolor == "" or textcolor == " ": + textcolor = getRandColor() + + f.write("simundump kenz /kinetics/" + trimPath(enz) + " " + str(int(0))+ " " + + str(concInit) + " " + + str(conc) + " " + + str(nInit) + " " + + str(n) + " " + + str(vol) + " " + + str(k1) + " " + + str(k2) + " " + + str(k3) + " " + + str(0) + " " + + str(isMichaelisMenten) + " " + + "\"\"" + " " + + str(textcolor) + " " + str(color) + " \"\"" + + " " + str(int(x)) + " " + str(int(y)) + " "+str(int(0))+"\n") + return enzList,error def nearestColorIndex(color, color_sequence): #Trying to find the index to closest color map from the rainbow pickle file for matching the Genesis color map @@ -272,34 +310,44 @@ def storeReacMsg(reacList,f): f.write( s) def writeReac(modelpath,f,sceneitems): + error = "" reacList = moose.wildcardFind(modelpath+'/##[ISA=ReacBase]') - for reac in reacList : - color = "" - textcolor = "" - kf = reac.numKf - kb = reac.numKb - # if sceneitems != None: - # value = sceneitems[reac] - # x = calPrime(value['x']) - # y = calPrime(value['y']) - x = sceneitems[reac]['x'] - y = sceneitems[reac]['y'] - rinfo = reac.path+'/info' - if moose.exists(rinfo): - color = moose.Annotator(rinfo).getField('color') - color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) - - textcolor = moose.Annotator(rinfo).getField('textColor') - textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) - - if color == "" or color == " ": - color = getRandColor() - if textcolor == "" or textcolor == " ": - textcolor = getRandColor() - f.write("simundump kreac /kinetics/" + trimPath(reac) + " " +str(0) +" "+ str(kf) + " " + str(kb) + " \"\" " + - str(color) + " " + str(textcolor) + " " + str(int(x)) + " " + str(int(y)) + " 0\n") - return reacList + for reac in reacList: + if findCompartment(reac) == moose.element('/'): + error = error + " \n "+reac.path+ " doesn't have compartment ignored to write to genesis" + else: + color = "" + textcolor = "" + kf = reac.numKf + kb = reac.numKb + # if sceneitems != None: + # value = sceneitems[reac] + # x = calPrime(value['x']) + # y = calPrime(value['y']) + + + rinfo = reac.path+'/info' + if moose.exists(rinfo): + x = sceneitems[reac]['x'] + y = sceneitems[reac]['y'] + + color = moose.Annotator(rinfo).getField('color') + color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) + textcolor = moose.Annotator(rinfo).getField('textColor') + textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) + else: + x = 0 + y = 0 + error = error + "\n x and y co-ordinates are not specified for `" + reac.name+ "` zero will be assigned \n " + if color == "" or color == " ": + color = getRandColor() + if textcolor == "" or textcolor == " ": + textcolor = getRandColor() + f.write("simundump kreac /kinetics/" + trimPath(reac) + " " +str(0) +" "+ str(kf) + " " + str(kb) + " \"\" " + + str(color) + " " + str(textcolor) + " " + str(int(x)) + " " + str(int(y)) + " "+ str(0)+"\n") + return reacList,error + def trimPath(mobj): original = mobj mobj = moose.element(mobj) @@ -308,7 +356,7 @@ def trimPath(mobj): mobj = moose.element(mobj.parent) found = True if mobj.path == "/": - print(("compartment is not found with the given path and the path has reached root ",original)) + print(original, " object doesn't have compartment as a parent ") return #other than the kinetics compartment, all the othername are converted to group in Genesis which are place under /kinetics # Any moose object comes under /kinetics then one level down the path is taken. @@ -324,14 +372,55 @@ def trimPath(mobj): s = splitpath.replace("_dash_",'-') return s +# def writeSumtotal( modelpath,f): +# funclist = moose.wildcardFind(modelpath+'/##[ISA=Function]') +# for func in funclist: +# funcInputs = moose.element(func.path+'/x[0]') +# s = "" +# for funcInput in funcInputs.neighbors["input"]: +# s = s+ "addmsg /kinetics/" + trimPath(funcInput)+ " /kinetics/" + trimPath(moose.element(func.parent)) + " SUMTOTAL n nInit\n" +# f.write(s) + def writeSumtotal( modelpath,f): + error = "" funclist = moose.wildcardFind(modelpath+'/##[ISA=Function]') + s = "" for func in funclist: + fInfound = True + fOutfound = True funcInputs = moose.element(func.path+'/x[0]') - s = "" - for funcInput in funcInputs.neighbors["input"]: - s = s+ "addmsg /kinetics/" + trimPath(funcInput)+ " /kinetics/" + trimPath(moose.element(func.parent)) + " SUMTOTAL n nInit\n" - f.write(s) + if not len(funcInputs.neighbors["input"]): + fInfound = False + error = error +' \n /'+ (moose.element(func)).parent.name+ '/'+moose.element(func).name + ' function doesn\'t have input which is not allowed in genesis. \n This function is not written down into genesis file\n' + + if not len(func.neighbors["valueOut"]): + error = error +'Function'+func.path+' has not been connected to any output, this function is not written to genesis file' + fOutfound = False + else: + for srcfunc in func.neighbors["valueOut"]: + if srcfunc.className in ["ZombiePool","ZombieBufPool","Pool","BufPool"]: + functionOut = moose.element(srcfunc) + else: + error = error +' \n Function output connected to '+srcfunc.name+ ' which is a '+ srcfunc.className+' which is not allowed in genesis, this function '+(moose.element(func)).path+' is not written to file' + fOutfound = False + + if fInfound and fOutfound: + srcPool = [] + for funcInput in funcInputs.neighbors["input"]: + if funcInput not in srcPool: + srcPool.append(funcInput) + if trimPath(funcInput) != None and trimPath(functionOut) != None: + s = "addmsg /kinetics/" + trimPath(funcInput)+ " /kinetics/"+ trimPath(functionOut)+ " SUMTOTAL n nInit\n" + f.write(s) + else: + error = error + '\n Genesis doesn\'t allow same moluecule connect to function mutiple times. \n Pool \''+ moose.element(funcInput).name + '\' connected to '+ (moose.element(func)).path + return error + +def writeStimulus(modelpath,f): + error = "" + if len(moose.wildcardFind(modelpath+'/##[ISA=StimulusTable]')): + error = error +'\n StimulusTable is not written into genesis. This is in Todo List' + return error def storePlotMsgs( tgraphs,f): s = "" @@ -407,55 +496,70 @@ def writeplot( tgraphs,f ): return first,second def writePool(modelpath,f,volIndex,sceneitems): + error = "" color = "" textcolor = "" for p in moose.wildcardFind(modelpath+'/##[ISA=PoolBase]'): - slave_enable = 0 - if (p.className == "BufPool" or p.className == "ZombieBufPool"): - pool_children = p.children - if pool_children== 0: - slave_enable = 4 - else: - for pchild in pool_children: - if not(pchild.className == "ZombieFunction") and not(pchild.className == "Function"): - slave_enable = 4 - else: - slave_enable = 0 - break - if (p.parent.className != "Enz" and p.parent.className !='ZombieEnz'): - #Assuming "p.parent.className !=Enzyme is cplx which is not written to genesis" - x = sceneitems[p]['x'] - y = sceneitems[p]['y'] - # if sceneitems != None: - # value = sceneitems[p] - # x = calPrime(value['x']) - # y = calPrime(value['y']) - - pinfo = p.path+'/info' - if moose.exists(pinfo): - color = moose.Annotator(pinfo).getField('color') - color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) - textcolor = moose.Annotator(pinfo).getField('textColor') - textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) - - geometryName = volIndex[p.volume] - volume = p.volume * NA * 1e-3 - if color == "" or color == " ": - color = getRandColor() - if textcolor == "" or textcolor == " ": - textcolor = getRandColor() - #print " trimPath",trimPath(p) - f.write("simundump kpool /kinetics/" + trimPath(p) + " 0 " + - str(p.diffConst) + " " + - str(0) + " " + - str(0) + " " + - str(0) + " " + - str(p.nInit) + " " + - str(0) + " " + str(0) + " " + - str(volume)+ " " + - str(slave_enable) + - " /kinetics"+ geometryName + " " + - str(color) +" " + str(textcolor) + " " + str(int(x)) + " " + str(int(y)) + " "+ str(0)+"\n") + if findCompartment(p) == moose.element('/'): + error = error + " \n "+p.path+ " doesn't have compartment ignored to write to genesis" + else: + slave_enable = 0 + if (p.className == "BufPool" or p.className == "ZombieBufPool"): + pool_children = p.children + if pool_children== 0: + slave_enable = 4 + else: + for pchild in pool_children: + if not(pchild.className == "ZombieFunction") and not(pchild.className == "Function"): + slave_enable = 4 + else: + slave_enable = 0 + break + if (p.parent.className != "Enz" and p.parent.className !='ZombieEnz'): + #Assuming "p.parent.className !=Enzyme is cplx which is not written to genesis" + + # x = sceneitems[p]['x'] + # y = sceneitems[p]['y'] + # if sceneitems != None: + # value = sceneitems[p] + # x = calPrime(value['x']) + # y = calPrime(value['y']) + + pinfo = p.path+'/info' + if moose.exists(pinfo): + x = sceneitems[p]['x'] + y = sceneitems[p]['y'] + else: + x = 0 + y = 0 + error = error + " \n x and y co-ordinates are not specified for `" + p.name+ "` zero will be assigned"+ " \n " + + if moose.exists(pinfo): + color = moose.Annotator(pinfo).getField('color') + color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) + textcolor = moose.Annotator(pinfo).getField('textColor') + textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) + + + geometryName = volIndex[p.volume] + volume = p.volume * NA * 1e-3 + if color == "" or color == " ": + color = getRandColor() + if textcolor == "" or textcolor == " ": + textcolor = getRandColor() + + f.write("simundump kpool /kinetics/" + trimPath(p) + " 0 " + + str(p.diffConst) + " " + + str(0) + " " + + str(0) + " " + + str(0) + " " + + str(p.nInit) + " " + + str(0) + " " + str(0) + " " + + str(volume)+ " " + + str(slave_enable) + + " /kinetics"+ geometryName + " " + + str(color) +" " + str(textcolor) + " " + str(int(x)) + " " + str(int(y)) + " "+ str(0)+"\n") + return error def getColorCheck(color,GENESIS_COLOR_SEQUENCE): if isinstance(color, str): @@ -486,12 +590,12 @@ def getColorCheck(color,GENESIS_COLOR_SEQUENCE): else: raise Exception("Invalid Color Value!") -ignoreColor= ["mistyrose","antiquewhite","aliceblue","azure","bisque","black","blanchedalmond","blue","cornsilk","darkolivegreen","darkslategray","dimgray","floralwhite","gainsboro","ghostwhite","honeydew","ivory","lavender","lavenderblush","lemonchiffon","lightcyan","lightgoldenrodyellow","lightgray","lightyellow","linen","mediumblue","mintcream","navy","oldlace","papayawhip","saddlebrown","seashell","snow","wheat","white","whitesmoke","aquamarine","lightsalmon","moccasin","limegreen","snow","sienna","beige","dimgrey","lightsage"] -matplotcolor = {} -for name,hexno in matplotlib.colors.cnames.items(): - matplotcolor[name]=hexno - def getRandColor(): + ignoreColor= ["mistyrose","antiquewhite","aliceblue","azure","bisque","black","blanchedalmond","blue","cornsilk","darkolivegreen","darkslategray","dimgray","floralwhite","gainsboro","ghostwhite","honeydew","ivory","lavender","lavenderblush","lemonchiffon","lightcyan","lightgoldenrodyellow","lightgray","lightyellow","linen","mediumblue","mintcream","navy","oldlace","papayawhip","saddlebrown","seashell","snow","wheat","white","whitesmoke","aquamarine","lightsalmon","moccasin","limegreen","snow","sienna","beige","dimgrey","lightsage"] + matplotcolor = {} + for name,hexno in matplotlib.colors.cnames.items(): + matplotcolor[name]=hexno + k = random.choice(list(matplotcolor.keys())) if k in ignoreColor: return getRandColor() @@ -536,6 +640,7 @@ def writeGroup(modelpath,f): f.write("simundump group /kinetics/" + trimPath(g) + " 0 " + "blue" + " " + "green" + " x 0 0 \"\" defaultfile \\\n") f.write(" defaultfile.g 0 0 0 " + str(int(x)) + " " + str(int(y)) + " 0\n") + def writeHeader(f,maxVol): simdt = 0.001 plotdt = 0.1 @@ -619,12 +724,18 @@ def writeFooter2(f): if __name__ == "__main__": import sys - + import os filename = sys.argv[1] - modelpath = filename[0:filename.find('.')] + filepath, filenameWithext = os.path.split(filename) + if filenameWithext.find('.') != -1: + modelpath = filenameWithext[:filenameWithext.find('.')] + else: + modelpath = filenameWithext + moose.loadModel(filename,'/'+modelpath,"gsl") output = modelpath+"_.g" - written = write('/'+modelpath,output) + written = mooseWriteKkit('/'+modelpath,output) + if written: print(" file written to ",output) else: diff --git a/moose-core/python/moose/moose.py b/moose-core/python/moose/moose.py index 3fa32f5c371d8b81f8dd34211f00625d5ad6ac59..7d045d3793fc3c711c75ad23224fda085f6dece3 100644 --- a/moose-core/python/moose/moose.py +++ b/moose-core/python/moose/moose.py @@ -9,9 +9,23 @@ import warnings import pydoc from io import StringIO -import moose.SBML.readSBML as _readSBML -import moose.SBML.writeSBML as _writeSBML -import moose.chemUtil as _chemUtil +import moose + +sbmlImport_, sbmlError_ = True, '' + +try: + import moose.SBML.readSBML as _readSBML + import moose.SBML.writeSBML as _writeSBML +except Exception as e: + sbmlImport_ = False + sbmlError_ = '%s' % e + +chemImport_, chemError_ = True, '' +try: + import moose.chemUtil as _chemUtil +except Exception as e: + chemImport_ = False + chemError_ = '%s' % e kkitImport_, kkitImport_error_ = True, '' try: @@ -20,6 +34,14 @@ except ImportError as e: kkitImport_ = False kkitImport_err_ = '%s' % e +mergechemImport_, mergechemError_ = True, '' +try: + import moose.chemMerge as _chemMerge +except Exception as e: + mergechemImport_ = False + mergechemError_ = '%s' % e + + # Import function from C++ module into moose namespace. from moose._moose import * @@ -63,7 +85,12 @@ def mooseReadSBML(filepath, loadpath, solver='ee'): solver -- Solver to use (default 'ee' ) \n """ - return _readSBML.mooseReadSBML( filepath, loadpath, solver ) + global sbmlImport_ + if sbmlImport_: + return _readSBML.mooseReadSBML( filepath, loadpath, solver ) + else: + print( sbmlError_ ) + return False def mooseWriteSBML(modelpath, filepath, sceneitems={}): @@ -82,7 +109,11 @@ def mooseWriteSBML(modelpath, filepath, sceneitems={}): --- else, auto-coordinates is used for layout position and passed """ - return _writeSBML.mooseWriteSBML(modelpath, filepath, sceneitems) + if sbmlImport_: + return _writeSBML.mooseWriteSBML(modelpath, filepath, sceneitems) + else: + print( sbmlError_ ) + return False def mooseWriteKkit(modelpath, filepath,sceneitems={}): @@ -107,7 +138,11 @@ def moosedeleteChemSolver(modelpath): this should be followed by mooseaddChemSolver for add solvers on to compartment to simulate else default is Exponential Euler (ee) """ - return _chemUtil.add_Delete_ChemicalSolver.moosedeleteChemSolver(modelpath) + if chemImport_: + return _chemUtil.add_Delete_ChemicalSolver.moosedeleteChemSolver(modelpath) + else: + print( chemError_ ) + return False def mooseaddChemSolver(modelpath, solver): @@ -121,7 +156,22 @@ def mooseaddChemSolver(modelpath, solver): "Runge Kutta" ("gsl") """ - return _chemUtil.add_Delete_ChemicalSolver.mooseaddChemSolver(modelpath, solver) + if chemImport_: + return _chemUtil.add_Delete_ChemicalSolver.mooseaddChemSolver(modelpath, solver) + else: + print( chemError_ ) + return False + +def mergeChemModel(src, des): + """ Merges two chemical model, \n + File or filepath can be passed + source is merged to destination + """ + #global mergechemImport_ + if mergechemImport_: + return _chemMerge.merge.mergeChemModel(src,des) + else: + return False ################################################################ # Wrappers for global functions @@ -134,7 +184,7 @@ def pwe(): the path, use moose.getCwe() """ - pwe_ = _moose.getCwe() + pwe_ = moose.getCwe() print(pwe_.getPath()) return pwe_ @@ -191,10 +241,10 @@ def syncDataHandler(target): raise NotImplementedError('The implementation is not working for IntFire - goes to invalid objects. \ First fix that issue with SynBase or something in that line.') if isinstance(target, str): - if not _moose.exists(target): + if not moose.exists(target): raise ValueError('%s: element does not exist.' % (target)) target = vec(target) - _moose.syncDataHandler(target) + moose.syncDataHandler(target) def showfield(el, field='*', showtype=False): @@ -349,7 +399,7 @@ def getfielddoc(tokens, indent=''): fieldname = tokens[1] while True: try: - classelement = _moose.element('/classes/' + classname) + classelement = moose.element('/classes/' + classname) for finfo in classelement.children: for fieldelement in finfo: baseinfo = '' @@ -412,7 +462,7 @@ def getmoosedoc(tokens, inherited=False): if not tokens: return "" try: - class_element = _moose.element('/classes/%s' % (tokens[0])) + class_element = moose.element('/classes/%s' % (tokens[0])) except ValueError: raise NameError('name \'%s\' not defined.' % (tokens[0])) if len(tokens) > 1: @@ -421,14 +471,14 @@ def getmoosedoc(tokens, inherited=False): docstring.write(toUnicode('%s\n' % (class_element.docs))) append_finfodocs(tokens[0], docstring, indent) if inherited: - mro = eval('_moose.%s' % (tokens[0])).mro() + mro = eval('moose.%s' % (tokens[0])).mro() for class_ in mro[1:]: - if class_ == _moose.melement: + if class_ == moose.melement: break docstring.write(toUnicode( '\n\n#Inherited from %s#\n' % (class_.__name__))) append_finfodocs(class_.__name__, docstring, indent) - if class_ == _moose.Neutral: # Neutral is the toplevel moose class + if class_ == moose.Neutral: # Neutral is the toplevel moose class break return docstring.getvalue() @@ -436,13 +486,13 @@ def getmoosedoc(tokens, inherited=False): def append_finfodocs(classname, docstring, indent): """Append list of finfos in class name to docstring""" try: - class_element = _moose.element('/classes/%s' % (classname)) + class_element = moose.element('/classes/%s' % (classname)) except ValueError: raise NameError('class \'%s\' not defined.' % (classname)) for ftype, rname in finfotypes: docstring.write(toUnicode('\n*%s*\n' % (rname.capitalize()))) try: - finfo = _moose.element('%s/%s' % (class_element.path, ftype)) + finfo = moose.element('%s/%s' % (class_element.path, ftype)) for field in finfo.vec: docstring.write(toUnicode( '%s%s: %s\n' % (indent, field.fieldName, field.type))) diff --git a/moose-core/python/moose/plot_utils.py b/moose-core/python/moose/plot_utils.py index d3f10c92b7e4241f08ef47700966159efafe6d77..9b309e46e92fcee319b4084ae019b27b19b3dc22 100644 --- a/moose-core/python/moose/plot_utils.py +++ b/moose-core/python/moose/plot_utils.py @@ -1,9 +1,8 @@ # -*- coding: utf-8 -*- -"""plot_utils.py: Some utility function for plotting data in moose. +# +# plot_utils.py: Some utility function for plotting data in moose. -Last modified: Sun Jan 10, 2016 04:04PM - -""" +from __future__ import print_function, division, absolute_import __author__ = "Dilawar Singh" __copyright__ = "Copyright 2013, NCBS Bangalore" @@ -15,10 +14,14 @@ __email__ = "dilawars@ncbs.res.in" __status__ = "Development" import numpy as np +import matplotlib.pyplot as plt import moose -import print_utils as pu +import moose.print_utils as pu +import re -import matplotlib.pyplot as plt +# To support python 2.6. On Python3.6+, dictionaries are ordered by default. +# This is here to fill this gap. +from moose.OrderedDict import OrderedDict def plotAscii(yvec, xvec = None, file=None): """Plot two list-like object in terminal using gnuplot. @@ -235,7 +238,7 @@ def saveRecords(records, xvec = None, **kwargs): def plotRecords(records, xvec = None, **kwargs): """Wrapper """ - dataDict = {} + dataDict = OrderedDict( ) try: for k in sorted(records.keys(), key=str.lower): dataDict[k] = records[k] @@ -284,12 +287,21 @@ def plotRecords(records, xvec = None, **kwargs): plt.close( ) -def plotTables( records, xvec = None, **kwargs ): - """Plot dictionary of moose.Table/moose.Table2 +def plotTables( regex = '.*', **kwargs ): + """plotTables Plot all moose.Table/moose.Table2 matching given regex. By + default plot all tables. Table names must be unique. Table name are used as + legend. - :param records: A dictionary of moose.Table. All tables must have same - length. - :param xvec: If None, moose.Clock is used to generate time-vector. + :param regex: Python regular expression to be matched. :param **kwargs: + - subplot = True/False; if True, each Table is plotted in a subplot. + - outfile = filepath; If given, plot will be saved to this path. """ - return plotRecords( records, xvec, **kwargs ) + + tables = moose.wildcardFind( '/##[TYPE=Table]' ) + tables += moose.wildcardFind( '/##[TYPE=Table2]' ) + toPlot = OrderedDict( ) + for t in sorted(tables, key = lambda x: x.name): + if re.search( regex, t.name ): + toPlot[ t.name ] = t + return plotRecords( toPlot, None, **kwargs ) diff --git a/moose-core/python/moose/print_utils.py b/moose-core/python/moose/print_utils.py index b3acd81bcd9b8dba28f95b842897107aa3b1a2f6..5cc08416a1c7784219272c2a779522d08f953748 100644 --- a/moose-core/python/moose/print_utils.py +++ b/moose-core/python/moose/print_utils.py @@ -1,10 +1,9 @@ # -*- coding: utf-8 -*- -"""print_utils.py: A library with some print functions. Very useful during -development. +# print_utils.py: +# +# A library with some print functions. Very useful during development. -Last modified: Sat Jan 18, 2014 05:01PM - -""" +from __future__ import print_function, division, absolute_import __author__ = "Dilawar Singh" __copyright__ = "Copyright 2013, NCBS Bangalore" @@ -15,9 +14,7 @@ __maintainer__ = "Dilawar Singh" __email__ = "dilawars@ncbs.res.in" __status__ = "Development" -import inspect -import sys -from . import _moose as moose +import moose HEADER = '\033[95m' OKBLUE = '\033[94m' @@ -79,7 +76,9 @@ def dump(label, msg, frame=None, exception=None): prefix = '[{0}] '.format(label) - ''' Enable it if you want indented messages + ''' + # Enable it if you want indented messages + import inspect stackLength = len(inspect.stack()) - 1 if stackLength == 1: prefix = '\n[{}] '.format(label) @@ -102,8 +101,11 @@ def dump(label, msg, frame=None, exception=None): if exception: print((" [Expcetion] {0}".format(exception))) -def info(msg): dump("INFO", msg) -def warn(msg): dump("WARN", msg) +def info(msg): + dump("INFO", msg) + +def warn(msg): + dump("WARN", msg) def error(msg): dump("ERROR", msg) diff --git a/moose-core/python/moose/utils.py b/moose-core/python/moose/utils.py index 5a288cc04f90b52e477e5e7eb7f2dcc9e14cf48b..d4fcf16e06632ea20e16b02933cf12c9bc57b3aa 100644 --- a/moose-core/python/moose/utils.py +++ b/moose-core/python/moose/utils.py @@ -1,15 +1,11 @@ # -*- coding: utf-8 -*- -""" -utils.py: - -Utility functions for moose. +# utils.py: +# +# Utility functions for moose. +# +# NOTE: Some function might break because unicode is default string in python3. -NOTE: Some function might break because unicode is default string in python3. - -""" - -from __future__ import print_function, division -from __future__ import absolute_import +from __future__ import print_function, division, absolute_import __author__ = 'Subhasis Ray, Aditya Gilra, Dilawar Singh' __copyright__ = "Copyright 2013, NCBS Bangalore" @@ -30,17 +26,9 @@ import re from moose.moose_constants import * -# Make these import non-important. -try: - from moose.plot_utils import * -except Exception as e: - pass - -try: - from moose.print_utils import * -except Exception as e: - pass - +# Print and Plot utilities. +from moose.plot_utils import * +from moose.print_utils import * def create_table_path(model, graph, element, field): diff --git a/moose-core/python/rdesigneur/rdesigneur.py b/moose-core/python/rdesigneur/rdesigneur.py index e5704f2a04e0d0787224e18b2075e92c793e161d..4220f3763ecc8e5ba84062b1d82cedcbe775148d 100644 --- a/moose-core/python/rdesigneur/rdesigneur.py +++ b/moose-core/python/rdesigneur/rdesigneur.py @@ -17,18 +17,15 @@ ## channel conductances, between them. ########################################################################## from __future__ import print_function -from __future__ import absolute_import - import imp import os import moose import numpy as np import pylab import math -import rdesigneur.rmoogli - -from rdesigneur.rdesigneurProtos import * - +import rmoogli +#import rdesigneurProtos +from rdesigneurProtos import * from moose.neuroml.NeuroML import NeuroML from moose.neuroml.ChannelML import ChannelML @@ -84,6 +81,7 @@ class rdesigneur: useGssa = True, combineSegments = True, stealCellFromLibrary = False, + verbose = True, diffusionLength= 2e-6, meshLambda = -1.0, #This is a backward compatibility hack temperature = 32, @@ -115,6 +113,7 @@ class rdesigneur: self.useGssa = useGssa self.combineSegments = combineSegments self.stealCellFromLibrary = stealCellFromLibrary + self.verbose = verbose self.diffusionLength= diffusionLength if meshLambda > 0.0: print("Warning: meshLambda argument is deprecated. Please use 'diffusionLength' instead.\nFor now rdesigneur will accept this argument.") @@ -205,7 +204,8 @@ class rdesigneur: self._buildMoogli() self._buildStims() self._configureClocks() - self._printModelStats() + if self.verbose: + self._printModelStats() self._savePlots() except BuildError as msg: @@ -319,7 +319,8 @@ class rdesigneur: return True if moose.exists( '/library/' + protoVec[0] ): #moose.copy('/library/' + protoVec[0], '/library/', protoVec[1]) - print('renaming /library/' + protoVec[0] + ' to ' + protoVec[1]) + if self.verbose: + print('renaming /library/' + protoVec[0] + ' to ' + protoVec[1]) moose.element( '/library/' + protoVec[0]).name = protoVec[1] #moose.le( '/library' ) return True @@ -458,6 +459,7 @@ class rdesigneur: "buildChemDistrib: No elec compartments found in path: '" \ + pair + "'" ) self.spineComptElist = self.elecid.spinesFromExpression[ pair ] + #print( 'LEN SPINECOMPTELIST =' + str( pair ) + ", " str( len( self.spineComptElist ) ) ) ''' if len( self.spineComptElist ) == 0: raise BuildError( \ @@ -522,7 +524,7 @@ class rdesigneur: objList = self._collapseElistToPathAndClass( comptList, plotSpec[2], kf[0] ) # print ("objList: ", len(objList), kf[1]) return objList, kf[1] - elif (field == 'n' or field == 'conc' ): + elif (field == 'n' or field == 'conc' or field == 'volume' ): path = plotSpec[2] pos = path.find( '/' ) if pos == -1: # Assume it is in the dend compartment. @@ -558,15 +560,18 @@ class rdesigneur: def _buildPlots( self ): knownFields = { 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ), + 'spikeTime':('CompartmentBase', 'getVm', 1, 'Spike Times (s)'), 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ), 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ), 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ), + 'modulation':('ChanBase', 'getModulation', 1, 'chan modulation (unitless)' ), 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)' ), 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)' ), 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)' ), 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ), 'n':('PoolBase', 'getN', 1, '# of molecules'), - 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ) + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ), + 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ) } graphs = moose.Neutral( self.modelPath + '/graphs' ) dummy = moose.element( '/' ) @@ -579,17 +584,21 @@ class rdesigneur: plotObj2, plotField2 = self._parseComptField( spineCompts, i, knownFields ) assert( plotField == plotField2 ) plotObj3 = plotObj + plotObj2 - numPlots = sum( i != dummy for i in plotObj3 ) + numPlots = sum( q != dummy for q in plotObj3 ) if numPlots > 0: tabname = graphs.path + '/plot' + str(k) scale = knownFields[i[3]][2] units = knownFields[i[3]][3] - self.plotNames.append( ( tabname, i[4], k, scale, units ) ) + self.plotNames.append( ( tabname, i[4], k, scale, units, i[3] ) ) k += 1 - if i[3] == 'n' or i[3] == 'conc': + if i[3] == 'n' or i[3] == 'conc' or i[3] == 'volume' or i[3] == 'Gbar': tabs = moose.Table2( tabname, numPlots ) else: tabs = moose.Table( tabname, numPlots ) + if i[3] == 'spikeTime': + tabs.vec.threshold = -0.02 # Threshold for classifying Vm as a spike. + tabs.vec.useSpikeMode = True # spike detect mode on + vtabs = moose.vec( tabs ) q = 0 for p in [ x for x in plotObj3 if x != dummy ]: @@ -599,6 +608,7 @@ class rdesigneur: def _buildMoogli( self ): knownFields = { 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)', -80.0, 40.0 ), + 'initVm':('CompartmentBase', 'getInitVm', 1000, 'Init. Memb. Potl (mV)', -80.0, 40.0 ), 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)', -10.0, 10.0 ), 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)', -10.0, 10.0 ), 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)', 0.0, 1.0 ), @@ -607,7 +617,8 @@ class rdesigneur: 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)', -10.0, 10.0 ), 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)', 0.0, 10.0 ), 'n':('PoolBase', 'getN', 1, '# of molecules', 0.0, 200.0 ), - 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)', 0.0, 2.0 ) + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)', 0.0, 2.0 ), + 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ) } moogliBase = moose.Neutral( self.modelPath + '/moogli' ) k = 0 @@ -639,14 +650,24 @@ class rdesigneur: def display( self ): for i in self.plotNames: + pylab.figure( i[2] ) pylab.title( i[1] ) pylab.xlabel( "Time (s)" ) pylab.ylabel( i[4] ) vtab = moose.vec( i[0] ) - t = np.arange( 0, vtab[0].vector.size, 1 ) * vtab[0].dt - for j in vtab: - pylab.plot( t, j.vector * i[3] ) + if i[5] == 'spikeTime': + k = 0 + tmax = moose.element( '/clock' ).currentTime + for j in vtab: # Plot a raster + y = [k] * len( j.vector ) + pylab.plot( j.vector * i[3], y, linestyle = 'None', marker = '.', markersize = 10 ) + pylab.xlim( 0, tmax ) + + else: + t = np.arange( 0, vtab[0].vector.size, 1 ) * vtab[0].dt + for j in vtab: + pylab.plot( t, j.vector * i[3] ) if len( self.moogList ) > 0: pylab.ion() pylab.show(block=True) @@ -666,15 +687,18 @@ class rdesigneur: knownFields = { 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ), + 'spikeTime':('CompartmentBase', 'getVm', 1, 'Spike Times (s)'), 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ), 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ), 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ), + 'modulation':('ChanBase', 'getModulation', 1, 'chan modulation (unitless)' ), 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)' ), 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)' ), 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)' ), 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ), 'n':('PoolBase', 'getN', 1, '# of molecules'), - 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ) + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ), + 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ) } save_graphs = moose.Neutral( self.modelPath + '/save_graphs' ) @@ -696,11 +720,15 @@ class rdesigneur: units = knownFields[i[3]][3] self.saveNames.append( ( save_tabname, i[4], k, scale, units ) ) k += 1 - if i[3] == 'n' or i[3] == 'conc': + if i[3] == 'n' or i[3] == 'conc' or i[3] == 'volume' or i[3] == 'Gbar': save_tabs = moose.Table2( save_tabname, numPlots ) + save_vtabs = moose.vec( save_tabs ) else: save_tabs = moose.Table( save_tabname, numPlots ) - save_vtabs = moose.vec( save_tabs ) + save_vtabs = moose.vec( save_tabs ) + if i[3] == 'spikeTime': + save_vtabs.threshold = -0.02 # Threshold for classifying Vm as a spike. + save_vtabs.useSpikeMode = True # spike detect mode on q = 0 for p in [ x for x in plotObj3 if x != dummy ]: moose.connect( save_vtabs[q], 'requestOut', p, plotField ) @@ -717,6 +745,7 @@ class rdesigneur: knownFields = { 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ), + 'spikeTime':('CompartmentBase', 'getVm', 1, 'Spike Times (s)'), 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ), 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ), 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ), @@ -725,7 +754,8 @@ class rdesigneur: 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)' ), 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ), 'n':('PoolBase', 'getN', 1, '# of molecules'), - 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ) + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ), + 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ) } ''' @@ -893,18 +923,17 @@ class rdesigneur: if self.turnOffElec: elecDt = 1e6 elecPlotDt = 1e6 - diffDt = 0.1 # Slow it down again because no multiscaling - chemDt = 0.1 # Slow it down again because no multiscaling else: elecDt = self.elecDt - diffDt = self.diffDt - chemDt = self.chemDt + elecPlotDt = self.elecPlotDt + diffDt = self.diffDt + chemDt = self.chemDt for i in range( 0, 9 ): moose.setClock( i, elecDt ) + moose.setClock( 8, elecPlotDt ) moose.setClock( 10, diffDt ) for i in range( 11, 18 ): moose.setClock( i, chemDt ) - moose.setClock( 8, self.elecPlotDt ) moose.setClock( 18, self.chemPlotDt ) hsolve = moose.HSolve( self.elecid.path + '/hsolve' ) hsolve.dt = elecDt @@ -1231,7 +1260,14 @@ class rdesigneur: mesh = moose.element( '/model/chem/' + meshName ) #elecComptList = mesh.elecComptList if elecRelPath == 'spine': - elecComptList = moose.vec( mesh.elecComptList[0].path + '/../spine' ) + # This is nasty. The spine indexing is different from + # the compartment indexing and the mesh indexing and the + # chem indexing. Need to fix at some time. + #elecComptList = moose.vec( mesh.elecComptList[0].path + '/../spine' ) + elecComptList = moose.element( '/model/elec').spineIdsFromCompartmentIds[ mesh.elecComptList ] + #print( len( mesh.elecComptList ) ) + # for i,j in zip( elecComptList, mesh.elecComptList ): + # print( "Lookup: {} {} {}; orig: {} {} {}".format( i.name, i.index, i.fieldIndex, j.name, j.index, j.fieldIndex )) else: elecComptList = mesh.elecComptList @@ -1270,18 +1306,23 @@ class rdesigneur: #print 'building ', len( elecComptList ), 'adaptors ', adName, ' for: ', mesh.name, elecRelPath, elecField, chemRelPath av = ad.vec chemVec = moose.element( mesh.path + '/' + chemRelPath ).vec + root = moose.element( '/' ) for i in zip( elecComptList, startVoxelInCompt, endVoxelInCompt, av ): i[3].inputOffset = 0.0 i[3].outputOffset = offset i[3].scale = scale if elecRelPath == 'spine': + # Check needed in case there were unmapped entries in + # spineIdsFromCompartmentIds elObj = i[0] + if elObj == root: + continue else: ePath = i[0].path + '/' + elecRelPath if not( moose.exists( ePath ) ): - raise BuildError( \ - "Error: buildAdaptor: no elec obj in " + ePath ) + continue + #raise BuildError( "Error: buildAdaptor: no elec obj in " + ePath ) elObj = moose.element( i[0].path + '/' + elecRelPath ) if ( isElecToChem ): elecFieldSrc = 'get' + capField diff --git a/moose-core/python/setup.py b/moose-core/python/setup.py index f8d537d0129e531d60288284535a50d0a6c649ab..93c432cd4afaeda4153e7a2d0e8909f27764da71 100644 --- a/moose-core/python/setup.py +++ b/moose-core/python/setup.py @@ -64,7 +64,7 @@ setup( , 'moose.neuroml' , 'moose.genesis' , 'moose.chemUtil' - , 'moose.merge' + , 'moose.chemMerge' ], package_dir = { 'moose' : 'moose' diff --git a/moose-core/tests/python/test_gc.py b/moose-core/tests/python/deprecated_test_gc.py similarity index 100% rename from moose-core/tests/python/test_gc.py rename to moose-core/tests/python/deprecated_test_gc.py diff --git a/moose-core/tests/python/test_connectionLists.py b/moose-core/tests/python/test_connectionLists.py new file mode 100644 index 0000000000000000000000000000000000000000..057cd5a47749d474f6fb05e1c2e2d987f6b69367 --- /dev/null +++ b/moose-core/tests/python/test_connectionLists.py @@ -0,0 +1,158 @@ +"""test_connectionLists.py: + +Test connectionList in SparseMsg. + +""" + +from __future__ import print_function + +__author__ = "Dilawar Singh" +__copyright__ = "Copyright 2016, Dilawar Singh" +__credits__ = ["NCBS Bangalore"] +__license__ = "GNU GPL" +__version__ = "1.0.0" +__maintainer__ = "Dilawar Singh" +__email__ = "dilawars@ncbs.res.in" +__status__ = "Development" + +import sys +import os +import moose +import numpy as np + +params = { + 'Rm': 1e10, + 'Cm': 1e-11, + 'randInputRate': 10, + 'randRefractTime': 0.005, + 'numInputs': 100, + 'numInhib': 10, + 'inhibRefractTime': 0.005, + 'inhibThresh': -0.04, + 'numOutput': 10, + 'outputRefractTime': 0.005, + 'outputThresh': 5.00e6, + 'stimToInhProb': 0.01, + 'stimToInhSeed': 11234, + 'inhVreset': -0.06, + 'stimToOutProb': 0.05, + 'stimToOutSeed': 1234, + 'inhToOutProb': 0.5, + 'inhToOutSeed': 1234, + 'runtime': 2.0, + 'wtStimToOut': 0.001, + 'wtStimToInh': 0.002, + 'wtInhToOut': -0.002, +} + + +def makeGlobalBalanceNetwork(): + stim = moose.RandSpike( '/model/stim', params['numInputs'] ) + inhib = moose.LIF( '/model/inhib', params['numInhib'] ) + insyn = moose.SimpleSynHandler(inhib.path + '/syns', params['numInhib']) + moose.connect( insyn, 'activationOut', inhib, 'activation', 'OneToOne' ) + output = moose.LIF( '/model/output', params['numOutput'] ) + outsyn = moose.SimpleSynHandler(output.path+'/syns',params['numOutput']) + moose.connect(outsyn, 'activationOut', output, 'activation', 'OneToOne') + outInhSyn = moose.SimpleSynHandler(output.path+'/inhsyns',params['numOutput']) + moose.connect(outInhSyn, 'activationOut', output, 'activation', 'OneToOne') + + iv = moose.vec( insyn.path + '/synapse' ) + ov = moose.vec( outsyn.path + '/synapse' ) + oiv = moose.vec( outInhSyn.path + '/synapse' ) + + temp = moose.connect( stim, 'spikeOut', iv, 'addSpike', 'Sparse' ) + inhibMatrix = moose.element( temp ) + inhibMatrix.setRandomConnectivity( + params['stimToInhProb'], params['stimToInhSeed'] ) + cl = inhibMatrix.connectionList + expectedCl = [ 1,4,13,13,26,42,52,56,80,82,95,97,4,9,0,9,4,8,0,6,1,6,6,7] + assert list(cl) == expectedCl, "Expected %s, got %s" % (expectedCl, cl ) + + temp = moose.connect( stim, 'spikeOut', ov, 'addSpike', 'Sparse' ) + excMatrix = moose.element( temp ) + excMatrix.setRandomConnectivity( + params['stimToOutProb'], params['stimToOutSeed'] ) + + temp = moose.connect( inhib, 'spikeOut', oiv, 'addSpike', 'Sparse' ) + negFFMatrix = moose.element( temp ) + negFFMatrix.setRandomConnectivity( + params['inhToOutProb'], params['inhToOutSeed'] ) + + # print("ConnMtxEntries: ", inhibMatrix.numEntries, excMatrix.numEntries, negFFMatrix.numEntries) + assert (12, 57, 48) == (inhibMatrix.numEntries, excMatrix.numEntries, negFFMatrix.numEntries) + + cl = negFFMatrix.connectionList + numInhSyns = [ ] + niv = 0 + nov = 0 + noiv = 0 + for i in moose.vec( insyn ): + niv += i.synapse.num + numInhSyns.append( i.synapse.num ) + if i.synapse.num > 0: + i.synapse.weight = params['wtStimToInh'] + + assert numInhSyns == [2,1,0,0,2,0,3,1,1,2] + + for i in moose.vec( outsyn ): + nov += i.synapse.num + if i.synapse.num > 0: + i.synapse.weight = params['wtStimToOut'] + for i in moose.vec( outInhSyn ): + noiv += i.synapse.num + #print i.synapse.num + if i.synapse.num > 0: + i.synapse.weight = params['wtInhToOut'] + + # print("SUMS: ", sum( iv.numField ), sum( ov.numField ), sum( oiv.numField )) + assert [4,49,9] == [sum( iv.numField ), sum( ov.numField ), sum( oiv.numField )] + # print("SUMS2: ", niv, nov, noiv) + assert [12,57,48] == [ niv, nov, noiv ] + # print("SUMS3: ", sum( insyn.vec.numSynapses ), sum( outsyn.vec.numSynapses ), sum( outInhSyn.vec.numSynapses )) + assert [12, 57, 48 ] == [ sum( insyn.vec.numSynapses ), sum( outsyn.vec.numSynapses ), sum( outInhSyn.vec.numSynapses ) ] + + # print(oiv.numField) + # print(insyn.vec[1].synapse.num) + # print(insyn.vec.numSynapses) + # print(sum( insyn.vec.numSynapses )) + # niv = iv.numSynapses + # ov = iv.numSynapses + + sv = moose.vec( stim ) + sv.rate = params['randInputRate'] + sv.refractT = params['randRefractTime'] + #moose.showfield( sv[7] ) + + inhib.vec.thresh = params['inhibThresh'] + inhib.vec.Rm = params['Rm'] + inhib.vec.Cm = params['Cm'] + inhib.vec.vReset = params['inhVreset'] + inhib.vec.refractoryPeriod = params['inhibRefractTime'] + output.vec.thresh = params['outputThresh'] + output.vec.Rm = params['Rm'] + output.vec.Cm = params['Cm'] + output.vec.refractoryPeriod = params['outputRefractTime'] + otab = moose.Table( '/model/otab', params['numOutput'] ) + moose.connect( otab, 'requestOut', output, 'getVm', 'OneToOne' ) + itab = moose.Table( '/model/itab', params['numInhib'] ) + moose.connect( itab, 'requestOut', inhib, 'getVm', 'OneToOne' ) + return inhib, output + +def makeRaster( network ): + raster = moose.Table(network.path + '/raster', network.numData) + moose.connect(network, 'spikeOut', raster, 'spike', 'OneToOne') + return raster + +def main(): + mod = moose.Neutral( '/model' ) + inhib, output = makeGlobalBalanceNetwork() + iraster = makeRaster( inhib ) + oraster = makeRaster( output ) + moose.reinit() + moose.start( params['runtime'] ) + for i in range( len(iraster.vec) ): + print( iraster.vec[i].vector ) + +if __name__ == '__main__': + main() diff --git a/moose-core/tests/python/test_docs.py b/moose-core/tests/python/test_docs.py new file mode 100644 index 0000000000000000000000000000000000000000..cd877d9941b99e04cb78911bbc89925fbacbe6e8 --- /dev/null +++ b/moose-core/tests/python/test_docs.py @@ -0,0 +1,25 @@ +# -*- coding: utf-8 -*- +"""test_docs.py: + +Test if moose.doc is working. + +""" + +__author__ = "Dilawar Singh" +__copyright__ = "Copyright 2017-, Dilawar Singh" +__version__ = "1.0.0" +__maintainer__ = "Dilawar Singh" +__email__ = "dilawars@ncbs.res.in" +__status__ = "Development" + +import sys +import os +import moose +print( '[INFO] Using moose from %s' % moose.__file__ ) +print( '[INFO] Version : %s' % moose.version( ) ) + +def main( ): + moose.doc( 'Compartment' ) + +if __name__ == '__main__': + main() diff --git a/moose-core/tests/python/test_import.py b/moose-core/tests/python/test_import.py new file mode 100644 index 0000000000000000000000000000000000000000..ec0913fb0e1f839af0c2fe9ea181dbade2579d3c --- /dev/null +++ b/moose-core/tests/python/test_import.py @@ -0,0 +1,8 @@ +# -*- coding: utf-8 -*- +# Script to test all import modules. + +import moose +from moose.genesis import * +from moose.SBML import * +import moose.chemMerge +import moose.utils as mu diff --git a/moose-core/tests/python/test_moogli.py b/moose-core/tests/python/test_moogli.py deleted file mode 100644 index 8533a395748a6a840aa4d4606b1ab35de776dd7a..0000000000000000000000000000000000000000 --- a/moose-core/tests/python/test_moogli.py +++ /dev/null @@ -1,17 +0,0 @@ -# -*- coding: utf-8 -*- -"""test_moogli.py: - - Some sanity tests on moogli - -""" - -__author__ = "Dilawar Singh" -__copyright__ = "Copyright 2015, Dilawar Singh and NCBS Bangalore" -__credits__ = ["NCBS Bangalore"] -__license__ = "GNU GPL" -__version__ = "1.0.0" -__maintainer__ = "Dilawar Singh" -__email__ = "dilawars@ncbs.res.in" -__status__ = "Development" - -import moogli diff --git a/moose-core/tests/python/test_mumbl.py b/moose-core/tests/python/test_mumbl.py deleted file mode 100644 index 702d0c2386cf4a114335aa2d13f17aab11609664..0000000000000000000000000000000000000000 --- a/moose-core/tests/python/test_mumbl.py +++ /dev/null @@ -1,64 +0,0 @@ -# -*- coding: utf-8 -*- -# This file is part of MOOSE simulator: http://moose.ncbs.res.in. - -# MOOSE is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. - -# MOOSE is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# You should have received a copy of the GNU General Public License -# along with MOOSE. If not, see <http://www.gnu.org/licenses/>. - - -"""test_mumbl.py: - - A test script to test MUMBL support in MOOSE. - -Last modified: Fri Jun 13, 2014 06:30PM - -""" - -__author__ = "Dilawar Singh" -__copyright__ = "Copyright 2013, NCBS Bangalore" -__credits__ = ["NCBS Bangalore", "Bhalla Lab"] -__license__ = "GNU GPL" -__version__ = "1.0.0" -__maintainer__ = "Dilawar Singh" -__email__ = "dilawars@ncbs.res.in" -__status__ = "Development" - -import sys -import os - -import moose -import moose.neuroml as nml -import moose.utils as utils -import moose.backend.graphviz as graphviz - -# the model lives in the same directory as the test script -modeldir = os.path.dirname(__file__) - -def main(): - utils.parser - p = os.path.join(modeldir, 'two_cells_nml_1.8/two_cells.nml') - nml.loadNeuroML_L123(p) - #mumbl.loadMumbl("./two_cells_nml_1.8/mumbl.xml") - table1 = utils.recordTarget('/tableA', '/cells/purkinjeGroup_0/Dend_37_41', 'vm') - table2 = utils.recordTarget('/tableB', '/cells/granuleGroup_0/Soma_0', 'vm') - moose.setClock(0, 5e-6) - moose.useClock(0, '/##', 'process') - moose.useClock(0, '/##', 'init') - moose.reinit() - utils.run(0.1, verify=True) - graphviz.writeGraphviz('test_mumble.dot', ignore='/library') - utils.plotRecords({ 'Dend 37' : table1, 'Soma 0' : table2 } - , outfile = '%s.png' % sys.argv[0] - , subplot = True - ) - -if __name__ == '__main__': - main() diff --git a/moose-core/tests/python/test_nsdf.py b/moose-core/tests/python/test_nsdf.py index 3a4da9fcc449d159c5376eb507642420ba0ec8bf..bdec9641e6dab406d585246dc353754c84acf256 100644 --- a/moose-core/tests/python/test_nsdf.py +++ b/moose-core/tests/python/test_nsdf.py @@ -90,4 +90,4 @@ if __name__ == '__main__': raise Exception("Test failed. No file : %s" % nsdfFile) data = np.loadtxt( nsdfFile ) - assert len(data) == 60001, "Expected 60001 entries" + assert len(data) == 4, "Expected 4 entries" diff --git a/moose-core/tests/python/test_utils.py b/moose-core/tests/python/test_utils.py new file mode 100644 index 0000000000000000000000000000000000000000..a57633b84eb83766bbe848312947c472b0b7ee5a --- /dev/null +++ b/moose-core/tests/python/test_utils.py @@ -0,0 +1,10 @@ +# -*- coding: utf-8 -*- +import moose +import moose.utils as mu + +def main( ): + print( dir( mu ) ) + mu.info( 'Hellow' ) + +if __name__ == '__main__': + main() diff --git a/moose-core/tests/python/test_warnigns.py b/moose-core/tests/python/test_warnigns.py deleted file mode 100644 index a68ed5e83614d183ddaf7ea176d6ac9952e0bd6a..0000000000000000000000000000000000000000 --- a/moose-core/tests/python/test_warnigns.py +++ /dev/null @@ -1,14 +0,0 @@ -# -*- coding: utf-8 -*- -import moose -print( 'Using moose from %s' % moose.__file__ ) -a = moose.Neutral( '/a' ) -b = moose.Ksolve( '/a/k' ) - -try: - c = moose.Ksolve( '/a/k' ) - d = moose.Neutral( '/a' ) - d = moose.Neutral( '/a' ) - quit(-1) -except Exception as e: - c = moose.element( '/a/k' ) - quit( 0 )