Implement back end of faster multivariate integration #3262

Closed
wants to merge 6 commits into
from

Conversation

Projects
None yet
10 participants
@BrianNewsom
Contributor

BrianNewsom commented Jan 31, 2014

This PR implements my recent work on the backend of scipy.integrate, adding:
A lower level (C) method of handling multivariate functions so that they may be passed to the quadpack routines
C Wrappers for all double fortran quadpack routines supporting this function type
A simple additional wrapper that allows the wrappers to be incorporated into the existing SciPy integration, for testing purposes.

1.) Currently Scipy’s integration is inhibited in that it’s multivariate integration requires nested calls between python and fortran that are slow. The intent of my work (which is not fully completed) is to allow the movement of the innermost-loop integration entirely to compiled languages, which should be much faster.
The first step of my solution is the cwrapper.h file, which converts a multivariate function of the form f(int nargs, double args[nargs]) to a function of a single variable - the form expected by quadpack. This is handled nearly exactly as it is currently done in Python, storing the original function and any additional arguments to global variables, then using them from inside a generic “call” function passed into quadpack.

An example function of this type would be:

    double f(int nargs, double args[nargs]){
        args[0] * args[1] + args[2] * args[2]
     }

Which would be equivalent to f(x) = x_y+z_z but now allows us to pass the number of extra arguments (2) and their values in the args array so that they may be integrated in fortran.

2.) This file then additionally includes C wrappers for each of quadpack’s double routines that handle this evaluation. I added wrappers for every routine, for future use if necessary, though the current schema seems to only use about 6 of them. These can be removed if superfluous. Between these two elements the foundation is laid for handling the n-dimensional integration at this lower level.

3.) Additionally, so that these wrappers may be tested and used, an additional wrapper has been added. This wrapper allows the currently handled python and ctypes functions to be evaluated in the f(nargs,args) form with another simple initialization and call. These functions are funcwrapper_init and funcwrapper, and these have been applied in the __quadpack.h file inside each routine so that the cwrapper is being applied as will be necessary with future additions. These functions should be removed, once SciPy’s integration module is adapted to use the new-style wrappers, as described earlier.

What does this mean for the end user?
As of now, nothing. No functionality has changed, quad works exactly the same way as before, passing all existing unit tests going through my wrappers. It can still handle ctypes or python functions of any number of variables, and this handling is still quite slow. Because of this, no additional documentation has been included.

Why does this PR matter, then?
This code adds a framework for doing exactly what I described earlier on, moving the handling of multiple variables down into a compiled language, which should provide a significant speed advantage to the power user, once a framework is written that allows ctypes or cython (or both) functions of this new format to be passed directly into the wrappers, without changing anything for the traditional user.

I am a first time contributor to SciPy but believe this will be an important contribution and am open to any advice, criticism, or comments on the code or process.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Feb 1, 2014

Member

This seems to have a merge issue. Could you rebase it onto latest master?

Member

rgommers commented Feb 1, 2014

This seems to have a merge issue. Could you rebase it onto latest master?

@rgommers

This comment has been minimized.

Show comment
Hide comment
Member

rgommers commented Feb 1, 2014

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Feb 1, 2014

Member

If there's no user-visible changes yet, I'd prefer to merge this after branching 0.14.x (planned for ~ last week of Feb)

Member

rgommers commented Feb 1, 2014

If there's no user-visible changes yet, I'd prefer to merge this after branching 0.14.x (planned for ~ last week of Feb)

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Feb 1, 2014

Member

I think it will be necessary to have a more complete plan for the rest of the feature, before we can consider merging this. The important questions:

  • how does the user pass in low-level functions to quad()?
  • how are these represented on the Python side? ctypes?

Notes on the present code:

  • use of global variables like this is not reentrant or threadsafe
  • that non-reentrancy doesn't break the test suite is due to the fact that funcwrapper_init is always called with the same function quad_function. The following test program crashes now:
from scipy import integrate
import ctypes
lib = ctypes.CDLL('libm.so')
lib.sin.restype = ctypes.c_double
lib.sin.argtypes = (ctypes.c_double,)
def b(y):
    return y + integrate.quad(lib.sin, 0, 1)[0]
print(integrate.quad(b, 0, 1))
  • C coding style needs refinement. Please run the cwrapper.h through GNU indent with indent -kr -bad -sc -nce cwrapper.h to make it more readable. Also, prefer foo_bar_quux over foobarquux in variable naming.
  • The functions dq* are not directly callable by users from C code, so it's not clear what is gained by adding another wrapper layer to this set of functions.
Member

pv commented Feb 1, 2014

I think it will be necessary to have a more complete plan for the rest of the feature, before we can consider merging this. The important questions:

  • how does the user pass in low-level functions to quad()?
  • how are these represented on the Python side? ctypes?

Notes on the present code:

  • use of global variables like this is not reentrant or threadsafe
  • that non-reentrancy doesn't break the test suite is due to the fact that funcwrapper_init is always called with the same function quad_function. The following test program crashes now:
from scipy import integrate
import ctypes
lib = ctypes.CDLL('libm.so')
lib.sin.restype = ctypes.c_double
lib.sin.argtypes = (ctypes.c_double,)
def b(y):
    return y + integrate.quad(lib.sin, 0, 1)[0]
print(integrate.quad(b, 0, 1))
  • C coding style needs refinement. Please run the cwrapper.h through GNU indent with indent -kr -bad -sc -nce cwrapper.h to make it more readable. Also, prefer foo_bar_quux over foobarquux in variable naming.
  • The functions dq* are not directly callable by users from C code, so it's not clear what is gained by adding another wrapper layer to this set of functions.
@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Feb 1, 2014

Contributor

I'm working with Brian on this, and I think I can answer some of the questions. We suspected the reentrancy problem, but neither of us is really clear on how to resolve it. We'll be looking at Oliphant's original implementation to see if we can figure it out, but any direction we can get on the issue would be welcome.

As for the plan going forward. We think the preferred end-goal is to allow users to code their functions in cython, since that's the direction SciPy seems to be going, however getting there from the status quo isn't trivial, and will likely take a few steps. Our immediate next-step would be a very small change in end-user functionality: getting a multivariate ctypes function of the form f(nargs,args[nargs]) to pass directly to quadpack, in much the same way as is possible now. This builds on already-implemented concepts and ideas, and is a small enough project to be manageable for us. Additionally, it provides a good entry point for any future improvements, since any new code will have to provide a pointer to a function f(nargs args[nargs]), anyway, if it's going to use this backend.

This is deep water for us already, so any advice and/or help would be welcome.

Contributor

woodscn commented Feb 1, 2014

I'm working with Brian on this, and I think I can answer some of the questions. We suspected the reentrancy problem, but neither of us is really clear on how to resolve it. We'll be looking at Oliphant's original implementation to see if we can figure it out, but any direction we can get on the issue would be welcome.

As for the plan going forward. We think the preferred end-goal is to allow users to code their functions in cython, since that's the direction SciPy seems to be going, however getting there from the status quo isn't trivial, and will likely take a few steps. Our immediate next-step would be a very small change in end-user functionality: getting a multivariate ctypes function of the form f(nargs,args[nargs]) to pass directly to quadpack, in much the same way as is possible now. This builds on already-implemented concepts and ideas, and is a small enough project to be manageable for us. Additionally, it provides a good entry point for any future improvements, since any new code will have to provide a pointer to a function f(nargs args[nargs]), anyway, if it's going to use this backend.

This is deep water for us already, so any advice and/or help would be welcome.

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Feb 3, 2014

Member

The trick in making it re-entrant is to save a copy of the global state before calling the function, and restore it afterward (this is what the code currently in quad does). This is most convenient to do by putting all of the state in a struct (or, better, making the global a pointer to a struct). To make it threadsafe, one should use thread-local storage instead, but I'm not sure what's the portable way to do that.

The simplest way to proceed for the rest is probably (i) write a class ScipyLowLevelFunction in Python, that does the checking and preprocessing currently done in get_func_type in the quad() wrappers, plus extend it for functions accepting extra arguments. (ii) make a clone of the mechanism currently in quad() and modify it so that it makes use of the checked results computed by ScipyLowLevelFunction, in addition to normal Python functions. For dealing with extra arguments, a number of trampoline functions similar to your call are necessary, one for each combination of possible additional argument types (plus one for usual Python functions).

Member

pv commented Feb 3, 2014

The trick in making it re-entrant is to save a copy of the global state before calling the function, and restore it afterward (this is what the code currently in quad does). This is most convenient to do by putting all of the state in a struct (or, better, making the global a pointer to a struct). To make it threadsafe, one should use thread-local storage instead, but I'm not sure what's the portable way to do that.

The simplest way to proceed for the rest is probably (i) write a class ScipyLowLevelFunction in Python, that does the checking and preprocessing currently done in get_func_type in the quad() wrappers, plus extend it for functions accepting extra arguments. (ii) make a clone of the mechanism currently in quad() and modify it so that it makes use of the checked results computed by ScipyLowLevelFunction, in addition to normal Python functions. For dealing with extra arguments, a number of trampoline functions similar to your call are necessary, one for each combination of possible additional argument types (plus one for usual Python functions).

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 11, 2014

Coverage Status

Changes Unknown when pulling d1fc3c3 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling d1fc3c3 on BrianNewsom:cwrapping into * on scipy:master*.

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Feb 11, 2014

Contributor

Sorry about the messed up history, everyone involved. Very strange things are happening and I am working on getting them fixed.

Brian

Contributor

BrianNewsom commented Feb 11, 2014

Sorry about the messed up history, everyone involved. Very strange things are happening and I am working on getting them fixed.

Brian

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 11, 2014

Coverage Status

Changes Unknown when pulling 1bb717d on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 1bb717d on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 11, 2014

Coverage Status

Changes Unknown when pulling 1bb717d on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 1bb717d on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 12, 2014

Coverage Status

Changes Unknown when pulling 2803f3e on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 2803f3e on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 12, 2014

Coverage Status

Changes Unknown when pulling 2803f3e on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 2803f3e on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 12, 2014

Coverage Status

Changes Unknown when pulling 2803f3e on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 2803f3e on BrianNewsom:cwrapping into * on scipy:master*.

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Feb 12, 2014

Contributor

Hello all,
I believe that I have taken in all of your comments and remedied the situations. I rebased and fixed the ugly commits.

I have also fixed the re-entrancy issue as described by @pv . The cwrapper.h file has been removed, as further examination shows that it creates overhead and the solution doesn't necessitate those wrappers.

So now I just have the extra two wrappers both allowing re-entrancy so that functions of form f(n, args[n]) can be accepted. I agree that steps should be taken towards getting that function pointer in cython, but as @thezealite described, immediate next steps should be allowing ctypes of this form as is done now for f(*x) functions.

Brian

Contributor

BrianNewsom commented Feb 12, 2014

Hello all,
I believe that I have taken in all of your comments and remedied the situations. I rebased and fixed the ugly commits.

I have also fixed the re-entrancy issue as described by @pv . The cwrapper.h file has been removed, as further examination shows that it creates overhead and the solution doesn't necessitate those wrappers.

So now I just have the extra two wrappers both allowing re-entrancy so that functions of form f(n, args[n]) can be accepted. I agree that steps should be taken towards getting that function pointer in cython, but as @thezealite described, immediate next steps should be allowing ctypes of this form as is done now for f(*x) functions.

Brian

@pv pv added the PR label Feb 19, 2014

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Feb 23, 2014

Contributor

Hello,

This PR now is functional, as it detects the ctypes function signature and handles it appropriately, producing a 2x increase in speed in the integrations I've tested so far.

Use:
Function written in c:

double new(int n, double args[n]){
    return args[0]*args[0] + args[1]*args[1] + args[0] * args[1] * args[2] + args[3]*args[3]*args[3]      +args[4];
}

Compile to library
Set restype and argtypes as in ctypes 1d version in Python:

lib = ctypes.CDLL('lib.so')
lib.new.restype = ctypes.c_double
lib.new.argtypes = (ctypes.c_int,ctypes.c_double*4)

As in 1d ctypes, run nquad on the function:

integrate.nquad(testlib.new,[[-1,1],[-1,1],[-1,1],[-1,1],[1,1]])

For this specific example the integration occurs at about a 1.9x faster speed.

I will try and get some test cases into the integrate unit testing when I get a chance. In the meantime, feel free to try it out and/or take a look at the code, and please bring up any concerns or deficiencies found.

Thanks,
Brian

Contributor

BrianNewsom commented Feb 23, 2014

Hello,

This PR now is functional, as it detects the ctypes function signature and handles it appropriately, producing a 2x increase in speed in the integrations I've tested so far.

Use:
Function written in c:

double new(int n, double args[n]){
    return args[0]*args[0] + args[1]*args[1] + args[0] * args[1] * args[2] + args[3]*args[3]*args[3]      +args[4];
}

Compile to library
Set restype and argtypes as in ctypes 1d version in Python:

lib = ctypes.CDLL('lib.so')
lib.new.restype = ctypes.c_double
lib.new.argtypes = (ctypes.c_int,ctypes.c_double*4)

As in 1d ctypes, run nquad on the function:

integrate.nquad(testlib.new,[[-1,1],[-1,1],[-1,1],[-1,1],[1,1]])

For this specific example the integration occurs at about a 1.9x faster speed.

I will try and get some test cases into the integrate unit testing when I get a chance. In the meantime, feel free to try it out and/or take a look at the code, and please bring up any concerns or deficiencies found.

Thanks,
Brian

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 23, 2014

Coverage Status

Coverage remained the same when pulling 0be6bc0 on BrianNewsom:cwrapping into dfd2e2b on scipy:master.

Coverage Status

Coverage remained the same when pulling 0be6bc0 on BrianNewsom:cwrapping into dfd2e2b on scipy:master.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Feb 23, 2014

Coverage Status

Coverage remained the same when pulling 0be6bc0 on BrianNewsom:cwrapping into dfd2e2b on scipy:master.

Coverage Status

Coverage remained the same when pulling 0be6bc0 on BrianNewsom:cwrapping into dfd2e2b on scipy:master.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Feb 28, 2014

Contributor

Brian, I notice that the Travis CI build failed. You should try to get that to finish.

Contributor

woodscn commented Feb 28, 2014

Brian, I notice that the Travis CI build failed. You should try to get that to finish.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Mar 4, 2014

Coverage Status

Coverage remained the same when pulling e3704e5 on BrianNewsom:cwrapping into c81c978 on scipy:master.

Coverage Status

Coverage remained the same when pulling e3704e5 on BrianNewsom:cwrapping into c81c978 on scipy:master.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Mar 18, 2014

Coverage Status

Changes Unknown when pulling 6eaefc7 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 6eaefc7 on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Mar 19, 2014

Coverage Status

Changes Unknown when pulling fb42e11 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling fb42e11 on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Apr 1, 2014

Coverage Status

Changes Unknown when pulling dd8d3d9 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling dd8d3d9 on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Apr 9, 2014

Coverage Status

Changes Unknown when pulling 7b1e4e6 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 7b1e4e6 on BrianNewsom:cwrapping into * on scipy:master*.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Apr 9, 2014

Contributor

Building a unit test for this feature seems difficult, since it would require a C library that could be loaded by ctypes. How would such a library be distributed? The need for a cython version of this is becoming more and more clear.

Any further comments? Has anyone had a chance to look over this more closely?

Contributor

woodscn commented Apr 9, 2014

Building a unit test for this feature seems difficult, since it would require a C library that could be loaded by ctypes. How would such a library be distributed? The need for a cython version of this is becoming more and more clear.

Any further comments? Has anyone had a chance to look over this more closely?

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Apr 10, 2014

Coverage Status

Changes Unknown when pulling 5e06ea6 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 5e06ea6 on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Apr 12, 2014

Coverage Status

Changes Unknown when pulling d8fa178 on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling d8fa178 on BrianNewsom:cwrapping into * on scipy:master*.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Apr 29, 2014

Coverage Status

Changes Unknown when pulling 611615d on BrianNewsom:cwrapping into * on scipy:master*.

Coverage Status

Changes Unknown when pulling 611615d on BrianNewsom:cwrapping into * on scipy:master*.

@rkern

This comment has been minimized.

Show comment
Hide comment
@rkern

rkern Apr 29, 2014

Member

Please follow the stdlib's ctypes test more closely. Make _test_multivariate an actual Python extension module. Import this module in the test suite. Use its __file__ attribute to pass to CDLL(); you don't need to do all of that filename munging to guess the name. You must make sure to EXPORT the functions that you are exposing for this to work on Windows; again look at the stdlib's ctypes test module for the appropriate macros.

Member

rkern commented Apr 29, 2014

Please follow the stdlib's ctypes test more closely. Make _test_multivariate an actual Python extension module. Import this module in the test suite. Use its __file__ attribute to pass to CDLL(); you don't need to do all of that filename munging to guess the name. You must make sure to EXPORT the functions that you are exposing for this to work on Windows; again look at the stdlib's ctypes test module for the appropriate macros.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Apr 29, 2014

Contributor

Thanks for the comments. I was very unclear on what EXPORT was doing, and I suppose that’s because I’m not running Windows.

As for the filename work, I was really hoping to avoid mucking around in the C API any more than absolutely necessary. I guess we’ll have to add that to our list of things to learn now.

Nathan Woods

On Apr 29, 2014, at 4:00 PM, Robert Kern notifications@github.com wrote:

Please follow the stdlib's ctypes test more closely. Make _test_multivariate an actual Python extension module. Import this module in the test suite. Use its file attribute to pass to CDLL(); you don't need to do all of that filename munging to guess the name. You must make sure to EXPORT the functions that you are exposing for this to work on Windows; again look at the stdlib's ctypes test module for the appropriate macros.


Reply to this email directly or view it on GitHub.

Contributor

woodscn commented Apr 29, 2014

Thanks for the comments. I was very unclear on what EXPORT was doing, and I suppose that’s because I’m not running Windows.

As for the filename work, I was really hoping to avoid mucking around in the C API any more than absolutely necessary. I guess we’ll have to add that to our list of things to learn now.

Nathan Woods

On Apr 29, 2014, at 4:00 PM, Robert Kern notifications@github.com wrote:

Please follow the stdlib's ctypes test more closely. Make _test_multivariate an actual Python extension module. Import this module in the test suite. Use its file attribute to pass to CDLL(); you don't need to do all of that filename munging to guess the name. You must make sure to EXPORT the functions that you are exposing for this to work on Windows; again look at the stdlib's ctypes test module for the appropriate macros.


Reply to this email directly or view it on GitHub.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Apr 30, 2014

Contributor

All right, first things first. I've been feeling somewhat ungrateful, especially to @rkern, who has been extremely helpful the last few days with this. Once I figured out what was going on (and why), your comments seem spot on, so thanks again for being patient with me.

Second, I noticed that the way modules are implemented in the C API changed with Python 3. I think I've figured it out for Python 2, which you can see in a PR to Brian's fork. I'm not sure if/how that might work with Python 3, though, or how that problem should be handled.

Nathan Woods

Contributor

woodscn commented Apr 30, 2014

All right, first things first. I've been feeling somewhat ungrateful, especially to @rkern, who has been extremely helpful the last few days with this. Once I figured out what was going on (and why), your comments seem spot on, so thanks again for being patient with me.

Second, I noticed that the way modules are implemented in the C API changed with Python 3. I think I've figured it out for Python 2, which you can see in a PR to Brian's fork. I'm not sure if/how that might work with Python 3, though, or how that problem should be handled.

Nathan Woods

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Apr 30, 2014

Coverage Status

Coverage remained the same when pulling fcea3af on BrianNewsom:cwrapping into d9a8c21 on scipy:master.

Coverage Status

Coverage remained the same when pulling fcea3af on BrianNewsom:cwrapping into d9a8c21 on scipy:master.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls May 1, 2014

Coverage Status

Coverage remained the same when pulling c68abf7 on BrianNewsom:cwrapping into d9a8c21 on scipy:master.

Coverage Status

Coverage remained the same when pulling c68abf7 on BrianNewsom:cwrapping into d9a8c21 on scipy:master.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls May 14, 2014

Coverage Status

Coverage remained the same when pulling 0847c38 on BrianNewsom:cwrapping into d9a8c21 on scipy:master.

Coverage Status

Coverage remained the same when pulling 0847c38 on BrianNewsom:cwrapping into d9a8c21 on scipy:master.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers May 28, 2014

Member

I've done some testing, and it works for me on 32-bit Linux (after the build fix I sent BrianNewsom#3 for). Except for this spurious failure, which happened only once in ~100 test runs:

======================================================================
FAIL: test_improvement (test_quadpack.TestMultivariateCtypesQuad)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/rgommers/.local/lib/python2.7/site-packages/numpy/testing/decorators.py", line 146, in skipper_func
    return f(*args, **kwargs)
  File "/home/rgommers/Code/scipy/scipy/integrate/tests/test_quadpack.py", line 121, in test_improvement
    assert_(fast < .666 * slow, (fast, slow))
  File "/home/rgommers/.local/lib/python2.7/site-packages/numpy/testing/utils.py", line 34, in assert_
    raise AssertionError(msg)
AssertionError: (0.004848003387451172, 0.005922079086303711)

It looks the same as gh-2209. You see that fast is still faster than slow, just not fast enough. This kind of test should probably simply use inputs that cause longer runtimes; 5 ms isn't a lot compared to timing variations due to other processes. Or time.time is just not that reliable. I suggest to change the test to run 10x slower, then merge and see if there are still failures reported.

Member

rgommers commented May 28, 2014

I've done some testing, and it works for me on 32-bit Linux (after the build fix I sent BrianNewsom#3 for). Except for this spurious failure, which happened only once in ~100 test runs:

======================================================================
FAIL: test_improvement (test_quadpack.TestMultivariateCtypesQuad)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/rgommers/.local/lib/python2.7/site-packages/numpy/testing/decorators.py", line 146, in skipper_func
    return f(*args, **kwargs)
  File "/home/rgommers/Code/scipy/scipy/integrate/tests/test_quadpack.py", line 121, in test_improvement
    assert_(fast < .666 * slow, (fast, slow))
  File "/home/rgommers/.local/lib/python2.7/site-packages/numpy/testing/utils.py", line 34, in assert_
    raise AssertionError(msg)
AssertionError: (0.004848003387451172, 0.005922079086303711)

It looks the same as gh-2209. You see that fast is still faster than slow, just not fast enough. This kind of test should probably simply use inputs that cause longer runtimes; 5 ms isn't a lot compared to timing variations due to other processes. Or time.time is just not that reliable. I suggest to change the test to run 10x slower, then merge and see if there are still failures reported.

@rgommers rgommers added this to the 0.15.0 milestone May 28, 2014

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers May 28, 2014

Member

Can you rebase this? Then it'll be mergeable again, and all the merge commits will be gone. It would be good to also remove c5b5db9 and 5d648a0

Member

rgommers commented May 28, 2014

Can you rebase this? Then it'll be mergeable again, and all the merge commits will be gone. It would be good to also remove c5b5db9 and 5d648a0

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers May 31, 2014

Member

@woodscn if you want to do it, that would be helpful. I was thinking maybe 3 commits:

  • new functionality
  • tests
  • THANKS and 0.15.0 release notes (the latter is not in here yet, a couple of sentences are justified though)

You should use git commit --author=... to ensure the right author name is kept. Looks like that should be Brian for functionality and Nathan for tests.

Member

rgommers commented May 31, 2014

@woodscn if you want to do it, that would be helpful. I was thinking maybe 3 commits:

  • new functionality
  • tests
  • THANKS and 0.15.0 release notes (the latter is not in here yet, a couple of sentences are justified though)

You should use git commit --author=... to ensure the right author name is kept. Looks like that should be Brian for functionality and Nathan for tests.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn May 31, 2014

Contributor

I've put things together and pushed it to https://github.com/woodscn/scipy/tree/fix_cwrapping. I've asked Brian to look over how I split up author attribution before it's final, but I'd like to make sure that I got the right idea, if anyone wants to look it over quickly.

Contributor

woodscn commented May 31, 2014

I've put things together and pushed it to https://github.com/woodscn/scipy/tree/fix_cwrapping. I've asked Brian to look over how I split up author attribution before it's final, but I'd like to make sure that I got the right idea, if anyone wants to look it over quickly.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 1, 2014

Member

That branch looks good to me.

Member

rgommers commented Jun 1, 2014

That branch looks good to me.

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Jun 1, 2014

Contributor

Looks great to me!

Contributor

BrianNewsom commented Jun 1, 2014

Looks great to me!

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 2, 2014

Member

@BrianNewsom could you push that branch to this PR?

Member

rgommers commented Jun 2, 2014

@BrianNewsom could you push that branch to this PR?

BrianNewsom and others added some commits May 31, 2014

ENH: Implement faster multivariate integration in scipy.integrate.
The numerical integration routines in SciPy use a well-tested, fast, Fortran
library. The computational cost of integration with these routines can be
greatly reduced if the function to be integrated is compiled and passed to
the library directly, avoiding callbacks to the Python interpreter.
@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Jun 2, 2014

Contributor

@rgommers I believe it's done now!

Contributor

BrianNewsom commented Jun 2, 2014

@rgommers I believe it's done now!

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Jun 2, 2014

Coverage Status

Coverage increased (+0.0%) when pulling 51bbf4b on BrianNewsom:cwrapping into 9f89371 on scipy:master.

Coverage Status

Coverage increased (+0.0%) when pulling 51bbf4b on BrianNewsom:cwrapping into 9f89371 on scipy:master.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Jun 2, 2014

Coverage Status

Coverage increased (+0.0%) when pulling 51bbf4b on BrianNewsom:cwrapping into 9f89371 on scipy:master.

Coverage Status

Coverage increased (+0.0%) when pulling 51bbf4b on BrianNewsom:cwrapping into 9f89371 on scipy:master.

@ev-br

This comment has been minimized.

Show comment
Hide comment
@ev-br

ev-br Jun 2, 2014

Member

A comment from the peanut gallery for this PR: it would be good to have some tutorial / user-manual type docs for this feature. For example, I as a user have an expensive function to integrate, what do I do? If I already have a function in cython (or Fortran), I can always wrap it up as a python callable and feed it to quad, is there a better way now? What about multivariate functions, what is the status of dblquad, tplquad, nquad?

Member

ev-br commented Jun 2, 2014

A comment from the peanut gallery for this PR: it would be good to have some tutorial / user-manual type docs for this feature. For example, I as a user have an expensive function to integrate, what do I do? If I already have a function in cython (or Fortran), I can always wrap it up as a python callable and feed it to quad, is there a better way now? What about multivariate functions, what is the status of dblquad, tplquad, nquad?

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Jun 2, 2014

Coverage Status

Coverage increased (+0.0%) when pulling 51bbf4b on BrianNewsom:cwrapping into 9f89371 on scipy:master.

Coverage Status

Coverage increased (+0.0%) when pulling 51bbf4b on BrianNewsom:cwrapping into 9f89371 on scipy:master.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Jun 3, 2014

Contributor

As a quick answer to Evgeni, the docs for quad and nquad include give some quick instructions. From nquad:

  •    If speed is desired, this function may be a ctypes function of
    
    •    the form::
      
    •        f(int n, double args[n])
      
    •    where `n` is the number of extra parameters and args is an array
      
    •    of doubles of the additional parameters.  This function may then
      
    •    be compiled to a dynamic/shared library then imported through
      
    •    `ctypes`, setting the function's argtypes to `(c_int, c_double)`,
      
    •    and the function's restype to `(c_double)`.  Its pointer may then be
      
    •    passed into `nquad` normally.
      

Unfortunately, ctypes is the only interface that works right now, so if you have a Cython or Fortran function, you would have to find some way to wrap it with ctypes first. With Fortran that's not so much of a problem. We've been looking at Cython, but haven't had any luck with that yet.

This was explicitly designed for use with multivariate functions, especially with nquad. Since the changes are the level of the interface between quad and quadpack, they should work just fine with dblquad and tplquad too.

The performance increases here come from two main sources. First, faster evaluation of compiled vs interpreted functions. Second, reduced overhead because the quadpack library doesn't have to callback to Python to evaluate the integrand function. If those are bottlenecks in a particular application, then it may be worthwhile for people to use the extra ctypes path. In my case, I need to evaluate volume integrals of messy analytic functions, and the speed improvements from compiling the integrand and cutting overhead in the innermost loop are significant.

As for including a tutorial of sorts, we could try to put something together, but my first thought is that there are a lot of ways for someone to get a ctypes pointer, so it would either be incredibly general, or else it would encourage specific design choices.

What are everyone's thoughts? If there's nothing else, is this ready to merge?

Contributor

woodscn commented Jun 3, 2014

As a quick answer to Evgeni, the docs for quad and nquad include give some quick instructions. From nquad:

  •    If speed is desired, this function may be a ctypes function of
    
    •    the form::
      
    •        f(int n, double args[n])
      
    •    where `n` is the number of extra parameters and args is an array
      
    •    of doubles of the additional parameters.  This function may then
      
    •    be compiled to a dynamic/shared library then imported through
      
    •    `ctypes`, setting the function's argtypes to `(c_int, c_double)`,
      
    •    and the function's restype to `(c_double)`.  Its pointer may then be
      
    •    passed into `nquad` normally.
      

Unfortunately, ctypes is the only interface that works right now, so if you have a Cython or Fortran function, you would have to find some way to wrap it with ctypes first. With Fortran that's not so much of a problem. We've been looking at Cython, but haven't had any luck with that yet.

This was explicitly designed for use with multivariate functions, especially with nquad. Since the changes are the level of the interface between quad and quadpack, they should work just fine with dblquad and tplquad too.

The performance increases here come from two main sources. First, faster evaluation of compiled vs interpreted functions. Second, reduced overhead because the quadpack library doesn't have to callback to Python to evaluate the integrand function. If those are bottlenecks in a particular application, then it may be worthwhile for people to use the extra ctypes path. In my case, I need to evaluate volume integrals of messy analytic functions, and the speed improvements from compiling the integrand and cutting overhead in the innermost loop are significant.

As for including a tutorial of sorts, we could try to put something together, but my first thought is that there are a lot of ways for someone to get a ctypes pointer, so it would either be incredibly general, or else it would encourage specific design choices.

What are everyone's thoughts? If there's nothing else, is this ready to merge?

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Jun 3, 2014

Contributor

@ev-br What @woodscn is saying is correct. Although I would like to expand on the cython question.

I believe this route provides two (primary) advantages over just compiling to Cython and passing it into quad - as this will work already.

1.) Removal of Python callbacks. I'm pretty sure that although Cython compiles to C, the way it's handled requires calls to Python for evaluation at each inner iteration of the integration loop. The compilation will provide speed increases, but the time necessary for callbacks between langauges remains because of the way the quadpack code is written as of now. This may be possible to change, but that is another topic.

2.) Ease of use. For users that have not worked with Cython much (e.g. me), the ctypes method presented here consists of a fairly concise set of steps that can be repeated without delving too much into what is going on below. Maybe that is just my familiarity with Ctypes, but it allows a user to avoid the Cython building infrastructure which can be confusing.

Everything I wrote applying to Cython also is true for Fortran using, say, F2PY. I have run timing tests and such on this and this method DOES provide a ~2x increase on simply compiled fortran, where this is provided by the removal of python/c callbacks.

Contributor

BrianNewsom commented Jun 3, 2014

@ev-br What @woodscn is saying is correct. Although I would like to expand on the cython question.

I believe this route provides two (primary) advantages over just compiling to Cython and passing it into quad - as this will work already.

1.) Removal of Python callbacks. I'm pretty sure that although Cython compiles to C, the way it's handled requires calls to Python for evaluation at each inner iteration of the integration loop. The compilation will provide speed increases, but the time necessary for callbacks between langauges remains because of the way the quadpack code is written as of now. This may be possible to change, but that is another topic.

2.) Ease of use. For users that have not worked with Cython much (e.g. me), the ctypes method presented here consists of a fairly concise set of steps that can be repeated without delving too much into what is going on below. Maybe that is just my familiarity with Ctypes, but it allows a user to avoid the Cython building infrastructure which can be confusing.

Everything I wrote applying to Cython also is true for Fortran using, say, F2PY. I have run timing tests and such on this and this method DOES provide a ~2x increase on simply compiled fortran, where this is provided by the removal of python/c callbacks.

@ev-br

This comment has been minimized.

Show comment
Hide comment
@ev-br

ev-br Jun 3, 2014

Member

@BrianNewsom @woodscn You seem to be implying that I criticize your design decisions or question the choice of ctypes or advocate for switching over to cython --- I am not.

All I am saying is that I think it would make sense to add an example of using this feature to eg this tutorial: http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html
This addition need not be large. A slightly expanded example from the docstrings of quad and nquad would already be very helpful for us mere mortals with nasty functions to integrate and without extensive experience with ctypes pointers.

Note that I am not saying this needs to be done, and I do not mean to block or postpone merging this PR. I just think this feature should be advertised more, and that now seems to be a good time to write things up.

Member

ev-br commented Jun 3, 2014

@BrianNewsom @woodscn You seem to be implying that I criticize your design decisions or question the choice of ctypes or advocate for switching over to cython --- I am not.

All I am saying is that I think it would make sense to add an example of using this feature to eg this tutorial: http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html
This addition need not be large. A slightly expanded example from the docstrings of quad and nquad would already be very helpful for us mere mortals with nasty functions to integrate and without extensive experience with ctypes pointers.

Note that I am not saying this needs to be done, and I do not mean to block or postpone merging this PR. I just think this feature should be advertised more, and that now seems to be a good time to write things up.

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Jun 3, 2014

Contributor

@ev-br I apologize if that came off as defensive, I was just trying to explain why this could be an improvement over cython - why it is a viable new option.

I absolutely agree that we could add something into the tutorial, it is definitely a slightly confusing feature meant for adventurous users and it would be great to throw a more explicit and elaborate explanation into the documentation. Thanks for the suggestion!

Contributor

BrianNewsom commented Jun 3, 2014

@ev-br I apologize if that came off as defensive, I was just trying to explain why this could be an improvement over cython - why it is a viable new option.

I absolutely agree that we could add something into the tutorial, it is definitely a slightly confusing feature meant for adventurous users and it would be great to throw a more explicit and elaborate explanation into the documentation. Thanks for the suggestion!

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Jun 4, 2014

Contributor

Added a more explicit explanation in the tutorials section. Take a look, let me know if this is helpful.

Contributor

BrianNewsom commented Jun 4, 2014

Added a more explicit explanation in the tutorials section. Take a look, let me know if this is helpful.

@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Jun 4, 2014

Coverage Status

Coverage increased (+0.0%) when pulling c8f742d on BrianNewsom:cwrapping into 9f89371 on scipy:master.

Coverage Status

Coverage increased (+0.0%) when pulling c8f742d on BrianNewsom:cwrapping into 9f89371 on scipy:master.

@ev-br

This comment has been minimized.

Show comment
Hide comment
@ev-br

ev-br Jun 4, 2014

It's definitely helpful, thanks! I intend to try it out --- when I've time (tm), and it is unlikely to happen until next week at the earliest. At any rate, it shouldn't postpone merging this PR. The doc tweaks can be done separately.

ev-br commented on c8f742d Jun 4, 2014

It's definitely helpful, thanks! I intend to try it out --- when I've time (tm), and it is unlikely to happen until next week at the earliest. At any rate, it shouldn't postpone merging this PR. The doc tweaks can be done separately.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Jun 8, 2014

Contributor

Is this just waiting for merge now?

Contributor

woodscn commented Jun 8, 2014

Is this just waiting for merge now?

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 8, 2014

Member

@woodscn I think so. I'll try to make some time this week (but I'm traveling, may take a few days).

Member

rgommers commented Jun 8, 2014

@woodscn I think so. I'll try to make some time this week (but I'm traveling, may take a few days).

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Jun 17, 2014

Contributor

Pinging again.

Contributor

woodscn commented Jun 17, 2014

Pinging again.

@BrianNewsom

This comment has been minimized.

Show comment
Hide comment
@BrianNewsom

BrianNewsom Jun 17, 2014

Contributor

If you'd like me to rollback the most recent commit (it was suggested to have it done separately) just let me know. Otherwise I'll leave it.

Contributor

BrianNewsom commented Jun 17, 2014

If you'd like me to rollback the most recent commit (it was suggested to have it done separately) just let me know. Otherwise I'll leave it.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 17, 2014

Member

Leave it in I'd say - doc updates are an integral part of PRs like this, and it helps me to play with the code.

Member

rgommers commented Jun 17, 2014

Leave it in I'd say - doc updates are an integral part of PRs like this, and it helps me to play with the code.

@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Jul 2, 2014

Contributor

Ping!

Contributor

woodscn commented Jul 2, 2014

Ping!

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jul 5, 2014

Member

Tutorial looks good to me. I verified the statements on speedup: example given is ~1.5x while np.sin(args[0] - np.exp(-args[1])) is more than 10x faster.

Member

rgommers commented Jul 5, 2014

Tutorial looks good to me. I verified the statements on speedup: example given is ~1.5x while np.sin(args[0] - np.exp(-args[1])) is more than 10x faster.

rgommers added a commit that referenced this pull request Jul 5, 2014

Merge branch 'pr/3262' into master.
Adds support for direct callbacks to C functions via ctypes
for multivariate integration.

Reviewed at #3262
@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jul 5, 2014

Member

Squashed the two commits that touched the release notes, fixed tutorial formatting and merged in 04cc74f.

Thanks a lot all!

Member

rgommers commented Jul 5, 2014

Squashed the two commits that touched the release notes, fixed tutorial formatting and merged in 04cc74f.

Thanks a lot all!

@rgommers rgommers closed this Jul 5, 2014

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jul 5, 2014

Member

Build fails on Python 3.4, which has been enabled on TravisCI after running the last tests on this PR.

I'll try to fix it tomorrow (football will get in the way tonight). Should be straightforward:

In file included from scipy/integrate/_quadpackmodule.c:4:0:
scipy/integrate/quadpack.h: In function ‘c_array_from_tuple’:
scipy/integrate/quadpack.h:131:5: error: ISO C90 forbids mixed declarations and
code [-Werror=declaration-after-statement]
Member

rgommers commented Jul 5, 2014

Build fails on Python 3.4, which has been enabled on TravisCI after running the last tests on this PR.

I'll try to fix it tomorrow (football will get in the way tonight). Should be straightforward:

In file included from scipy/integrate/_quadpackmodule.c:4:0:
scipy/integrate/quadpack.h: In function ‘c_array_from_tuple’:
scipy/integrate/quadpack.h:131:5: error: ISO C90 forbids mixed declarations and
code [-Werror=declaration-after-statement]
@woodscn

This comment has been minimized.

Show comment
Hide comment
@woodscn

woodscn Jul 6, 2014

Contributor

I tried to see if I could fix it, but building on my machine didn't reproduce the error, and I didn't want to mess without being able to check the results.

Contributor

woodscn commented Jul 6, 2014

I tried to see if I could fix it, but building on my machine didn't reproduce the error, and I didn't want to mess without being able to check the results.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jul 6, 2014

Member

Looks like this is due to http://bugs.python.org/issue21121; -Werror=declaration-after-statement was by accident added to flags that distutils passes on for 3.4. Is already reversed. I can't reproduce it either, but I'll still send a fix.

Member

rgommers commented Jul 6, 2014

Looks like this is due to http://bugs.python.org/issue21121; -Werror=declaration-after-statement was by accident added to flags that distutils passes on for 3.4. Is already reversed. I can't reproduce it either, but I'll still send a fix.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jul 6, 2014

Member

fixed in gh-3774

Member

rgommers commented Jul 6, 2014

fixed in gh-3774

+ return NPY_FAIL;
+ }
+ int n_args_ref = PyTuple_Size(args);
+ global_n_args = &n_args_ref;

This comment has been minimized.

@pitrou

pitrou Apr 14, 2016

This stores a global pointer to a local variable that will soon go out of scope. Clearly the value dereferenced later will be non-sensical.

@pitrou

pitrou Apr 14, 2016

This stores a global pointer to a local variable that will soon go out of scope. Clearly the value dereferenced later will be non-sensical.

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