Merge pull request #1460 from brian-team/spatialneuron_without_scipy

Remove numpy-fallback for SpatialNeuron

brian2/codegen/generators/base.py

@@ -120,7 +120,7 @@ def determine_keywords(self):
         used for example by the `CPPCodeGenerator` to set up all the supporting
         code
         """
-        raise NotImplementedError
+        return {}
 
     def translate_one_statement_sequence(self, statements, scalar=False):
         raise NotImplementedError

brian2/codegen/generators/numpy_generator.py

@@ -330,16 +330,6 @@ def translate_one_statement_sequence(self, statements, scalar=False):
 
         return lines
 
-    def determine_keywords(self):
-        try:
-            import scipy  # noqa: F401
-
-            scipy_available = True
-        except ImportError:
-            scipy_available = False
-
-        return {"_scipy_available": scipy_available}
-
 
 ################################################################################
 # Implement functions

brian2/codegen/runtime/numpy_rt/templates/spatialstateupdate.py_

@@ -51,6 +51,10 @@ for _counter, (_first, _last) in enumerate(zip({{_starts}},
 {% block maincode %}
 from numpy import pi
 from numpy import zeros
+try:
+    from scipy.linalg import solve_banded
+except ImportError:
+    raise ImportError("Install 'scipy' to run multi-compartmental neurons with numpy")
 
 # scalar code
 _vectorisation_idx = 1
@@ -61,8 +65,6 @@ _vectorisation_idx = LazyArange(N)
 {{vector_code|autoindent}}
 
 {{_v_previous}}[:] = {{v}}
-{%if _scipy_available %}
-from scipy.linalg import solve_banded
 
 # Particular solution
 _b=-({{Cm}}/{{dt}}*{{v}})-_I0
@@ -82,39 +84,6 @@ _ab[0,:] = {{_ab_star0}}
 _ab[1,:] = {{_ab_star1}} - _gtot
 _ab[2,:] = {{_ab_star2}}
 {{_u_minus}}[:] = solve_banded((1,1),_ab,_b,overwrite_ab=True,overwrite_b=True)
-{% else %}
-# Pure numpy solution, very slow because it needs forward and backward passes
-# along the array that can't be vectorized
-_b = -({{Cm}} / {{dt}}*{{v}}) - _I0
-
-# Tridiagonal solving
-_forward = range(1, N)
-_backward = range(N-2, -1, -1)
-
-_bi = {{_ab_star1}} - _gtot  # main diagonal
-
-# First compartment
-{{_c}}[0] = {{_ab_star0}}[1] / _bi[0]
-{{_v_star}}[0] = _b[0] / _bi[0]
-{{_u_plus}}[0] = {{_b_plus}}[0]/_bi[0]
-{{_u_minus}}[0] = {{_b_minus}}[0]/_bi[0]
-
-# Other compartments
-{{_c}}[1:-1] = {{_ab_star0}}[2:]
-for _i in _forward:
-    _ai = {{_ab_star2}}[_i-1]  # subdiagonal
-    _m = 1.0/(_bi[_i]-_ai*{{_c}}[_i-1])
-    {{_c}}[_i] *= _m
-    {{_v_star}}[_i] = (_b[_i] - _ai*{{_v_star}}[_i-1])*_m
-    {{_u_plus}}[_i] = ({{_b_plus}}[_i] - _ai*{{_u_plus}}[_i-1])*_m
-    {{_u_minus}}[_i] = ({{_b_minus}}[_i] - _ai*{{_u_minus}}[_i-1])*_m
-
-for _i in _backward:
-    {{_v_star}}[_i] -= {{_c}}[_i]*{{_v_star}}[_i+1]
-    {{_u_plus}}[_i] -= {{_c}}[_i]*{{_u_plus}}[_i+1]
-    {{_u_minus}}[_i] -= {{_c}}[_i]*{{_u_minus}}[_i+1]
-
-{% endif %}
 
 # indexing for _P_children which contains the elements above the diagonal of the coupling matrix _P
 children_rowlength = len({{_morph_children}})//len({{_morph_children_num}})

brian2/spatialneuron/spatialneuron.py

@@ -706,21 +706,3 @@ def __init__(self, group, method, clock, order=0):
         self._morph_idxchild = group.flat_morphology.morph_idxchild
         self._starts = group.flat_morphology.starts
         self._ends = group.flat_morphology.ends
-
-    def before_run(self, run_namespace):
-        super().before_run(run_namespace)
-        # Raise a warning if the slow pure numpy version is used
-        from brian2.codegen.runtime.numpy_rt.numpy_rt import NumpyCodeObject
-
-        if type(self.code_objects[0]) == NumpyCodeObject:
-            # If numpy is used, raise a warning if scipy is not present
-            try:
-                import scipy  # noqa: F401
-            except ImportError:
-                logger.info(
-                    "SpatialNeuron will use numpy to do the numerical "
-                    "integration -- this will be very slow. Either "
-                    "switch to a different code generation target "
-                    "(e.g. cython) or install scipy.",
-                    once=True,
-                )

brian2/tests/test_spatialneuron.py

@@ -8,8 +8,22 @@
 from brian2.devices.device import reinit_and_delete
 from brian2.tests.utils import assert_allclose
 
+try:
+    import scipy
+except ImportError:
+    scipy = None
+
+
+numpy_needs_scipy = pytest.mark.skipif(
+    # Using condition string, since we cannot yet know
+    # prefs.codegen.target at module import time
+    "prefs.codegen.target == 'numpy' and not scipy",
+    reason="multi-compartmental models need scipy to run with numpy",
+)
+
 
 @pytest.mark.codegen_independent
+@numpy_needs_scipy
 def test_custom_events():
     # Set (could be moved in a setup)
     EL = -65 * mV
@@ -166,6 +180,7 @@ def test_construction_coordinates():
 
 
 @pytest.mark.long
+@numpy_needs_scipy
 def test_infinitecable():
     """
     Test simulation of an infinite cable vs. theory for current pulse (Green function)
@@ -218,6 +233,7 @@ def test_infinitecable():
 
 
 @pytest.mark.standalone_compatible
+@numpy_needs_scipy
 def test_finitecable():
     """
     Test simulation of short cylinder vs. theory for constant current.
@@ -259,6 +275,7 @@ def test_finitecable():
 
 
 @pytest.mark.standalone_compatible
+@numpy_needs_scipy
 def test_rallpack1():
     """
     Rallpack 1
@@ -318,6 +335,7 @@ def test_rallpack1():
 
 
 @pytest.mark.standalone_compatible
+@numpy_needs_scipy
 def test_rallpack2():
     """
     Rallpack 2
@@ -401,6 +419,7 @@ def test_rallpack2():
 
 @pytest.mark.standalone_compatible
 @pytest.mark.long
+@numpy_needs_scipy
 def test_rallpack3():
     """
     Rallpack 3
@@ -482,6 +501,7 @@ def test_rallpack3():
 
 
 @pytest.mark.standalone_compatible
+@numpy_needs_scipy
 def test_rall():
     """
     Test simulation of a cylinder plus two branches, with diameters according to Rall's formula
@@ -555,6 +575,7 @@ def test_rall():
 
 
 @pytest.mark.standalone_compatible
+@numpy_needs_scipy
 def test_basic_diffusion():
     # A very basic test that shows that propagation is working in a very basic
     # sense, testing all morphological classes
@@ -815,6 +836,7 @@ def test_spatialneuron_morphology_assignment():
 
 @pytest.mark.standalone_compatible
 @pytest.mark.multiple_runs
+@numpy_needs_scipy
 def test_spatialneuron_capacitive_currents():
     if prefs.core.default_float_dtype is np.float32:
         pytest.skip("Need double precision for this test")
@@ -878,6 +900,7 @@ def test_point_current():
 
 @pytest.mark.standalone_compatible
 @pytest.mark.multiple_runs
+@numpy_needs_scipy
 def test_spatialneuron_threshold_location():
     morpho = Soma(10 * um)
     morpho.axon = Cylinder(1 * um, n=2, length=20 * um)
@@ -934,6 +957,7 @@ def test_spatialneuron_threshold_location():
 
 
 @pytest.mark.standalone_compatible
+@numpy_needs_scipy
 def test_spatialneuron_timedarray():
     # See GitHub issue 1427
     ta = TimedArray([0, 1] * nA, dt=1 * ms)