Adding scipy.special modified Bessel function ops.

theano/scalar/basic_scipy.py

@@ -429,3 +429,61 @@ class J0(UnaryScalarOp):
                 j0(%(x)s);""" % locals()
         raise NotImplementedError('only floating point is implemented')
 j0 = J0(upgrade_to_float, name='j0')
+
+
+class I1(UnaryScalarOp):
+    """
+    Modified Bessel function of order 1.
+    """
+
+    @staticmethod
+    def st_impl(x):
+        return scipy.special.i1(x)
+
+    def impl(self, x):
+        if imported_scipy_special:
+            return self.st_impl(x)
+        else:
+            super(I1, self).impl(x)
+
+    def grad(self, inp, grads):
+        raise NotImplementedError()
+
+    def c_code(self, node, name, inp, out, sub):
+        x, = inp
+        z, = out
+        if node.inputs[0].type in float_types:
+            return """%(z)s =
+                i1(%(x)s);""" % locals()
+        raise NotImplementedError('only floating point is implemented')
+i1 = I1(upgrade_to_float, name='i1')
+
+
+class I0(UnaryScalarOp):
+    """
+    Modified Bessel function of order 0.
+    """
+
+    @staticmethod
+    def st_impl(x):
+        return scipy.special.i0(x)
+
+    def impl(self, x):
+        if imported_scipy_special:
+            return self.st_impl(x)
+        else:
+            super(I0, self).impl(x)
+
+    def grad(self, inp, grads):
+        x, = inp
+        gz, = grads
+        return [gz * i1(x)]
+
+    def c_code(self, node, name, inp, out, sub):
+        x, = inp
+        z, = out
+        if node.inputs[0].type in float_types:
+            return """%(z)s =
+                i0(%(x)s);""" % locals()
+        raise NotImplementedError('only floating point is implemented')
+i0 = I0(upgrade_to_float, name='i0')

theano/tensor/basic.py

@@ -2309,6 +2309,16 @@ def j1(a):
     """Bessel function of the 1'th kind"""
 
 
+@scal_elemwise
+def i0(a):
+    """Modified Bessel function of order 0."""
+
+
+@scal_elemwise
+def i1(a):
+    "Modified Bessel function of order 1."
+
+
 @_scal_elemwise
 def real(z):
     """Return real component of complex-valued tensor `z`"""

theano/tensor/inplace.py

@@ -290,6 +290,16 @@ def j1_inplace(a):
     """Bessel function of the 0'th kind"""
 
 
+@_scal_inplace
+def i0_inplace(a):
+    """Modified Bessel function of order 0."""
+
+
+@_scal_inplace
+def i1_inplace(a):
+    "Modified Bessel function of order 1."
+
+
 @_scal_inplace
 def second_inplace(a):
     """Fill `a` with `b`"""

theano/tensor/tests/test_basic.py

@@ -1721,6 +1721,8 @@ if imported_scipy_special:
     expected_chi2sf = scipy.stats.chi2.sf
     expected_j0 = scipy.special.j0
     expected_j1 = scipy.special.j1
+    expected_i0 = scipy.special.i0
+    expected_i1 = scipy.special.i1
     skip_scipy = False
     expected_erfcx = scipy.special.erfcx
 else:
@@ -1735,6 +1737,8 @@ else:
     expected_chi2sf = []
     expected_j0 = []
     expected_j1 = []
+    expected_i0 = []
+    expected_i1 = []
     skip_scipy = "scipy is not present"
 
 ErfTester = makeBroadcastTester(