Merge pull request #470 from yarikoptic/bf-sympy

MRG: avoid relying on Pythonic casting of bool into numeric for sympy expressions

sympy 1.6 introduced breaking change (see
https://github.com/sympy/sympy/wiki/Release-Notes-for-1.6 for full log):

    Relational is no longer a subclass of Expr and does not produce nonsensical
    results in arithmetic operations. This affects all Relational subclasses (Eq,
    Ne, Gt, Ge, Lt, Le). It is no longer possible to call meaningless Expr methods
    like as_coeff_Mul on Relational instances. (#18053 by @oscarbenjamin)

which broke all Pythonish expressions like `(x > 0) * (x < 1)` to define a square
wave.  Thanks to guidance from @oscarbenjamin in
sympy/sympy#20981 I RFed all such uses to use
sympy.Piecewise instead to make it more explicit.

While at it, I also adjusted comments claiming intervals like `[0, 1]` (which
includes the endpoints), whenever in reality they were `(0, 1)` (endpoints not
included).

This is closes #466

nipy/modalities/fmri/hrf.py

@@ -103,7 +103,8 @@ def gamma_params(peak_location, peak_fwhm):
 def gamma_expr(peak_location, peak_fwhm):
     shape, scale, coef = gamma_params(peak_location, peak_fwhm)
     return (
-        coef * ((T >= 0) * (T+1.0e-14))**(shape-1)
+        coef
+        * sympy.Piecewise((T + 1e-14, T >= 0), (0, True))**(shape-1)
         * sympy.exp(-(T+1.0e-14)/scale)
         )
 
@@ -140,7 +141,7 @@ def _get_num_int(lf, dt=0.02, t=50):
 del(_gexpr); del(_dpos); del(_dgexpr)
 
 # AFNI's HRF
-_aexpr = ((T >= 0) * T)**8.6 * sympy.exp(-T/0.547)
+_aexpr = sympy.Piecewise((T, T >= 0), (0, True))**8.6 * sympy.exp(-T/0.547)
 _aexpr = _aexpr / _get_sym_int(_aexpr)
 # Numerical function
 afnit = lambdify_t(_aexpr)

nipy/modalities/fmri/tests/test_utils.py

@@ -170,8 +170,8 @@ def numerical_convolve(func1, func2, interval, dt):
 
 def test_convolve_functions():
     # replicate convolution
-    # This is a square wave on [0,1]
-    f1 = (t > 0) * (t < 1)
+    # This is a square wave on (0,1)
+    f1 = sympy.Piecewise((0, t <= 0), (1, t < 1), (0, True))
     # ff1 is the numerical implementation of same
     ff1 = lambdify(t, f1)
     # Time delta
@@ -205,7 +205,7 @@ def kern_conv2(f1, f2, f1_interval, f2_interval, dt, fill=0, name=None):
             y = ftri(time)
             assert_array_almost_equal(value, y)
         # offset square wave by 1 - offset triangle by 1
-        f2 = (t > 1) * (t < 2)
+        f2 = sympy.Piecewise((0, t <= 1), (1, t < 2), (0, True))
         tri = cfunc(f1, f2, [0, 3], [0, 3], dt)
         ftri = lambdify(t, tri)
         o1_time = np.arange(0, 3, dt)
@@ -221,7 +221,7 @@ def kern_conv2(f1, f2, f1_interval, f2_interval, dt, fill=0, name=None):
         o2_time = np.arange(0, 4, dt)
         assert_array_almost_equal(ftri(o2_time), np.r_[z1s, z1s, value])
         # offset by -0.5 - offset triangle by -0.5
-        f3 = (t > -0.5) * (t < 0.5)
+        f3 = sympy.Piecewise((0, t <= -0.5), (1, t < 0.5), (0, True))
         tri = cfunc(f1, f3, [0, 2], [-0.5, 1.5], dt)
         ftri = lambdify(t, tri)
         o1_time = np.arange(-0.5, 1.5, dt)

nipy/modalities/fmri/utils.py

@@ -528,9 +528,9 @@ def convolve_functions(f, g, f_interval, g_interval, dt,
     >>> from nipy.algorithms.statistics.formula.formulae import Term
     >>> t = Term('t')
 
-    This is a square wave on [0,1]
+    This is a square wave on (0,1)
 
-    >>> f1 = (t > 0) * (t < 1)
+    >>> f1 = sympy.Piecewise((0, t <= 0), (1, t < 1), (0, True))
 
     The convolution of ``f1`` with itself is a triangular wave on [0, 2],
     peaking at 1 with height 1