This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Other format: | [Raw text] |
At the same time I was a bit disappointed in its performance on the hypersphere integration problem and so I instrumented the code with some debugging stuff so I could see roughly what it is doing.
The library is NOT (yet) in a form suitable for direct inclusion in the GSL -- in my opinion the source file is too monolithic and the code really needs to be sanely partitioned and perhaps decrufted and otherwise cleaned up before even trying. Some parts (as Steve noted) are even FROM the GSL and need to be replaced with suitable GSL calls where I've left them alone to avoid "breaking" the code while porting it. This just takes more time than I've got to spend on it right now, but eventually I'll try to get to it, if somebody else doesn't do it first.
I do think cubature should be in the GSL and think that it should be there in an extensible way, where more than one rule is available (an interface similar to the rng interface, perhaps, where one can simply change back ends by altering a single parameter).
http://www.phy.duke.edu/~rgb/libcubature-0.1.0-1.i386.rpm http://www.phy.duke.edu/~rgb/libcubature-0.1.0-1.src.rpm
Tarball: http://www.phy.duke.edu/~rgb/libcubature-0.1.0.tgz
The binary rpm will install the library and include files and libcubature(3) and adapt_integrate(3) man pages as well as /usr/share/doc/libcubature-0.1.0/, which contains both a README (containing among other things Steven's original announcement to this list) and test_cubature.c, which is a VERY simple program you should be able to compile and run once you've installed the rpm (read the comment at the head of the source for very explicit instructions).
One hopes that the src rpm will rebuild for people who want to make a native x86_64 version. I can also do this on request -- I've got several Athlon 64's and Opterons handy to serve as build hosts -- so just ask.
The tarball is probably the most "useful" to anybody wanting to play, though, as it has both the library sources and a "test" subdirectory where all three of Steve's examples are programmed into an example shell called "multidim". Again, I instrumented the code a bit (and made it more readable for C novices as Steve likes the ? operator which can obfuscate just a bit:-). The most useful addition I made is I added a point dump feature -- e.g. multidim -i 0 -d 3 -m 100000 -p will do the hypersphere example problem in 3 dimensions with at most 100000 function evals AND will dump the coordinates of those function evals to stdout.
If you then put them into a file and do e.g. an xy projective scatterplot, you can see some of the "quirks" of the Genz and Malik algorithm when faced with a discontinuous integrand -- whole chunklets of the 3d sphere are missed even in 2d projection (where in 2d the routine does quite well on finding the circle boundary and refining).
On the continuous objectives (gaussian and product of cosines) the routine does really, really well, very, very quickly as Steve noted and is also evident from the projective point coverage. It really is lovely code for continuous integrands.
For a curved discontinuous surface with smooth angle relative to the cartesian directions, though, this means that some directions always win at the expense of other directions (because the "gradient" is really "binary" across the surface discontinuity) -- it doesn't form a "true" gradient at right angles to the surface and subdivide both directions. Apparently holes can appear when the points selected have the right modulus relative to the surface-spanning dimensions. Since the algorithm doesn't appear to be x-y-z symmetric, the resulting point distribution is both asymmetric and highly uneven in 2d projection, indicating that large parts of the surface just aren't sampled.
Or, of course, I may have just plain broken some of Steven's code turning it into a library, or he may have broken some of HintLib's code in some subtle way. But it works so WELL on smooth problems...
Attachment:
pgp00000.pgp
Description: PGP signature
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |