Steven Johnson's Cubature routines...

["Robert G. Brown" <rgb at phy dot duke dot edu>]

Dear GSLers,

OK, as I noted over the summer, a user here really needs/wants cubature
(or else they're going to have to buy NAG:-(.  To make them (maybe)
happy, I've worked pretty hard taking Steven G. Johnson's port of the
HintLib rule75genzmalik adaptive integration routine and turning it into
an rpm-ified or at least rpm-ifiable library.  

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).

Anybody who wants to try out the result can grab any of the following:

FC4 builds:

 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.

Here is where having multiple rules might be very handy.  I >>think<<
that the problem with Genz and Malik is that the routine does a local
gradient search for the "best dimension" to refine in via hypercube
bisection -- effectively a projection of the error gradient onto the
problem unit vectors.  

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...

Anyway, enjoy (if this is your sort of thing...:-).

   rgb

Attachment: pgp00000.pgp
Description: PGP signature