Merge pull request #686 from jsalvatier/master

theano/scalar/basic_scipy.py

@@ -61,3 +61,99 @@ class Erfc(UnaryScalarOp):
 
 # scipy.special.erfc don't support complex. Why?
 erfc = Erfc(upgrade_to_float_no_complex, name = 'erfc')
+
+
+
+class GammaLn(UnaryScalarOp):
+    """
+    Log gamma function.
+    """
+    @staticmethod
+    def st_impl(x):
+        return 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')
+    def __eq__(self, other):
+        return type(self) == type(other)
+    def __hash__(self):
+        return hash(type(self))
+    
+gammaln  = GammaLn(upgrade_to_float, name='gammaln')
+
+
+class Psi(UnaryScalarOp):
+    """
+    Derivative of log gamma function.
+    """
+    @staticmethod
+    def st_impl(x):
+        return special.psi(x)
+    def impl(self, x):
+        return Psi.st_impl(x)
+    
+    #def grad()  no gradient now 
+    
+    def c_support_code(self):
+        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. 
+    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  ;   
+
+    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')
+    
+    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,6 +2639,14 @@ 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"""
+
 
 @_scal_elemwise_with_nfunc('real', 1, -1)
 def real(z):

theano/tensor/inplace.py

@@ -202,6 +202,14 @@ 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"""
 
 @_scal_inplace
 def second_inplace(a):

theano/tensor/tests/test_basic.py

@@ -1237,6 +1237,9 @@ 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
 else:
     expected_erf = []
@@ -1275,6 +1278,38 @@ ErfcInplaceTester = makeBroadcastTester(op = inplace.erfc_inplace,
                                         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)
+
 ZerosLikeTester = makeBroadcastTester(
         op=tensor.zeros_like,
         expected=numpy.zeros_like,