[Scipy-svn] r4171 - in trunk/scipy/special: . c_misc

scipy-svn@scip... scipy-svn@scip...
Thu Apr 24 00:03:13 CDT 2008


Author: cookedm
Date: 2008-04-24 00:03:09 -0500 (Thu, 24 Apr 2008)
New Revision: 4171

Added:
   trunk/scipy/special/c_misc/fsolve.c
   trunk/scipy/special/c_misc/gammaincinv.c
Modified:
   trunk/scipy/special/_cephesmodule.c
   trunk/scipy/special/basic.py
   trunk/scipy/special/c_misc/misc.h
   trunk/scipy/special/setup.py
Log:
scipy.special: Fix for #299 (gammaincinv gives incorrect answers)
- implement a gammaincinv(a,y) routine in C that inverts gammainc when
  y < 0.25 instead of using gammainccinv(a,1-y).
- add a generic robust root-finding routine using false position for the above


Modified: trunk/scipy/special/_cephesmodule.c
===================================================================
--- trunk/scipy/special/_cephesmodule.c	2008-04-24 01:19:57 UTC (rev 4170)
+++ trunk/scipy/special/_cephesmodule.c	2008-04-24 05:03:09 UTC (rev 4171)
@@ -1,16 +1,15 @@
-
 /* Cephes module version 1.5
  *  This module defines the functions in the cephes and amos libraries as
- *   Numerical python ufunc objects so that they can operate on arbitrary 
+ *   Numerical python ufunc objects so that they can operate on arbitrary
  *   NumPy arrays with broadcasting and typecasting rules implemented.
- *  
+ *
  *  Copyright 1999  Travis E. Oliphant
  * Revisions 2002 (added functions from cdflib)
  */
 
 #include "Python.h"
 #include "numpy/arrayobject.h"
-#include "numpy/ufuncobject.h" 
+#include "numpy/ufuncobject.h"
 #include "ufunc_extras.h"
 #include "abstract.h"
 #include "cephes.h"
@@ -22,10 +21,9 @@
 
 /* Defined in mtherr in the cephes library */
 extern int scipy_special_print_error_messages;
- 
+
 #include "cephes_doc.h"
 
-
 static PyUFuncGenericFunction cephes1_functions[] = { NULL, NULL, };
 static PyUFuncGenericFunction cephes1rc_functions[] = { NULL, NULL, NULL, NULL};
 static PyUFuncGenericFunction cephes1_2_functions[] = { NULL, NULL, NULL, NULL,};
@@ -145,6 +143,8 @@
 static void * igamc_data[] = { (void *)igamc, (void *)igamc, };
 static void * igam_data[] = { (void *)igam, (void *)igam, };
 static void * igami_data[] = { (void *)igami, (void *)igami, };
+static void * gammaincinv_data[] = { (void *)gammaincinv,
+                                     (void *)gammaincinv, };
 
 static void * iv_data[] = { (void *)iv, (void *)iv, (void *)cbesi_wrap, (void *)cbesi_wrap,};
 static void * ive_data[] = { (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, };
@@ -306,7 +306,7 @@
 static char cephes_1c_types[] = { PyArray_CFLOAT, PyArray_CFLOAT, PyArray_CDOUBLE, PyArray_CDOUBLE, };
 
 
-/* Some functions needed from ufunc object, so that Py_complex's aren't being returned 
+/* Some functions needed from ufunc object, so that Py_complex's aren't being returned
 between code possibly compiled with different compilers.
 */
 
@@ -397,7 +397,7 @@
         cephes4a_2_functions[1] = PyUFunc_dddd_dd_As_dddi_dd;
         cephes5_2_functions[0] = PyUFunc_fffff_ff_As_ddddd_dd;
         cephes5_2_functions[1] = PyUFunc_ddddd_dd;
-	
+
 	/* Create function objects for each function call and insert
 	   them in the dictionary */
 	f = PyUFunc_FromFuncAndData(cephes3a_functions, bdtrc_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "bdtrc", bdtrc_doc, 0);
@@ -424,7 +424,7 @@
 	Py_DECREF(f);
 	f = PyUFunc_FromFuncAndData(cephes3_functions, fdtri_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "fdtri", fdtri_doc, 0);
 	PyDict_SetItemString(dictionary, "fdtri", f);
- 	Py_DECREF(f);
+	Py_DECREF(f);
 
 	f = PyUFunc_FromFuncAndData(cephes3_functions, gdtrc_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "gdtrc", gdtrc_doc, 0);
 	PyDict_SetItemString(dictionary, "gdtrc", f);
@@ -545,7 +545,13 @@
 	f = PyUFunc_FromFuncAndData(cephes2_functions, igami_data, cephes_3_types, 2, 2, 1, PyUFunc_None, "gammainccinv", gammainccinv_doc, 0);
 	PyDict_SetItemString(dictionary, "gammainccinv", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes2c_functions, iv_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "iv", iv_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes2_functions, gammaincinv_data,
+                                    cephes_3_types, 2, 2, 1, PyUFunc_None,
+                                    "gammaincinv", gammaincinv_doc, 0);
+	PyDict_SetItemString(dictionary, "gammaincinv", f);
+	Py_DECREF(f);
+
+        f = PyUFunc_FromFuncAndData(cephes2c_functions, iv_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "iv", iv_doc, 0);
 	PyDict_SetItemString(dictionary, "iv", f);
 	Py_DECREF(f);
 	f = PyUFunc_FromFuncAndData(cephes2cp_functions, ive_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "ive", ive_doc, 0);
@@ -580,7 +586,7 @@
 	PyDict_SetItemString(dictionary, "pdtri", f);
 	Py_DECREF(f);
         /*  Use the student t library from cdflib (it supports doubles for
-              degrees of freedom 
+              degrees of freedom
 	f = PyUFunc_FromFuncAndData(cephes2a_functions, stdtr_data, cephes_3_types, 2, 2, 1, PyUFunc_None, "stdtr", stdtr_doc, 0);
 	PyDict_SetItemString(dictionary, "stdtr", f);
 	Py_DECREF(f);
@@ -812,7 +818,7 @@
 	PyDict_SetItemString(dictionary, "kolmogi", f);
 	Py_DECREF(f);
 
-	f = PyUFunc_FromFuncAndData(cephes1c_functions, wofz_data, cephes_1c_types, 2, 1, 1, PyUFunc_None, "wofz", wofz_doc, 0); 
+	f = PyUFunc_FromFuncAndData(cephes1c_functions, wofz_data, cephes_1c_types, 2, 1, 1, PyUFunc_None, "wofz", wofz_doc, 0);
 	PyDict_SetItemString(dictionary, "wofz", f);
 	Py_DECREF(f);
 
@@ -858,14 +864,14 @@
 	Py_DECREF(f);
 	f = PyUFunc_FromFuncAndData(cephes3_functions, cdff2_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "fdtrix", fdtri_doc, 0);
 	PyDict_SetItemString(dictionary, "fdtrix", f);
- 	Py_DECREF(f);
+	Py_DECREF(f);
         */
-        
+
         /*  The Fortran code for this one seems not to be working properly.
 	f = PyUFunc_FromFuncAndData(cephes3_functions, cdff3_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "fdtridfn", "", 0);
 	PyDict_SetItemString(dictionary, "fdtridfn", f);
 	Py_DECREF(f);
-        */ 
+        */
 	f = PyUFunc_FromFuncAndData(cephes3_functions, cdff4_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "fdtridfd", "", 0);
 	PyDict_SetItemString(dictionary, "fdtridfd", f);
 	Py_DECREF(f);
@@ -947,7 +953,7 @@
 	PyDict_SetItemString(dictionary, "tklmbda", f);
 	Py_DECREF(f);
 
-        
+
 	f = PyUFunc_FromFuncAndData(cephes2_functions, mathieu_a_data, cephes_3_types, 2, 2, 1, PyUFunc_None, "mathieu_a", mathieu_a_doc, 0);
 	PyDict_SetItemString(dictionary, "mathieu_a", f);
 	Py_DECREF(f);
@@ -1059,14 +1065,14 @@
   int oldflag = 0;
   if (!PyArg_ParseTuple ( args, "|i;cephes.errprint", &inflag)) return NULL;
 
-  oldflag = scipy_special_print_error_messages;  
+  oldflag = scipy_special_print_error_messages;
   if (inflag != -37) {
     scipy_special_print_error_messages = (inflag != 0);
   }
   return PyInt_FromLong((long) oldflag);
 }
 
-  
+
 static struct PyMethodDef methods[] = {
   {"errprint", errprint_func, METH_VARARGS, errprint_doc},
   {NULL,		NULL, 0}		/* sentinel */
@@ -1075,9 +1081,9 @@
 
 PyMODINIT_FUNC init_cephes(void) {
   PyObject *m, *d, *s;
-  
+
   /* Create the module and add the functions */
-  m = Py_InitModule("_cephes", methods); 
+  m = Py_InitModule("_cephes", methods);
 
   /* Import the ufunc objects */
   import_array();
@@ -1094,10 +1100,9 @@
   /*  No, instead acessible through errprint */
 
   /* Load the cephes operators into the array module's namespace */
-  Cephes_InitOperators(d); 
-  
+  Cephes_InitOperators(d);
+
   /* Check for errors */
   if (PyErr_Occurred())
     Py_FatalError("can't initialize module _cephes");
 }
-

Modified: trunk/scipy/special/basic.py
===================================================================
--- trunk/scipy/special/basic.py	2008-04-24 01:19:57 UTC (rev 4170)
+++ trunk/scipy/special/basic.py	2008-04-24 05:03:09 UTC (rev 4171)
@@ -394,12 +394,6 @@
         raise ValueError, "Argument must be positive scalar integer."
     return specfun.fcszo(2,nt), specfun.fcszo(1,nt)
 
-def gammaincinv(a,y):
-    """returns the inverse of the incomplete gamma integral in that it
-    finds x such that gammainc(a,x)=y
-    """
-    return gammainccinv(a,1-y)
-
 def hyp0f1(v,z):
     """Confluent hypergeometric limit function 0F1.
     Limit as q->infinity of 1F1(q;a;z/q)

Added: trunk/scipy/special/c_misc/fsolve.c
===================================================================
--- trunk/scipy/special/c_misc/fsolve.c	2008-04-24 01:19:57 UTC (rev 4170)
+++ trunk/scipy/special/c_misc/fsolve.c	2008-04-24 05:03:09 UTC (rev 4171)
@@ -0,0 +1,159 @@
+#include "misc.h"
+#include <math.h>
+
+static inline double
+max(double a, double b)
+{
+    return (a > b ? a : b);
+}
+
+/*
+   Use a combination of bisection and false position to find a root
+   of a function within a given interval. This is guaranteed to converge,
+   and always keeps a bounding interval, unlike Newton's method.
+
+   The false position steps are either unmodified, or modified with
+   the Anderson-Bjorck method as appropiate. Theoretically, this has
+   a "speed of convergence" of 1.7 (bisection is 1, Newton is 2).
+
+   Input
+   -----
+    a, b:   initial bounding interval
+    fa, fb: value of f() at a and b
+    f, f_extra: function to find root of is f(x, f_extra)
+    abserr, relerr: absolute and relative errors on the bounding interval
+    bisect_til: If > 0.0, perform bisection until the width of the
+                bounding interval is less than this.
+
+   Output
+   ------
+    a, b, fa, fb: Final bounding interval and function values
+    best_x, best_f: Best root approximation and the function value there
+
+   Returns
+   -------
+    FSOLVE_CONVERGED: Bounding interval is smaller than required error.
+    FSOLVE_NOT_BRACKET: Initial interval is not a bounding interval.
+    FSOLVE_EXACT: An exact root was found (best_f = 0)
+
+
+   Note that this routine was designed initially to work with gammaincinv, so
+   it may not be tuned right for other problems. Don't use it blindly.
+ */
+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 x1=*a, f1=*fa, x2=*b, f2=*fb;
+    fsolve_result_t r = FSOLVE_CONVERGED;
+    double gamma = 1.0;
+    enum {bisect, falsep} state = bisect;
+    int n_falsep = 0;
+    double x3, f3;
+    double w, last_bisect_width;
+    double tol;
+
+    if (f1*f2 >= 0.0) {
+        return FSOLVE_NOT_BRACKET;
+    }
+    if (bisect_til > 0.0) {
+        state = bisect;
+    } else {
+        state = falsep;
+    }
+    w = fabs(x2 - x1);
+    last_bisect_width = w;
+    while (1) {
+        switch (state) {
+        case bisect: {
+            x3 = 0.5 * (x1 + x2);
+            if (x3 == x1 || x3 == x2) {
+                /* i.e., x1 and x2 are successive floating-point numbers. */
+                *best_x = x3;
+                *best_f = (x3==x1) ? f1 : f2;
+                goto finish;
+            }
+            f3 = f(x3, f_extra);
+            if (f3 == 0.0) {
+                goto exact_soln;
+            }
+            if (f3*f2 < 0.0) {
+                x1 = x2; f1 = f2;
+            }
+            x2 = x3; f2 = f3;
+            w = fabs(x2 - x1);
+            last_bisect_width = w;
+            if (bisect_til > 0.0) {
+                if (w < bisect_til) {
+                    bisect_til = -1.0;
+                    gamma = 1.0;
+                    n_falsep = 0;
+                    state = falsep;
+                }
+            } else {
+                gamma = 1.0;
+                n_falsep = 0;
+                state = falsep;
+            }
+            break;
+        }
+        case falsep: {
+            double s12 = (f2 - gamma*f1) / (x2 - x1);
+            x3 = x2 - f2/s12;
+            f3 = f(x3, f_extra);
+            if (f3 == 0.0) {
+                goto exact_soln;
+            }
+            n_falsep += 1;
+            if (f3*f2 < 0.0) {
+                gamma = 1.0;
+                x1 = x2; f1 = f2;
+            } else {
+                /* Anderson-Bjorck method */
+                double g = 1.0 - f3 / f2;
+                if (g <= 0.0) { g = 0.5; }
+                /* It's not really clear from the sources I've looked at,
+                   but I believe this is *= instead of =. */
+                gamma *= g;
+            }
+            x2 = x3; f2 = f3;
+            w = fabs(x2 - x1);
+            /* Sanity check. For every 4 false position checks, see if we
+               really are decreasing the interval by comparing to what
+               bisection would have achieved (or, rather, a bit more lenient
+               than that -- interval decreased by 4 instead of by 16, as
+               the fp could be decreasing gamma for a bit).
+
+               Note that this should guarantee convergence, as it makes
+               sure that we always end up decreasing the interval width
+               with a bisection.
+             */
+            if (n_falsep > 4) {
+                if (w*4 > last_bisect_width) {
+                    state = bisect;
+                }
+                n_falsep = 0;
+                last_bisect_width = w;
+            }
+            break;
+        }
+        }
+        tol = abserr + relerr*max(max(fabs(x1), fabs(x2)), 1.0);
+        if (w <= tol) {
+            if (fabs(f1) < fabs(f2)) {
+                *best_x = x1; *best_f = f1;
+            } else {
+                *best_x = x2; *best_f = f2;
+            }
+            goto finish;
+        }
+    }
+exact_soln:
+    *best_x = x3; *best_f = 0.0;
+    r = FSOLVE_EXACT;
+finish:
+    *a = x1; *fa = f1; *b = x2; *fb = f2;
+    return r;
+}

Added: trunk/scipy/special/c_misc/gammaincinv.c
===================================================================
--- trunk/scipy/special/c_misc/gammaincinv.c	2008-04-24 01:19:57 UTC (rev 4170)
+++ trunk/scipy/special/c_misc/gammaincinv.c	2008-04-24 05:03:09 UTC (rev 4171)
@@ -0,0 +1,52 @@
+#include <stdio.h>
+#include <math.h>
+#include "../cephes.h"
+#undef fabs
+#include "misc.h"
+
+/*
+  Inverse of the (regularised) incomplete Gamma integral.
+
+  Given a, find x such that igam(a, x) = y.
+  For y not small, we just use igami(a, 1-y) (inverse of the complemented
+  incomplete Gamma integral). For y small, however, 1-y is about 1, and we
+  lose digits.
+
+*/
+
+extern double MACHEP;
+
+static double
+gammainc(double x, double params[2])
+{
+    return cephes_igam(params[0], x) - params[1];
+}
+
+double
+gammaincinv(double a, double y)
+{
+    if (a <= 0.0 || y <= 0.0 || y > 0.25) {
+        return cephes_igami(a, 1-y);
+    }
+
+    /* I found Newton to be unreliable. Also, after we generate a small
+       interval by bisection above, false position will do a large step
+       from an interval of width ~1e-4 to ~1e-14 in one step (a=10, x=0.05,
+       but similiar for other values).
+     */
+
+    double lo = 0.0, hi = cephes_igami(a, 0.75);
+    double flo = -y, fhi = 0.25 - y;
+    double params[2] = {a, y};
+    double best_x, best_f;
+    fsolve_result_t r;
+
+    r = false_position(&lo, &flo, &hi, &fhi,
+                       (objective_function)gammainc, params,
+                       MACHEP, MACHEP, 1e-2*a,
+                       &best_x, &best_f);
+    if (r == FSOLVE_NOT_BRACKET) {
+        best_x = 0.0;
+    }
+    return best_x;
+}

Modified: trunk/scipy/special/c_misc/misc.h
===================================================================
--- trunk/scipy/special/c_misc/misc.h	2008-04-24 01:19:57 UTC (rev 4170)
+++ trunk/scipy/special/c_misc/misc.h	2008-04-24 05:03:09 UTC (rev 4171)
@@ -1,6 +1,26 @@
+#ifndef C_MISC_MISC_H
+#define C_MISC_MISC_H
 
-double besselpoly(double a, double lambda, double nu);
+typedef enum {
+  /* An exact solution was found, in which case the first point
+     on the interval is the value */
+  FSOLVE_EXACT,
+  /* Interval width is less than the tolerance */
+  FSOLVE_CONVERGED,
+  /* Not a bracket */
+  FSOLVE_NOT_BRACKET,
+} fsolve_result_t;
 
+typedef double (*objective_function)(double, void *);
 
+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 besselpoly(double a, double lambda, double nu);
+double gammaincinv(double a, double x);
 
+#define gammaincinv_doc """gammaincinv(a, y) returns x such that gammainc(a, x) = y."""
+
+#endif /* C_MISC_MISC_H */

Modified: trunk/scipy/special/setup.py
===================================================================
--- trunk/scipy/special/setup.py	2008-04-24 01:19:57 UTC (rev 4170)
+++ trunk/scipy/special/setup.py	2008-04-24 05:03:09 UTC (rev 4171)
@@ -42,8 +42,8 @@
                                   "cephes/mconf.h", "cephes/cephes_names.h"],
                          define_macros = define_macros
                          )
+
     # Extension specfun
-
     config.add_extension('specfun',
                          sources=['specfun.pyf'],
                          f2py_options=['--no-wrap-functions'],



More information about the Scipy-svn mailing list