Adding real-order Bessel / Modified Bessel function ops.

theano/scalar/basic_scipy.py

 from __future__ import absolute_import, print_function, division
-# definition theano.scalar op that have their python implementation taked from scipy
-# as scipy is not always available, we treat them separatly
+# Definitions of theano.scalar ops that have their python implementation taken
+# from SciPy. As SciPy is not always available, we treat them separately.
+
 import numpy
 
 import theano
+from theano.gradient import grad_not_implemented
 from theano.scalar.basic import (UnaryScalarOp, BinaryScalarOp,
                                  exp, upgrade_to_float,
                                  upgrade_to_float64,
@@ -373,9 +375,36 @@ class Chi2SF(BinaryScalarOp):
 chi2sf = Chi2SF(upgrade_to_float64, name='chi2sf')
 
 
+class Jv(BinaryScalarOp):
+    """
+    Bessel function of the first kind of order v (real).
+    """
+
+    @staticmethod
+    def st_impl(v, x):
+        return scipy.special.jv(v, x)
+
+    def impl(self, v, x):
+        if imported_scipy_special:
+            return self.st_impl(v, x)
+        else:
+            super(Jv, self).impl(v, x)
+
+    def grad(self, inputs, grads):
+        if imported_scipy_special:
+            v, x = inputs
+            gz, = grads
+            return [grad_not_implemented(self, 0, v),
+                    gz * (jv(v - 1, x) - jv(v + 1, x)) / 2.]
+        else:
+            super(Jv, self).grad(v, x)
+
+jv = Jv(upgrade_to_float, name='jv')
+
+
 class J1(UnaryScalarOp):
     """
-    Bessel function of the 1'th kind
+    Bessel function of the first kind of order 1.
     """
 
     @staticmethod
@@ -388,8 +417,13 @@ class J1(UnaryScalarOp):
         else:
             super(J1, self).impl(x)
 
-    def grad(self, inp, grads):
-        raise NotImplementedError()
+    def grad(self, inputs, grads):
+        if imported_scipy_special:
+            x, = inputs
+            gz, = grads
+            return [gz * (j0(x) - jv(2, x)) / 2.]
+        else:
+            super(J1, self).grad(v, x)
 
     def c_code(self, node, name, inp, out, sub):
         x, = inp
@@ -398,12 +432,13 @@ class J1(UnaryScalarOp):
             return """%(z)s =
                 j1(%(x)s);""" % locals()
         raise NotImplementedError('only floating point is implemented')
+
 j1 = J1(upgrade_to_float, name='j1')
 
 
 class J0(UnaryScalarOp):
     """
-    Bessel function of the 0'th kind
+    Bessel function of the first kind of order 0.
     """
 
     @staticmethod
@@ -428,12 +463,40 @@ class J0(UnaryScalarOp):
             return """%(z)s =
                 j0(%(x)s);""" % locals()
         raise NotImplementedError('only floating point is implemented')
+
 j0 = J0(upgrade_to_float, name='j0')
 
 
+class Iv(BinaryScalarOp):
+    """
+    Modified Bessel function of the first kind of order v (real).
+    """
+
+    @staticmethod
+    def st_impl(v, x):
+        return scipy.special.iv(v, x)
+
+    def impl(self, v, x):
+        if imported_scipy_special:
+            return self.st_impl(v, x)
+        else:
+            super(Iv, self).impl(v, x)
+
+    def grad(self, inputs, grads):
+        if imported_scipy_special:
+            v, x = inputs
+            gz, = grads
+            return [grad_not_implemented(self, 0, v),
+                    gz * (iv(v - 1, x) + iv(v + 1, x)) / 2.]
+        else:
+            super(Iv, self).grad(v, x)
+
+iv = Iv(upgrade_to_float, name='iv')
+
+
 class I1(UnaryScalarOp):
     """
-    Modified Bessel function of order 1.
+    Modified Bessel function of the first kind of order 1.
     """
 
     @staticmethod
@@ -446,15 +509,20 @@ class I1(UnaryScalarOp):
         else:
             super(I1, self).impl(x)
 
-    def grad(self, inp, grads):
-        raise NotImplementedError()
+    def grad(self, inputs, grads):
+        if imported_scipy_special:
+            x, = inputs
+            gz, = grads
+            return [gz * (i0(x) + iv(2, x)) / 2.]
+        else:
+            super(I1, self).grad(v, x)
 
 i1 = I1(upgrade_to_float, name='i1')
 
 
 class I0(UnaryScalarOp):
     """
-    Modified Bessel function of order 0.
+    Modified Bessel function of the first kind of order 0.
     """
 
     @staticmethod

theano/tensor/basic.py

@@ -2300,23 +2300,33 @@ def chi2sf(x, k):
 
 
 @_scal_elemwise
-def j0(a):
-    """Bessel function of the 0'th kind"""
+def j0(x):
+    """Bessel function of the first kind of order 0."""
 
 
 @_scal_elemwise
-def j1(a):
-    """Bessel function of the 1'th kind"""
+def j1(x):
+    """Bessel function of the first kind of order 1."""
 
 
 @_scal_elemwise
-def i0(a):
-    """Modified Bessel function of order 0."""
+def jv(v, x):
+    """Bessel function of the first kind of order v (real)."""
 
 
 @_scal_elemwise
-def i1(a):
-    "Modified Bessel function of order 1."
+def i0(x):
+    """Modified Bessel function of the first kind of order 0."""
+
+
+@_scal_elemwise
+def i1(x):
+    """Modified Bessel function of the first kind of order 1."""
+
+
+@_scal_elemwise
+def iv(v, x):
+    """Modified Bessel function of the first kind of order v (real)."""
 
 
 @_scal_elemwise

theano/tensor/inplace.py

@@ -281,23 +281,33 @@ def chi2sf_inplace(x, k):
 
 
 @_scal_inplace
-def j0_inplace(a):
-    """Bessel function of the 0'th kind"""
+def j0_inplace(x):
+    """Bessel function of the first kind of order 0."""
 
 
 @_scal_inplace
-def j1_inplace(a):
-    """Bessel function of the 0'th kind"""
+def j1_inplace(x):
+    """Bessel function of the first kind of order 1."""
 
 
 @_scal_inplace
-def i0_inplace(a):
-    """Modified Bessel function of order 0."""
+def jv_inplace(v, x):
+    """Bessel function of the first kind of order v (real)."""
 
 
 @_scal_inplace
-def i1_inplace(a):
-    "Modified Bessel function of order 1."
+def i0_inplace(x):
+    """Modified Bessel function of the first kind of order 0."""
+
+
+@_scal_inplace
+def i1_inplace(x):
+    """Modified Bessel function of the first kind of order 1."""
+
+
+@_scal_inplace
+def iv_inplace(v, x):
+    """Modified Bessel function of the first kind of order v (real)."""
 
 
 @_scal_inplace

theano/tensor/tests/test_basic.py

@@ -1721,8 +1721,10 @@ if imported_scipy_special:
     expected_chi2sf = scipy.stats.chi2.sf
     expected_j0 = scipy.special.j0
     expected_j1 = scipy.special.j1
+    expected_jv = scipy.special.jv
     expected_i0 = scipy.special.i0
     expected_i1 = scipy.special.i1
+    expected_iv = scipy.special.iv
     skip_scipy = False
     expected_erfcx = scipy.special.erfcx
 else:
@@ -1737,8 +1739,10 @@ else:
     expected_chi2sf = []
     expected_j0 = []
     expected_j1 = []
+    expected_jv = []
     expected_i0 = []
     expected_i1 = []
+    expected_iv = []
     skip_scipy = "scipy is not present"
 
 ErfTester = makeBroadcastTester(
@@ -1928,6 +1932,7 @@ J0Tester = makeBroadcastTester(
     eps=2e-10,
     mode=mode_no_scipy,
     skip=skip_scipy)
+
 J0InplaceTester = makeBroadcastTester(
     op=inplace.j0_inplace,
     expected=expected_j0,
@@ -1945,6 +1950,7 @@ J1Tester = makeBroadcastTester(
     eps=2e-10,
     mode=mode_no_scipy,
     skip=skip_scipy)
+
 J1InplaceTester = makeBroadcastTester(
     op=inplace.j1_inplace,
     expected=expected_j1,
@@ -1954,13 +1960,41 @@ J1InplaceTester = makeBroadcastTester(
     inplace=True,
     skip=skip_scipy)
 
+_good_broadcast_binary_jv = dict(
+    normal=(rand_ranged(-5, 5, (2, 3)),
+            rand_ranged(0.1, 8, (2, 3))),
+    empty=(numpy.asarray([], dtype=config.floatX),
+           numpy.asarray([], dtype=config.floatX)),
+    integers=(randint_ranged(-5, 5, (2, 3)),
+              randint_ranged(1, 8, (2, 3))),
+    uint8=(randint_ranged(-5, 5, (2, 3)).astype('uint8'),
+           randint_ranged(1, 8, (2, 3)).astype('uint8')),
+    uint16=(randint_ranged(-5, 5, (2, 3)).astype('uint16'),
+            randint_ranged(1, 8, (2, 3)).astype('uint16')))
+
+JvTester = makeBroadcastTester(
+    op=tensor.jv,
+    expected=expected_jv,
+    good=_good_broadcast_binary_jv,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy)
+
+JvInplaceTester = makeBroadcastTester(
+    op=inplace.jv_inplace,
+    expected=expected_jv,
+    good=_good_broadcast_binary_jv,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    inplace=True,
+    skip=skip_scipy)
+
 _good_broadcast_unary_i = dict(
     normal=(rand_ranged(-10, 10, (2, 3)),),
     empty=(numpy.asarray([], dtype=config.floatX),),
     int=(randint_ranged(-10, 10, (2, 3)),),
     uint8=(randint_ranged(0, 10, (2, 3)).astype('uint8'),),
-    uint16=(randint_ranged(0, 10, (2, 3)).astype('uint16'),)
-)
+    uint16=(randint_ranged(0, 10, (2, 3)).astype('uint16'),))
 
 I0Tester = makeBroadcastTester(
     op=tensor.i0,
@@ -1973,7 +2007,7 @@ I0Tester = makeBroadcastTester(
 
 I0InplaceTester = makeBroadcastTester(
     op=inplace.i0_inplace,
-    expected=expected_j0,
+    expected=expected_i0,
     good=_good_broadcast_unary_i,
     grad=_good_broadcast_unary_i,
     eps=2e-10,
@@ -1998,6 +2032,35 @@ I1InplaceTester = makeBroadcastTester(
     inplace=True,
     skip=skip_scipy)
 
+_good_broadcast_binary_iv = dict(
+    normal=(rand_ranged(-5, 5, (2, 3)),
+            rand_ranged(-10, 10, (2, 3))),
+    empty=(numpy.asarray([], dtype=config.floatX),
+           numpy.asarray([], dtype=config.floatX)),
+    integers=(randint_ranged(-5, 5, (2, 3)),
+              randint_ranged(-10, 10, (2, 3))),
+    uint8=(randint_ranged(-5, 5, (2, 3)).astype('uint8'),
+           randint_ranged(-10, 10, (2, 3)).astype('uint8')),
+    uint16=(randint_ranged(-5, 5, (2, 3)).astype('uint16'),
+            randint_ranged(-10, 10, (2, 3)).astype('uint16')))
+
+IvTester = makeBroadcastTester(
+    op=tensor.iv,
+    expected=expected_iv,
+    good=_good_broadcast_binary_iv,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy)
+
+IvInplaceTester = makeBroadcastTester(
+    op=inplace.iv_inplace,
+    expected=expected_iv,
+    good=_good_broadcast_binary_iv,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    inplace=True,
+    skip=skip_scipy)
+
 ZerosLikeTester = makeBroadcastTester(
         op=tensor.zeros_like,
         expected=numpy.zeros_like,