Merge pull request #984 from brian-team/division

Establish consistent division semantics

brian2/codegen/generators/cpp_generator.py

@@ -88,33 +88,47 @@ def c_data_type(dtype):
     )
 
 
-mod_support_code = ''
 typestrs = ['int', 'long', 'long long', 'float', 'double', 'long double']
 floattypestrs = ['float', 'double', 'long double']
+hightype_support_code = 'template < typename T1, typename T2 > struct _higher_type;\n'
 for ix, xtype in enumerate(typestrs):
     for iy, ytype in enumerate(typestrs):
         hightype = typestrs[max(ix, iy)]
-        if xtype in floattypestrs or ytype in floattypestrs:
-            expr = 'fmod(fmod(x, y)+y, y)'
-        else:
-            expr = '((x%y)+y)%y'
-        mod_support_code += '''
-        inline {hightype} _brian_mod({xtype} ux, {ytype} uy)
-        {{
-            const {hightype} x = ({hightype})ux;
-            const {hightype} y = ({hightype})uy;
-            return {expr};
-        }}
-        '''.format(hightype=hightype, xtype=xtype, ytype=ytype, expr=expr)
-
-_universal_support_code = deindent(mod_support_code)+'''
+        hightype_support_code += '''
+template < > struct _higher_type<{xtype},{ytype}> {{ typedef {hightype} type; }};
+        '''.format(hightype=hightype, xtype=xtype, ytype=ytype)
+
+mod_support_code = '''
+
+template < typename T1, typename T2 >
+static inline typename _higher_type<T1,T2>::type
+_brian_mod(T1 x, T2 y)
+{{
+    return x-y*floor(1.0*x/y);
+}}
+'''
+
+floordiv_support_code = '''
+
+template < typename T1, typename T2 >
+static inline typename _higher_type<T1,T2>::type
+_brian_floordiv(T1 x, T2 y)
+{{
+    return floor(1.0*x/y);
+}}
+'''
+
+pow_support_code = '''
 #ifdef _MSC_VER
 #define _brian_pow(x, y) (pow((double)(x), (y)))
 #else
 #define _brian_pow(x, y) (pow((x), (y)))
 #endif
 '''
 
+_universal_support_code = (hightype_support_code + mod_support_code +
+                           floordiv_support_code + pow_support_code)
+
 
 class CPPCodeGenerator(CodeGenerator):
     '''

brian2/codegen/optimisation.py

@@ -248,6 +248,17 @@ def render_BinOp(self, node):
                 # only simplify this if the type wouldn't be cast by the operation
                 if dtype_hierarchy[right.dtype] <= dtype_hierarchy[left.dtype]:
                     return left
+        elif op.__class__.__name__ == 'FloorDiv':
+            if left.__class__.__name__ == 'Num' and left.n == 0:  # 0//x
+                if node.stateless:
+                    # Do not remove stateful functions
+                    return _replace_with_zero(left, node)
+            # Only optimise floor division by 1 if both numbers are integers,
+            # for floating point values, floor division by 1 changes the value,
+            # and division by 1.0 can change the type for an integer value
+            if (left.dtype == right.dtype == 'integer' and
+                    right.__class__.__name__ == 'Num' and right.n == 1):  # x//1
+                    return left
         # Handle addition of 0
         elif op.__class__.__name__ == 'Add':
             for operand, other in [(left, right),

brian2/codegen/runtime/cython_rt/templates/common.pyx

@@ -1,7 +1,8 @@
 #cython: boundscheck=False
 #cython: wraparound=False
-#cython: cdivision=True
+#cython: cdivision=False
 #cython: infer_types=True
+from __future__ import division
 
 import numpy as _numpy
 cimport numpy as _numpy

brian2/codegen/runtime/cython_rt/templates/spatialstateupdate.pyx

@@ -80,7 +80,7 @@
     # STEP 3: solve the coupling system
 
     # indexing for _P_children which contains the elements above the diagonal of the coupling matrix _P
-    _children_rowlength = _num{{_morph_children}}/_num{{_morph_children_num}}
+    _children_rowlength = _num{{_morph_children}}//_num{{_morph_children_num}}
 
     # STEP 3a: construct the coupling system with matrix _P in sparse form. s.t.
     # _P_diag contains the diagonal elements

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

@@ -1,4 +1,5 @@
 {# USES_VARIABLES { N } #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

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

@@ -1,8 +1,8 @@
 {# Note that we use this template only for subexpressions -- for normal arrays
    we do not generate any code but simply access the data in the underlying
    array directly. See RuntimeDevice.get_with_array #}
-
 {# USES_VARIABLES { _group_idx } #}
+from __future__ import division
 import numpy as _numpy
 # scalar code
 _vectorisation_idx = 1

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

@@ -1,5 +1,6 @@
 {# ITERATE_ALL { _idx } #}
 {# USES_VARIABLES { N } #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

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

@@ -1,5 +1,6 @@
 {# USES_VARIABLES { _group_idx, N } #}
 {# ITERATE_ALL { _target_idx } #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

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

@@ -1,5 +1,6 @@
 {# USES_VARIABLES { N } #}
 {# ALLOWS_SCALAR_WRITE #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

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

@@ -1,5 +1,6 @@
 {# USES_VARIABLES { rate, t, _spikespace, _num_source_neurons,
                     _clock_t, _clock_dt, _source_start, _source_stop, N } #}
+from __future__ import division
 import numpy as _numpy
 _spikes = {{_spikespace}}[:{{_spikespace}}[-1]]
 # Take subgroups into account

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

@@ -1,6 +1,7 @@
 {# Note that we don't write to scalar variables conditionally. The scalar code
    should therefore only include the calculation of scalar expressions
    that are used below for writing to a vector variable #}
+from __future__ import division
 import numpy as _numpy
 
 # scalar code

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

@@ -1,6 +1,7 @@
 {# USES_VARIABLES { _invr, Ri, Cm, dt, area, r_length_1, r_length_2,
                     _ab_star0, _ab_star1, _ab_star2,
                     _starts, _ends, _invr0, _invrn, _b_plus, _b_minus } #}
+from __future__ import division
 import numpy as _numpy
 from numpy import pi
 

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

@@ -9,7 +9,7 @@
                   _morph_children, _morph_children_num, _morph_idxchild,
                   _invr0, _invrn} #}
 {# ITERATE_ALL { _idx } #}
-
+from __future__ import division
 '''
 Solves the cable equation (spatial diffusion of currents).
 This is where most time-consuming time computations are done.

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

@@ -1,4 +1,5 @@
 {# USES_VARIABLES { _spikespace, t, dt, neuron_index, spike_time, period, _lastindex } #}
+from __future__ import division
 import numpy as _numpy
 
 # TODO: We don't deal with more than one spike per neuron yet

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

@@ -1,4 +1,5 @@
 {# USES_VARIABLES {N, count, _clock_t, _source_start, _source_stop, _source_N} #}
+from __future__ import division
 import numpy as _numpy
 
 {#  Get the name of the array that stores these events (e.g. the spikespace array) #}

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

@@ -1,4 +1,5 @@
 {# USES_VARIABLES { t, _clock_t, _indices } #}
+from __future__ import division
 import numpy as _numpy
 
 # Resize dynamic arrays

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

@@ -1,6 +1,7 @@
 {# ITERATE_ALL { _idx } #}
 {# USES_VARIABLES { N } #}
 {# ALLOWS_SCALAR_WRITE #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

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

@@ -1,5 +1,6 @@
 {# USES_VARIABLES { N } #}
 {# ITERATE_ALL { _idx } #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

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

@@ -1,4 +1,5 @@
 {# USES_VARIABLES { _queue } #}
+from __future__ import division
 import numpy as _numpy
 
 _spiking_synapses = _queue.peek()

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

@@ -6,12 +6,12 @@ USES_VARIABLES { _synaptic_pre, _synaptic_post, sources, targets, N,
 #}
 {# This is to show that we don't need to index the sources/targets #}
 {# ITERATE_ALL { _idx } #}
+from __future__ import division
+import numpy as _numpy
 
 {# After this code has been executed, the arrays _real_sources and
-   _real_variables contain the final indices. Having any code here it all is
+   _real_variables contain the final indices. Having any code here at all is
    only necessary for supporting subgroups #}
-import numpy as _numpy
-
 {{vector_code|autoindent}}
 
 _old_num_synapses = {{N}}

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

@@ -5,6 +5,7 @@ USES_VARIABLES { _synaptic_pre, _synaptic_post, N,
 {# WRITES_TO_READ_ONLY_VARIABLES { _synaptic_pre, _synaptic_post, N}
 #}
 {# ITERATE_ALL { _idx } #}
+from __future__ import division
 # Only need to do the initial import once, and trying is not expensive
 # only catching exceptions is expensive (which we just do the first time)
 try:

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

@@ -5,6 +5,7 @@
   Same goes for "eventspace" (e.g. spikespace) which depends on the type of
   event #}
 {# ITERATE_ALL { _idx } #}
+from __future__ import division
 import numpy as _numpy
 
 from brian2.codegen.runtime.numpy_rt.numpy_rt import LazyArange

brian2/parsing/bast.py

@@ -5,12 +5,19 @@
 '''
 
 import ast
+import weakref
+
 import numpy
 from __builtin__ import all as logical_all # defensive programming against numpy import
 
+from brian2.parsing.rendering import NodeRenderer
+from brian2.utils.logger import get_logger
+
 __all__ = ['brian_ast', 'BrianASTRenderer', 'dtype_hierarchy']
 
 
+logger = get_logger(__name__)
+
 # This codifies the idea that operations involving e.g. boolean and integer will end up
 # as integer. In general the output type will be the max of the hierarchy values here.
 dtype_hierarchy = {'boolean': 0,
@@ -163,7 +170,12 @@ def render_Call(self, node):
             raise ValueError("Variable number of arguments not supported")
         elif getattr(node, 'kwargs', None) is not None:
             raise ValueError("Keyword arguments not supported")
-        node.args = [self.render_node(subnode) for subnode in node.args]
+        args = []
+        for subnode in node.args:
+            subnode.parent = weakref.proxy(node)
+            subnode = self.render_node(subnode)
+            args.append(subnode)
+        node.args = args
         node.dtype = 'float' # default dtype
         # Condition for scalarity of function call: stateless and arguments are scalar
         node.scalar = False
@@ -192,19 +204,53 @@ def render_Call(self, node):
         return node
 
     def render_BinOp(self, node):
+        node.left.parent = weakref.proxy(node)
+        node.right.parent = weakref.proxy(node)
         node.left = self.render_node(node.left)
         node.right = self.render_node(node.right)
         # TODO: we could capture some syntax errors here, e.g. bool+bool
         # captures, e.g. int+float->float
         newdtype = dtype_hierarchy[max(dtype_hierarchy[subnode.dtype] for subnode in [node.left, node.right])]
+        if node.op.__class__.__name__ == 'Div':
+            # Division turns integers into floating point values
+            newdtype = 'float'
+            # Give a warning if the code uses floating point division where it
+            # previously might have used floor division
+            if node.left.dtype == node.right.dtype == 'integer':
+                # This would have led to floor division in earlier versions of
+                # Brian (except for the numpy target on Python 3)
+                # Ignore cases where the user already took care of this by
+                # wrapping the result of the division in int(...) or
+                # floor(...)
+                if not (hasattr(node, 'parent') and
+                        node.parent.__class__.__name__ == 'Call' and
+                        node.parent.func.id in ['int', 'floor']):
+                    rendered_expr = NodeRenderer().render_node(node)
+                    msg = ('The expression "{}" divides two integer values. '
+                           'In previous versions of Brian, this would have '
+                           'used either an integer ("flooring") or a floating '
+                           'point division, depending on the Python version '
+                           'and the code generation target. In the current '
+                           'version, it always uses a floating point '
+                           'division. Explicitly ask for an  integer division '
+                           '("//"), or turn one of the operands into a '
+                           'floating point value (e.g. replace "1/2" by '
+                           '"1.0/2") to no longer receive this '
+                           'warning.'.format(rendered_expr))
+                    logger.warn(msg, 'floating_point_division', once=True)
         node.dtype = newdtype
         node.scalar = node.left.scalar and node.right.scalar
         node.complexity = 1+node.left.complexity+node.right.complexity
         node.stateless = node.left.stateless and node.right.stateless
         return node
 
     def render_BoolOp(self, node):
-        node.values = [self.render_node(subnode) for subnode in node.values]
+        values = []
+        for subnode in node.values:
+            subnode.parent = node
+            subnode = self.render_node(subnode)
+            values.append(subnode)
+        node.values = values
         node.dtype = 'boolean'
         for subnode in node.values:
             if subnode.dtype!='boolean':
@@ -217,7 +263,12 @@ def render_BoolOp(self, node):
 
     def render_Compare(self, node):
         node.left = self.render_node(node.left)
-        node.comparators = [self.render_node(subnode) for subnode in node.comparators]
+        comparators = []
+        for subnode in node.comparators:
+            subnode.parent = node
+            subnode = self.render_node(subnode)
+            comparators.append(subnode)
+        node.comparators = comparators
         node.dtype = 'boolean'
         comparators = [node.left]+node.comparators
         node.scalar = logical_all(subnode.scalar for subnode in comparators)
@@ -227,6 +278,7 @@ def render_Compare(self, node):
         return node
 
     def render_UnaryOp(self, node):
+        node.operand.parent = node
         node.operand = self.render_node(node.operand)
         node.dtype = node.operand.dtype
         if node.dtype=='boolean' and node.op.__class__.__name__ != 'Not':