提交 e503b3f0 authored 作者: Frederic's avatar Frederic

pep8

上级 08d16c0a
...@@ -289,20 +289,23 @@ DEVICE double _psi(double x){ ...@@ -289,20 +289,23 @@ DEVICE double _psi(double x){
return hash(type(self)) return hash(type(self))
psi = Psi(upgrade_to_float, name='psi') psi = Psi(upgrade_to_float, name='psi')
class Chi2SF(BinaryScalarOp): class Chi2SF(BinaryScalarOp):
""" """
Compute (1 - chi2_cdf(x)) Compute (1 - chi2_cdf(x))
ie. chi2 pvalue (chi2 'survival function') ie. chi2 pvalue (chi2 'survival function')
""" """
@staticmethod @staticmethod
def st_impl(x, k): def st_impl(x, k):
return scipy.stats.chi2.sf(x, k) return scipy.stats.chi2.sf(x, k)
def impl(self, x, k): def impl(self, x, k):
if imported_scipy_special: if imported_scipy_special:
return Chi2SF.st_impl(x, k) return Chi2SF.st_impl(x, k)
else: else:
super(Chi2SF, self).impl(x, k) super(Chi2SF, self).impl(x, k)
def c_support_code(self): def c_support_code(self):
return( return(
""" """
...@@ -312,10 +315,10 @@ class Chi2SF(BinaryScalarOp): ...@@ -312,10 +315,10 @@ class Chi2SF(BinaryScalarOp):
#else #else
#define DEVICE #define DEVICE
#endif #endif
#ifndef _CHI2FUNCDEFINED #ifndef _CHI2FUNCDEFINED
#define _CHI2FUNCDEFINED #define _CHI2FUNCDEFINED
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
File : gamma.c File : gamma.c
Contents: computation of the (incomplete/regularized) gamma function Contents: computation of the (incomplete/regularized) gamma function
...@@ -328,8 +331,8 @@ class Chi2SF(BinaryScalarOp): ...@@ -328,8 +331,8 @@ class Chi2SF(BinaryScalarOp):
----------------------------------------------------------------------*/ ----------------------------------------------------------------------*/
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <assert.h> #include <assert.h>
#include <float.h> #include <float.h>
#include <math.h> #include <math.h>
...@@ -346,7 +349,7 @@ class Chi2SF(BinaryScalarOp): ...@@ -346,7 +349,7 @@ class Chi2SF(BinaryScalarOp):
#define MAXFACT 170 #define MAXFACT 170
#define MAXITER 1024 #define MAXITER 1024
#define TINY (EPSILON *EPSILON *EPSILON) #define TINY (EPSILON *EPSILON *EPSILON)
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
Table of Factorials/Gamma Values Table of Factorials/Gamma Values
----------------------------------------------------------------------*/ ----------------------------------------------------------------------*/
...@@ -354,7 +357,7 @@ class Chi2SF(BinaryScalarOp): ...@@ -354,7 +357,7 @@ class Chi2SF(BinaryScalarOp):
static double _logfs[MAXFACT+1]; static double _logfs[MAXFACT+1];
static double _halfs[MAXFACT+1]; static double _halfs[MAXFACT+1];
static double _loghs[MAXFACT+1]; static double _loghs[MAXFACT+1];
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
Functions Functions
----------------------------------------------------------------------*/ ----------------------------------------------------------------------*/
...@@ -362,7 +365,7 @@ class Chi2SF(BinaryScalarOp): ...@@ -362,7 +365,7 @@ class Chi2SF(BinaryScalarOp):
{ /* --- init. factorial tables */ { /* --- init. factorial tables */
int i; /* loop variable */ int i; /* loop variable */
double x = 1; /* factorial */ double x = 1; /* factorial */
_facts[0] = _facts[1] = 1; /* store factorials for 0 and 1 */ _facts[0] = _facts[1] = 1; /* store factorials for 0 and 1 */
_logfs[0] = _logfs[1] = 0; /* and their logarithms */ _logfs[0] = _logfs[1] = 0; /* and their logarithms */
for (i = 1; ++i <= MAXFACT; ) { for (i = 1; ++i <= MAXFACT; ) {
...@@ -376,14 +379,14 @@ class Chi2SF(BinaryScalarOp): ...@@ -376,14 +379,14 @@ class Chi2SF(BinaryScalarOp):
_loghs[i] = log(x); /* the Gamma function of half numbers */ _loghs[i] = log(x); /* the Gamma function of half numbers */
} /* and the table of their logarithms */ } /* and the table of their logarithms */
} /* _init() */ } /* _init() */
/*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/
#if 0 #if 0
double logGamma (double n) double logGamma (double n)
{ /* --- compute ln(Gamma(n)) */ { /* --- compute ln(Gamma(n)) */
double s; /* = ln((n-1)!), n \in IN */ double s; /* = ln((n-1)!), n \in IN */
assert(n > 0); /* check the function argument */ assert(n > 0); /* check the function argument */
if (_facts[0] <= 0) _init(); /* initialize the tables */ if (_facts[0] <= 0) _init(); /* initialize the tables */
if (n < MAXFACT +1 +4 *EPSILON) { if (n < MAXFACT +1 +4 *EPSILON) {
...@@ -401,13 +404,13 @@ class Chi2SF(BinaryScalarOp): ...@@ -401,13 +404,13 @@ class Chi2SF(BinaryScalarOp):
- 0.5395239384953e-5 /(n+6); - 0.5395239384953e-5 /(n+6);
return (n+0.5) *log((n+5.5)/LN_BASE) +(LN_SQRT_2PI +log(s/n) -5.0); return (n+0.5) *log((n+5.5)/LN_BASE) +(LN_SQRT_2PI +log(s/n) -5.0);
} /* logGamma() */ } /* logGamma() */
#else /*--------------------------------------------------------------*/ #else /*--------------------------------------------------------------*/
double logGamma (double n) double logGamma (double n)
{ /* --- compute ln(Gamma(n)) */ { /* --- compute ln(Gamma(n)) */
double s; /* = ln((n-1)!), n \in IN */ double s; /* = ln((n-1)!), n \in IN */
assert(n > 0); /* check the function argument */ assert(n > 0); /* check the function argument */
if (_facts[0] <= 0) _init(); /* initialize the tables */ if (_facts[0] <= 0) _init(); /* initialize the tables */
if (n < MAXFACT +1 +4 *EPSILON) { if (n < MAXFACT +1 +4 *EPSILON) {
...@@ -427,7 +430,7 @@ class Chi2SF(BinaryScalarOp): ...@@ -427,7 +430,7 @@ class Chi2SF(BinaryScalarOp):
+ 1.50563273514931155834e-7 /(n+8); + 1.50563273514931155834e-7 /(n+8);
return (n+0.5) *log((n+7.5)/LN_BASE) +(LN_SQRT_2PI +log(s/n) -7.0); return (n+0.5) *log((n+7.5)/LN_BASE) +(LN_SQRT_2PI +log(s/n) -7.0);
} /* logGamma() */ } /* logGamma() */
#endif #endif
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
Use Lanczos' approximation Use Lanczos' approximation
...@@ -437,19 +440,19 @@ class Chi2SF(BinaryScalarOp): ...@@ -437,19 +440,19 @@ class Chi2SF(BinaryScalarOp):
* (c_0 +c_1/(n+1) +c_2/(n+2) +...+c_n/(n+k) +\epsilon) * (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, and exploit the recursion \Gamma(n+1) = n *\Gamma(n) once,
i.e., compute \Gamma(n) as \Gamma(n+1) /n. 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 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. 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 Source: W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery
Numerical Recipes in C - The Art of Scientific Computing Numerical Recipes in C - The Art of Scientific Computing
Cambridge University Press, Cambridge, United Kingdom 1992 Cambridge University Press, Cambridge, United Kingdom 1992
pp. 213-214 pp. 213-214
For the choices gamma = 7, k = 8, and c_0 to c_8 as defined 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. in the second version, the value is slightly more accurate.
----------------------------------------------------------------------*/ ----------------------------------------------------------------------*/
double Gamma (double n) double Gamma (double n)
{ /* --- compute Gamma(n) = (n-1)! */ { /* --- compute Gamma(n) = (n-1)! */
assert(n > 0); /* check the function argument */ assert(n > 0); /* check the function argument */
...@@ -462,14 +465,14 @@ class Chi2SF(BinaryScalarOp): ...@@ -462,14 +465,14 @@ class Chi2SF(BinaryScalarOp):
} /* try to get the value from a table */ } /* try to get the value from a table */
return exp(logGamma(n)); /* compute through natural logarithm */ return exp(logGamma(n)); /* compute through natural logarithm */
} /* Gamma() */ } /* Gamma() */
/*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/
static double _series (double n, double x) static double _series (double n, double x)
{ /* --- series approximation */ { /* --- series approximation */
int i; /* loop variable */ int i; /* loop variable */
double t, sum; /* buffers */ double t, sum; /* buffers */
sum = t = 1/n; /* compute initial values */ sum = t = 1/n; /* compute initial values */
for (i = MAXITER; --i >= 0; ) { for (i = MAXITER; --i >= 0; ) {
sum += t *= x/++n; /* add one term of the series */ sum += t *= x/++n; /* add one term of the series */
...@@ -477,25 +480,25 @@ class Chi2SF(BinaryScalarOp): ...@@ -477,25 +480,25 @@ class Chi2SF(BinaryScalarOp):
} /* if term is small enough, abort */ } /* if term is small enough, abort */
return sum; /* return the computed factor */ return sum; /* return the computed factor */
} /* _series() */ } /* _series() */
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
series approximation: series approximation:
P(a,x) = \gamma(a,x)/\Gamma(a) 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 \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 Source: W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery
Numerical Recipes in C - The Art of Scientific Computing Numerical Recipes in C - The Art of Scientific Computing
Cambridge University Press, Cambridge, United Kingdom 1992 Cambridge University Press, Cambridge, United Kingdom 1992
formula: pp. 216-219 formula: pp. 216-219
The factor exp(n *log(x) -x) is added in the functions below. The factor exp(n *log(x) -x) is added in the functions below.
----------------------------------------------------------------------*/ ----------------------------------------------------------------------*/
static double _cfrac (double n, double x) static double _cfrac (double n, double x)
{ /* --- continued fraction approx. */ { /* --- continued fraction approx. */
int i; /* loop variable */ int i; /* loop variable */
double a, b, c, d, e, f; /* buffers */ double a, b, c, d, e, f; /* buffers */
b = x+1-n; c = 1/TINY; f = d = 1/b; b = x+1-n; c = 1/TINY; f = d = 1/b;
for (i = 1; i < MAXITER; i++) { for (i = 1; i < MAXITER; i++) {
a = i*(n-i); /* use Lentz's algorithm to compute */ a = i*(n-i); /* use Lentz's algorithm to compute */
...@@ -508,38 +511,38 @@ class Chi2SF(BinaryScalarOp): ...@@ -508,38 +511,38 @@ class Chi2SF(BinaryScalarOp):
} /* if factor is small enough, abort */ } /* if factor is small enough, abort */
return f; /* return the computed factor */ return f; /* return the computed factor */
} /* _cfrac() */ } /* _cfrac() */
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
continued fraction approximation: continued fraction approximation:
P(a,x) = 1 -\Gamma(a,x)/\Gamma(a) 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- ...)))) \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 Source: W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery
Numerical Recipes in C - The Art of Scientific Computing Numerical Recipes in C - The Art of Scientific Computing
Cambridge University Press, Cambridge, United Kingdom 1992 Cambridge University Press, Cambridge, United Kingdom 1992
formula: pp. 216-219 formula: pp. 216-219
Lentz's algorithm: p. 171 Lentz's algorithm: p. 171
The factor exp(n *log(x) -x) is added in the functions below. The factor exp(n *log(x) -x) is added in the functions below.
----------------------------------------------------------------------*/ ----------------------------------------------------------------------*/
double lowerGamma (double n, double x) double lowerGamma (double n, double x)
{ /* --- lower incomplete Gamma fn. */ { /* --- lower incomplete Gamma fn. */
assert((n > 0) && (x > 0)); /* check the function arguments */ assert((n > 0) && (x > 0)); /* check the function arguments */
return _series(n, x) *exp(n *log(x) -x); return _series(n, x) *exp(n *log(x) -x);
} /* lowerGamma() */ } /* lowerGamma() */
/*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/
double upperGamma (double n, double x) double upperGamma (double n, double x)
{ /* --- upper incomplete Gamma fn. */ { /* --- upper incomplete Gamma fn. */
assert((n > 0) && (x > 0)); /* check the function arguments */ assert((n > 0) && (x > 0)); /* check the function arguments */
return _cfrac(n, x) *exp(n *log(x) -x); return _cfrac(n, x) *exp(n *log(x) -x);
} /* upperGamma() */ } /* upperGamma() */
/*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/
double GammaP (double n, double x) double GammaP (double n, double x)
{ /* --- regularized Gamma function P */ { /* --- regularized Gamma function P */
assert((n > 0) && (x >= 0)); /* check the function arguments */ assert((n > 0) && (x >= 0)); /* check the function arguments */
...@@ -547,8 +550,8 @@ class Chi2SF(BinaryScalarOp): ...@@ -547,8 +550,8 @@ class Chi2SF(BinaryScalarOp):
if (x < n+1) return _series(n, x) *exp(n *log(x) -x -logGamma(n)); 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)); return 1 -_cfrac(n, x) *exp(n *log(x) -x -logGamma(n));
} /* GammaP() */ } /* GammaP() */
//ebuchman: this function is equivalent to scipy.stats.chi2.sf //ebuchman: this function is equivalent to scipy.stats.chi2.sf
//it's the pvalue (survival function) of a chi2 distribution //it's the pvalue (survival function) of a chi2 distribution
DEVICE double Chi2SF (double k, double x) DEVICE double Chi2SF (double k, double x)
...@@ -556,10 +559,8 @@ class Chi2SF(BinaryScalarOp): ...@@ -556,10 +559,8 @@ class Chi2SF(BinaryScalarOp):
return 1 - GammaP(k/2., x/2.); return 1 - GammaP(k/2., x/2.);
} }
""") """)
def c_code(self, node, name, inp, out, sub): def c_code(self, node, name, inp, out, sub):
x, k = inp x, k = inp
z, = out z, = out
if node.inputs[0].type in float_types: if node.inputs[0].type in float_types:
...@@ -567,9 +568,10 @@ class Chi2SF(BinaryScalarOp): ...@@ -567,9 +568,10 @@ class Chi2SF(BinaryScalarOp):
return """%(z)s = return """%(z)s =
(%(dtype)s)Chi2SF(%(k)s, %(x)s);""" % locals() (%(dtype)s)Chi2SF(%(k)s, %(x)s);""" % locals()
raise NotImplementedError('only floatingpoint is implemented') raise NotImplementedError('only floatingpoint is implemented')
def __eq__(self, other): def __eq__(self, other):
return type(self) == type(other) return type(self) == type(other)
def __hash__(self): def __hash__(self):
return hash(type(self)) return hash(type(self))
......
Markdown 格式
0%
您添加了 0 到此讨论。请谨慎行事。
请先完成此评论的编辑!
注册 或者 后发表评论