Skip to content

P
pytensor
提交 426686fe authored 作者: ebuchman's avatar ebuchman
浏览文件
操作

Merge pull request #1 from nouiz/ebuchman-chi2sf/master

Finish chi2 tests and GPU code.
上级 08d16c0a 914eb1be
隐藏空白字符变更
内嵌 并排
正在显示 包含 107 行增加 和 93 行删除
theano/scalar/basic_scipy.py
浏览文件 @ 426686fe
......@@ -289,20 +289,23 @@ DEVICE double _psi(double x){
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(
"""
......@@ -312,10 +315,10 @@ class Chi2SF(BinaryScalarOp):
#else
#define DEVICE
#endif
#ifndef _CHI2FUNCDEFINED
#define _CHI2FUNCDEFINED
/*----------------------------------------------------------------------
File : gamma.c
Contents: computation of the (incomplete/regularized) gamma function
......@@ -328,8 +331,8 @@ class Chi2SF(BinaryScalarOp):
----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <float.h>
#include <math.h>
......@@ -346,23 +349,23 @@ class Chi2SF(BinaryScalarOp):
#define MAXFACT 170
#define MAXITER 1024
#define TINY (EPSILON *EPSILON *EPSILON)
/*----------------------------------------------------------------------
Table of Factorials/Gamma Values
----------------------------------------------------------------------*/
static double _facts[MAXFACT+1] = { 0 };
static double _logfs[MAXFACT+1];
static double _halfs[MAXFACT+1];
static double _loghs[MAXFACT+1];
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
----------------------------------------------------------------------*/
static void _init (void)
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; ) {
......@@ -376,14 +379,14 @@ class Chi2SF(BinaryScalarOp):
_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) {
......@@ -401,13 +404,13 @@ class Chi2SF(BinaryScalarOp):
- 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 /*--------------------------------------------------------------*/
double logGamma (double n)
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) {
......@@ -427,7 +430,7 @@ class Chi2SF(BinaryScalarOp):
+ 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
......@@ -437,20 +440,20 @@ class Chi2SF(BinaryScalarOp):
* (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.
----------------------------------------------------------------------*/
double Gamma (double n)
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 */
......@@ -462,14 +465,14 @@ class Chi2SF(BinaryScalarOp):
} /* try to get the value from a table */
return exp(logGamma(n)); /* compute through natural logarithm */
} /* Gamma() */
/*--------------------------------------------------------------------*/
static double _series (double n, double x)
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 */
......@@ -477,25 +480,25 @@ class Chi2SF(BinaryScalarOp):
} /* 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.
----------------------------------------------------------------------*/
static double _cfrac (double n, double x)
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 */
......@@ -508,68 +511,67 @@ class Chi2SF(BinaryScalarOp):
} /* 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.
----------------------------------------------------------------------*/
double lowerGamma (double n, double x)
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() */
/*--------------------------------------------------------------------*/
double upperGamma (double n, double x)
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() */
/*--------------------------------------------------------------------*/
double GammaP (double n, double x)
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 = z.dtype
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))
......
theano/tensor/tests/test_basic.py
浏览文件 @ 426686fe
......@@ -69,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):
......@@ -292,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)
......@@ -382,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)
......@@ -536,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()
......@@ -740,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.
......@@ -749,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))
)
......@@ -842,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 = (
......@@ -889,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,
......@@ -909,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(
......@@ -1123,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)),),)
......@@ -1182,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':
......@@ -1231,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,
......@@ -1269,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)),))
......@@ -1304,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(
......@@ -1338,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,
......@@ -1386,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)),),)
......@@ -1420,7 +1430,7 @@ if imported_scipy_special:
expected_gamma = scipy.special.gamma
expected_gammaln = scipy.special.gammaln
expected_psi = scipy.special.psi
expected_chi2sf = scipy.stats.chi2.sf
expected_chi2sf = lambda x, df: scipy.stats.chi2.sf(x, df).astype(x.dtype)
skip_scipy = False
else:
expected_erf = []
......@@ -1489,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)),),)
......@@ -1532,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,
......@@ -1551,13 +1561,13 @@ PsiInplaceTester = makeBroadcastTester(
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)),),
empty=(numpy.asarray([]),),)
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,
......@@ -1565,16 +1575,17 @@ Chi2SFTester = makeBroadcastTester(
good=_good_broadcast_unary_chi2sf,
eps=2e-10,
mode=mode_no_scipy,
skip=skip_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)
'''
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,
......@@ -1598,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)),
......@@ -1637,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)),
),
......
Markdown 格式
0% 或
您添加了 0 人 到此讨论。请谨慎行事。
请先完成此评论的编辑！
请 注册 或者 后发表评论