Skip to content

Commit

Permalink
Merge 3b3b72d into b9d9be8
Browse files Browse the repository at this point in the history
  • Loading branch information
mstimberg committed Mar 14, 2023
2 parents b9d9be8 + 3b3b72d commit 5b485e6
Show file tree
Hide file tree
Showing 5 changed files with 31 additions and 64 deletions.
2 changes: 1 addition & 1 deletion brian2/codegen/generators/base.py
Expand Up @@ -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
Expand Down
10 changes: 0 additions & 10 deletions brian2/codegen/generators/numpy_generator.py
Expand Up @@ -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
Expand Down
41 changes: 6 additions & 35 deletions brian2/codegen/runtime/numpy_rt/templates/spatialstateupdate.py_
Expand Up @@ -51,6 +51,12 @@ 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("Simulating multi-compartmental neurons with the "
"numpy code generation target requires the 'scipy' "
"package.")

# scalar code
_vectorisation_idx = 1
Expand All @@ -61,8 +67,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
Expand All @@ -82,39 +86,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}})
Expand Down
18 changes: 0 additions & 18 deletions brian2/spatialneuron/spatialneuron.py
Expand Up @@ -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,
)
24 changes: 24 additions & 0 deletions brian2/tests/test_spatialneuron.py
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -259,6 +275,7 @@ def test_finitecable():


@pytest.mark.standalone_compatible
@numpy_needs_scipy
def test_rallpack1():
"""
Rallpack 1
Expand Down Expand Up @@ -318,6 +335,7 @@ def test_rallpack1():


@pytest.mark.standalone_compatible
@numpy_needs_scipy
def test_rallpack2():
"""
Rallpack 2
Expand Down Expand Up @@ -401,6 +419,7 @@ def test_rallpack2():

@pytest.mark.standalone_compatible
@pytest.mark.long
@numpy_needs_scipy
def test_rallpack3():
"""
Rallpack 3
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 5b485e6

Please sign in to comment.