Merge pull request #1410 from ebuchman/chi2sf/master

theano/scalar/basic_scipy.py

 #definition theano.scalar op that have their python implementation taked from scipy
-#as scipy is not always available, we put threat them separatly
+#as scipy is not always available, we treat them separatly
 import numpy
 
-from theano.scalar.basic import (UnaryScalarOp,
+from theano.scalar.basic import (UnaryScalarOp, BinaryScalarOp,
                                  exp, upgrade_to_float,
                                  float_types)
 from theano.scalar.basic import (upgrade_to_float_no_complex,
@@ -12,6 +12,7 @@ from theano.scalar.basic import (upgrade_to_float_no_complex,
 imported_scipy_special = False
 try:
     import scipy.special
+    import scipy.stats
     imported_scipy_special = True
 # Importing scipy.special may raise ValueError.
 # See http://projects.scipy.org/scipy/ticket/1739
@@ -152,7 +153,10 @@ class Gamma(UnaryScalarOp):
         return scipy.special.gamma(x)
 
     def impl(self, x):
+        if imported_scipy_special:
             return Gamma.st_impl(x)
+        else:
+            super(Gamma, self).impl(x)
 
     def grad(self, (x, ), (gz, )):
         return gz * gamma(x) * psi(x),
@@ -179,8 +183,10 @@ class GammaLn(UnaryScalarOp):
         return scipy.special.gammaln(x)
 
     def impl(self, x):
+        if imported_scipy_special:
             return GammaLn.st_impl(x)
-
+        else:
+            super(GammaLn, self).impl(x)
     def grad(self, inp, grads):
         x, = inp
         gz, = grads
@@ -211,7 +217,10 @@ class Psi(UnaryScalarOp):
         return scipy.special.psi(x)
 
     def impl(self, x):
+        if imported_scipy_special:
             return Psi.st_impl(x)
+        else:
+            super(Psi, self).impl(x)
 
     def grad(self, inputs, outputs_gradients):
         raise NotImplementedError()
@@ -279,3 +288,291 @@ DEVICE double _psi(double x){
     def __hash__(self):
         return hash(type(self))
 psi = Psi(upgrade_to_float, name='psi')
+
+
+class Chi2SF(BinaryScalarOp):
+    """
+    Compute (1 - chi2_cdf(x))
+        ie. chi2 pvalue (chi2 'survival function')
+    """
+
+    @staticmethod
+    def st_impl(x, k):
+        return scipy.stats.chi2.sf(x, k)
+
+    def impl(self, x, k):
+        if imported_scipy_special:
+            return Chi2SF.st_impl(x, k)
+        else:
+            super(Chi2SF, self).impl(x, k)
+
+    def c_support_code(self):
+        return(
+               """
+                   //For GPU support
+                   #ifdef __CUDACC__
+                   #define DEVICE __device__
+                   #else
+                   #define DEVICE
+                   #endif
+
+                   #ifndef _CHI2FUNCDEFINED
+                   #define _CHI2FUNCDEFINED
+
+                   /*----------------------------------------------------------------------
+                   File    : gamma.c
+                   Contents: computation of the (incomplete/regularized) gamma function
+                   Author  : Christian Borgelt
+                   History : 2002.07.04 file created
+                   2003.05.19 incomplete Gamma function added
+                   2008.03.14 more incomplete Gamma functions added
+                   2008.03.15 table of factorials and logarithms added
+                   2008.03.17 gamma distribution functions added
+                   ----------------------------------------------------------------------*/
+                   #include <stdio.h>
+                   #include <stdlib.h>
+
+
+                   #include <assert.h>
+                   #include <float.h>
+                   #include <math.h>
+                   /*----------------------------------------------------------------------
+                   Preprocessor Definitions
+                   ----------------------------------------------------------------------*/
+                   #define LN_BASE      2.71828182845904523536028747135  /* e */
+                   #define SQRT_PI      1.77245385090551602729816748334  /* \sqrt(\pi) */
+                   #define LN_PI        1.14472988584940017414342735135  /* \ln(\pi) */
+                   #define LN_SQRT_2PI  0.918938533204672741780329736406
+                   /* \ln(\sqrt(2\pi)) */
+                   #define EPSILON      2.2204460492503131e-16
+                   #define EPS_QTL      1.4901161193847656e-08
+                   #define MAXFACT      170
+                   #define MAXITER      1024
+                   #define TINY         (EPSILON *EPSILON *EPSILON)
+
+                   /*----------------------------------------------------------------------
+                   Table of Factorials/Gamma Values
+                   ----------------------------------------------------------------------*/
+                   DEVICE static double _facts[MAXFACT+1] = { 0 };
+                   DEVICE static double _logfs[MAXFACT+1];
+                   DEVICE static double _halfs[MAXFACT+1];
+                   DEVICE static double _loghs[MAXFACT+1];
+
+                   /*----------------------------------------------------------------------
+                   Functions
+                   ----------------------------------------------------------------------*/
+                   DEVICE static void _init (void)
+                   {                               /* --- init. factorial tables */
+                   int    i;                     /* loop variable */
+                   double x = 1;                 /* factorial */
+
+                   _facts[0] = _facts[1] = 1;    /* store factorials for 0 and 1 */
+                   _logfs[0] = _logfs[1] = 0;    /* and their logarithms */
+                   for (i = 1; ++i <= MAXFACT; ) {
+                   _facts[i] = x *= i;         /* initialize the factorial table */
+                   _logfs[i] = log(x);         /* and the table of their logarithms */
+                   }
+                   _halfs[0] = x = SQRT_PI;      /* store Gamma(0.5) */
+                   _loghs[0] = 0.5*LN_PI;        /* and its logarithm */
+                   for (i = 0; ++i < MAXFACT; ) {
+                   _halfs[i] = x *= i-0.5;     /* initialize the table for */
+                   _loghs[i] = log(x);         /* the Gamma function of half numbers */
+                   }                             /* and the table of their logarithms */
+                   }  /* _init() */
+
+                   /*--------------------------------------------------------------------*/
+                   #if 0
+
+                   double logGamma (double n)
+                   {                               /* --- compute ln(Gamma(n))         */
+                   double s;                     /*           = ln((n-1)!), n \in IN */
+
+                   assert(n > 0);                /* check the function argument */
+                   if (_facts[0] <= 0) _init();  /* initialize the tables */
+                   if (n < MAXFACT +1 +4 *EPSILON) {
+                   if (fabs(  n -floor(  n)) < 4 *EPSILON)
+                   return _logfs[(int)floor(n)-1];
+                   if (fabs(2*n -floor(2*n)) < 4 *EPSILON)
+                   return _loghs[(int)floor(n)];
+                   }                             /* try to get the value from a table */
+                   s =  1.000000000190015        /* otherwise compute it */
+                   + 76.18009172947146      /(n+1)
+                   - 86.50532032941677      /(n+2)
+                   + 24.01409824083091      /(n+3)
+                   -  1.231739572450155     /(n+4)
+                   +  0.1208650972866179e-2 /(n+5)
+                   -  0.5395239384953e-5    /(n+6);
+                   return (n+0.5) *log((n+5.5)/LN_BASE) +(LN_SQRT_2PI +log(s/n) -5.0);
+                   }  /* logGamma() */
+
+                   #else /*--------------------------------------------------------------*/
+
+                   DEVICE double logGamma (double n)
+                   {                               /* --- compute ln(Gamma(n))         */
+                   double s;                     /*           = ln((n-1)!), n \in IN */
+
+                   assert(n > 0);                /* check the function argument */
+                   if (_facts[0] <= 0) _init();  /* initialize the tables */
+                   if (n < MAXFACT +1 +4 *EPSILON) {
+                   if (fabs(  n -floor(  n)) < 4 *EPSILON)
+                   return _logfs[(int)floor(n)-1];
+                   if (fabs(2*n -floor(2*n)) < 4 *EPSILON)
+                   return _loghs[(int)floor(n)];
+                   }                             /* try to get the value from a table */
+                   s =    0.99999999999980993227684700473478  /* otherwise compute it */
+                   +  676.520368121885098567009190444019 /(n+1)
+                   - 1259.13921672240287047156078755283  /(n+2)
+                   +  771.3234287776530788486528258894   /(n+3)
+                   -  176.61502916214059906584551354     /(n+4)
+                   +   12.507343278686904814458936853    /(n+5)
+                   -    0.13857109526572011689554707     /(n+6)
+                   +    9.984369578019570859563e-6       /(n+7)
+                   +    1.50563273514931155834e-7        /(n+8);
+                   return (n+0.5) *log((n+7.5)/LN_BASE) +(LN_SQRT_2PI +log(s/n) -7.0);
+                   }  /* logGamma() */
+
+                   #endif
+                   /*----------------------------------------------------------------------
+                   Use Lanczos' approximation
+                   \Gamma(n+1) = (n+\gamma+0.5)^(n+0.5)
+                   * e^{-(n+\gamma+0.5)}
+                   * \sqrt{2\pi}
+                   * (c_0 +c_1/(n+1) +c_2/(n+2) +...+c_n/(n+k) +\epsilon)
+                   and exploit the recursion \Gamma(n+1) = n *\Gamma(n) once,
+                   i.e., compute \Gamma(n) as \Gamma(n+1) /n.
+
+                   For the choices \gamma = 5, k = 6, and c_0 to c_6 as defined
+                   in the first version, it is |\epsilon| < 2e-10 for all n > 0.
+
+                   Source: W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery
+                   Numerical Recipes in C - The Art of Scientific Computing
+                   Cambridge University Press, Cambridge, United Kingdom 1992
+                   pp. 213-214
+
+                   For the choices gamma = 7, k = 8, and c_0 to c_8 as defined
+                   in the second version, the value is slightly more accurate.
+                   ----------------------------------------------------------------------*/
+
+                   DEVICE double Gamma (double n)
+                   {                               /* --- compute Gamma(n) = (n-1)! */
+                   assert(n > 0);                /* check the function argument */
+                   if (_facts[0] <= 0) _init();  /* initialize the tables */
+                   if (n < MAXFACT +1 +4 *EPSILON) {
+                   if (fabs(  n -floor(  n)) < 4 *EPSILON)
+                   return _facts[(int)floor(n)-1];
+                   if (fabs(2*n -floor(2*n)) < 4 *EPSILON)
+                   return _halfs[(int)floor(n)];
+                   }                             /* try to get the value from a table */
+                   return exp(logGamma(n));      /* compute through natural logarithm */
+                   }  /* Gamma() */
+
+                   /*--------------------------------------------------------------------*/
+
+                   DEVICE static double _series (double n, double x)
+                   {                               /* --- series approximation */
+                   int    i;                     /* loop variable */
+                   double t, sum;                /* buffers */
+
+                   sum = t = 1/n;                /* compute initial values */
+                   for (i = MAXITER; --i >= 0; ) {
+                   sum += t *= x/++n;          /* add one term of the series */
+                   if (fabs(t) < fabs(sum) *EPSILON) break;
+                   }                             /* if term is small enough, abort */
+                   return sum;                   /* return the computed factor */
+                   }  /* _series() */
+
+                   /*----------------------------------------------------------------------
+                   series approximation:
+                   P(a,x) =    \gamma(a,x)/\Gamma(a)
+                   \gamma(a,x) = e^-x x^a \sum_{n=0}^\infty (\Gamma(a)/\Gamma(a+1+n)) x^n
+
+                   Source: W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery
+                   Numerical Recipes in C - The Art of Scientific Computing
+                   Cambridge University Press, Cambridge, United Kingdom 1992
+                   formula: pp. 216-219
+
+                   The factor exp(n *log(x) -x) is added in the functions below.
+                   ----------------------------------------------------------------------*/
+
+                   DEVICE static double _cfrac (double n, double x)
+                   {                               /* --- continued fraction approx. */
+                   int    i;                     /* loop variable */
+                   double a, b, c, d, e, f;      /* buffers */
+
+                   b = x+1-n; c = 1/TINY; f = d = 1/b;
+                   for (i = 1; i < MAXITER; i++) {
+                   a = i*(n-i);                /* use Lentz's algorithm to compute */
+                   d = a *d +(b += 2);         /* consecutive approximations */
+                   if (fabs(d) < TINY) d = TINY;
+                   c = b +a/c;
+                   if (fabs(c) < TINY) c = TINY;
+                   d = 1/d; f *= e = d *c;
+                   if (fabs(e-1) < EPSILON) break;
+                   }                             /* if factor is small enough, abort */
+                   return f;                     /* return the computed factor */
+                   }  /* _cfrac() */
+
+                   /*----------------------------------------------------------------------
+                   continued fraction approximation:
+                   P(a,x) = 1 -\Gamma(a,x)/\Gamma(a)
+                   \Gamma(a,x) = e^-x x^a (1/(x+1-a- 1(1-a)/(x+3-a- 2*(2-a)/(x+5-a- ...))))
+
+                   Source: W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery
+                   Numerical Recipes in C - The Art of Scientific Computing
+                   Cambridge University Press, Cambridge, United Kingdom 1992
+                   formula:           pp. 216-219
+                   Lentz's algorithm: p.  171
+
+                   The factor exp(n *log(x) -x) is added in the functions below.
+                   ----------------------------------------------------------------------*/
+
+                   DEVICE double lowerGamma (double n, double x)
+                   {                               /* --- lower incomplete Gamma fn. */
+                   assert((n > 0) && (x > 0));   /* check the function arguments */
+                   return _series(n, x) *exp(n *log(x) -x);
+                   }  /* lowerGamma() */
+
+                   /*--------------------------------------------------------------------*/
+
+                   DEVICE double upperGamma (double n, double x)
+                   {                               /* --- upper incomplete Gamma fn. */
+                   assert((n > 0) && (x > 0));   /* check the function arguments */
+                   return _cfrac(n, x) *exp(n *log(x) -x);
+                   }  /* upperGamma() */
+
+                   /*--------------------------------------------------------------------*/
+
+                   DEVICE double GammaP (double n, double x)
+                   {                               /* --- regularized Gamma function P */
+                   assert((n > 0) && (x >= 0));  /* check the function arguments */
+                   if (x <=  0) return 0;        /* treat x = 0 as a special case */
+                   if (x < n+1) return _series(n, x) *exp(n *log(x) -x -logGamma(n));
+                   return 1 -_cfrac(n, x) *exp(n *log(x) -x -logGamma(n));
+                   }  /* GammaP() */
+
+
+                   //ebuchman: this function is equivalent to scipy.stats.chi2.sf
+                   //it's the pvalue (survival function) of a chi2 distribution
+                   DEVICE double Chi2SF (double k, double x)
+                   {
+                   return 1 - GammaP(k/2., x/2.);
+                   }
+                   #endif
+                   """)
+
+    def c_code(self, node, name, inp, out, sub):
+        x, k = inp
+        z, = out
+        if node.inputs[0].type in float_types:
+            dtype = 'npy_' + node.outputs[0].dtype
+            return """%(z)s =
+                (%(dtype)s)Chi2SF(%(k)s, %(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))
+
+chi2sf = Chi2SF(upgrade_to_float, name='chi2sf')

theano/tensor/basic.py

@@ -3260,6 +3260,11 @@ def gammaln(a):
 def psi(a):
     """derivative of log gamma function"""
 
+@_scal_elemwise
+def chi2sf(x, k):
+    """chi squared survival function"""
+
+
 
 @_scal_elemwise_with_nfunc('real', 1, -1)
 def real(z):

theano/tensor/inplace.py

@@ -229,6 +229,10 @@ def gammaln_inplace(a):
 def psi_inplace(a):
     """derivative of log gamma function"""
 
+@_scal_inplace
+def chi2sf_inplace(x, k):
+    """chi squared survival function"""
+
 @_scal_inplace
 def second_inplace(a):
     """Fill `a` with `b`"""

theano/tensor/tests/test_basic.py

@@ -52,6 +52,7 @@ imported_scipy_special = False
 mode_no_scipy = get_default_mode()
 try:
     import scipy.special
+    import scipy.stats
     imported_scipy_special = True
 except ImportError:
     if config.mode == "FAST_COMPILE":
@@ -68,14 +69,15 @@ utt.seed_rng()
 
 
 def inplace_func(inputs, outputs, mode=None, allow_input_downcast=False,
-        on_unused_input='raise'):
+                 on_unused_input='raise', name=None):
     if mode is None:
         mode = get_default_mode()
     return function(inputs, outputs,
             mode=mode,
             allow_input_downcast=allow_input_downcast,
             accept_inplace=True,
-            on_unused_input=on_unused_input)
+            on_unused_input=on_unused_input,
+            name=name)
 
 
 def eval_outputs(outputs):
@@ -291,7 +293,7 @@ def makeTester(name, op, expected, checks=None, good=None, bad_build=None,
                     raise
 
                 try:
-                    f = inplace_func(inputrs, node.outputs, mode=mode)
+                    f = inplace_func(inputrs, node.outputs, mode=mode, name='test_good')
                 except Exception, exc:
                     err_msg = ("Test %s::%s: Error occurred while"
                         " trying to make a Function") % (self.op, testname)
@@ -381,7 +383,7 @@ def makeTester(name, op, expected, checks=None, good=None, bad_build=None,
                     raise
 
                 try:
-                    f = inplace_func([], node.outputs, mode=mode)
+                    f = inplace_func([], node.outputs, mode=mode, name="test_bad_runtime")
                 except Exception, exc:
                     err_msg = ("Test %s::%s: Error occurred while trying"
                         " to make a Function") % (self.op, testname)
@@ -535,7 +537,8 @@ _good_broadcast_binary_normal = dict(
     # Disabled as we test the case where we reuse the same output as the
     # first inputs.
     # complex3=(rand(2,3),randcomplex(2,3)),
-    empty=(numpy.asarray([]), numpy.asarray([1])),
+    empty=(numpy.asarray([], dtype=config.floatX),
+           numpy.asarray([1], dtype=config.floatX)),
     )
 
 _bad_build_broadcast_binary_normal = dict()
@@ -739,7 +742,8 @@ if PY3:
 else:
     _good_broadcast_div_mod_normal_float_inplace = copymod(
         _good_broadcast_div_mod_normal_float_no_complex,
-        empty1=(numpy.asarray([]), numpy.asarray([1])),
+        empty1=(numpy.asarray([], dtype=config.floatX),
+                numpy.asarray([1], dtype=config.floatX)),
         complex1=(randcomplex(2, 3), randcomplex_nonzero((2, 3))),
         complex2=(randcomplex(2, 3), rand_nonzero((2, 3))),
         # Inplace on the first element. Must have the same type.
@@ -748,7 +752,8 @@ else:
 
 _good_broadcast_div_mod_normal_float = copymod(
     _good_broadcast_div_mod_normal_float_inplace,
-    empty2=(numpy.asarray([0]), numpy.asarray([]))
+    empty2=(numpy.asarray([0], dtype=config.floatX),
+            numpy.asarray([], dtype=config.floatX))
     )
 
 
@@ -841,8 +846,13 @@ _good_broadcast_pow_normal_float = dict(same_shapes = (rand_ranged(1, 5, (2, 3))
                                         complex1 = (randcomplex(2,3),randcomplex(2,3)),
                                         complex2 = (randcomplex(2,3),rand(2,3)),
                                         #complex3 = (rand(2,3),randcomplex(2,3)), # Inplace on the first element.
-                                        empty1 = (numpy.asarray([]), numpy.asarray([1])),
-                                        empty2 = (numpy.asarray([0]), numpy.asarray([])),)
+                                        empty1 = (numpy.asarray([], dtype=config.floatX),
+                                                  numpy.asarray([1], dtype=config.floatX)),
+                                        empty2 = (numpy.asarray([0], dtype=config.floatX),
+                                                  numpy.asarray([], dtype=config.floatX)),
+                                        empty3 = (numpy.asarray([], dtype=config.floatX),
+                                                  numpy.asarray([], dtype=config.floatX)),
+)
 _grad_broadcast_pow_normal = dict(same_shapes = (rand_ranged(1, 5, (2, 3)), rand_ranged(-3, 3, (2, 3))),
                                   scalar = (rand_ranged(1, 5, (2, 3)), rand_ranged(-3, 3, (1, 1))),
                                   row = (
@@ -888,7 +898,7 @@ _good_broadcast_unary_normal_float = dict(
         normal=[rand_ranged(-5, 5, (2, 3))],
         corner_case=[corner_case],
         complex=[randcomplex(2, 3)],
-        empty=[numpy.asarray([])])
+        empty=[numpy.asarray([], dtype=config.floatX)])
 
 _good_broadcast_unary_normal_float_no_empty = copymod(
         _good_broadcast_unary_normal_float,
@@ -908,14 +918,14 @@ _good_broadcast_unary_normal = dict(
         integers=[randint_ranged(-5, 5, (2, 3))],
         corner_case=[corner_case],
         complex=[randcomplex(2, 3)],
-        empty=[numpy.asarray([])],
+        empty=[numpy.asarray([], dtype=config.floatX)],
         )
 
 _good_broadcast_unary_normal_no_complex = dict(
         normal=[numpy.asarray(rand_ranged(-5, 5, (2, 3)), dtype=floatX)],
         integers=[randint_ranged(-5, 5, (2, 3))],
         corner_case=[corner_case],
-        empty=[numpy.asarray([])],
+        empty=[numpy.asarray([], dtype=config.floatX)],
         )
 
 _grad_broadcast_unary_normal_no_complex = dict(
@@ -1122,7 +1132,7 @@ Expm1InplaceTester = makeBroadcastTester(op=inplace.expm1_inplace,
 _good_broadcast_unary_positive = dict(normal=(rand_ranged(0.001, 5, (2, 3)),),
                                       integers=(randint_ranged(1, 5, (2, 3)),),
                                       complex=(randc128_ranged(1, 5, (2, 3)),),
-                                      empty=(numpy.asarray([]),),
+                                      empty=(numpy.asarray([], dtype=config.floatX),),
                                       )
 
 _grad_broadcast_unary_positive = dict(normal=(rand_ranged(0.001, 5, (2, 3)),),)
@@ -1181,7 +1191,7 @@ _good_broadcast_unary_wide = dict(
     normal=(rand_ranged(-1000, 1000, (2, 3)),),
     integers=(randint_ranged(-1000, 1000, (2, 3)),),
     complex=(randc128_ranged(-1000, 1000, (2, 3)),),
-    empty=(numpy.asarray([]),),)
+    empty=(numpy.asarray([], dtype=config.floatX),),)
 _grad_broadcast_unary_wide = dict(normal=(rand_ranged(-1000, 1000, (2, 3)),),)
 
 if theano.config.floatX == 'float32':
@@ -1230,7 +1240,7 @@ SinInplaceTester = makeBroadcastTester(op=inplace.sin_inplace,
 _good_broadcast_unary_arcsin = dict(normal=(rand_ranged(-1, 1, (2, 3)),),
                                     integers=(randint_ranged(-1, 1, (2, 3)),),
                                     complex=(randc128_ranged(-1, 1, (2, 3)),),
-                                    empty=(numpy.asarray([]),),)
+                                    empty=(numpy.asarray([], dtype=config.floatX),),)
 _grad_broadcast_unary_arcsin = dict(normal=(rand_ranged(-1, 1, (2, 3)),),)
 
 ArcsinTester = makeBroadcastTester(op=tensor.arcsin,
@@ -1268,7 +1278,7 @@ _good_broadcast_unary_tan = dict(
     shifted=(rand_ranged(3.15, 6.28, (2, 3)),),
     integers=(randint_ranged(-3, 3, (2, 3)),),
     complex=(randc128_ranged(-3.14, 3.14, (2, 3)),),
-    empty=(numpy.asarray([]),),)
+    empty=(numpy.asarray([], dtype=config.floatX),),)
 #We do not want to test around the discontinuity.
 _grad_broadcast_unary_tan = dict(normal=(rand_ranged(-1.5, 1.5, (2, 3)),),
                                  shifted=(rand_ranged(1.6, 4.6, (2, 3)),))
@@ -1303,7 +1313,8 @@ _good_broadcast_binary_arctan2 = dict(
     integers=(randint(2, 3), randint(2, 3)),
     dtype_mixup_1=(rand(2, 3), randint(2, 3)),
     dtype_mixup_2=(randint(2, 3), rand(2, 3)),
-    empty=(numpy.asarray([]), numpy.asarray([1])),
+    empty=(numpy.asarray([], dtype=config.floatX),
+           numpy.asarray([1], dtype=config.floatX)),
     )
 
 _grad_broadcast_binary_arctan2 = dict(
@@ -1337,7 +1348,7 @@ _good_broadcast_unary_arccosh = dict(
     normal=(rand_ranged(1, 1000, (2, 3)),),
     integers=(randint_ranged(1, 1000, (2, 3)),),
     complex=(randc128_ranged(1, 1000, (2, 3)),),
-    empty=(numpy.asarray([]),),)
+    empty=(numpy.asarray([], dtype=config.floatX),),)
 _grad_broadcast_unary_arccosh = dict(normal=(rand_ranged(1, 1000, (2, 3)),),)
 
 ArccoshTester = makeBroadcastTester(op=tensor.arccosh,
@@ -1385,7 +1396,7 @@ _good_broadcast_unary_arctanh = dict(
     normal=(rand_ranged(-1 + _eps, 1 - _eps, (2, 3)),),
     integers=(randint_ranged(-1 + _eps, 1 - _eps, (2, 3)),),
     complex=(randc128_ranged(-1 + _eps, 1 - _eps, (2, 3)),),
-    empty=(numpy.asarray([]),),)
+    empty=(numpy.asarray([], dtype=config.floatX),),)
 _grad_broadcast_unary_arctanh = dict(
     normal=(rand_ranged(-1 + _eps, 1 - _eps, (2, 3)),),)
 
@@ -1419,6 +1430,7 @@ if imported_scipy_special:
     expected_gamma = scipy.special.gamma
     expected_gammaln = scipy.special.gammaln
     expected_psi = scipy.special.psi
+    expected_chi2sf = lambda x, df: scipy.stats.chi2.sf(x, df).astype(x.dtype)
     skip_scipy = False
 else:
     expected_erf = []
@@ -1428,6 +1440,7 @@ else:
     expected_gamma = []
     expected_gammaln = []
     expected_psi = []
+    expected_chi2sf = []
     skip_scipy = "scipy is not present"
 
 ErfTester = makeBroadcastTester(
@@ -1486,7 +1499,7 @@ ErfcinvTester = makeBroadcastTester(
 
 _good_broadcast_unary_gammaln = dict(
     normal=(rand_ranged(-1 + 1e-2, 10, (2, 3)),),
-    empty=(numpy.asarray([]),),)
+    empty=(numpy.asarray([], dtype=config.floatX),),)
 _grad_broadcast_unary_gammaln = dict(
     # smaller range as our grad method does not estimate it well enough.
     normal=(rand_ranged(1e-8, 8, (2, 3)),),)
@@ -1529,7 +1542,7 @@ GammalnInplaceTester = makeBroadcastTester(
 
 _good_broadcast_unary_psi = dict(
     normal=(rand_ranged(1, 10, (2, 3)),),
-    empty=(numpy.asarray([]),),)
+    empty=(numpy.asarray([], dtype=config.floatX),),)
 
 PsiTester = makeBroadcastTester(
     op=tensor.psi,
@@ -1547,6 +1560,33 @@ PsiInplaceTester = makeBroadcastTester(
     inplace=True,
     skip=skip_scipy)
 
+
+#chi2sf takes two inputs, a value (x) and a degrees of freedom (k).
+# not sure how to deal with that here...
+
+_good_broadcast_unary_chi2sf = dict(
+    normal=(rand_ranged(1, 10, (2, 3)), numpy.asarray(1, dtype=config.floatX)),
+    empty=(numpy.asarray([], dtype=config.floatX),
+           numpy.asarray(1, dtype=config.floatX)))
+
+Chi2SFTester = makeBroadcastTester(
+    op=tensor.chi2sf,
+    expected=expected_chi2sf,
+    good=_good_broadcast_unary_chi2sf,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    skip=skip_scipy,
+    name='Chi2SF')
+Chi2SFInplaceTester = makeBroadcastTester(
+    op=inplace.chi2sf_inplace,
+    expected=expected_chi2sf,
+    good=_good_broadcast_unary_chi2sf,
+    eps=2e-10,
+    mode=mode_no_scipy,
+    inplace=True,
+    skip=skip_scipy,
+    name='Chi2SF')
+
 ZerosLikeTester = makeBroadcastTester(
         op=tensor.zeros_like,
         expected=numpy.zeros_like,
@@ -1569,7 +1609,8 @@ _good_complex_from_polar = dict(
     row=(abs(rand(2, 3)), rand(1, 3)),
     column=(abs(rand(2, 3)), rand(2, 1)),
     integers=(abs(randint(2, 3)), randint(2, 3)),
-    empty=(numpy.asarray([]), numpy.asarray([1])),)
+    empty=(numpy.asarray([], dtype=config.floatX),
+           numpy.asarray([1], dtype=config.floatX)),)
 _grad_complex_from_polar = dict(
     same_shapes=(abs(rand(2, 3)), rand(2, 3)),
     scalar=(abs(rand(2, 3)), rand(1, 1)),
@@ -1608,8 +1649,8 @@ DotTester = makeTester(name='DotTester',
                                               randcomplex(7)),
                                     complex2=(rand(5, 7), randcomplex(7)),
                                     complex3=(randcomplex(5, 7), rand(7)),
-                                    empty1=(numpy.asarray([]),
-                                            numpy.asarray([])),
+                                    empty1=(numpy.asarray([], dtype=config.floatX),
+                                            numpy.asarray([], dtype=config.floatX)),
                                     empty2=(rand(5, 0), rand(0, 2)),
                                     empty3=(rand(0, 5), rand(5, 0)),
                                     ),