BhallaLab · Sep 18, 2018
diff --git a/‎moose-core/CMakeLists.txt
+1 b/‎moose-core/CMakeLists.txt
+1
diff --git a/‎moose-core/biophysics/Neuron.cpp
+21-9 b/‎moose-core/biophysics/Neuron.cpp
+21-9
diff --git a/‎moose-core/ksolve/Ksolve.cpp
+6-10 b/‎moose-core/ksolve/Ksolve.cpp
+6-10
diff --git a/‎moose-core/ksolve/OdeSystem.h
+5-4 b/‎moose-core/ksolve/OdeSystem.h
+5-4
diff --git a/‎moose-core/ksolve/SteadyStateBoost.cpp
+48-73 b/‎moose-core/ksolve/SteadyStateBoost.cpp
+48-73
@@ -130,6 +130,7 @@ include_directories( ${CMAKE_CURRENT_SOURCE_DIR} )
 
 if(WITH_BOOST)
     set(WITH_BOOST_ODE ON)
+    set(WITH_GSL       OFF)
 endif(WITH_BOOST)
 
 # If using BOOST ODE2 library to solve ODE system, then don't use GSL.
 
@@ -315,11 +315,17 @@ const Cinfo* Neuron::initCinfo()
         "The expression for the *spacing* field must evaluate to > 0 for "
         "the spine to be installed. For example, if the expresssion is\n"
         "		H(1 - L) \n"
-        "then the systemwill only put spines closer than "
+        "then the system will only put spines closer than "
         "one length constant from the soma, and zero elsewhere. \n"
         "Available spine parameters are: \n"
         "spacing, minSpacing, size, sizeDistrib "
-        "angle, angleDistrib \n",
+        "angle, angleDistrib \n"
+		"minSpacing sets the granularity of sampling (typically about 0.1*"
+	    "spacing) for the usual case where spines are spaced randomly. "
+		"If minSpacing < 0 then the spines are spaced equally at "
+		"'spacing', unless the dendritic segment length is smaller than "
+		"'spacing'. In that case it falls back to the regular random "
+		"placement method.",
         &Neuron::setSpineDistribution,
         &Neuron::getSpineDistribution
     );
@@ -1788,6 +1794,19 @@ static void addPos( unsigned int segIndex, unsigned int eIndex,
                     vector< unsigned int >& elistIndex,
                     vector< double >& pos )
 {
+	if ( minSpacing < 0.0 ) {
+		// Use uniform spacing
+		for ( double position = spacing * 0.5;
+				position < dendLength; position += spacing ) {
+			seglistIndex.push_back( segIndex );
+			elistIndex.push_back( eIndex );
+			pos.push_back( position );
+		}
+		if ( dendLength > spacing * 0.5 )
+			return;
+		// If the dend length is too small for regular placement, 
+		// fall back to using probability to decide if segment gets spine
+	}
     if ( minSpacing < spacing * 0.1 && minSpacing < 1e-7 )
         minSpacing = spacing * 0.1;
     if ( minSpacing > spacing * 0.5 )
@@ -1860,13 +1879,6 @@ void Neuron::makeSpacingDistrib( const vector< ObjId >& elist,
             {
                 double spacing = val[ j + nuParser::EXPR ];
                 double spacingDistrib = parser.eval( val.begin() + j );
-                if ( spacingDistrib > spacing || spacingDistrib < 0 )
-                {
-                    cout << "Warning: Neuron::makeSpacingDistrib: " <<
-                         "0 < " << spacingDistrib << " < " << spacing <<
-                         " fails on " << elist[i].path() << ". Using 0.\n";
-                    spacingDistrib = 0.0;
-                }
                 map< Id, unsigned int>::const_iterator
                 lookupDend = segIndex_.find( elist[i] );
                 if ( lookupDend != segIndex_.end() )
 
@@ -235,7 +235,7 @@ static const Cinfo* ksolveCinfo = Ksolve::initCinfo();
 
 Ksolve::Ksolve()
     :
-#if USE_GSL
+#ifdef USE_GSL
     method_( "rk5" ),
 #elif USE_BOOST_ODE
     method_( "rk5a" ),
@@ -248,19 +248,12 @@ Ksolve::Ksolve()
     dsolve_(),
     dsolvePtr_( 0 )
 {
+    ;
 }
 
 Ksolve::~Ksolve()
 {
-#if 0
-    char* p = getenv( "MOOSE_SHOW_SOLVER_PERF" );
-    if( p != NULL )
-    {
-        cout << "Info: Ksolve (+Dsolve) took " << totalTime_ << " seconds and took " << numSteps_
-             << " steps." << endl;
-
-    }
-#endif
+    ;
 }
 
 //////////////////////////////////////////////////////////////
@@ -378,6 +371,7 @@ void Ksolve::setStoich( Id stoich )
     assert( stoich.element()->cinfo()->isA( "Stoich" ) );
     stoich_ = stoich;
     stoichPtr_ = reinterpret_cast< Stoich* >( stoich.eref().data() );
+
     if ( !isBuilt_ )
     {
         OdeSystem ode;
@@ -386,6 +380,7 @@ void Ksolve::setStoich( Id stoich )
         // ode.initStepSize = getEstimatedDt();
         ode.initStepSize = 0.01; // This will be overridden at reinit.
         ode.method = method_;
+
 #ifdef USE_GSL
         ode.gslSys.dimension = stoichPtr_->getNumAllPools();
         if ( ode.gslSys.dimension == 0 )
@@ -416,6 +411,7 @@ void Ksolve::setStoich( Id stoich )
 #endif
         isBuilt_ = true;
     }
+
 }
 
 Id Ksolve::getDsolve() const
 
@@ -26,15 +26,16 @@ class OdeSystem {
     {;}
 
         std::string method;
-        // GSL stuff
+
+        double initStepSize;
+        double epsAbs; // Absolute error
+        double epsRel; // Relative error
 
 #ifdef USE_GSL
+        // GSL stuff
         gsl_odeiv2_system gslSys;
         const gsl_odeiv2_step_type* gslStep;
 #endif
-        double initStepSize;
-        double epsAbs; // Absolute error
-        double epsRel; // Relative error
 
 #if  USE_BOOST_ODE
         size_t dimension;
 
@@ -85,27 +85,6 @@ struct reac_info
 
 const Cinfo* SteadyState::initCinfo()
 {
-    /**
-     * This picks up the entire Stoich data structure
-    static Finfo* gslShared[] =
-    {
-    	new SrcFinfo( "reinitSrc", Ftype0() ),
-    	new DestFinfo( "assignStoich",
-    		Ftype1< void* >(),
-    		RFCAST( &SteadyState::assignStoichFunc )
-    		),
-    	new DestFinfo( "setMolN",
-    		Ftype2< double, unsigned int >(),
-    		RFCAST( &SteadyState::setMolN )
-    		),
-    	new SrcFinfo( "requestYsrc", Ftype0() ),
-    	new DestFinfo( "assignY",
-    		Ftype1< double* >(),
-    		RFCAST( &SteadyState::assignY )
-    		),
-    };
-     */
-
     /**
      * These are the fields of the SteadyState class
      */
@@ -207,59 +186,59 @@ const Cinfo* SteadyState::initCinfo()
     // MsgDest definitions
     ///////////////////////////////////////////////////////
     static DestFinfo setupMatrix( "setupMatrix",
-                                  "This function initializes and rebuilds the matrices used "
-                                  "in the calculation.",
-                                  new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
-                                );
+            "This function initializes and rebuilds the matrices used "
+            "in the calculation.",
+            new OpFunc0< SteadyState >(&SteadyState::setupMatrix)
+            );
 
     static DestFinfo settle( "settle",
-                             "Finds the nearest steady state to the current initial "
-                             "conditions. This function rebuilds the entire calculation "
-                             "only if the object has not yet been initialized.",
-                             new OpFunc0< SteadyState >( &SteadyState::settleFunc )
-                           );
+            "Finds the nearest steady state to the current initial "
+            "conditions. This function rebuilds the entire calculation "
+            "only if the object has not yet been initialized.",
+            new OpFunc0< SteadyState >( &SteadyState::settleFunc )
+            );
     static DestFinfo resettle( "resettle",
-                               "Finds the nearest steady state to the current initial "
-                               "conditions. This function rebuilds the entire calculation ",
-                               new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
-                             );
+            "Finds the nearest steady state to the current initial "
+            "conditions. This function rebuilds the entire calculation ",
+            new OpFunc0< SteadyState >( &SteadyState::resettleFunc )
+            );
     static DestFinfo showMatrices( "showMatrices",
-                                   "Utility function to show the matrices derived for the calculations on the reaction system. Shows the Nr, gamma, and total matrices",
-                                   new OpFunc0< SteadyState >( &SteadyState::showMatrices )
-                                 );
+            "Utility function to show the matrices derived for the calculations on the reaction system. Shows the Nr, gamma, and total matrices",
+            new OpFunc0< SteadyState >( &SteadyState::showMatrices )
+            );
     static DestFinfo randomInit( "randomInit",
-                                 "Generate random initial conditions consistent with the mass"
-                                 "conservation rules. Typically invoked in order to scan"
-                                 "states",
-                                 new EpFunc0< SteadyState >(
-                                     &SteadyState::randomizeInitialCondition )
-                               );
+            "Generate random initial conditions consistent with the mass"
+            "conservation rules. Typically invoked in order to scan"
+            "states",
+            new EpFunc0< SteadyState >(
+                &SteadyState::randomizeInitialCondition )
+            );
     ///////////////////////////////////////////////////////
     // Shared definitions
     ///////////////////////////////////////////////////////
 
     static Finfo * steadyStateFinfos[] =
     {
-        &stoich,				// Value
-        &badStoichiometry,		// ReadOnlyValue
-        &isInitialized,			// ReadOnlyValue
-        &nIter,					// ReadOnlyValue
-        &status,				// ReadOnlyValue
-        &maxIter,				// Value
-        &convergenceCriterion,	// ReadOnlyValue
-        &numVarPools,			// ReadOnlyValue
-        &rank,					// ReadOnlyValue
-        &stateType,				// ReadOnlyValue
-        &nNegEigenvalues,		// ReadOnlyValue
-        &nPosEigenvalues,		// ReadOnlyValue
-        &solutionStatus,		// ReadOnlyValue
-        &total,					// LookupValue
-        &eigenvalues,			// ReadOnlyLookupValue
-        &setupMatrix,			// DestFinfo
-        &settle,				// DestFinfo
-        &resettle,				// DestFinfo
-        &showMatrices,			// DestFinfo
-        &randomInit,			// DestFinfo
+        &stoich,                // Value
+        &badStoichiometry,        // ReadOnlyValue
+        &isInitialized,            // ReadOnlyValue
+        &nIter,                    // ReadOnlyValue
+        &status,                // ReadOnlyValue
+        &maxIter,                // Value
+        &convergenceCriterion,    // ReadOnlyValue
+        &numVarPools,            // ReadOnlyValue
+        &rank,                    // ReadOnlyValue
+        &stateType,                // ReadOnlyValue
+        &nNegEigenvalues,        // ReadOnlyValue
+        &nPosEigenvalues,        // ReadOnlyValue
+        &solutionStatus,        // ReadOnlyValue
+        &total,                    // LookupValue
+        &eigenvalues,            // ReadOnlyLookupValue
+        &setupMatrix,            // DestFinfo
+        &settle,                // DestFinfo
+        &resettle,                // DestFinfo
+        &showMatrices,            // DestFinfo
+        &randomInit,            // DestFinfo
 
 
     };
@@ -278,7 +257,7 @@ const Cinfo* SteadyState::initCinfo()
         "Note that the method finds unstable as well as stable fixed "
         "points.\n "
         "The SteadyState class also provides a utility function "
-        "*randomInit()*	to "
+        "*randomInit()*    to "
         "randomly initialize the concentrations, within the constraints "
         "of stoichiometry. This is useful if you are trying to find "
         "the major fixed points of the system. Note that this is "
@@ -368,12 +347,10 @@ void SteadyState::setStoich( Id value )
     numVarPools_ = Field< unsigned int >::get( stoich_, "numVarPools" );
     nReacs_ = Field< unsigned int >::get( stoich_, "numRates" );
     setupSSmatrix();
-    double vol = LookupField< unsigned int, double >::get(
-                     stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
+    double vol = LookupField< unsigned int, double >::get(stoichPtr->getCompartment(), "oneVoxelVolume", 0 );
     pool_.setVolume( vol );
     pool_.setStoich( stoichPtr, 0 );
-    pool_.updateAllRateTerms( stoichPtr->getRateTerms(),
-                              stoichPtr->getNumCoreRates() );
+    pool_.updateAllRateTerms( stoichPtr->getRateTerms(), stoichPtr->getNumCoreRates() );
     isInitialized_ = 1;
 }
 
@@ -513,9 +490,9 @@ void SteadyState::showMatrices()
         return;
     }
     int numConsv = numVarPools_ - rank_;
-    cout << "Totals:	";
+    cout << "Totals:    ";
     for ( int i = 0; i < numConsv; ++i )
-        cout << total_[i] << "	";
+        cout << total_[i] << "    ";
     cout << endl;
     cout << "gamma " << gamma_ << endl;
     cout << "Nr " << Nr_ << endl;
@@ -549,9 +526,7 @@ void SteadyState::setupSSmatrix()
         {
             double x = 0;
             if ( j == colIndex[k] && k < rowStart[i+1] )
-            {
                 x = entry[k++];
-            }
             N(i,j) = x;
             LU_(i,j) = x;
         }
@@ -714,7 +689,7 @@ static bool isSolutionValid( const vector< double >& x )
     for( size_t i = 0; i < x.size(); i++ )
     {
         double v = x[i];
-        if ( std::isnan( v ) or std::isinf( v ) )
+        if ( std::isnan( v ) || std::isinf( v ) )
         {
             cout << "Warning: SteadyState iteration gave nan/inf concs\n";
             return false;