Skip to content
This repository

Ctypes callback #190

Merged
merged 1 commit into from almost 2 years ago

4 participants

Travis E. Oliphant Ralf Gommers Pauli Virtanen Charles Harris
Travis E. Oliphant
Owner

This Pull Request allows integrate.quad to by-pass the Python-layer when given a Ctypes function pointer of the correct signature. When the callable passed to integrate.quad is a ctypes function pointer, it is checked for the correct signature, and then the raw C-function is handed to the underlying integration routine which does not call-back into Python.

Ralf Gommers
Owner

On the C code changes I can't really comment (except that it looks useful).

The test is going to fail on some systems, because ctypes is not always available. So the import should be conditional, and the test should be decorated with @skipif.

About the timing test, doesn't that belong under benchmarks instead of tests? Each module has a bench() function for this sort of thing.

scipy/integrate/__quadpack.h
((62 lines not shown))
  128
+}
  129
+
  130
+#define CHECK_FUNC(func_type, fcn) { \
  131
+  func_type = get_func_type(fcn); \
  132
+  if (func_type == ERROR) return NULL; \
  133
+  if (func_type == NOTCALLABLE) { \
  134
+      PyErr_SetString(PyExc_TypeError, "quad: first argument is not callable"); \
  135
+      return NULL; \
  136
+  } \
  137
+  if (func_type == INVALID_CTYPE) { \
  138
+      PyErr_SetString(PyExc_TypeError, \
  139
+	   "quad: first argument is a ctypes function pointer with incorrect signature"); \
  140
+      return NULL; \
  141
+  }}
  142
+
  143
+/* Just the first part of the ctypes CFuncPtrObject */
8
Pauli Virtanen Owner
pv added a note April 07, 2012

This hack should be replaced with a call to ctypes.addressof via Python space -- I'd like to avoid mucking around with CPython internals.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

This part is not a hack --- look at get_func_type again. It's just making calls to the ctypes module to identify what fcn is (a callable, a ctype-function with the right signature, or whatnot).

It sounds like you might be talking about the _GET_FUNC macro which is assuming something about the structure of ctypes objects. We could use ctypes.addressof via Python to get the function pointer, but that is a lot of C-code and reference count manipulation when it's a simple pointer de-reference in C-space. It would be very surprising if ctypes ever changed it's ABI so that the pointer immediately following PyObject_HEAD is not the function pointer we are after.

I just don't see it as a major problem in this case.

Pauli Virtanen Owner
pv added a note April 09, 2012

Yep, it's about the _GET_FUNC below. I agree with you that it's unlikely that the ctypes API will change. However, this call is not on a hot path, so I would rather prefer the code use public APIs, if the performance gain from the hack is not significant.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

This is going to come-up again on the NumPy list. I'm not excited about being dependent on a slow Python-only API to get this information out of CTypes whenever we need it --- so we are going to do something like this in NumPy.

One thought is to do this in NumPy and then have SciPy rely on the NumPy API. This seems like a reasonable path forward.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

To that end, I will create a function for getting the raw function pointer using the Python code: ctypes.cast(func, ctypes.c_void_p).value. This is a bit un-pleasant in C and I would rather replace the code with something more direct (like an API that is supported somewhere else --- perhaps in NumPy).

Pauli Virtanen Owner
pv added a note April 09, 2012

Ok, sounds good. Doing Python on C level is indeed messy (a similar situation with no C-level APIs occurs also with I/O on Python 3). It's probably best to move grabbing the pointer to a helper routine in a separate utilities header file. Also using the direct access hack could be considered if it is better encapsulated this way. But anyway, this is a minor issue, so it's best not to waste too much time on it.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

I wrote a C-function to do the casting and my timing test is now 4x slower instead of 2x faster. Of course this is in a loop of 100 calls to the integration routine for a simple function.

This is really unfortunate. Let me see if I can play with it some more.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

I changed the code back to using the direct function-pointer access. It was much, much faster. It makes me also wonder how much it is costing to do the ctypes type checks as well.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/integrate/__quadpack.h
... ...
@@ -67,6 +67,94 @@
67 67
   already_printed_python_error = 0;\
68 68
 }
69 69
 
  70
+#define VALID_CTYPE  1
2
Pauli Virtanen Owner
pv added a note April 07, 2012

Maybe enum?

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

Sure.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/integrate/__quadpack.h
((28 lines not shown))
162 258
 
163  
-  if (PyErr_Occurred()) {
164  
-    ier = 80;             /* Python error */
165  
-    PyErr_Clear();
  259
+    if (PyErr_Occurred()) {
  260
+      ier = 80;             /* Python error */
  261
+      PyErr_Clear();
  262
+    }
  263
+  }
  264
+  else { /* func_type == VALID_CTYPE */
  265
+      /* Can't allow another thread to run because of the global variables
  266
+	 quadpack_raw_function and quad_function2 being used */
2
Pauli Virtanen Owner
pv added a note April 07, 2012

Could also this be made re-entrant by storing the quadpack_raw_function in store_quadpack_globals, like is done for the Python function? Granted, it's less likely that someone calls quad again from C code, but it could happen.

A thunk library would make life easier...

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

Yes, I thought about that. It would be easy enough to make it re-entrant and probably worth it, even with the low-probability event.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Travis E. Oliphant
Owner

Ralf and Pauli thanks for the feedback. I've commented on Pauli's notes in-line.

Ralf, thanks for the skipif suggestion. The timing test is not really a bench-mark. It's more of a "is this code calling the C-function pointer directly, or is it interpreting the Ctypes object as a Python-callback". So, I'm thinking of it as a test to see whether or not the C-function pointer is being called directly by QUADPACK. But, other ways to test this would be helpful. The speed test is a bit brittle.

scipy/integrate/__quadpack.h
... ...
@@ -20,35 +20,35 @@
20 20
  */
21 21
 
22 22
 #if defined(NO_APPEND_FORTRAN)
23  
-	#if defined(UPPERCASE_FORTRAN)
24  
-	/* nothing to do here */
25  
-	#else
26  
-		#define DQAGSE dqagse
27  
-		#define DQAGIE dqagie
28  
-		#define DQAGPE dqagpe
29  
-		#define DQAWOE dqawoe
30  
-		#define DQAWFE dqawfe
31  
-		#define DQAWSE dqawse
32  
-		#define DQAWCE dqawce
33  
-	#endif
  23
+        #if defined(UPPERCASE_FORTRAN)
  24
+        /* nothing to do here */
  25
+        #else
  26
+                #define DQAGSE dqagse
2
Charles Harris Collaborator
charris added a note April 09, 2012

Tabs?

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

I'm not really sure what is going on here. Perhaps there are tabs in the original file. These have been replaced with spaces now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/integrate/__quadpack.h
((13 lines not shown))
  79
+     Returns INVALID_CTYPE if it is a CType Function but of the incorrect type
  80
+     Returns NOTCALLABLE if it is not a Python callable
  81
+     Returns ERROR if other error occurs.
  82
+*/
  83
+
  84
+static int
  85
+get_func_type(PyObject *func) {
  86
+    PyObject *ctypes_module, *CFuncPtr, *check, *c_double;
  87
+    int result;
  88
+
  89
+    if (!PyCallable_Check(func)) return NOTCALLABLE;
  90
+    ctypes_module = PyImport_ImportModule("ctypes");
  91
+    if (ctypes_module == NULL) return ERROR;
  92
+    CFuncPtr = PyObject_GetAttrString(ctypes_module, "_CFuncPtr");
  93
+    if (CFuncPtr == NULL) {
  94
+	Py_DECREF(ctypes_module);
1
Charles Harris Collaborator
charris added a note April 09, 2012

Indentation

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/integrate/__quadpack.h
((54 lines not shown))
  120
+	goto fail;
  121
+
  122
+    return VALID_CTYPE;
  123
+
  124
+fail:
  125
+    Py_DECREF(check);
  126
+    Py_XDECREF(c_double);
  127
+    return INVALID_CTYPE;
  128
+}
  129
+
  130
+#define CHECK_FUNC(func_type, fcn) { \
  131
+  func_type = get_func_type(fcn); \
  132
+  if (func_type == ERROR) return NULL; \
  133
+  if (func_type == NOTCALLABLE) { \
  134
+      PyErr_SetString(PyExc_TypeError, "quad: first argument is not callable"); \
  135
+      return NULL; \
8
Charles Harris Collaborator
charris added a note April 09, 2012

Putting control transfer in a macro is just sinful. I cry every time I see this in Numpy, and I'm sure that somewhere a puppy dies.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

That is a bit melo-dramatic. However, this is not general "control transfer" it's early exit from the function because of errors. It's semantically contained and easily identified upon reviewing the code. The point of any software engineering "ethics" should be to make code easier to write, read, and maintain. You always have to balance many things in that process. I find that zealous over-adherence to any one principle typically causes unintended effects that make things more difficult.

In this particular case, do you have an alternative suggestion that does not require removing the benefit of the macro in encapsulating code that is identical?

Charles Harris Collaborator
charris added a note April 09, 2012

What benefits? Modern compilers will inline small functions without prompting, plus you get type checking on the arguments and your eye can follow the flow of control without reference to some macro off in the woods. Really, I'm trying to break you of your macro habit, they simply aren't necessary for efficiency and when used like this obfuscate the code.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

I told you the benefit. It prevents writing a lot of similar code over and over again. I'm not convinced that being able to "follow the flow of control" actually matters when the flow on error is out-of-the-function.

Any desire to "break" a habit requires first seeing the need to. You do a lot of complaining about style, code re-writing, but not a lot of demonstrating a better approach. This rarely works and generally alienates the people you are purportedly trying to "help". It also takes a lot of time and resources -- resources that are already in short supply. You have to also look at how often someone will care about a particular piece of code. The fact that this bit of code (its supposedly ugly macros and all) has not been really touched in 10 years says something about the need to spend time agonizing over how it gets written.

I'm left to do the work of inferring what you would rather see. Are you preferring to see something like:

if (!check_func(fcn)) return NULL;

where check_func is a regular function? How else would you do things?

-Travis

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

I've got code working that does the check_func and replaces that macro. I do not know how to replace the other macros and keep the benefit of storage variables and error-returning from within the other macros without major code re-design. Eventually, this C-code might be replaced with Cython anyway, and I'd rather not spend more time on it.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

The presumed check_func function and the get_func_type function have been merged into a single function (no macros). There was the need to create a macro to do the temporary function storage trick which is similar to the python-function code-path. If a better approach can be written, I would be happy to see it.

Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

Chuck: Sorry to be a bit prickly before on your observations about the use of Macros. My tendency is towards getting as quickly as possible to working code. Coupling that with most of the code being written when I knew less about the real reasons for object-oriented style and when it matters and when it doesn't often leaves me copying style instead of trying to improve it. I'm sure I'll continue to be more focused on getting something working rather than style, but I don't want to discourage you (or really anyone else as I know you are unflappable) from voicing your stylistic opinions.

Charles Harris Collaborator
charris added a note April 09, 2012

Not to worry. I figure it's my job to play the critic in these situations as your standing and reputation could easily intimidate some folks. Now I'm waiting for Mark to instruct me in the finer points...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/integrate/quadpack.h
... ...
@@ -60,6 +59,18 @@
60 59
   quadpack_extra_arguments = store_quadpack_globals[1]; \
61 60
   memcpy(&quadpack_jmpbuf, &store_jmp, sizeof(jmp_buf));
62 61
 
  62
+#define INIT_CTYPES_FUNC(fcn) \
2
Charles Harris Collaborator
charris added a note April 09, 2012

In C++ this would be a class with a constructor and a restore method. I'd suggest a C implementation along the lines:

struct cfunc_store {
    void *qfunc;
    void *cfunc;
}

NPY_INLINE int
cfunc_store_init(struct cfunc_store *store, void *cfunc, void *qfunc)
{
    store->qfunc = qfunc;
    if ((store->cfunc = get_ctypes_function_pointer(cfunc)) == NULL) {
        return -1;
    }
    return 0;
}

NPY_INLINE void
cfunc_store_restore(struct cfunc_store *store)
{
    store->cfunc = store->qfunc;
}
Travis E. Oliphant Owner
teoliphant added a note April 09, 2012

Thanks. This is helpful. In just trying to get something to work, I was just mirroring the code that was there (from 1999). I'm not sure if it's the best use of time or not, but I agree it's better to define a structure like this. I'm changing all the code to reflect this style.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Travis E. Oliphant ENH: Add ability to use ctypes function pointers in quadpack and clea…
…n-up quadpack to use Object-style C instead of Macros.
9ba5fcc
Travis E. Oliphant
Owner

O.K. I've redone the PR to address the major concerns. The ctypes implementation dependent stuff is reduced to a single function which could be updated to use NumPy's function call when it is created.

Pauli Virtanen
Owner
pv commented May 19, 2012

Is there something additional that needs to be done to this? As far as I see, we don't need to wait here for the Cython/Numba/Numpy fast callback discussion to settle down. The quadpack wrapper code would call for some more reorganization, but that can maybe wait.

Travis E. Oliphant
Owner

I think this is ready to merge. I agree, we don't need to wait for the fast callback discussion.

Ralf Gommers rgommers merged commit b2afa94 into from June 03, 2012
Ralf Gommers rgommers closed this June 03, 2012
Ralf Gommers
Owner

OK, merged to get it in for 0.11.

Ralf Gommers
Owner

Ctypes is used in quadpack, this doesn't work. See http://projects.scipy.org/scipy/ticket/1670.

Travis E. Oliphant
Owner

I will fix this up. It's an import in the C-code. I thought the only reference to ctypes itself was commented out. But, there is one more place. I will submit a new PR.

Travis E. Oliphant teoliphant referenced this pull request June 07, 2012
Closed

Fix Issue #1670 #246

jnothman jnothman referenced this pull request from a commit May 08, 2013
Commit has since been removed from the repository and is no longer available.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Showing 1 unique commit by 1 author.

Apr 09, 2012
Travis E. Oliphant ENH: Add ability to use ctypes function pointers in quadpack and clea…
…n-up quadpack to use Object-style C instead of Macros.
9ba5fcc
This page is out of date. Refresh to see the latest.
450  scipy/integrate/__quadpack.h
@@ -20,35 +20,35 @@
20 20
  */
21 21
 
22 22
 #if defined(NO_APPEND_FORTRAN)
23  
-	#if defined(UPPERCASE_FORTRAN)
24  
-	/* nothing to do here */
25  
-	#else
26  
-		#define DQAGSE dqagse
27  
-		#define DQAGIE dqagie
28  
-		#define DQAGPE dqagpe
29  
-		#define DQAWOE dqawoe
30  
-		#define DQAWFE dqawfe
31  
-		#define DQAWSE dqawse
32  
-		#define DQAWCE dqawce
33  
-	#endif
  23
+  #if defined(UPPERCASE_FORTRAN)
  24
+  /* nothing to do here */
  25
+  #else
  26
+    #define DQAGSE dqagse
  27
+    #define DQAGIE dqagie
  28
+    #define DQAGPE dqagpe
  29
+    #define DQAWOE dqawoe
  30
+    #define DQAWFE dqawfe
  31
+    #define DQAWSE dqawse
  32
+    #define DQAWCE dqawce
  33
+  #endif
34 34
 #else
35  
-	#if defined(UPPERCASE_FORTRAN)
36  
-		#define DQAGSE DQAGSE_
37  
-		#define DQAGIE DQAGIE_
38  
-		#define DQAGPE DQAGPE_
39  
-		#define DQAWOE DQAWOE_
40  
-		#define DQAWFE DQAWFE_
41  
-		#define DQAWSE DQAWSE_
42  
-		#define DQAWCE DQAWCE_
43  
-	#else
44  
-		#define DQAGSE dqagse_
45  
-		#define DQAGIE dqagie_
46  
-		#define DQAGPE dqagpe_
47  
-		#define DQAWOE dqawoe_
48  
-		#define DQAWFE dqawfe_
49  
-		#define DQAWSE dqawse_
50  
-		#define DQAWCE dqawce_
51  
-	#endif
  35
+  #if defined(UPPERCASE_FORTRAN)
  36
+    #define DQAGSE DQAGSE_
  37
+    #define DQAGIE DQAGIE_
  38
+    #define DQAGPE DQAGPE_
  39
+    #define DQAWOE DQAWOE_
  40
+    #define DQAWFE DQAWFE_
  41
+    #define DQAWSE DQAWSE_
  42
+    #define DQAWCE DQAWCE_
  43
+#else
  44
+    #define DQAGSE dqagse_
  45
+    #define DQAGIE dqagie_
  46
+    #define DQAGPE dqagpe_
  47
+    #define DQAWOE dqawoe_
  48
+    #define DQAWFE dqawfe_
  49
+    #define DQAWSE dqawse_
  50
+    #define DQAWCE dqawce_ 
  51
+  #endif
52 52
 #endif
53 53
 
54 54
 void DQAGSE();
@@ -60,11 +60,123 @@ void DQAWSE();
60 60
 void DQAWCE();
61 61
 
62 62
 
63  
-static int already_printed_python_error = 0;
  63
+typedef enum {
  64
+  Error=-3, 
  65
+  Not_Callable=-2,
  66
+  Invalid_Ctype=-1,
  67
+  /* Acceptable returns below */
  68
+  Callable=1,
  69
+  Valid_Ctype=2
  70
+} FuncType;
  71
+
  72
+
  73
+/* Checks a callable object:
  74
+   Returns Valid_Ctype if the Python Object is a CType Function of the type (double) -> (double)
  75
+   Return  Callable if it is a Python callable (including a Ctypes function with no stored types)
  76
+   Returns Invalid_Ctype if it is a CType Function but of the incorrect type
  77
+   Returns Not_Callable if it is not a Python callable
  78
+   Returns Error if other error occurs.
  79
+*/
  80
+
  81
+static FuncType
  82
+get_func_type(PyObject *func) {
  83
+  PyObject *ctypes_module, *CFuncPtr, *check, *c_double;
  84
+  int result;
  85
+  
  86
+  if (!PyCallable_Check(func)) {
  87
+    PyErr_SetString(quadpack_error, "quad: first argument is not callable");
  88
+    return Not_Callable;
  89
+  }
  90
+  ctypes_module = PyImport_ImportModule("ctypes");
  91
+  if (ctypes_module == NULL) return Error;
  92
+  CFuncPtr = PyObject_GetAttrString(ctypes_module, "_CFuncPtr");
  93
+  if (CFuncPtr == NULL) {
  94
+    Py_DECREF(ctypes_module);
  95
+    return Error;
  96
+  }
  97
+  result = PyObject_TypeCheck(func, (PyTypeObject *)CFuncPtr);
  98
+  Py_DECREF(CFuncPtr);
  99
+  if (!result) {
  100
+    Py_DECREF(ctypes_module);
  101
+    return Callable;
  102
+  }
  103
+  
  104
+  /* We are a ctypes function */
  105
+  /* Check for restype and argtypes */
  106
+  if (!PyObject_HasAttrString(func, "restype") ||
  107
+      !PyObject_HasAttrString(func, "argtypes")) {
  108
+    Py_DECREF(ctypes_module);
  109
+    return Callable;
  110
+  }
  111
+  
  112
+  c_double = PyObject_GetAttrString(ctypes_module, "c_double");
  113
+  Py_DECREF(ctypes_module);
  114
+  
  115
+  check = PyObject_GetAttrString(func, "restype");
  116
+  if (check != c_double) {  /* Checking for pointer identity */
  117
+    goto fail;
  118
+  }
  119
+  Py_DECREF(check);
  120
+  check = PyObject_GetAttrString(func, "argtypes");
  121
+  if (!PyTuple_Check(check) || (PyTuple_GET_SIZE(check) != 1) || 
  122
+      (PyTuple_GET_ITEM(check, 0) != c_double)) 
  123
+    goto fail;
  124
+  
  125
+  Py_DECREF(check);
  126
+  Py_DECREF(c_double);
  127
+  return Valid_Ctype;
  128
+  
  129
+fail:
  130
+  Py_DECREF(check);
  131
+  Py_XDECREF(c_double);
  132
+  PyErr_SetString(quadpack_error,
  133
+                  "quad: first argument is a ctypes function pointer with incorrect signature");
  134
+  return Invalid_Ctype;
  135
+}
64 136
 
65  
-#define QUAD_INIT_FUNC(fun,arg) {\
66  
-  INIT_FUNC(fun,arg,quadpack_error); \
67  
-  already_printed_python_error = 0;\
  137
+/*
  138
+  This is the equivalent to the above which only uses ctypes calls... However it is about 8x slower in a tight loop
  139
+  than the above code: 
  140
+ 
  141
+static _sp_double_func
  142
+get_ctypes_function_pointer(PyObject *obj) {
  143
+  PyObject *ctypes_module=NULL, *c_void_p=NULL, *castfunc=NULL, *value=NULL, *result=NULL;
  144
+  void *final;
  145
+
  146
+  ctypes_module = PyImport_ImportModule("ctypes");
  147
+  if (ctypes_module == NULL) return NULL;
  148
+  
  149
+  castfunc = PyObject_GetAttrString(ctypes_module, "cast");
  150
+  if (castfunc == NULL) goto fail;
  151
+  c_void_p = PyObject_GetAttrString(ctypes_module, "c_void_p");
  152
+  if (c_void_p == NULL) goto fail;
  153
+  result = PyObject_CallFunctionObjArgs(castfunc, obj, c_void_p);
  154
+  if (result == NULL) goto fail;
  155
+  value = PyObject_GetAttrString(result, "value");
  156
+  if (value == NULL) goto fail;
  157
+  final = PyLong_AsVoidPtr(value);
  158
+  if (PyErr_Occurred()) goto fail;
  159
+  
  160
+  Py_DECREF(ctypes_module);
  161
+  Py_DECREF(castfunc);
  162
+  Py_DECREF(c_void_p);
  163
+  Py_DECREF(result);
  164
+  Py_DECREF(value);
  165
+  
  166
+  return (_sp_double_func) final;
  167
+  
  168
+fail:
  169
+  Py_XDECREF(ctypes_module);
  170
+  Py_XDECREF(castfunc);
  171
+  Py_XDECREF(c_void_p);
  172
+  Py_XDECREF(result);
  173
+  Py_XDECREF(value);
  174
+  return NULL;    
  175
+}
  176
+*/
  177
+
  178
+double quad_function2(double *x) {
  179
+  return quadpack_ctypes_function(*x);
68 180
 }
69 181
 
70 182
 double quad_function(double *x) {
@@ -87,7 +199,6 @@ double quad_function(double *x) {
87 199
 
88 200
   /* Have to do own error checking because PyFloat_AsDouble returns -1 on
89 201
      error -- making that return value from the function unusable.
90  
-
91 202
      No; Solution is to test for Python Error Occurrence if -1 is return of PyFloat_AsDouble.
92 203
   */
93 204
 
@@ -126,8 +237,8 @@ static PyObject *quadpack_qagse(PyObject *dummy, PyObject *args) {
126 237
   int      neval=0, ier=6, last=0, *iord;
127 238
   double   result=0.0, abserr=0.0;
128 239
   double   *alist, *blist, *rlist, *elist;
129  
-
130  
-  STORE_VARS();
  240
+  FuncType func_type;
  241
+  QStorage storevar;
131 242
 
132 243
   if (!PyArg_ParseTuple(args, "Odd|Oiddi", &fcn, &a, &b, &extra_args, &full_output, &epsabs, &epsrel, &limit)) return NULL;
133 244
   limit_shape[0] = limit;
@@ -136,7 +247,8 @@ static PyObject *quadpack_qagse(PyObject *dummy, PyObject *args) {
136 247
   if (limit < 1)
137 248
     return Py_BuildValue("ddi",result,abserr,ier);
138 249
 
139  
-  QUAD_INIT_FUNC(fcn,extra_args)
  250
+  if ((func_type = get_func_type(fcn)) < Callable) 
  251
+      return NULL;
140 252
 
141 253
   /* Setup iwork and work arrays */
142 254
   ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_INT);
@@ -151,20 +263,30 @@ static PyObject *quadpack_qagse(PyObject *dummy, PyObject *args) {
151 263
   rlist = (double *)ap_rlist->data;
152 264
   elist = (double *)ap_elist->data;
153 265
 
154  
-  if (setjmp(quadpack_jmpbuf)) {
155  
-    goto fail;
156  
-  }
157  
-  else {
158  
-    DQAGSE(quad_function, &a, &b, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
159  
-  }
  266
+  if (func_type == Callable) {
  267
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  268
+      goto fail;
160 269
 
161  
-  RESTORE_FUNC();
  270
+    if (setjmp(quadpack_jmpbuf)) {
  271
+      quad_restore_func(&storevar, NULL);
  272
+      goto fail;
  273
+    }
  274
+    else {
  275
+      DQAGSE(quad_function, &a, &b, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, 
  276
+    blist, rlist, elist, iord, &last);
  277
+    }
162 278
 
163  
-  if (PyErr_Occurred()) {
164  
-    ier = 80;             /* Python error */
165  
-    PyErr_Clear();
  279
+    quad_restore_func(&storevar, &ier);
  280
+  }
  281
+  else { /* func_type == VALID_CTYPE */
  282
+      /* Can't allow another thread to run because of the global variables
  283
+         quadpack_raw_function and quad_function2 being used */
  284
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL)
  285
+      goto fail;
  286
+    DQAGSE(quad_function2, &a, &b, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist,
  287
+           blist, rlist, elist, iord, &last);
  288
+    restore_ctypes_func(&storevar);
166 289
   }
167  
-  Py_DECREF(extra_args);
168 290
 
169 291
   if (full_output) {
170 292
     return Py_BuildValue("dd{s:i,s:i,s:N,s:N,s:N,s:N,s:N}i", result, abserr, "neval", neval, "last", last, "iord", PyArray_Return(ap_iord), "alist", PyArray_Return(ap_alist), "blist", PyArray_Return(ap_blist), "rlist", PyArray_Return(ap_rlist), "elist", PyArray_Return(ap_elist),ier);
@@ -179,8 +301,6 @@ static PyObject *quadpack_qagse(PyObject *dummy, PyObject *args) {
179 301
   }
180 302
 
181 303
  fail:
182  
-  RESTORE_FUNC();
183  
-  Py_XDECREF(extra_args);
184 304
   Py_XDECREF(ap_alist);
185 305
   Py_XDECREF(ap_blist);
186 306
   Py_XDECREF(ap_rlist);
@@ -207,17 +327,20 @@ static PyObject *quadpack_qagie(PyObject *dummy, PyObject *args) {
207 327
   int      inf, neval=0, ier=6, last=0, *iord;
208 328
   double   result=0.0, abserr=0.0;
209 329
   double   *alist, *blist, *rlist, *elist;
  330
+  FuncType func_type;
  331
+  QStorage storevar; 
210 332
 
211  
-  STORE_VARS();
212  
-
213  
-  if (!PyArg_ParseTuple(args, "Odi|Oiddi", &fcn, &bound, &inf, &extra_args, &full_output, &epsabs, &epsrel, &limit)) return NULL;
  333
+  if (!PyArg_ParseTuple(args, "Odi|Oiddi", &fcn, &bound, &inf, &extra_args, 
  334
+                        &full_output, &epsabs, &epsrel, &limit)) 
  335
+    return NULL;
214 336
   limit_shape[0] = limit;
215 337
 
216 338
   /* Need to check that limit is bigger than 1 */
217 339
   if (limit < 1)
218 340
     return Py_BuildValue("ddi",result,abserr,ier);
219 341
 
220  
-  QUAD_INIT_FUNC(fcn,extra_args);
  342
+  if ((func_type = get_func_type(fcn)) < Callable) 
  343
+      return NULL;
221 344
 
222 345
   /* Setup iwork and work arrays */
223 346
   ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_INT);
@@ -225,28 +348,34 @@ static PyObject *quadpack_qagie(PyObject *dummy, PyObject *args) {
225 348
   ap_blist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE);
226 349
   ap_rlist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE);
227 350
   ap_elist = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_DOUBLE);
228  
-  if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL || ap_elist == NULL) goto fail;
  351
+  if (ap_iord == NULL || ap_alist == NULL || ap_blist == NULL || ap_rlist == NULL
  352
+      || ap_elist == NULL) goto fail;
229 353
   iord = (int *)ap_iord->data;
230 354
   alist = (double *)ap_alist->data;
231 355
   blist = (double *)ap_blist->data;
232 356
   rlist = (double *)ap_rlist->data;
233 357
   elist = (double *)ap_elist->data;
234 358
 
235  
-  if (setjmp(quadpack_jmpbuf)) {
236  
-    goto fail;
237  
-  }
238  
-  else {
239  
-    DQAGIE(quad_function, &bound, &inf, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
240  
-  }
241  
-
242  
-  RESTORE_FUNC();
  359
+  if (func_type == Callable) {
  360
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  361
+      goto fail;
243 362
 
244  
-  if (PyErr_Occurred()) {
245  
-    ier = 80;             /* Python error */
246  
-    PyErr_Clear();
  363
+    if (setjmp(quadpack_jmpbuf)) {
  364
+      quad_restore_func(&storevar, NULL);
  365
+      goto fail;
  366
+    }
  367
+    else {
  368
+      DQAGIE(quad_function, &bound, &inf, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  369
+    }
  370
+    quad_restore_func(&storevar, &ier);
  371
+  }
  372
+  else { /* func_type == VALID_CTYPE */
  373
+      /* Can't allow another thread to run because of the global variables
  374
+         quadpack_raw_function and quad_function2 being used */
  375
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL) goto fail;
  376
+    DQAGIE(quad_function2, &bound, &inf, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  377
+    restore_ctypes_func(&storevar);
247 378
   }
248  
-
249  
-  Py_DECREF(extra_args);
250 379
 
251 380
   if (full_output) {
252 381
     return Py_BuildValue("dd{s:i,s:i,s:N,s:N,s:N,s:N,s:N}i", result, abserr, "neval", neval, "last", last, "iord", PyArray_Return(ap_iord), "alist", PyArray_Return(ap_alist), "blist", PyArray_Return(ap_blist), "rlist", PyArray_Return(ap_rlist), "elist", PyArray_Return(ap_elist),ier);
@@ -261,8 +390,6 @@ static PyObject *quadpack_qagie(PyObject *dummy, PyObject *args) {
261 390
   }
262 391
 
263 392
  fail:
264  
-  RESTORE_FUNC();
265  
-  Py_XDECREF(extra_args);
266 393
   Py_XDECREF(ap_alist);
267 394
   Py_XDECREF(ap_blist);
268 395
   Py_XDECREF(ap_rlist);
@@ -294,8 +421,8 @@ static PyObject *quadpack_qagpe(PyObject *dummy, PyObject *args) {
294 421
   double   result=0.0, abserr=0.0;
295 422
   double   *alist, *blist, *rlist, *elist;
296 423
   double   *pts, *points;
297  
-
298  
-  STORE_VARS();
  424
+  FuncType func_type;
  425
+  QStorage storevar;    
299 426
 
300 427
   if (!PyArg_ParseTuple(args, "OddO|Oiddi", &fcn, &a, &b, &o_points, &extra_args, &full_output, &epsabs, &epsrel, &limit)) return NULL;
301 428
   limit_shape[0] = limit;
@@ -304,7 +431,8 @@ static PyObject *quadpack_qagpe(PyObject *dummy, PyObject *args) {
304 431
   if (limit < 1)
305 432
     return Py_BuildValue("ddi",result,abserr,ier);
306 433
 
307  
-  QUAD_INIT_FUNC(fcn,extra_args)
  434
+  if ((func_type = get_func_type(fcn)) < Callable)
  435
+      return NULL;
308 436
 
309 437
   ap_points = (PyArrayObject *)PyArray_ContiguousFromObject(o_points, NPY_DOUBLE, 1, 1);
310 438
   if (ap_points == NULL) goto fail;
@@ -331,20 +459,26 @@ static PyObject *quadpack_qagpe(PyObject *dummy, PyObject *args) {
331 459
   level = (int *)ap_level->data;
332 460
   ndin = (int *)ap_level->data;
333 461
 
334  
-  if (setjmp(quadpack_jmpbuf)) {
335  
-    goto fail;
  462
+  if (func_type == Callable) {
  463
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  464
+      goto fail;
  465
+
  466
+    if (setjmp(quadpack_jmpbuf)) {
  467
+      quad_restore_func(&storevar, NULL);
  468
+      goto fail;
  469
+    }
  470
+    else {
  471
+      DQAGPE(quad_function, &a, &b, &npts2, points, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, pts, iord, level, ndin, &last);
  472
+    }
  473
+    
  474
+    quad_restore_func(&storevar, &ier);
336 475
   }
337 476
   else {
338  
-    DQAGPE(quad_function, &a, &b, &npts2, points, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, pts, iord, level, ndin, &last);
  477
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL) goto fail;
  478
+    DQAGPE(quad_function2, &a, &b, &npts2, points, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, pts, iord, level, ndin, &last);
  479
+    restore_ctypes_func(&storevar);
339 480
   }
340 481
 
341  
-  RESTORE_FUNC()
342  
-
343  
-  if (PyErr_Occurred()) {
344  
-    ier = 80;             /* Python error */
345  
-    PyErr_Clear();
346  
-  }
347  
-  Py_DECREF(extra_args);
348 482
   Py_DECREF(ap_points);
349 483
 
350 484
   if (full_output) {
@@ -363,8 +497,6 @@ static PyObject *quadpack_qagpe(PyObject *dummy, PyObject *args) {
363 497
   }
364 498
 
365 499
  fail:
366  
-  RESTORE_FUNC();
367  
-  Py_XDECREF(extra_args);
368 500
   Py_XDECREF(ap_alist);
369 501
   Py_XDECREF(ap_blist);
370 502
   Py_XDECREF(ap_rlist);
@@ -399,8 +531,8 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) {
399 531
   double   result=0.0, abserr=0.0, omega=0.0;
400 532
   double   *chebmo;
401 533
   double   *alist, *blist, *rlist, *elist;
402  
-
403  
-  STORE_VARS();
  534
+  FuncType func_type;  
  535
+  QStorage storevar;  
404 536
 
405 537
   if (!PyArg_ParseTuple(args, "Odddi|OiddiiiiO", &fcn, &a, &b, &omega, &integr, &extra_args, &full_output, &epsabs, &epsrel, &limit, &maxp1, &icall, &momcom, &o_chebmo)) return NULL;
406 538
   limit_shape[0] = limit;
@@ -409,7 +541,8 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) {
409 541
   if (limit < 1)
410 542
     return Py_BuildValue("ddi",result,abserr,ier);
411 543
 
412  
-  QUAD_INIT_FUNC(fcn,extra_args)
  544
+  if ((func_type = get_func_type(fcn)) < Callable)
  545
+      return NULL;
413 546
 
414 547
   if (o_chebmo != NULL) {
415 548
     ap_chebmo = (PyArrayObject *)PyArray_ContiguousFromObject(o_chebmo, NPY_DOUBLE, 2, 2);
@@ -440,20 +573,26 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) {
440 573
   rlist = (double *)ap_rlist->data;
441 574
   elist = (double *)ap_elist->data;
442 575
 
443  
-  if (setjmp(quadpack_jmpbuf)) {
444  
-    goto fail;
445  
-  }
446  
-  else {
447  
-    DQAWOE(quad_function, &a, &b, &omega, &integr, &epsabs, &epsrel, &limit, &icall, &maxp1, &result, &abserr, &neval, &ier, &last, alist, blist, rlist, elist, iord, nnlog, &momcom, chebmo);
448  
-  }
  576
+  if (func_type == Callable) {
  577
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  578
+      goto fail;
  579
+    
  580
+    if (setjmp(quadpack_jmpbuf)) {
  581
+      quad_restore_func(&storevar, NULL);
  582
+      goto fail;
  583
+    }
  584
+    else {
  585
+      DQAWOE(quad_function, &a, &b, &omega, &integr, &epsabs, &epsrel, &limit, &icall, &maxp1, &result, &abserr, &neval, &ier, &last, alist, blist, rlist, elist, iord, nnlog, &momcom, chebmo);
  586
+    }
449 587
 
450  
-  RESTORE_FUNC();
  588
+    quad_restore_func(&storevar, &ier);
451 589
 
452  
-  if (PyErr_Occurred()) {
453  
-    ier = 80;             /* Python error */
454  
-    PyErr_Clear();
455 590
   }
456  
-  Py_DECREF(extra_args);
  591
+  else {
  592
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL) goto fail;
  593
+    DQAWOE(quad_function2, &a, &b, &omega, &integr, &epsabs, &epsrel, &limit, &icall, &maxp1, &result, &abserr, &neval, &ier, &last, alist, blist, rlist, elist, iord, nnlog, &momcom, chebmo);
  594
+    restore_ctypes_func(&storevar);
  595
+  }
457 596
 
458 597
   if (full_output) {
459 598
     return Py_BuildValue("dd{s:i,s:i,s:N,s:N,s:N,s:N,s:N,s:N,s:i,s:N}i", result, abserr, "neval", neval, "last", last, "iord", PyArray_Return(ap_iord), "alist", PyArray_Return(ap_alist), "blist", PyArray_Return(ap_blist), "rlist", PyArray_Return(ap_rlist), "elist", PyArray_Return(ap_elist), "nnlog", PyArray_Return(ap_nnlog), "momcom", momcom, "chebmo", PyArray_Return(ap_chebmo),ier);
@@ -470,8 +609,6 @@ static PyObject *quadpack_qawoe(PyObject *dummy, PyObject *args) {
470 609
   }
471 610
 
472 611
  fail:
473  
-  RESTORE_FUNC();
474  
-  Py_XDECREF(extra_args);
475 612
   Py_XDECREF(ap_alist);
476 613
   Py_XDECREF(ap_blist);
477 614
   Py_XDECREF(ap_rlist);
@@ -505,8 +642,8 @@ static PyObject *quadpack_qawfe(PyObject *dummy, PyObject *args) {
505 642
   double   *chebmo, *rslst, *erlst;
506 643
   double   result=0.0, abserr=0.0, omega=0.0;
507 644
   double   *alist, *blist, *rlist, *elist;
508  
-
509  
-  STORE_VARS();
  645
+  FuncType func_type;  
  646
+  QStorage storevar;  
510 647
 
511 648
   if (!PyArg_ParseTuple(args, "Oddi|Oidiii", &fcn, &a, &omega, &integr, &extra_args, &full_output, &epsabs, &limlst, &limit, &maxp1)) return NULL;
512 649
   limit_shape[0] = limit;
@@ -516,7 +653,8 @@ static PyObject *quadpack_qawfe(PyObject *dummy, PyObject *args) {
516 653
   if (limit < 1)
517 654
     return Py_BuildValue("ddi",result,abserr,ier);
518 655
 
519  
-  QUAD_INIT_FUNC(fcn,extra_args)
  656
+  if ((func_type = get_func_type(fcn)) < Callable)
  657
+      return NULL;
520 658
 
521 659
   sz[0] = 25;
522 660
   sz[1] = maxp1;
@@ -545,20 +683,26 @@ static PyObject *quadpack_qawfe(PyObject *dummy, PyObject *args) {
545 683
   erlst = (double *)ap_erlst->data;
546 684
   ierlst = (int *)ap_ierlst->data;
547 685
 
548  
-  if (setjmp(quadpack_jmpbuf)) {
549  
-    goto fail;
  686
+  if (func_type == Callable) {
  687
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  688
+      goto fail;
  689
+
  690
+    if (setjmp(quadpack_jmpbuf)) {
  691
+      quad_restore_func(&storevar, NULL);
  692
+      goto fail;
  693
+    }
  694
+    else {
  695
+      DQAWFE(quad_function, &a, &omega, &integr, &epsabs, &limlst, &limit, &maxp1, &result, &abserr, &neval, &ier, rslst, erlst, ierlst, &lst, alist, blist, rlist, elist, iord, nnlog, chebmo);
  696
+    }
  697
+
  698
+    quad_restore_func(&storevar, &ier);
550 699
   }
551 700
   else {
552  
-    DQAWFE(quad_function, &a, &omega, &integr, &epsabs, &limlst, &limit, &maxp1, &result, &abserr, &neval, &ier, rslst, erlst, ierlst, &lst, alist, blist, rlist, elist, iord, nnlog, chebmo);
  701
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL) goto fail;
  702
+    DQAWFE(quad_function2, &a, &omega, &integr, &epsabs, &limlst, &limit, &maxp1, &result, &abserr, &neval, &ier, rslst, erlst, ierlst, &lst, alist, blist, rlist, elist, iord, nnlog, chebmo);
  703
+    restore_ctypes_func(&storevar);
553 704
   }
554 705
 
555  
-  RESTORE_FUNC();
556  
-
557  
-  if (PyErr_Occurred()) {
558  
-    ier = 80;             /* Python error */
559  
-    PyErr_Clear();
560  
-  }
561  
-  Py_DECREF(extra_args);
562 706
   Py_DECREF(ap_nnlog);
563 707
   Py_DECREF(ap_alist);
564 708
   Py_DECREF(ap_blist);
@@ -578,8 +722,6 @@ static PyObject *quadpack_qawfe(PyObject *dummy, PyObject *args) {
578 722
   }
579 723
 
580 724
  fail:
581  
-  RESTORE_FUNC();
582  
-  Py_XDECREF(extra_args);
583 725
   Py_XDECREF(ap_alist);
584 726
   Py_XDECREF(ap_blist);
585 727
   Py_XDECREF(ap_rlist);
@@ -612,8 +754,8 @@ static PyObject *quadpack_qawce(PyObject *dummy, PyObject *args) {
612 754
   int      neval=0, ier=6, last=0, *iord;
613 755
   double   result=0.0, abserr=0.0;
614 756
   double   *alist, *blist, *rlist, *elist;
615  
-
616  
-  STORE_VARS();
  757
+  FuncType func_type;
  758
+  QStorage storevar;  
617 759
 
618 760
   if (!PyArg_ParseTuple(args, "Oddd|Oiddi", &fcn, &a, &b, &c, &extra_args, &full_output, &epsabs, &epsrel, &limit)) return NULL;
619 761
   limit_shape[0] = limit;
@@ -622,7 +764,8 @@ static PyObject *quadpack_qawce(PyObject *dummy, PyObject *args) {
622 764
   if (limit < 1)
623 765
     return Py_BuildValue("ddi",result,abserr,ier);
624 766
 
625  
-  QUAD_INIT_FUNC(fcn,extra_args)
  767
+  if ((func_type = get_func_type(fcn)) < Callable)
  768
+      return NULL;
626 769
 
627 770
   /* Setup iwork and work arrays */
628 771
   ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_INT);
@@ -637,20 +780,26 @@ static PyObject *quadpack_qawce(PyObject *dummy, PyObject *args) {
637 780
   rlist = (double *)ap_rlist->data;
638 781
   elist = (double *)ap_elist->data;
639 782
 
640  
-  if (setjmp(quadpack_jmpbuf)) {
641  
-    goto fail;
642  
-  }
643  
-  else {
644  
-    DQAWCE(quad_function, &a, &b, &c, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
645  
-  }
  783
+  if (func_type == Callable) {
  784
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  785
+      goto fail;
  786
+
  787
+    if (setjmp(quadpack_jmpbuf)) {
  788
+      quad_restore_func(&storevar, NULL);
  789
+      goto fail;
  790
+    }
  791
+    else {
  792
+      DQAWCE(quad_function, &a, &b, &c, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  793
+    }
646 794
 
647  
-  RESTORE_FUNC();
  795
+    quad_restore_func(&storevar, &ier);
648 796
 
649  
-  if (PyErr_Occurred()) {
650  
-    ier = 80;             /* Python error */
651  
-    PyErr_Clear();
  797
+  } 
  798
+  else {
  799
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL) goto fail;
  800
+    DQAWCE(quad_function2, &a, &b, &c, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  801
+    restore_ctypes_func(&storevar);
652 802
   }
653  
-  Py_DECREF(extra_args);
654 803
 
655 804
   if (full_output) {
656 805
     return Py_BuildValue("dd{s:i,s:i,s:N,s:N,s:N,s:N,s:N}i", result, abserr, "neval", neval, "last", last, "iord", PyArray_Return(ap_iord), "alist", PyArray_Return(ap_alist), "blist", PyArray_Return(ap_blist), "rlist", PyArray_Return(ap_rlist), "elist", PyArray_Return(ap_elist),ier);
@@ -665,8 +814,6 @@ static PyObject *quadpack_qawce(PyObject *dummy, PyObject *args) {
665 814
   }
666 815
 
667 816
  fail:
668  
-  RESTORE_FUNC();
669  
-  Py_XDECREF(extra_args);
670 817
   Py_XDECREF(ap_alist);
671 818
   Py_XDECREF(ap_blist);
672 819
   Py_XDECREF(ap_rlist);
@@ -695,8 +842,8 @@ static PyObject *quadpack_qawse(PyObject *dummy, PyObject *args) {
695 842
   int      neval=0, ier=6, last=0, *iord;
696 843
   double   result=0.0, abserr=0.0;
697 844
   double   *alist, *blist, *rlist, *elist;
698  
-
699  
-  STORE_VARS();
  845
+  FuncType func_type;
  846
+  QStorage storevar;  
700 847
 
701 848
   if (!PyArg_ParseTuple(args, "Odd(dd)i|Oiddi", &fcn, &a, &b, &alfa, &beta, &integr, &extra_args, &full_output, &epsabs, &epsrel, &limit)) return NULL;
702 849
   limit_shape[0] = limit;
@@ -705,7 +852,8 @@ static PyObject *quadpack_qawse(PyObject *dummy, PyObject *args) {
705 852
   if (limit < 1)
706 853
     return Py_BuildValue("ddi",result,abserr,ier);
707 854
 
708  
-  QUAD_INIT_FUNC(fcn,extra_args)
  855
+  if ((func_type = get_func_type(fcn)) < Callable)
  856
+      return NULL;
709 857
 
710 858
   /* Setup iwork and work arrays */
711 859
   ap_iord = (PyArrayObject *)PyArray_SimpleNew(1,limit_shape,NPY_INT);
@@ -720,21 +868,26 @@ static PyObject *quadpack_qawse(PyObject *dummy, PyObject *args) {
720 868
   rlist = (double *)ap_rlist->data;
721 869
   elist = (double *)ap_elist->data;
722 870
 
723  
-  if (setjmp(quadpack_jmpbuf)) {
724  
-    goto fail;
  871
+  if (func_type == Callable) {
  872
+    if (quad_init_func(&storevar, fcn, extra_args) == NPY_FAIL)
  873
+      goto fail;
  874
+   
  875
+    if (setjmp(quadpack_jmpbuf)) {
  876
+      quad_restore_func(&storevar, NULL);
  877
+      goto fail;
  878
+    }
  879
+    else {
  880
+      DQAWSE(quad_function, &a, &b, &alfa, &beta, &integr, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  881
+    }
  882
+
  883
+    quad_restore_func(&storevar, &ier);
725 884
   }
726 885
   else {
727  
-    DQAWSE(quad_function, &a, &b, &alfa, &beta, &integr, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  886
+    if (init_ctypes_func(&storevar, fcn) == NPY_FAIL) goto fail;
  887
+    DQAWSE(quad_function2, &a, &b, &alfa, &beta, &integr, &epsabs, &epsrel, &limit, &result, &abserr, &neval, &ier, alist, blist, rlist, elist, iord, &last);
  888
+    restore_ctypes_func(&storevar);
728 889
   }
729  
-
730  
-  RESTORE_FUNC();
731  
-
732  
-  if (PyErr_Occurred()) {
733  
-    ier = 80;             /* Python error */
734  
-    PyErr_Clear();
735  
-  }
736  
-  Py_DECREF(extra_args);
737  
-
  890
+  
738 891
   if (full_output) {
739 892
     return Py_BuildValue("dd{s:i,s:i,s:N,s:N,s:N,s:N,s:N}i", result, abserr, "neval", neval, "last", last, "iord", PyArray_Return(ap_iord), "alist", PyArray_Return(ap_alist), "blist", PyArray_Return(ap_blist), "rlist", PyArray_Return(ap_rlist), "elist", PyArray_Return(ap_elist),ier);
740 893
   }
@@ -748,8 +901,6 @@ static PyObject *quadpack_qawse(PyObject *dummy, PyObject *args) {
748 901
   }
749 902
 
750 903
  fail:
751  
-  RESTORE_FUNC();
752  
-  Py_XDECREF(extra_args);
753 904
   Py_XDECREF(ap_alist);
754 905
   Py_XDECREF(ap_blist);
755 906
   Py_XDECREF(ap_rlist);
@@ -758,6 +909,3 @@ static PyObject *quadpack_qawse(PyObject *dummy, PyObject *args) {
758 909
   return NULL;
759 910
 }
760 911
 
761  
-
762  
-
763  
-
1  scipy/integrate/_quadpackmodule.c
@@ -2,7 +2,6 @@
2 2
   From Multipack project
3 3
  */
4 4
 #include "quadpack.h"
5  
-static PyObject *quadpack_error;
6 5
 #include "__quadpack.h"
7 6
 static struct PyMethodDef quadpack_module_methods[] = {
8 7
 {"_qagse", quadpack_qagse, METH_VARARGS, doc_qagse},
100  scipy/integrate/quadpack.h
@@ -27,44 +27,96 @@ results are placed in a dictionary and returned as an optional member of
27 27
 the result tuple when the full_output argument is non-zero.
28 28
 */
29 29
 
30  
-#include "Python.h"
31 30
 
  31
+#include "Python.h"
32 32
 #include "numpy/npy_3kcompat.h"
33 33
 #include "numpy/arrayobject.h"
34 34
 #include <setjmp.h>
35 35
 
  36
+
36 37
 #define PYERR(errobj,message) {PyErr_SetString(errobj,message); goto fail;}
37 38
 #define PYERR2(errobj,message) {PyErr_Print(); PyErr_SetString(errobj, message); goto fail;}
38 39
 #define ISCONTIGUOUS(m) ((m)->flags & CONTIGUOUS)
39 40
 
40  
-#define STORE_VARS() PyObject *store_quadpack_globals[2]; jmp_buf store_jmp;
41  
-
42  
-#define INIT_FUNC(fun,arg,errobj) { /* Get extra arguments or set to zero length tuple */ \
43  
-  store_quadpack_globals[0] = quadpack_python_function; \
44  
-  store_quadpack_globals[1] = quadpack_extra_arguments; \
45  
-  memcpy(&store_jmp,&quadpack_jmpbuf,sizeof(jmp_buf)); \
46  
-  if (arg == NULL) { \
47  
-    if ((arg = PyTuple_New(0)) == NULL) goto fail; \
48  
-  } \
49  
-  else \
50  
-    Py_INCREF(arg);   /* We decrement on exit. */ \
51  
-  if (!PyTuple_Check(arg))  \
52  
-    PYERR(errobj,"Extra Arguments must be in a tuple"); \
53  
-  /* Set up callback functions */ \
54  
-  if (!PyCallable_Check(fun)) \
55  
-    PYERR(errobj,"First argument must be a callable function."); \
56  
-  quadpack_python_function = fun; \
57  
-  quadpack_extra_arguments = arg;}
58  
-
59  
-#define RESTORE_FUNC() quadpack_python_function = store_quadpack_globals[0]; \
60  
-  quadpack_extra_arguments = store_quadpack_globals[1]; \
61  
-  memcpy(&quadpack_jmpbuf, &store_jmp, sizeof(jmp_buf));
62 41
 
63 42
 static PyObject *quadpack_python_function=NULL;
64 43
 static PyObject *quadpack_extra_arguments=NULL;    /* a tuple */
65 44
 static jmp_buf quadpack_jmpbuf;
66 45
 
67  
-
  46
+static double (*quadpack_ctypes_function)(double) = NULL;
  47
+
  48
+static PyObject *quadpack_error;
  49
+
  50
+/* Stack Storage for re-entrant capability */
  51
+typedef struct {
  52
+    void *global0;
  53
+    void *global1;
  54
+    jmp_buf jmp;    
  55
+    PyObject *arg;
  56
+} QStorage;
  57
+
  58
+typedef double (*_sp_double_func)(double);
  59
+
  60
+typedef struct {
  61
+    PyObject_HEAD
  62
+    char *b_ptr;
  63
+} _sp_cfuncptr_object;
  64
+
  65
+static _sp_double_func
  66
+get_ctypes_function_pointer(PyObject *obj) {
  67
+    return (*((void **)(((_sp_cfuncptr_object *)(obj))->b_ptr)));
  68
+}
  69
+
  70
+static int 
  71
+quad_init_func(QStorage *store, PyObject *fun, PyObject *arg) {
  72
+    store->global0 = (void *)quadpack_python_function;
  73
+    store->global1 = (void *)quadpack_extra_arguments;
  74
+    memcpy(&(store->jmp), &quadpack_jmpbuf, sizeof(jmp_buf));
  75
+    store->arg = arg;
  76
+    if (store->arg == NULL) {
  77
+        if ((store->arg = PyTuple_New(0)) == NULL) 
  78
+            return NPY_FAIL;
  79
+    }
  80
+    else {
  81
+        Py_INCREF(store->arg);  /* We decrement on restore */
  82
+    }
  83
+    if (!PyTuple_Check(store->arg)) {
  84
+        PyErr_SetString(quadpack_error, "Extra Arguments must be in a tuple");
  85
+        Py_XDECREF(store->arg);
  86
+        return NPY_FAIL;
  87
+    }
  88
+    quadpack_python_function = fun;
  89
+    quadpack_extra_arguments = store->arg;
  90
+    return NPY_SUCCEED;
  91
+}
  92
+
  93
+static void
  94
+quad_restore_func(QStorage *store, int *ierr) {
  95
+    quadpack_python_function = (PyObject *)store->global0;
  96
+    quadpack_extra_arguments = (PyObject *)store->global1;
  97
+    memcpy(&quadpack_jmpbuf, &(store->jmp), sizeof(jmp_buf));
  98
+    Py_XDECREF(store->arg);
  99
+    if (ierr != NULL) {
  100
+        if (PyErr_Occurred()) {
  101
+            *ierr = 80;             /* Python error */
  102
+            PyErr_Clear();
  103
+        }
  104
+    }
  105
+}
  106
+
  107
+static int
  108
+init_ctypes_func(QStorage *store, PyObject *fun) {
  109
+    store->global0 = quadpack_ctypes_function;
  110
+    store->global1 = get_ctypes_function_pointer(fun);
  111
+    if (store->global1 == NULL) return NPY_FAIL;
  112
+    quadpack_ctypes_function = store->global1;
  113
+    return NPY_SUCCEED;
  114
+}
  115
+
  116
+static void
  117
+restore_ctypes_func(QStorage *store) {
  118
+    quadpack_ctypes_function = store->global0;
  119
+}
68 120
 
69 121
 
70 122
 
47  scipy/integrate/tests/test_quadpack.py
... ...
@@ -1,12 +1,57 @@
1 1
 from numpy import sqrt, cos, sin, arctan, exp, log, pi, Inf
2  
-from numpy.testing import assert_, TestCase, run_module_suite
  2
+from numpy.testing import assert_, TestCase, run_module_suite, dec
3 3
 from scipy.integrate import quad, dblquad, tplquad
  4
+import sys
  5
+import math
  6
+
  7
+try:
  8
+    import ctypes
  9
+    _ctypes_missing = False
  10
+except ImportError:
  11
+    _ctypes_missing = True    
4 12
 
5 13
 def assert_quad((value, err), tabledValue, errTol=1.5e-8):
6 14
     assert_(abs(value-tabledValue) < err, (value, tabledValue, err))
7 15
     if errTol is not None:
8 16
         assert_(err < errTol, (err, errTol))
9 17
 
  18
+
  19
+class TestCtypesQuad(TestCase):
  20
+    @dec.skipif(_ctypes_missing, msg="Ctypes library could not be found")
  21
+    def setUp(self):
  22
+        if sys.platform == 'win32':
  23
+            file = 'msvcrt.dll'
  24
+        elif sys.platform == 'darwin':
  25
+            file = 'libm.dylib'
  26
+        else:
  27
+            file = 'libm.so'
  28
+        self.lib = ctypes.CDLL(file)
  29
+        restype = ctypes.c_double
  30
+        argtypes = (ctypes.c_double,)
  31
+        for name in ['sin', 'cos', 'tan']:
  32
+            func = getattr(self.lib, name)
  33
+            func.restype = restype
  34
+            func.argtypes = argtypes
  35
+
  36
+    @dec.skipif(_ctypes_missing, msg="Ctypes library could not be found")
  37
+    def test_typical(self):
  38
+        assert_quad(quad(self.lib.sin,0,5),quad(math.sin,0,5)[0])
  39
+        assert_quad(quad(self.lib.cos,0,5),quad(math.cos,0,5)[0])
  40
+        assert_quad(quad(self.lib.tan,0,1),quad(math.tan,0,1)[0])
  41
+
  42
+    @dec.skipif(_ctypes_missing, msg="Ctypes library could not be found")
  43
+    def test_improvement(self):
  44
+        import time
  45
+        start = time.time()
  46
+        for i in xrange(100):
  47
+            quad(self.lib.sin, 0, 100)
  48
+        fast = time.time() - start    
  49
+        start = time.time()
  50
+        for i in xrange(100):
  51
+            quad(math.sin, 0, 100)
  52
+        slow = time.time() - start
  53
+        assert_(fast < 0.5*slow, (fast, slow))
  54
+
10 55
 class TestQuad(TestCase):
11 56
     def test_typical(self):
12 57
         # 1) Typical function with two extra arguments:
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.