This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
contribution: code for multi-dimensional adaptive integration
- From: "Steven G. Johnson" <stevenj at fftw dot org>
- To: gsl-discuss at sources dot redhat dot com, help-gsl at gnu dot org
- Date: Thu, 9 Jun 2005 21:13:55 -0400 (EDT)
- Subject: contribution: code for multi-dimensional adaptive integration
- Reply-to: stevenj at math dot mit dot edu
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
/*
* 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, ¢er, 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(®ions, eval_region(make_region(h), f, fdata, r));
/* or:
if (initialRegions != 1)
partition(h, initialRegions, EQUIDISTANT, ®ions, 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(®ions); /* get worst region */
cut_region(&R, &R2);
heap_push(®ions, eval_region(R, f, fdata, r));
heap_push(®ions, 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(®ions.items[i]);
}
printf("regions.nalloc = %d\n", regions.nalloc);
heap_free(®ions);
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