contribution: code for multi-dimensional adaptive integration

contribution: code for multi-dimensional adaptive integration

[Steven G. Johnson]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 5397 bytes --]

Dear GSL folks,

For some time now, I've felt that GSL has cried out for a multidimensional 
cubature routine, similar to its 1d adaptive quadrature routines.  Now, 
I've written one (~700 lines in C) based on the well-known Genz-Malik 
embedded cubature rule (with help from a GPL'ed library called HIntLib, 
see below), and I hope that you're interested to have it included in GSL. 
I've attached the unmodified code; I'm willing to do the work to modify it 
for your coding style, etcetera.

Adaptive cubature is a huge win over competing methods in GSL for moderate 
dimensionality (dimensions 2-6), as I demonstrate with a couple benchmarks 
below, and these kinds of dimensionalities are obviously important in lots 
of scientific applications.

My Genz-Malik cubature is based in part on code from a GPL'ed library 
called HIntLib (http://www.cosy.sbg.ac.at/~rschuer/hintlib/).  HIntLib is 
a much more sophisticated, full-featured package that is focused on 
extremely high-dimensional integration, including many other cubature and 
Monte-Carlo methods.  Unfortunately, this means it is a large and complex 
library that is overkill for low-dimensional integrals, and is not 
suitable for inclusion in GSL (also, it is in C++).  My implementation 
strips out all the dozens of C++ classes defined in HIntLib and implements 
just the serial adaptive cubature with just the degree-7 Genz-Malik rule 
for dimensions > 1 (although it would be easy to add more embedded 
cubature rules), and I think it should be easy to include in GSL.

Actually, Genz and Malik wrote their own implementation of their adaptive 
cubature rule as a Fortran code adapt.f which can be found as part of the 
public-domain CMLIB distributed by NIST.  Interestingly, however, I found 
that HIntLib's (and my) implementation of the same cubature rule to often 
be more robust -- the Genz Malik code seems to be easily "fooled" 
(converges to the wrong result) for discontinuous and/or strongly peaked 
objectives, at least in my experience.  (These problems were the whole 
reason I switched from adapt.f to HIntLib in my own applications.)

The only other alternative that I know of is CUBPACK, which is non-free.

Cordially,
Steven G. Johnson

PS. Note that my attached code includes a 1d quadrature based on GSL's 
15-point QAGS.  This for completeness, since the Genz-Malik rule does not 
work in 1d.  Obviously, a GSL version of my routine would just call the 
your 1d integration routine when dimensions = 1.

------------------------------benchmarks-------------------------------- 
Currently, to perform multi-dimensional integration in GSL you have to 
either use Monte-Carlo integration or recursively call 1d adaptive 
integrations.  Both of these are extremely inefficient in many cases 
compared to adaptive cubature, and to show this I've benchmarked my 
routine against GSL for two simple test functions (whose integrals I know 
analytically):

cos: a simple smooth integrand = product of cos(x[i]) for dimensions i

gaussian: a gaussian integral that tests facility with strongly peaked
           integrands, with the axes rescaled via x -> (1-x)/x so that
           the integration goes from 0 to infinity.  The integrand
 	  (before coordinate rescaling) is:
 	        product of 2/sqrt(pi) * exp(-x[i]^2)
           (which integrates to 1).

The integration volume was from x[i] = 0 to 1+0.4*sin(i) for the cos 
integrand, and from 0 to 1 for the gaussian integrand (rescaled to 
0..infinity).

In each case, I integrated to a relative error tolerance of 1e-3, and 
counted the function evaluations. (For the Monte-Carlo algorithms, I 
repeatedly doubled the number of iterations until the estimated error 
achieved the given relative error.)  I used the Mersenne Twister RNG. 
For recursive quadrature, I used QAG with GSL_INTEG_GAUSS15.  The 
resulting counts for each method were as follows:

                ***** cos integrand (smooth) *****
#dimensions   cubature  recursive-QAG  MC-plain  MC-miser  MC-vegas
      2           17          225        163840    20462     2890
      3           33          3375       327680    40936     5120
      4           57         50625       327680    81885     12500
      5           93        759375       327680   163788     10240
      6          447       11390625         -     163794     10935
      7          1205     170859375      655360   163801     12800
      8          5213         -             -     327616     25600
      9         31185         -             -     327620     25600
     10        160605         -         1310720   655268     25600
     11        603693         -         1310720   655274     20480

(A "-" indicates that the "converged" value had an error bigger than the 
specified tolerance, or that, in the case of QAG, it was taking longer 
than I cared to wait.)

***** gaussian integrand (smooth, peaked, boundary singularity) *****
#dimensions   cubature  recursive-QAG  MC-plain  MC-miser  MC-vegas
      2          255        5085           -       40923      6250
      3         1353       337635       5242880    163759    49130
      4         8721     22314225      10485760    327564      -
      5        58869         -         20971520    655194    168070
      6       327055         -         20971520    1310458   390625
      7      2079589         -            -        2621013   409600

[-- Attachment #2: Type: TEXT/x-csrc, Size: 26177 bytes --]

/*
 * Copyright (c) 2005 Steven G. Johnson
 *
 * Portions (see comments) based on HIntLib (also distributed under
 * the GNU GPL), copyright (c) 2002-2005 Rudolf Schuerer.
 *
 * Portions (see comments) based on GNU GSL (also distributed under
 * the GNU GPL), copyright (c) 1996-2000 Brian Gough.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>

/* Adaptive multidimensional integration on hypercubes (or, really,
   hyper-rectangles) using cubature rules.

   A cubature rule takes a function and a hypercube and evaluates
   the function at a small number of points, returning an estimate
   of the integral as well as an estimate of the error, and also
   a suggested dimension of the hypercube to subdivide.

   Given such a rule, the adaptive integration is simple:

   1) Evaluate the cubature rule on the hypercube(s).
      Stop if converged.

   2) Pick the hypercube with the largest estimated error,
      and divide it in two along the suggested dimension.

   3) Goto (1).

*/

typedef double (*integrand) (unsigned ndim, const double *x, void *);

/* Integrate the function f from xmin[dim] to xmax[dim], with at
   most maxEval function evaluations (0 for no limit),
   until the given absolute is achieved relative error.  val returns
   the integral, and estimated_error returns the estimate for the
   absolute error in val.  The return value of the function is 0
   on success and non-zero if there was an error. */
int adapt_integrate(integrand f, void *fdata,
		    unsigned dim, const double *xmin, const double *xmax, 
		    unsigned maxEval, 
		    double reqAbsError, double reqRelError, 
		    double *val, double *estimated_error);

/***************************************************************************/
/* Basic datatypes */

typedef struct {
     double val, err;
} esterr;

static double relError(esterr ee)
{
     return (ee.val == 0.0 ? HUGE_VAL : fabs(ee.err / ee.val));
}

typedef struct {
     unsigned dim;
     double *data;	/* length 2*dim = center followed by half-width */
     double vol;	/* cache volume = product of widths */
} hypercube;

static double compute_vol(const hypercube *h)
{
     unsigned i;
     double vol = 1;
     for (i = 0; i < h->dim; ++i)
	  vol *= 2 * h->data[i + h->dim];
     return vol;
}

static hypercube make_hypercube(unsigned dim, const double *center, const double *width)
{
     unsigned i;
     hypercube h;
     h.dim = dim;
     h.data = (double *) malloc(sizeof(double) * dim * 2);
     for (i = 0; i < dim; ++i) {
	  h.data[i] = center[i];
	  h.data[i + dim] = width[i];
     }
     h.vol = compute_vol(&h);
     return h;
}

static hypercube make_hypercube_range(unsigned dim, const double *xmin, const double *xmax)
{
     hypercube h = make_hypercube(dim, xmin, xmax);
     unsigned i;
     for (i = 0; i < dim; ++i) {
	  h.data[i] = 0.5 * (xmin[i] + xmax[i]);
	  h.data[i + dim] = 0.5 * (xmax[i] - xmin[i]);
     }
     h.vol = compute_vol(&h);
     return h;
}

static void destroy_hypercube(hypercube *h)
{
     free(h->data);
     h->dim = 0;
}

typedef struct {
     hypercube h;
     esterr ee;
     unsigned splitDim;
} region;

static region make_region(const hypercube *h)
{
     region R;
     R.h = make_hypercube(h->dim, h->data, h->data + h->dim);
     R.splitDim = 0;
     return R;
}

static void destroy_region(region *R)
{
     destroy_hypercube(&R->h);
}

static void cut_region(region *R, region *R2)
{
     unsigned d = R->splitDim, dim = R->h.dim;
     *R2 = *R;
     R->h.data[d + dim] *= 0.5;
     R->h.vol *= 0.5;
     R2->h = make_hypercube(dim, R->h.data, R->h.data + dim);
     R->h.data[d] -= R->h.data[d + dim];
     R2->h.data[d] += R->h.data[d + dim];
}

typedef struct rule_s {
     unsigned dim;              /* the dimensionality */
     unsigned num_points;       /* number of evaluation points */
     unsigned (*evalError)(struct rule_s *r, integrand f, void *fdata,
			   const hypercube *h, esterr *ee);
     void (*destroy)(struct rule_s *r);
} rule;

static void destroy_rule(rule *r)
{
     if (r->destroy) r->destroy(r);
     free(r);
}

static region eval_region(region R, integrand f, void *fdata, rule *r)
{
     R.splitDim = r->evalError(r, f, fdata, &R.h, &R.ee);
     return R;
}

/***************************************************************************/
/* Functions to loop over points in a hypercube. */

/* Based on orbitrule.cpp in HIntLib-0.0.10 */

/* ls0 returns the least-significant 0 bit of n (e.g. it returns
   0 if the LSB is 0, it returns 1 if the 2 LSBs are 01, etcetera). */

#if (defined(__GNUC__) || defined(__ICC)) && (defined(__i386__) || defined (__x86_64__))
/* use x86 bit-scan instruction, based on count_trailing_zeros()
   macro in GNU GMP's longlong.h. */
static unsigned ls0(unsigned n)
{
     unsigned count;
     n = ~n;
   __asm__("bsfl %1,%0": "=r"(count):"rm"(n));
     return count;
}
#else
static unsigned ls0(unsigned n)
{
     const unsigned bits[] = {
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
	  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8,
     };
     unsigned bit;
     bit = 0;
     while ((n & 0xff) == 0xff) {
	  n >>= 8;
	  bit += 8;
     }
     return bit + bits[n & 0xff];
}
#endif

/**
 *  Evaluate the integral on all 2^n points (+/-r,...+/-r)
 *
 *  A Gray-code ordering is used to minimize the number of coordinate updates
 *  in p.
 */
static double evalR_Rfs(integrand f, void *fdata, unsigned dim, double *p, const double *c, const double *r)
{
     double sum = 0;
     unsigned i;
     unsigned signs = 0;	/* 0/1 bit = +/- for corresponding element of r[] */

     /* We start with the point where r is ADDed in every coordinate (This implies signs=0) */
     for (i = 0; i < dim; ++i)
	  p[i] = c[i] + r[i];

     /* Loop through the points in gray-code ordering */
     for (i = 0;; ++i) {
	  unsigned mask;

	  sum += f(dim, p, fdata);

	  unsigned d = ls0(i);	/* which coordinate to flip */
	  if (d >= dim)
	       break;

	  /* flip the d-th bit and add/subtract r[d] */
	  mask = 1U << d;
	  signs ^= mask;
	  p[d] = (signs & mask) ? c[d] - r[d] : c[d] + r[d];
     }
     return sum;
}

static double evalRR0_0fs(integrand f, void *fdata, unsigned dim, double *p, const double *c, const double *r)
{
     unsigned i, j;
     double sum = 0;

     for (i = 0; i < dim - 1; ++i) {
	  p[i] = c[i] - r[i];
	  for (j = i + 1; j < dim; ++j) {
	       p[j] = c[j] - r[j];
	       sum += f(dim, p, fdata);
	       p[i] = c[i] + r[i];
	       sum += f(dim, p, fdata);
	       p[j] = c[j] + r[j];
	       sum += f(dim, p, fdata);
	       p[i] = c[i] - r[i];
	       sum += f(dim, p, fdata);

	       p[j] = c[j];	/* Done with j -> Restore p[j] */
	  }
	  p[i] = c[i];		/* Done with i -> Restore p[i] */
     }
     return sum;
}

static double evalR0_0fs4d(integrand f, void *fdata, unsigned dim, double *p, const double *c, double *sum0_, const double *r1, double *sum1_, const double *r2, double *sum2_)
{
     double maxdiff = 0;
     unsigned i, dimDiffMax = 0;
     double sum0, sum1 = 0, sum2 = 0; /* copies for aliasing, performance */

     double ratio = r1[0] / r2[0];

     ratio *= ratio;
     sum0 = f(dim, p, fdata);

     for (i = 0; i < dim; i++) {
	  double f1a, f1b, f2a, f2b, diff;

	  p[i] = c[i] - r1[i];
	  sum1 += (f1a = f(dim, p, fdata));
	  p[i] = c[i] + r1[i];
	  sum1 += (f1b = f(dim, p, fdata));
	  p[i] = c[i] - r2[i];
	  sum2 += (f2a = f(dim, p, fdata));
	  p[i] = c[i] + r2[i];
	  sum2 += (f2b = f(dim, p, fdata));
	  p[i] = c[i];

	  diff = fabs(f1a + f1b - 2 * sum0 - ratio * (f2a + f2b - 2 * sum0));

	  if (diff > maxdiff) {
	       maxdiff = diff;
	       dimDiffMax = i;
	  }
     }

     *sum0_ += sum0;
     *sum1_ += sum1;
     *sum2_ += sum2;

     return dimDiffMax;
}

#define num0_0(dim) (1U)
#define numR0_0fs(dim) (2 * (dim))
#define numRR0_0fs(dim) (2 * (dim) * (dim-1))
#define numR_Rfs(dim) (1U << (dim))

/***************************************************************************/
/* Based on rule75genzmalik.cpp in HIntLib-0.0.10: An embedded
   cubature rule of degree 7 (embedded rule degree 5) due to A. C. Genz
   and A. A. Malik.  See:

         A. C. Genz and A. A. Malik, "An imbedded [sic] family of fully
         symmetric numerical integration rules," SIAM
         J. Numer. Anal. 20 (3), 580-588 (1983).
*/

typedef struct {
     rule parent;

     /* temporary arrays of length dim */
     double *widthLambda, *widthLambda2, *p;

     /* dimension-dependent constants */
     double weight1, weight3, weight5;
     double weightE1, weightE3;
} rule75genzmalik;

#define real(x) ((double)(x))
#define to_int(n) ((int)(n))

static int isqr(int x)
{
     return x * x;
}

static void destroy_rule75genzmalik(rule *r_)
{
     rule75genzmalik *r = (rule75genzmalik *) r_;
     free(r->p);
}

static unsigned rule75genzmalik_evalError(rule *r_, integrand f, void *fdata, const hypercube *h, esterr *ee)
{
     /* lambda2 = sqrt(9/70), lambda4 = sqrt(9/10), lambda5 = sqrt(9/19) */
     const double lambda2 = 0.3585685828003180919906451539079374954541;
     const double lambda4 = 0.9486832980505137995996680633298155601160;
     const double lambda5 = 0.6882472016116852977216287342936235251269;
     const double weight2 = 980. / 6561.;
     const double weight4 = 200. / 19683.;
     const double weightE2 = 245. / 486.;
     const double weightE4 = 25. / 729.;

     rule75genzmalik *r = (rule75genzmalik *) r_;
     unsigned i, dimDiffMax, dim = r_->dim;
     double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, sum4, result, res5th;
     const double *center = h->data;
     const double *width = h->data + dim;

     for (i = 0; i < dim; ++i)
	  r->p[i] = center[i];

     for (i = 0; i < dim; ++i)
	  r->widthLambda2[i] = width[i] * lambda2;
     for (i = 0; i < dim; ++i)
	  r->widthLambda[i] = width[i] * lambda4;

     /* Evaluate function in the center, in f(lambda2,0,...,0) and
        f(lambda3=lambda4, 0,...,0).  Estimate dimension with largest error */
     dimDiffMax = evalR0_0fs4d(f, fdata, dim, r->p, center, &sum1, r->widthLambda2, &sum2, r->widthLambda, &sum3);

     /* Calculate sum4 for f(lambda4, lambda4, 0, ...,0) */
     sum4 = evalRR0_0fs(f, fdata, dim, r->p, center, r->widthLambda);

     /* Calculate sum5 for f(lambda5, lambda5, ..., lambda5) */
     for (i = 0; i < dim; ++i)
	  r->widthLambda[i] = width[i] * lambda5;
     double sum5 = evalR_Rfs(f, fdata, dim, r->p, center, r->widthLambda);

     /* Calculate fifth and seventh order results */

     result = h->vol * (r->weight1 * sum1 + weight2 * sum2 + r->weight3 * sum3 + weight4 * sum4 + r->weight5 * sum5);
     res5th = h->vol * (r->weightE1 * sum1 + weightE2 * sum2 + r->weightE3 * sum3 + weightE4 * sum4);

     ee->val = result;
     ee->err = fabs(res5th - result);

     return dimDiffMax;
}

static rule *make_rule75genzmalik(unsigned dim)
{
     rule75genzmalik *r;

     if (dim < 2) return 0; /* this rule does not support 1d integrals */

     r = (rule75genzmalik *) malloc(sizeof(rule75genzmalik));
     r->parent.dim = dim;

     r->weight1 = (real(12824 - 9120 * to_int(dim) + 400 * isqr(to_int(dim)))
		   / real(19683));
     r->weight3 = real(1820 - 400 * to_int(dim)) / real(19683);
     r->weight5 = real(6859) / real(19683) / real(1U << dim);
     r->weightE1 = (real(729 - 950 * to_int(dim) + 50 * isqr(to_int(dim)))
		    / real(729));
     r->weightE3 = real(265 - 100 * to_int(dim)) / real(1458);

     r->p = (double *) malloc(sizeof(double) * dim * 3);
     r->widthLambda = r->p + dim;
     r->widthLambda2 = r->p + 2 * dim;

     r->parent.num_points = num0_0(dim) + 2 * numR0_0fs(dim)
	  + numRR0_0fs(dim) + numR_Rfs(dim);

     r->parent.evalError = rule75genzmalik_evalError;
     r->parent.destroy = destroy_rule75genzmalik;

     return (rule *) r;
}

/***************************************************************************/
/* 1d 15-point Gaussian quadrature rule, based on qk15.c and qk.c in
   GNU GSL (which in turn is based on QUADPACK). */

static unsigned rule15gauss_evalError(rule *r, integrand f, void *fdata,
				      const hypercube *h, esterr *ee)
{
     /* Gauss quadrature weights and kronrod quadrature abscissae and
	weights as evaluated with 80 decimal digit arithmetic by
	L. W. Fullerton, Bell Labs, Nov. 1981. */
     const unsigned n = 8;
     const double xgk[8] = {  /* abscissae of the 15-point kronrod rule */
	  0.991455371120812639206854697526329,
	  0.949107912342758524526189684047851,
	  0.864864423359769072789712788640926,
	  0.741531185599394439863864773280788,
	  0.586087235467691130294144838258730,
	  0.405845151377397166906606412076961,
	  0.207784955007898467600689403773245,
	  0.000000000000000000000000000000000
	  /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule. 
	     xgk[0], xgk[2], ... to optimally extend the 7-point gauss rule */
     };
     static const double wg[4] = {  /* weights of the 7-point gauss rule */
	  0.129484966168869693270611432679082,
	  0.279705391489276667901467771423780,
	  0.381830050505118944950369775488975,
	  0.417959183673469387755102040816327
     };
     static const double wgk[8] = { /* weights of the 15-point kronrod rule */
	  0.022935322010529224963732008058970,
	  0.063092092629978553290700663189204,
	  0.104790010322250183839876322541518,
	  0.140653259715525918745189590510238,
	  0.169004726639267902826583426598550,
	  0.190350578064785409913256402421014,
	  0.204432940075298892414161999234649,
	  0.209482141084727828012999174891714
     };

     const double center = h->data[0];
     const double width = h->data[1];
     double fv1[n - 1], fv2[n - 1];
     const double f_center = f(1, &center, fdata);
     double result_gauss = f_center * wg[n/2 - 1];
     double result_kronrod = f_center * wgk[n - 1];
     double result_abs = fabs(result_kronrod);
     double result_asc, mean, err;
     unsigned j;

     for (j = 0; j < (n - 1) / 2; ++j) {
	  int j2 = 2*j + 1;
	  double x, f1, f2, fsum, w = width * xgk[j2];
	  x = center - w; fv1[j2] = f1 = f(1, &x, fdata);
	  x = center + w; fv2[j2] = f2 = f(1, &x, fdata);
	  fsum = f1 + f2;
	  result_gauss += wg[j] * fsum;
	  result_kronrod += wgk[j2] * fsum;
	  result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
     }

     for (j = 0; j < n/2; ++j) {
	  int j2 = 2*j;
	  double x, f1, f2, w = width * xgk[j2];
	  x = center - w; fv1[j2] = f1 = f(1, &x, fdata);
          x = center + w; fv2[j2] = f2 = f(1, &x, fdata);
          result_kronrod += wgk[j2] * (f1 + f2);
          result_abs += wgk[j2] * (fabs(f1) + fabs(f2));

     }

     ee->val = result_kronrod * width;

     /* compute error estimate: */
     mean = result_kronrod * 0.5;
     result_asc = wgk[n - 1] * fabs(f_center - mean);
     for (j = 0; j < n - 1; ++j)
	  result_asc += wgk[j] * (fabs(fv1[j]-mean) + fabs(fv2[j]-mean));
     err = (result_kronrod - result_gauss) * width;
     result_abs *= width;
     result_asc *= width;
     if (result_asc != 0 && err != 0) {
	  double scale = pow((200 * err / result_asc), 1.5);
	  if (scale < 1)
	       err = result_asc * scale;
	  else
	       err = result_asc;
     }
     if (result_abs > DBL_MIN / (50 * DBL_EPSILON)) {
	  double min_err = 50 * DBL_EPSILON * result_abs;
	  if (min_err > err)
	       err = min_err;
     }
     ee->err = err;
     
     return 0; /* no choice but to divide 0th dimension */
}

static rule *make_rule15gauss(unsigned dim)
{
     rule *r;
     if (dim != 1) return 0; /* this rule is only for 1d integrals */
     r = (rule *) malloc(sizeof(rule));
     r->dim = dim;
     r->num_points = 15;
     r->evalError = rule15gauss_evalError;
     r->destroy = 0;
     return r;
}

/***************************************************************************/
/* binary heap implementation (ala _Introduction to Algorithms_ by
   Cormen, Leiserson, and Rivest), for use as a priority queue of
   regions to integrate. */

typedef region heap_item;
#define KEY(hi) ((hi).ee.err)

typedef struct {
     unsigned n, nalloc;
     heap_item *items;
     esterr ee;
} heap;

static void heap_resize(heap *h, unsigned nalloc)
{
     h->nalloc = nalloc;
     h->items = (heap_item *) realloc(h->items, sizeof(heap_item) * nalloc);
}

static heap heap_alloc(unsigned nalloc)
{
     heap h;
     h.n = 0;
     h.nalloc = 0;
     h.items = 0;
     h.ee.val = h.ee.err = 0;
     heap_resize(&h, nalloc);
     return h;
}

/* note that heap_free does not deallocate anything referenced by the items */
static void heap_free(heap *h)
{
     h->n = 0;
     heap_resize(h, 0);
}

static void heap_push(heap *h, heap_item hi)
{
     int insert;

     h->ee.val += hi.ee.val;
     h->ee.err += hi.ee.err;
     insert = h->n;
     if (++(h->n) > h->nalloc)
	  heap_resize(h, h->n * 2);

     while (insert) {
	  int parent = (insert - 1) / 2;
	  if (KEY(hi) <= KEY(h->items[parent]))
	       break;
	  h->items[insert] = h->items[parent];
	  insert = parent;
     }
     h->items[insert] = hi;
}

static heap_item heap_pop(heap *h)
{
     heap_item ret;
     int i, n, child;

     if (!(h->n)) {
	  fprintf(stderr, "attempted to pop an empty heap\n");
	  exit(EXIT_FAILURE);
     }

     ret = h->items[0];
     h->items[i = 0] = h->items[n = --(h->n)];
     while ((child = i * 2 + 1) < n) {
	  int largest;
	  heap_item swap;

	  if (KEY(h->items[child]) <= KEY(h->items[i]))
	       largest = i;
	  else
	       largest = child;
	  if (++child < n && KEY(h->items[largest]) < KEY(h->items[child]))
	       largest = child;
	  if (largest == i)
	       break;
	  swap = h->items[i];
	  h->items[i] = h->items[largest];
	  h->items[i = largest] = swap;
     }

     h->ee.val -= ret.ee.val;
     h->ee.err -= ret.ee.err;
     return ret;
}

/***************************************************************************/

/* adaptive integration, analogous to adaptintegrator.cpp in HIntLib */

static int ruleadapt_integrate(rule *r, integrand f, void *fdata, const hypercube *h, unsigned maxEval, double reqAbsError, double reqRelError, esterr *ee)
{
     unsigned initialRegions;	/* number of initial regions (non-adaptive) */
     unsigned minIter;		/* minimum number of adaptive subdivisions */
     unsigned maxIter;		/* maximum number of adaptive subdivisions */
     unsigned initialPoints;
     heap regions;
     unsigned i;
     int status = -1; /* = ERROR */

     initialRegions = 1; /* or: use some percentage of maxEval/r->num_points */
     initialPoints = initialRegions * r->num_points;
     minIter = 0;
     if (maxEval) {
	  if (initialPoints > maxEval) {
	       initialRegions = maxEval / r->num_points;
	       initialPoints = initialRegions * r->num_points;
	  }
	  maxEval -= initialPoints;
	  maxIter = maxEval / (2 * r->num_points);
     } else
	  maxIter = UINT_MAX;

     if (initialRegions == 0)
	  return status;	/* ERROR */

     regions = heap_alloc(initialRegions);

     heap_push(&regions, eval_region(make_region(h), f, fdata, r));
     /* or: 
	if (initialRegions != 1)
	   partition(h, initialRegions, EQUIDISTANT, &regions, f,fdata, r); */

     for (i = 0; i < maxIter; ++i) {
	  region R, R2;
	  if (i >= minIter && (regions.ee.err <= reqAbsError 
			       || relError(regions.ee) <= reqRelError)) {
	       status = 0; /* converged! */
	       break;
	  }
	  R = heap_pop(&regions); /* get worst region */
	  cut_region(&R, &R2);
	  heap_push(&regions, eval_region(R, f, fdata, r));
	  heap_push(&regions, eval_region(R2, f, fdata, r));
     }

     ee->val = ee->err = 0;  /* re-sum integral and errors */
     for (i = 0; i < regions.n; ++i) {
	  ee->val += regions.items[i].ee.val;
	  ee->err += regions.items[i].ee.err;
	  destroy_region(&regions.items[i]);
     }
     printf("regions.nalloc = %d\n", regions.nalloc);
     heap_free(&regions);

     return status;
}

int adapt_integrate(integrand f, void *fdata, 
		    unsigned dim, const double *xmin, const double *xmax, 
		    unsigned maxEval, double reqAbsError, double reqRelError, 
		    double *val, double *estimated_error)
{
     rule *r;
     hypercube h;
     esterr ee;
     
     if (dim == 0) { /* trivial integration */
	  ee.val = f(0, xmin, fdata);
	  ee.err = 0;
	  return 0;
     }
     r = dim == 1 ? make_rule15gauss(dim) : make_rule75genzmalik(dim);
     h = make_hypercube_range(dim, xmin, xmax);
     int status = ruleadapt_integrate(r, f, fdata, &h,
				      maxEval, reqAbsError, reqRelError,
				      &ee);
     *val = ee.val;
     *estimated_error = ee.err;
     destroy_hypercube(&h);
     destroy_rule(r);
     return status;
}

/***************************************************************************/

/* Compile with -DTEST_INTEGRATOR to generate this little test program.
   
   Usage: ./integrator <dim> <tol> <integrand> <maxeval>

   where <dim> = # dimensions, <tol> = relative tolerance,
   <integrand> is either 0/1/2 for the three test integrands (see below),
   and <maxeval> is the maximum # function evaluations (0 for none).
*/
   
#ifdef TEST_INTEGRATOR

int count = 0;
int which_integrand = 0;
const double radius = 0.50124145262344534123412; /* random */

double f_test(unsigned dim, const double *x, void *data)
{
     double val;
     unsigned i;
     ++count;
     switch (which_integrand) {
	 case 0: /* discontinuous objective: volume of hypersphere */
	      val = 0;
	      for (i = 0; i < dim; ++i)
		   val += x[i] * x[i];
	      val = val < radius * radius;
	      break;
	 case 1: /* simple smooth (separable) objective: prod. cos(x[i]). */
	      val = 1;
	      for (i = 0; i < dim; ++i)
		   val *= cos(x[i]);
	      break;
	 case 2: { /* integral of exp(-x^2), rescaled to (0,infinity) limits */
	      double scale = 1.0;
	      val = 0;
	      for (i = 0; i < dim; ++i) {
		   double z = (1 - x[i]) / x[i];
		   val += z * z;
		   scale *= M_2_SQRTPI / (x[i] * x[i]);
	      }
	      val = exp(-val) * scale;
	      break;
	 }
		   
	 default:
	      fprintf(stderr, "unknown integrand %d\n", which_integrand);
	      exit(EXIT_FAILURE);
     }
     /* if (count < 100) printf("%d: f(%g, ...) = %g\n", count, x[0], val); */
     return val;
}

/* surface area of n-dimensional unit hypersphere */
static double S(unsigned n)
{
     double val;
     int fact = 1;
     if (n % 2 == 0) { /* n even */
	  val = 2 * pow(M_PI, n * 0.5);
	  n = n / 2;
	  while (n > 1) fact *= (n -= 1);
	  val /= fact;
     }
     else { /* n odd */
	  val = (1 << (n/2 + 1)) * pow(M_PI, n/2);
	  while (n > 2) fact *= (n -= 2);
	  val /= fact;
     }
     return val;
}

static double exact_integral(unsigned dim, const double *xmax) {
     unsigned i;
     double val;
     switch(which_integrand) {
	 case 0:
	      val = dim == 0 ? 1 : S(dim) * pow(radius * 0.5, dim) / dim;
	      break;
	 case 1:
	      val = 1;
	      for (i = 0; i < dim; ++i)
		   val *= sin(xmax[i]);
	      break;
	 case 2:
	      val = 1;
	      break;
	 default:
	      fprintf(stderr, "unknown integrand %d\n", which_integrand);
              exit(EXIT_FAILURE);
     }
     return val;
}

int main(int argc, char **argv)
{
     double *xmin, *xmax;
     double tol, val, err;
     unsigned i, dim, maxEval;

     dim = argc > 1 ? atoi(argv[1]) : 2;
     tol = argc > 2 ? atof(argv[2]) : 1e-2;
     which_integrand = argc > 3 ? atoi(argv[3]) : 1;
     maxEval = argc > 4 ? atoi(argv[4]) : 0;

     xmin = (double *) malloc(dim * sizeof(double));
     xmax = (double *) malloc(dim * sizeof(double));
     for (i = 0; i < dim; ++i) {
	  xmin[i] = 0;
	  xmax[i] = 1 + (which_integrand == 2 ? 0 : 0.4 * sin(i));
     }

     printf("%u-dim integral, tolerance = %g, integrand = %d\n",
	    dim, tol, which_integrand);
     adapt_integrate(f_test, 0, dim, xmin, xmax, maxEval, 0, tol, &val, &err);
     printf("integration val = %g, est. err = %g, true err = %g\n", 
	    val, err, fabs(val - exact_integral(dim, xmax)));
     printf("#evals = %d\n", count);

     free(xmax);
     free(xmin);

     return 0;
}

#endif

Re: contribution: code for multi-dimensional adaptive integration

[Brian Gough]

Steven G. Johnson writes:
 > For some time now, I've felt that GSL has cried out for a multidimensional 
 > cubature routine, similar to its 1d adaptive quadrature routines.  Now, 
 > I've written one (~700 lines in C) based on the well-known Genz-Malik 
 > embedded cubature rule (with help from a GPL'ed library called HIntLib, 
 > see below), and I hope that you're interested to have it included in GSL. 
 > I've attached the unmodified code; I'm willing to do the work to modify it 
 > for your coding style, etcetera.

Thanks.. that does look useful.  You're right that cubature is an
obvious gap in GSL, and I would like to see that filled.

The first step would be make a package that builds outside of gsl,
there are some instructions at http://sources.redhat.com/gsl/devel.html
(notes on the coding conventions can also be found in the design document 
linked there)

-- 
Brian Gough

p.s. Sorry for the delay in replying.