[Scipy-svn] r6521 - in branches/0.8.x/scipy/special: . c_misc

[scipy-svn at scipy.org]

Author: ptvirtan
Date: 2010-06-17 16:40:16 -0500 (Thu, 17 Jun 2010)
New Revision: 6521

Modified:
   branches/0.8.x/scipy/special/c_misc/fsolve.c
   branches/0.8.x/scipy/special/c_misc/gammaincinv.c
   branches/0.8.x/scipy/special/c_misc/misc.h
   branches/0.8.x/scipy/special/setup.py
Log:
BUG: special: make gammaincinv behave more sanely when it thinks the iteration did not converge

(cherry picked from r6519)

Modified: branches/0.8.x/scipy/special/c_misc/fsolve.c
===================================================================
--- branches/0.8.x/scipy/special/c_misc/fsolve.c	2010-06-17 20:36:51 UTC (rev 6520)
+++ branches/0.8.x/scipy/special/c_misc/fsolve.c	2010-06-17 21:40:16 UTC (rev 6521)
@@ -48,7 +48,7 @@
 false_position(double *a, double *fa, double *b, double *fb,
                objective_function f, void *f_extra,
                double abserr, double relerr, double bisect_til,
-               double *best_x, double *best_f)
+               double *best_x, double *best_f, double *errest)
 {
     double x1=*a, f1=*fa, x2=*b, f2=*fb;
     fsolve_result_t r = FSOLVE_CONVERGED;
@@ -163,5 +163,6 @@
     r = FSOLVE_EXACT;
 finish:
     *a = x1; *fa = f1; *b = x2; *fb = f2;
+    *errest = w;
     return r;
 }

Modified: branches/0.8.x/scipy/special/c_misc/gammaincinv.c
===================================================================
--- branches/0.8.x/scipy/special/c_misc/gammaincinv.c	2010-06-17 20:36:51 UTC (rev 6520)
+++ branches/0.8.x/scipy/special/c_misc/gammaincinv.c	2010-06-17 21:40:16 UTC (rev 6521)
@@ -1,9 +1,17 @@
+#include <Python.h>
+#include <numpy/npy_math.h>
+
 #include <stdio.h>
 #include <math.h>
+
 #include "../cephes.h"
 #undef fabs
 #include "misc.h"
 
+/* Limits after which to issue warnings about non-convergence */
+#define ALLOWED_ATOL (1e-306)
+#define ALLOWED_RTOL (1e-9)
+
 void scipy_special_raise_warning(char *fmt, ...);
 
 /*
@@ -16,7 +24,7 @@
 
 */
 
-extern double MACHEP;
+extern double MACHEP, MAXNUM;
 
 static double
 gammainc(double x, double params[2])
@@ -30,7 +38,7 @@
     double lo = 0.0, hi;
     double flo = -y, fhi = 0.25 - y;
     double params[2];
-    double best_x, best_f;
+    double best_x, best_f, errest;
     fsolve_result_t r;
 
     if (a <= 0.0 || y <= 0.0 || y >= 0.25) {
@@ -53,12 +61,13 @@
     r = false_position(&lo, &flo, &hi, &fhi,
                        (objective_function)gammainc, params,
                        2*MACHEP, 2*MACHEP, 1e-2*a,
-                       &best_x, &best_f);
-    if (!(r == FSOLVE_CONVERGED || r == FSOLVE_EXACT)) {
+                       &best_x, &best_f, &errest);
+    if (!(r == FSOLVE_CONVERGED || r == FSOLVE_EXACT) &&
+            errest > ALLOWED_ATOL + ALLOWED_RTOL*fabs(best_x)) {
         scipy_special_raise_warning(
-            "gammaincinv: failed to converge at (a, y) = (%.20g, %.20g): %d\n",
-            a, y, r);
-        best_x = 0.0;
+            "gammaincinv: failed to converge at (a, y) = (%.20g, %.20g): got %g +- %g, code %d\n",
+            a, y, best_x, errest, r);
+        best_x = NPY_NAN;
     }
     return best_x;
 }

Modified: branches/0.8.x/scipy/special/c_misc/misc.h
===================================================================
--- branches/0.8.x/scipy/special/c_misc/misc.h	2010-06-17 20:36:51 UTC (rev 6520)
+++ branches/0.8.x/scipy/special/c_misc/misc.h	2010-06-17 21:40:16 UTC (rev 6521)
@@ -18,7 +18,7 @@
 fsolve_result_t false_position(double *a, double *fa, double *b, double *fb,
                        objective_function f, void *f_extra,
                        double abserr, double relerr, double bisect_til,
-                       double *best_x, double *best_f);
+                       double *best_x, double *best_f, double *errest);
 
 double besselpoly(double a, double lambda, double nu);
 double gammaincinv(double a, double x);

Modified: branches/0.8.x/scipy/special/setup.py
===================================================================
--- branches/0.8.x/scipy/special/setup.py	2010-06-17 20:36:51 UTC (rev 6520)
+++ branches/0.8.x/scipy/special/setup.py	2010-06-17 21:40:16 UTC (rev 6521)
@@ -24,7 +24,9 @@
         define_macros.append(('_USE_MATH_DEFINES',None))
 
     # C libraries
-    config.add_library('sc_c_misc',sources=[join('c_misc','*.c')])
+    config.add_library('sc_c_misc',sources=[join('c_misc','*.c')],
+                       include_dirs=[get_python_inc(), get_numpy_include_dirs()],
+                       macros=define_macros)
     config.add_library('sc_cephes',sources=[join('cephes','*.c')],
                        include_dirs=[get_python_inc(), get_numpy_include_dirs()],
                        macros=define_macros)