Merge pull request #691 from bouchnic/new_elemwise

theano/scalar/basic_scipy.py

@@ -2,8 +2,12 @@
 #as scipy is not always available, we put threat them separatly
 import numpy
 
-from theano.scalar.basic import UnaryScalarOp,exp,upgrade_to_float,float_types
-from theano.scalar.basic import upgrade_to_float_no_complex,complex_types,upcast
+from theano.scalar.basic import (UnaryScalarOp,
+                                 exp, upgrade_to_float,
+                                 float_types)
+from theano.scalar.basic import (upgrade_to_float_no_complex,
+                                 complex_types,
+                                 upcast)
 
 imported_scipy_special = False
 try:
@@ -12,46 +16,54 @@ try:
 except ImportError:
     pass
 
+
 class Erf(UnaryScalarOp):
     def impl(self, x):
         if imported_scipy_special:
             return scipy.special.erf(x)
         else:
-            super(Erf,self).impl(x)
+            super(Erf, self).impl(x)
+
     def grad(self, inp, grads):
         x, = inp
         gz, = grads
         if x.type in complex_types:
             raise NotImplementedError()
         elif x.type in float_types:
-            cst = numpy.asarray(2. / numpy.sqrt(numpy.pi),dtype=upcast(x.type.dtype,gz.type.dtype))
-            return gz * cst * exp(-x*x),
+            cst = numpy.asarray(2. / numpy.sqrt(numpy.pi),
+                                dtype=upcast(x.type.dtype, gz.type.dtype))
+            return gz * cst * exp(-x * x),
         else:
             return None,
+
     def c_code(self, node, name, inp, out, sub):
         x, = inp
         z, = out
         if node.inputs[0].type in complex_types:
             raise NotImplementedError('type not supported', type)
         return "%(z)s = erf(%(x)s);" % locals()
-erf = Erf(upgrade_to_float, name= 'erf')
+erf = Erf(upgrade_to_float, name='erf')
+
 
 class Erfc(UnaryScalarOp):
     def impl(self, x):
         if imported_scipy_special:
             return scipy.special.erfc(x)
         else:
-            super(Erfc,self).impl(x)
+            super(Erfc, self).impl(x)
+
     def grad(self, inp, grads):
         x, = inp
         gz, = grads
         if x.type in complex_types:
             raise NotImplementedError()
         elif x.type in float_types:
-            cst = numpy.asarray(2. / numpy.sqrt(numpy.pi),dtype=upcast(x.type.dtype,gz.type.dtype))
-            return - gz * cst * exp(-x*x),
+            cst = numpy.asarray(2. / numpy.sqrt(numpy.pi),
+                                dtype=upcast(x.type.dtype, gz.type.dtype))
+            return - gz * cst * exp(-x * x),
         else:
             return None,
+
     def c_code(self, node, name, inp, out, sub):
         x, = inp
         z, = out
@@ -60,8 +72,7 @@ class Erfc(UnaryScalarOp):
         return "%(z)s = erfc(%(x)s);" % locals()
 
 # scipy.special.erfc don't support complex. Why?
-erfc = Erfc(upgrade_to_float_no_complex, name = 'erfc')
-
+erfc = Erfc(upgrade_to_float_no_complex, name='erfc')
 
 
 class GammaLn(UnaryScalarOp):
@@ -70,26 +81,30 @@ class GammaLn(UnaryScalarOp):
     """
     @staticmethod
     def st_impl(x):
-        return special.gammaln(x)
+        return scipy.special.gammaln(x)
+
     def impl(self, x):
         return GammaLn.st_impl(x)
+
     def grad(self, inp, grads):
         x, = inp
         gz, = grads
         return [gz * psi(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 =
                 lgamma(%(x)s);""" % locals()
-        raise NotImplementedError('only floatingpoint is implemented')
+        raise NotImplementedError('only floating point is implemented')
+
     def __eq__(self, other):
         return type(self) == type(other)
+
     def __hash__(self):
         return hash(type(self))
-    
-gammaln  = GammaLn(upgrade_to_float, name='gammaln')
+gammaln = GammaLn(upgrade_to_float, name='gammaln')
 
 
 class Psi(UnaryScalarOp):
@@ -98,62 +113,67 @@ class Psi(UnaryScalarOp):
     """
     @staticmethod
     def st_impl(x):
-        return special.psi(x)
+        return scipy.special.psi(x)
+
     def impl(self, x):
         return Psi.st_impl(x)
-    
-    #def grad()  no gradient now 
-    
+
+    def grad(self, inputs, outputs_gradients):
+        raise NotImplementedError()
+        return [None]
+
     def c_support_code(self):
-        return ( 
-        """
+        return (
+"""
 #ifndef _PSIFUNCDEFINED
 #define _PSIFUNCDEFINED
 double _psi(double x){
 
-    /*taken from 
-    Bernardo, J. M. (1976). Algorithm AS 103: Psi (Digamma) Function. Applied Statistics. 25 (3), 315-317. 
+    /*taken from
+    Bernardo, J. M. (1976). Algorithm AS 103:
+    Psi (Digamma) Function. Applied Statistics. 25 (3), 315-317.
     http://www.uv.es/~bernardo/1976AppStatist.pdf */
-    
+
     double y, R, psi_ = 0;
     double S  = 1.0e-5;
     double C = 8.5;
     double S3 = 8.333333333e-2;
     double S4 = 8.333333333e-3;
     double S5 = 3.968253968e-3;
-    double D1 = -0.5772156649  ;   
+    double D1 = -0.5772156649;
 
     y = x;
-    
+
     if (y <= 0.0)
         return psi_;
-        
+
     if (y <= S )
         return D1 - 1.0/y;
-    
+
     while (y < C){
         psi_ = psi_ - 1.0 / y;
         y = y + 1;}
-    
+
     R = 1.0 / y;
     psi_ = psi_ + log(y) - .5 * R ;
     R= R*R;
     psi_ = psi_ - R * (S3 - R * (S4 - R * S5));
-    
+
     return psi_;}
     #endif
         """ )
+
     def c_code(self, node, name, inp, out, sub):
         x, = inp
         z, = out
         if node.inputs[0].type in float_types:
             return """%(z)s =
                 _psi(%(x)s);""" % locals()
-        raise NotImplementedError('only floatingpoint is implemented')
-    
+        raise NotImplementedError('only floating point is implemented')
+
     def __eq__(self, other):
         return type(self) == type(other)
+
     def __hash__(self):
         return hash(type(self))
-    
 psi = Psi(upgrade_to_float, name='psi')

theano/tensor/basic.py

@@ -2639,10 +2639,12 @@ def erf(a):
 def erfc(a):
     """complementary error function"""
 
+
 @_scal_elemwise
 def gammaln(a):
     """log gamma function"""
-    
+
+
 @_scal_elemwise
 def psi(a):
     """derivative of log gamma function"""

theano/tensor/inplace.py

@@ -202,11 +202,11 @@ def erf_inplace(a):
 @_scal_inplace
 def erfc_inplace(a):
     """complementary error function"""
-    
+
 @_scal_inplace
 def gammaln_inplace(a):
     """log gamma function"""
-    
+
 @_scal_inplace
 def psi_inplace(a):
     """derivative of log gamma function"""

theano/tensor/tests/test_basic.py

@@ -1237,7 +1237,7 @@ del _good_broadcast_unary_normal_no_int['integers']
 if imported_scipy_special:
     expected_erf = scipy.special.erf
     expected_erfc = scipy.special.erfc
-    
+
     expected_gammaln = scipy.special.gammaln
     expected_psi = scipy.special.psi
     skip_scipy = False
@@ -1246,69 +1246,85 @@ else:
     expected_erfc = []
     skip_scipy = "scipy is not present"
 
-ErfTester = makeBroadcastTester(op = tensor.erf,
-                                expected = expected_erf,
-                                good = _good_broadcast_unary_normal,
-                                grad = _grad_broadcast_unary_normal,
-                                eps = 2e-10,
-                                mode = mode_no_scipy,
-                                skip = skip_scipy)
-ErfInplaceTester = makeBroadcastTester(op = inplace.erf_inplace,
-                                       expected = expected_erf,
-                                       good = _good_broadcast_unary_normal_no_int,
-                                       grad = _grad_broadcast_unary_normal,
-                                       mode = mode_no_scipy,
-                                       eps = 2e-10,
-                                       inplace = True,
-                                       skip = skip_scipy)
-
-ErfcTester = makeBroadcastTester(op = tensor.erfc,
-                                 expected = expected_erfc,
-                                 good = _good_broadcast_unary_normal_no_int_no_complex,
-                                 grad = _grad_broadcast_unary_normal,
-                                 eps = 2e-10,
-                                 mode = mode_no_scipy,
-                                 skip = skip_scipy)
-ErfcInplaceTester = makeBroadcastTester(op = inplace.erfc_inplace,
-                                        expected = expected_erfc,
-                                        good = _good_broadcast_unary_normal_no_int_no_complex,
-                                        grad = _grad_broadcast_unary_normal,
-                                        eps = 2e-10,
-                                        mode = mode_no_scipy,
-                                        inplace = True,
-                                        skip = skip_scipy)
-
-GammaLnTester = makeBroadcastTester(op = tensor.gammaln,
-                                 expected = expected_gammaln,
-                                 good = _good_broadcast_unary_normal_no_int_no_complex,
-                                 grad = _grad_broadcast_unary_normal,
-                                 eps = 2e-10,
-                                 mode = mode_no_scipy,
-                                 skip = skip_scipy)
-GammaLnInplaceTester = makeBroadcastTester(op = inplace.gammaln_inplace,
-                                        expected = expected_erfc,
-                                        good = _good_broadcast_unary_normal_no_int_no_complex,
-                                        grad = _grad_broadcast_unary_normal,
-                                        eps = 2e-10,
-                                        mode = mode_no_scipy,
-                                        inplace = True,
-                                        skip = skip_scipy)
-
-PsiTester = makeBroadcastTester(op = tensor.psi,
-                                 expected = expected_psi,
-                                 good = _good_broadcast_unary_normal_no_int_no_complex,
-                                 grad = _grad_broadcast_unary_normal,
-                                 eps = 2e-10,
-                                 mode = mode_no_scipy,
-                                 skip = skip_scipy)
-PsiInplaceTester = makeBroadcastTester(op = inplace.psi_inplace,
-                                        expected = expected_psi,
-                                        good = _good_broadcast_unary_normal_no_int_no_complex,
-                                        grad = _grad_broadcast_unary_normal,
-                                        eps = 2e-10,
-                                        mode = mode_no_scipy,
-                                        inplace = True,
-                                        skip = skip_scipy)
+ErfTester = makeBroadcastTester(
+    op=tensor.erf,
+    expected=expected_erf,
+    good=_good_broadcast_unary_normal,
+    grad=_grad_broadcast_unary_normal,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy)
+ErfInplaceTester = makeBroadcastTester(
+    op=inplace.erf_inplace,
+    expected=expected_erf,
+    good=_good_broadcast_unary_normal_no_int,
+    grad=_grad_broadcast_unary_normal,
+    mode=mode_no_scipy,
+    eps=2e-10,
+    inplace=True,
+    skip=skip_scipy)
+
+ErfcTester = makeBroadcastTester(
+    op=tensor.erfc,
+    expected=expected_erfc,
+    good=_good_broadcast_unary_normal_no_int_no_complex,
+    grad=_grad_broadcast_unary_normal,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy)
+ErfcInplaceTester = makeBroadcastTester(
+    op=inplace.erfc_inplace,
+    expected=expected_erfc,
+    good=_good_broadcast_unary_normal_no_int_no_complex,
+    grad=_grad_broadcast_unary_normal,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    inplace=True,
+    skip=skip_scipy)
+
+_good_broadcast_unary_gammaln = dict(
+    normal=(rand_ranged(-1 + 1e-2, 10, (2, 3)),),
+    empty=(numpy.asarray([]),),)
+_grad_broadcast_unary_gammaln = dict(
+    normal=(rand_ranged(1e-8, 10, (2, 3)),),)
+
+GammaLnTester = makeBroadcastTester(
+    op=tensor.gammaln,
+    expected=expected_gammaln,
+    good=_good_broadcast_unary_gammaln,
+    grad=_grad_broadcast_unary_gammaln,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy)
+GammaLnInplaceTester = makeBroadcastTester(
+    op=inplace.gammaln_inplace,
+    expected=expected_gammaln,
+    good=_good_broadcast_unary_gammaln,
+    grad=_grad_broadcast_unary_gammaln,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    inplace=True,
+    skip=skip_scipy)
+
+_good_broadcast_unary_psi = dict(
+    normal=(rand_ranged(1, 10, (2, 3)),),
+    empty=(numpy.asarray([]),),)
+
+PsiTester = makeBroadcastTester(
+    op=tensor.psi,
+    expected=expected_psi,
+    good=_good_broadcast_unary_psi,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy)
+PsiInplaceTester = makeBroadcastTester(
+    op=inplace.psi_inplace,
+    expected=expected_psi,
+    good=_good_broadcast_unary_psi,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    inplace=True,
+    skip=skip_scipy)
 
 ZerosLikeTester = makeBroadcastTester(
         op=tensor.zeros_like,