From 7b3c355a04d22f279a88bfae4097e9cf3f62439b Mon Sep 17 00:00:00 2001
From: Dilawar Singh <dilawars@ncbs.res.in>
Date: Tue, 31 Oct 2017 18:39:15 +0530
Subject: [PATCH] Squashed 'moose-core/' changes from 2657e3345..268a28862

268a28862 Merge pull request #224 from upibhalla/master
3ab452e97 Fixed test_nsdf.py. It had relied on the old (and incorrect) function of the spike message to the Table. I have now updated the test_nsdf.py to work with the fixed version of the Table.
616f33a60 Fixed spike timing functionality in Table class. Incorporated into rdesigneur so that tables can record spike times for saving, and users can plot spike times as a dot raster.
7f372b6aa Merge pull request #223 from upibhalla/master
186e91d68 Merge branch 'master' into master
d91a681bf Fix on diffusion scaling in spines, and added Pool::volume field for plotting in rdesigneur.
41c71b147 Merge pull request #221 from upibhalla/master
e27ef8a45 Merge branch 'master' of https://github.com/BhallaLab/moose-core
b07046e15 New Neuron::getSpineIdsFromCompartmentIds function, needed for rdesineur to connect adaptors for spine size and other spine properties.
69d985877 Merge pull request #220 from bhanu151/master
4f67f3338 retracting seqActivation normalization by dt to restore time-step invariance
ca1f36ca6 Merge pull request #219 from hrani/master
be6dc0a28 Hopefully last cleanup
e03c69ddd Hopefully last cleanup
7332ff66b spell correction
d17d383dc absolute import with mtypes just for python3
fb28e1009 Hopefully last one, tested for filepath,moose path and moose object
b782b9552 In merge.py, 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
dc0c83027 clean up with in mergeChemModel file
13152346b cleanup in setup file
1f9100b7f cleanup in init file
50cd8feee removed unwanted lines and clean up
e9b88cee7 folder "merge" renamed to "chemMerge", moose.py: func added so that chemMerge can be imported directly,cleanup in test_import.py with the changes made to this commit,merge.py file Neutral object is copied into des path instead of creating and now user have to genesis/SBML format
072913a9c Merge pull request #218 from BhallaLab/missing_merge
fab362fa7 Couple of lines added to moose import test.
064fbea56 Fixed moose.merge import. Added a test which tests moose module import (incomplete).
4211da0e8 Merge pull request #216 from hrani/master
a45e9f0a1 Merge branch 'master' of http://github.com/hrani/moose-core
dacd64039 some more checks made for empty model, empty compartment and for display atleast one reaction should exist for displaying into widget
e0566d476 Merge branch 'master' of https://github.com/BhallaLab/moose-core
2db1b4dac Merge branch 'master' into master
a62961a27 Merge pull request #217 from upibhalla/master
15ee2cc7f Merge branch 'master' of https://github.com/BhallaLab/moose-core
aca164612 Two important fixes for numerics of stochastic solver with diffusion. Affects models using sumtotals and functions.
f5302b08f validator was switched off, which is not correct, so switching ON
87e66a18d if in reaction stoichiometry is not defined then assume 1 and if no unit base unit defined then assume 1 which is mole
38879fd5b Merge pull request #215 from hrani/master
2384657e0 readSBML: explicitly ratio is set to 4, so that e2,e1 is balanced, If units are defined in the rate law for the reaction then checks are made and if not in milli mole the base unit then recalculate, spell correction. writeSBML: instead of printing error msg, saving into string for later use writekkit : instead of printing error msg, saving into string for later use
798c0b6a7 Merge pull request #213 from hrani/master
c1c26f698 Merge branch 'master' into master
aa1cb43e2 Merge pull request #214 from bhanu151/master
b60b73b3a SeqSynHnadler fixes - normalized seqActivation by dt
3d36e7bd1 unwanted print statement is removed
67e026ab8 function definition are read from SBML files, rules which are not read to moose a errorFlag is set out, rateLaw without unit definition, are taken from the substance units and recalculated for all higher reaction, validator checks if model is valid and send out errormsg if not valid, if model doesn't have atleast one compartment, then errormsg is send out
91180c067 Merge pull request #210 from dilawar/master
1f37c10f0 Python3/Python2 import thread has been fixed.
012b5f3ca Preserve the order of plots. Use OrderedDict. Since OrderedDict is not available in python2.6; added implementation.
78f244b0c Utility function to plot tables matching a regular expression.
03dc5a8fc Fixed import error.
2b9e30565 try/except import in moose.py file.
574a12503 Added a test script.
f5202b0cd Merge branch 'master' of github.com:BhallaLab/moose-core
b792838ee Merge pull request #207 from hrani/master
e39c76df2 resolved conflicts
75e289296 all conflict a removed
1e7ebbe3c unwanted statements
98d625b62 Merge branch 'master' into master
cdbe6add9 Fixed clock dt assignment bug in rdesigneur, came up when turnOffElec flag is true (#208)
0da306f5a Fixed clock dt assignment bug in rdesigneur, came up when turnOffElec flag is true
e9526f7c3 Fixes to import plot_util library.
abd664d1e unwanted print statement removed
fd794de74 in print stm space clean up,added findcompartment function to one common place, some space char are allowed in name field, in kkit writer, moose obj which doesn't have compartment are not written, which then appended to error msg for user
5015b14cd Merge branch 'master' of github.com:BhallaLab/moose-core
b32a0bbe9 Fixes to issue204 (#205)
9f88d7c5f Merge pull request #203 from hrani/master
6dbf8bef2 Merge branch 'master' of http://github.com/hrani/moose-core
3c6d41ce4 Restoring contents after bad commit and bad merge
1d34b9b2c Restored some old commit, after the bad commit and after revert fail
f2d3e94ad Restored some old commit, after the bad commit and after revert fail
daa657802 Restoring contents after bad commit and bad merge
30b1c36b2 Merge pull request #202 from upibhalla/master
24c079561 Restoring updates to ksolve including isBuffered flag
d71f6eb53 Merge branch 'master' of https://github.com/BhallaLab/moose-core
eb43e5524 Fixed Rdesigneur verbose flag
5e42df95c Merge branch 'master' of github.com:BhallaLab/moose-core
3c46150a1 Merge branch 'master' of github.com:dilawar/moose-core
b04156b03 Fixes to issue #200. (#201)
1e3bbe04f Merge branch 'master' of github.com:BhallaLab/moose-core
f11ef0a3f Added test script but has not been enabled. This tests the SparseMsg::connectionLists.

git-subtree-dir: moose-core
git-subtree-split: 268a2886219290af163c9c95b755901fa68cfa5d
---
 .travis/travis_prepare_osx.sh                 |   4 +-
 CMakeLists.txt                                |  33 +-
 CheckCXXCompiler.cmake                        |  17 +-
 MooseTests.cmake                              | 154 ----
 biophysics/Neuron.cpp                         |  50 +-
 biophysics/Neuron.h                           |   2 +
 biophysics/Spine.cpp                          |   2 +-
 builtins/Table.cpp                            |  58 +-
 builtins/Table.h                              |   5 +
 cmake_modules/FindGSL.cmake                   |  49 +-
 kinetics/BufPool.cpp                          |   2 +-
 kinetics/Pool.cpp                             |  86 +-
 kinetics/Pool.h                               |  10 +-
 kinetics/PoolBase.cpp                         |  75 ++
 kinetics/PoolBase.h                           |  14 +
 ksolve/GssaVoxelPools.cpp                     |  20 +-
 ksolve/KinSparseMatrix.cpp                    |  15 -
 ksolve/KinSparseMatrix.h                      |   4 -
 ksolve/Stoich.cpp                             |  15 +-
 ksolve/Stoich.h                               |   9 +
 ksolve/VoxelPoolsBase.cpp                     |   7 +-
 ksolve/ZombieBufPool.cpp                      |   4 +
 ksolve/ZombieBufPool.h                        |   1 +
 ksolve/ZombiePool.cpp                         |   5 +
 ksolve/ZombiePool.h                           |   1 +
 python/moose/OrderedDict.py                   | 261 ++++++
 python/moose/SBML/readSBML.py                 | 750 ++++++++++++------
 python/moose/SBML/validation.py               |  22 +-
 python/moose/SBML/writeSBML.py                | 328 ++++++--
 python/moose/chemMerge/__init__.py            |   2 +
 python/moose/{merge => chemMerge}/merge.py    | 517 +++++++-----
 python/moose/{merge => chemMerge}/mtypes.py   |   0
 python/moose/chemUtil/chemConnectUtil.py      |  34 +-
 python/moose/chemUtil/graphUtils.py           |  24 +-
 python/moose/genesis/writeKkit.py             | 557 +++++++------
 python/moose/moose.py                         |  84 +-
 python/moose/plot_utils.py                    |  38 +-
 python/moose/print_utils.py                   |  24 +-
 python/moose/utils.py                         |  30 +-
 python/rdesigneur/rdesigneur.py               |  99 ++-
 python/setup.py                               |   2 +-
 .../{test_gc.py => deprecated_test_gc.py}     |   0
 tests/python/test_connectionLists.py          | 158 ++++
 tests/python/test_docs.py                     |  25 +
 tests/python/test_import.py                   |   8 +
 tests/python/test_moogli.py                   |  17 -
 tests/python/test_mumbl.py                    |  64 --
 tests/python/test_nsdf.py                     |   2 +-
 tests/python/test_utils.py                    |  10 +
 tests/python/test_warnigns.py                 |  14 -
 50 files changed, 2487 insertions(+), 1225 deletions(-)
 delete mode 100644 MooseTests.cmake
 create mode 100644 python/moose/OrderedDict.py
 create mode 100644 python/moose/chemMerge/__init__.py
 rename python/moose/{merge => chemMerge}/merge.py (55%)
 rename python/moose/{merge => chemMerge}/mtypes.py (100%)
 rename tests/python/{test_gc.py => deprecated_test_gc.py} (100%)
 create mode 100644 tests/python/test_connectionLists.py
 create mode 100644 tests/python/test_docs.py
 create mode 100644 tests/python/test_import.py
 delete mode 100644 tests/python/test_moogli.py
 delete mode 100644 tests/python/test_mumbl.py
 create mode 100644 tests/python/test_utils.py
 delete mode 100644 tests/python/test_warnigns.py

diff --git a/.travis/travis_prepare_osx.sh b/.travis/travis_prepare_osx.sh
index 8c67277a..61fdea5b 100755
--- a/.travis/travis_prepare_osx.sh
+++ b/.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/CMakeLists.txt b/CMakeLists.txt
index fa8df758..14bedd86 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -150,11 +150,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"
@@ -450,8 +450,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/CheckCXXCompiler.cmake b/CheckCXXCompiler.cmake
index 15a4d730..92fbfdfa 100644
--- a/CheckCXXCompiler.cmake
+++ b/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/MooseTests.cmake b/MooseTests.cmake
deleted file mode 100644
index 6026946d..00000000
--- a/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/biophysics/Neuron.cpp b/biophysics/Neuron.cpp
index 4ca30e10..4ae85145 100644
--- a/biophysics/Neuron.cpp
+++ b/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/biophysics/Neuron.h b/biophysics/Neuron.h
index 3cef4dab..777c5448 100644
--- a/biophysics/Neuron.h
+++ b/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/biophysics/Spine.cpp b/biophysics/Spine.cpp
index 6909b545..03ca6eac 100644
--- a/biophysics/Spine.cpp
+++ b/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/builtins/Table.cpp b/builtins/Table.cpp
index e07e4172..c82c6a8d 100644
--- a/builtins/Table.cpp
+++ b/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/builtins/Table.h b/builtins/Table.h
index 5a21b2d3..9fb78c8b 100644
--- a/builtins/Table.h
+++ b/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/cmake_modules/FindGSL.cmake b/cmake_modules/FindGSL.cmake
index 93c2686f..b45aca3a 100644
--- a/cmake_modules/FindGSL.cmake
+++ b/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/kinetics/BufPool.cpp b/kinetics/BufPool.cpp
index 230be5d6..54e18378 100644
--- a/kinetics/BufPool.cpp
+++ b/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/kinetics/Pool.cpp b/kinetics/Pool.cpp
index 86508cf4..8b0eb629 100644
--- a/kinetics/Pool.cpp
+++ b/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/kinetics/Pool.h b/kinetics/Pool.h
index d55696aa..bc4e8877 100644
--- a/kinetics/Pool.h
+++ b/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/kinetics/PoolBase.cpp b/kinetics/PoolBase.cpp
index b67e70cf..f1d7e0f4 100644
--- a/kinetics/PoolBase.cpp
+++ b/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/kinetics/PoolBase.h b/kinetics/PoolBase.h
index 3821569a..b5aaa1b8 100644
--- a/kinetics/PoolBase.h
+++ b/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/ksolve/GssaVoxelPools.cpp b/ksolve/GssaVoxelPools.cpp
index 0b09a418..438652d8 100644
--- a/ksolve/GssaVoxelPools.cpp
+++ b/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/ksolve/KinSparseMatrix.cpp b/ksolve/KinSparseMatrix.cpp
index 59783b29..1b1c1d31 100644
--- a/ksolve/KinSparseMatrix.cpp
+++ b/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/ksolve/KinSparseMatrix.h b/ksolve/KinSparseMatrix.h
index 4cc48aec..52132dcc 100644
--- a/ksolve/KinSparseMatrix.h
+++ b/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/ksolve/Stoich.cpp b/ksolve/Stoich.cpp
index 425c19dd..198f1978 100644
--- a/ksolve/Stoich.cpp
+++ b/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/ksolve/Stoich.h b/ksolve/Stoich.h
index 2951e365..6dbe51de 100644
--- a/ksolve/Stoich.h
+++ b/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/ksolve/VoxelPoolsBase.cpp b/ksolve/VoxelPoolsBase.cpp
index b9e4e306..22831c66 100644
--- a/ksolve/VoxelPoolsBase.cpp
+++ b/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/ksolve/ZombieBufPool.cpp b/ksolve/ZombieBufPool.cpp
index a00a6499..82665529 100644
--- a/ksolve/ZombieBufPool.cpp
+++ b/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/ksolve/ZombieBufPool.h b/ksolve/ZombieBufPool.h
index 8fda4df1..a00344e3 100644
--- a/ksolve/ZombieBufPool.h
+++ b/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/ksolve/ZombiePool.cpp b/ksolve/ZombiePool.cpp
index b87beb6e..6d612c83 100644
--- a/ksolve/ZombiePool.cpp
+++ b/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/ksolve/ZombiePool.h b/ksolve/ZombiePool.h
index d0cefe8f..147823fc 100644
--- a/ksolve/ZombiePool.h
+++ b/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/python/moose/OrderedDict.py b/python/moose/OrderedDict.py
new file mode 100644
index 00000000..cf7e1716
--- /dev/null
+++ b/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/python/moose/SBML/readSBML.py b/python/moose/SBML/readSBML.py
index aa321652..c5b00c52 100644
--- a/python/moose/SBML/readSBML.py
+++ b/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": "&#176", "-": "_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/python/moose/SBML/validation.py b/python/moose/SBML/validation.py
index 8b656922..ee2b2cae 100644
--- a/python/moose/SBML/validation.py
+++ b/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/python/moose/SBML/writeSBML.py b/python/moose/SBML/writeSBML.py
index 82ce209f..99fb3f21 100644
--- a/python/moose/SBML/writeSBML.py
+++ b/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": "&#176"}
     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"     : "&#176", 
+            "'"     : "_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": "&#176", "-": "_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/python/moose/chemMerge/__init__.py b/python/moose/chemMerge/__init__.py
new file mode 100644
index 00000000..0116e5db
--- /dev/null
+++ b/python/moose/chemMerge/__init__.py
@@ -0,0 +1,2 @@
+# -*- coding: utf-8 -*-
+from .merge import *
diff --git a/python/moose/merge/merge.py b/python/moose/chemMerge/merge.py
similarity index 55%
rename from python/moose/merge/merge.py
rename to python/moose/chemMerge/merge.py
index cb1ab62d..5b406a25 100644
--- a/python/moose/merge/merge.py
+++ b/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/python/moose/merge/mtypes.py b/python/moose/chemMerge/mtypes.py
similarity index 100%
rename from python/moose/merge/mtypes.py
rename to python/moose/chemMerge/mtypes.py
diff --git a/python/moose/chemUtil/chemConnectUtil.py b/python/moose/chemUtil/chemConnectUtil.py
index be0c4ba6..48e06e94 100644
--- a/python/moose/chemUtil/chemConnectUtil.py
+++ b/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/python/moose/chemUtil/graphUtils.py b/python/moose/chemUtil/graphUtils.py
index 061b8ea0..2ca109f3 100644
--- a/python/moose/chemUtil/graphUtils.py
+++ b/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/python/moose/genesis/writeKkit.py b/python/moose/genesis/writeKkit.py
index 2e25bc7e..c94ca3da 100644
--- a/python/moose/genesis/writeKkit.py
+++ b/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/python/moose/moose.py b/python/moose/moose.py
index 3fa32f5c..7d045d37 100644
--- a/python/moose/moose.py
+++ b/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/python/moose/plot_utils.py b/python/moose/plot_utils.py
index d3f10c92..9b309e46 100644
--- a/python/moose/plot_utils.py
+++ b/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/python/moose/print_utils.py b/python/moose/print_utils.py
index b3acd81b..5cc08416 100644
--- a/python/moose/print_utils.py
+++ b/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/python/moose/utils.py b/python/moose/utils.py
index 5a288cc0..d4fcf16e 100644
--- a/python/moose/utils.py
+++ b/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/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py
index e5704f2a..4220f376 100644
--- a/python/rdesigneur/rdesigneur.py
+++ b/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/python/setup.py b/python/setup.py
index f8d537d0..93c432cd 100644
--- a/python/setup.py
+++ b/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/tests/python/test_gc.py b/tests/python/deprecated_test_gc.py
similarity index 100%
rename from tests/python/test_gc.py
rename to tests/python/deprecated_test_gc.py
diff --git a/tests/python/test_connectionLists.py b/tests/python/test_connectionLists.py
new file mode 100644
index 00000000..057cd5a4
--- /dev/null
+++ b/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/tests/python/test_docs.py b/tests/python/test_docs.py
new file mode 100644
index 00000000..cd877d99
--- /dev/null
+++ b/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/tests/python/test_import.py b/tests/python/test_import.py
new file mode 100644
index 00000000..ec0913fb
--- /dev/null
+++ b/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/tests/python/test_moogli.py b/tests/python/test_moogli.py
deleted file mode 100644
index 8533a395..00000000
--- a/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/tests/python/test_mumbl.py b/tests/python/test_mumbl.py
deleted file mode 100644
index 702d0c23..00000000
--- a/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/tests/python/test_nsdf.py b/tests/python/test_nsdf.py
index 3a4da9fc..bdec9641 100644
--- a/tests/python/test_nsdf.py
+++ b/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/tests/python/test_utils.py b/tests/python/test_utils.py
new file mode 100644
index 00000000..a57633b8
--- /dev/null
+++ b/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/tests/python/test_warnigns.py b/tests/python/test_warnigns.py
deleted file mode 100644
index a68ed5e8..00000000
--- a/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 )
-- 
GitLab