https://github.com/stevengj/cubature
Tip revision: d1955d864f53ecb2b95a802c4b9cceeeccae7b71 authored by Steven G. Johnson on 31 May 2023, 14:27:06 UTC
added citation info
added citation info
Tip revision: d1955d8
pcubature.c
/* Adaptive multidimensional integration of a vector of integrands.
*
* Copyright (c) 2005-2013 Steven G. Johnson
*
* 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
*
*/
/* p-adaptive cubature (adaptive by increasing the degree of the
cubature rule rather than subdividing the domain), using products
of Clenshaw-Curtis rules. This algorithm may be superior to
Genz-Malik for smooth integrands lacking strongly-localized
features, in moderate dimensions. */
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "cubature.h"
/* error return codes */
#define SUCCESS 0
#define FAILURE 1
/* pre-generated Clenshaw-Curtis rules and weights */
#include "clencurt.h"
/* no point in supporting very high dimensional integrals here */
#define MAXDIM (20U)
/***************************************************************************/
/* For adaptive cubature, thanks to the nesting of the C-C rules, we
can re-use the values from coarser grids for finer grids, and the
coarser grids are also used for error estimation.
A grid is determined by an m[dim] array, where m[i] denotes
2^(m[i]+1)+1 points in the i-th dimension.
*/
/* cache of the values for the m[dim] grid. If mi < dim, then we only
store the values corresponding to the difference between the m grid
and the grid with m[mi] -> m[mi]-1. (m[mi]-1 == -1 corresponds to
the trivial grid of one point in the center.) */
typedef struct cacheval_s {
unsigned m[MAXDIM];
unsigned mi;
double *val;
} cacheval;
/* array of ncache cachevals c[i] */
typedef struct valcache_s {
size_t ncache;
cacheval *c;
} valcache;
static void free_cachevals(valcache *v)
{
if (!v) return;
if (v->c) {
size_t i;
for (i = 0; i < v->ncache; ++i)
free(v->c[i].val);
free(v->c);
v->c = NULL;
}
v->ncache = 0;
}
/***************************************************************************/
/* recursive loop over all cubature points for the given (m,mi) cache entry:
add each point to the buffer buf, evaluating all at once whenever the
buffer is full or when we are done */
static int compute_cacheval(const unsigned *m, unsigned mi,
double *val, size_t *vali,
unsigned fdim, integrand_v f, void *fdata,
unsigned dim, unsigned id, double *p,
const double *xmin, const double *xmax,
double *buf, size_t nbuf, size_t *ibuf)
{
if (id == dim) { /* add point to buffer of points */
memcpy(buf + (*ibuf)++ * dim, p, sizeof(double) * dim);
if (*ibuf == nbuf) { /* flush buffer */
if (f(dim, nbuf, buf, fdata, fdim, val + *vali))
return FAILURE;
*vali += *ibuf * fdim;
*ibuf = 0;
}
}
else {
double c = (xmin[id] + xmax[id]) * 0.5;
double r = (xmax[id] - xmin[id]) * 0.5;
const double *x = clencurt_x
+ ((id == mi) ? (m[id] ? (1 << (m[id] - 1)) : 0) : 0);
unsigned i, nx = (id == mi ? (m[id] ? (1 << (m[id] - 1)) : 1)
: (1 << (m[id])));
if (id != mi) {
p[id] = c;
if (compute_cacheval(m, mi, val, vali, fdim, f, fdata,
dim, id + 1, p,
xmin, xmax, buf, nbuf, ibuf))
return FAILURE;
}
for (i = 0; i < nx; ++i) {
p[id] = c + r * x[i];
if (compute_cacheval(m, mi, val, vali, fdim, f, fdata,
dim, id + 1, p,
xmin, xmax, buf, nbuf, ibuf))
return FAILURE;
p[id] = c - r * x[i];
if (compute_cacheval(m, mi, val, vali, fdim, f, fdata,
dim, id + 1, p,
xmin, xmax, buf, nbuf, ibuf))
return FAILURE;
}
}
return SUCCESS;
}
static size_t num_cacheval(const unsigned *m, unsigned mi, unsigned dim)
{
unsigned i;
size_t nval = 1;
for (i = 0; i < dim; ++i) {
if (i == mi)
nval *= m[i] == 0 ? 2 : (1 << (m[i]));
else
nval *= (1 << (m[i] + 1)) + 1;
}
return nval;
}
static int add_cacheval(valcache *vc,
const unsigned *m, unsigned mi,
unsigned fdim, integrand_v f, void *fdata,
unsigned dim, const double *xmin, const double *xmax,
double *buf, size_t nbuf)
{
size_t ic = vc->ncache;
size_t nval, vali = 0, ibuf = 0;
double p[MAXDIM];
vc->c = (cacheval *) realloc(vc->c, sizeof(cacheval) * ++(vc->ncache));
if (!vc->c) return -1;
vc->c[ic].mi = mi;
memcpy(vc->c[ic].m, m, sizeof(unsigned) * dim);
nval = fdim * num_cacheval(m, mi, dim);
vc->c[ic].val = (double *) malloc(sizeof(double) * nval);
if (!vc->c[ic].val) return FAILURE;
if (compute_cacheval(m, mi, vc->c[ic].val, &vali,
fdim, f, fdata,
dim, 0, p, xmin, xmax,
buf, nbuf, &ibuf))
return FAILURE;
if (ibuf > 0) /* flush remaining buffer */
return f(dim, ibuf, buf, fdata, fdim, vc->c[ic].val + vali);
return SUCCESS;
}
/***************************************************************************/
/* recursive loop to evaluate the integral contribution from the cache
entry c, accumulating in val, for the given m[] except with m[md]
-> m[md] - 1 if md < dim, using the cached values (cm,cmi,cval). id is the
current loop dimension (from 0 to dim-1). */
static unsigned eval(const unsigned *cm, unsigned cmi, double *cval,
const unsigned *m, unsigned md,
unsigned fdim, unsigned dim, unsigned id,
double weight, double *val)
{
size_t voff = 0; /* amount caller should offset cval array afterwards */
if (id == dim) {
unsigned i;
for (i = 0; i < fdim; ++i) val[i] += cval[i] * weight;
voff = fdim;
}
else if (m[id] == 0 && id == md) /* using trivial rule for this dim */ {
voff = eval(cm, cmi, cval, m, md, fdim, dim, id+1, weight*2, val);
voff += fdim * (1 << cm[id]) * 2
* num_cacheval(cm + id+1, cmi - (id+1), dim - (id+1));
}
else {
unsigned i;
unsigned mid = m[id] - (id == md); /* order of C-C rule */
const double *w = clencurt_w + mid + (1 << mid) - 1
+ (id == cmi ? (cm[id] ? 1 + (1 << (cm[id]-1)) : 1) : 0);
unsigned cnx = (id == cmi ? (cm[id] ? (1 << (cm[id]-1)) : 1)
: (1 << (cm[id])));
unsigned nx = cm[id] <= mid ? cnx : (1 << mid);
if (id != cmi) {
voff = eval(cm, cmi, cval, m, md, fdim, dim, id + 1,
weight * w[0], val);
++w;
}
for (i = 0; i < nx; ++i) {
voff += eval(cm, cmi, cval + voff, m, md, fdim, dim, id + 1,
weight * w[i], val);
voff += eval(cm, cmi, cval + voff, m, md, fdim, dim, id + 1,
weight * w[i], val);
}
voff += (cnx - nx) * fdim * 2
* num_cacheval(cm + id+1, cmi - (id+1), dim - (id+1));
}
return voff;
}
/* loop over all cache entries that contribute to the integral,
(with m[md] decremented by 1) */
static void evals(valcache vc, const unsigned *m, unsigned md,
unsigned fdim, unsigned dim,
double V, double *val)
{
size_t i;
memset(val, 0, sizeof(double) * fdim);
for (i = 0; i < vc.ncache; ++i) {
if (vc.c[i].mi >= dim ||
vc.c[i].m[vc.c[i].mi] + (vc.c[i].mi == md) <= m[vc.c[i].mi])
eval(vc.c[i].m, vc.c[i].mi, vc.c[i].val,
m, md, fdim, dim, 0, V, val);
}
}
/* evaluate the integrals for the given m[] using the cached values in vc,
storing the integrals in val[], the error estimate in err[], and the
dimension to subdivide next (the largest error contribution) in *mi */
static void eval_integral(valcache vc, const unsigned *m,
unsigned fdim, unsigned dim, double V,
unsigned *mi, double *val, double *err, double *val1)
{
double maxerr = 0;
unsigned i, j;
evals(vc, m, dim, fdim, dim, V, val);
/* error estimates along each dimension by comparing val with
lower-order rule in that dimension; overall (conservative)
error estimate from maximum error of lower-order rules. */
memset(err, 0, sizeof(double) * fdim);
*mi = 0;
for (i = 0; i < dim; ++i) {
double emax = 0;
evals(vc, m, i, fdim, dim, V, val1);
for (j = 0; j < fdim; ++j) {
double e = fabs(val[j] - val1[j]);
if (e > emax) emax = e;
if (e > err[j]) err[j] = e;
}
if (emax > maxerr) {
maxerr = emax;
*mi = i;
}
}
/* printf("eval: %g +/- %g (dim %u)\n", val[0], err[0], *mi); */
}
/***************************************************************************/
static int converged(unsigned fdim, const double *vals, const double *errs,
double reqAbsError, double reqRelError, error_norm norm)
#define ERR(j) errs[j]
#define VAL(j) vals[j]
#include "converged.h"
/***************************************************************************/
/* Vectorized version with user-supplied buffer to store points and values.
The buffer *buf should be of length *nbuf * dim on entry (these parameters
are changed upon return to the final buffer and length that was used).
The buffer length will be kept <= max(max_nbuf, 1) * dim.
Also allows the caller to specify an array m[dim] of starting degrees
for the rule, which upon return will hold the final degrees. The
number of points in each dimension i is 2^(m[i]+1) + 1. */
int pcubature_v_buf(unsigned fdim, integrand_v f, void *fdata,
unsigned dim, const double *xmin, const double *xmax,
size_t maxEval,
double reqAbsError, double reqRelError,
error_norm norm,
unsigned *m,
double **buf, size_t *nbuf, size_t max_nbuf,
double *val, double *err)
{
int ret = FAILURE;
double V = 1;
size_t numEval = 0, new_nbuf;
unsigned i;
valcache vc = {0, NULL};
double *val1 = NULL;
if (fdim <= 1) norm = ERROR_INDIVIDUAL; /* norm is irrelevant */
if (norm < 0 || norm > ERROR_LINF) return FAILURE; /* invalid norm */
if (fdim == 0) return SUCCESS; /* nothing to do */
if (dim > MAXDIM) return FAILURE; /* unsupported */
if (dim == 0) { /* trivial case */
if (f(0, 1, xmin, fdata, fdim, val)) return FAILURE;
for (i = 0; i < fdim; ++i) err[i] = 0;
return SUCCESS;
}
for (i = 0; i < fdim; ++i) {
val[i] = 0;
err[i] = HUGE_VAL;
}
for (i = 0; i < dim; ++i)
V *= (xmax[i] - xmin[i]) * 0.5; /* scale factor for C-C volume */
new_nbuf = num_cacheval(m, dim, dim);
if (max_nbuf < 1) max_nbuf = 1;
if (new_nbuf > max_nbuf) new_nbuf = max_nbuf;
if (*nbuf < new_nbuf) {
free(*buf);
*buf = (double *) malloc(sizeof(double)
* (*nbuf = new_nbuf) * dim);
if (!*buf) goto done;
}
/* start by evaluating the m=0 cubature rule */
if (add_cacheval(&vc, m, dim, fdim, f, fdata, dim, xmin, xmax,
*buf, *nbuf) != SUCCESS)
goto done;
val1 = (double *) malloc(sizeof(double) * fdim);
while (1) {
unsigned mi;
eval_integral(vc, m, fdim, dim, V, &mi, val, err, val1);
if (converged(fdim, val, err, reqAbsError, reqRelError, norm)
|| (numEval > maxEval && maxEval)) {
ret = SUCCESS;
goto done;
}
m[mi] += 1;
if (m[mi] > clencurt_M) goto done; /* FAILURE */
new_nbuf = num_cacheval(m, mi, dim);
if (new_nbuf > *nbuf && *nbuf < max_nbuf) {
*nbuf = new_nbuf;
if (*nbuf > max_nbuf) *nbuf = max_nbuf;
free(*buf);
*buf = (double *) malloc(sizeof(double) * *nbuf * dim);
if (!*buf) goto done; /* FAILURE */
}
if (add_cacheval(&vc, m, mi, fdim, f, fdata,
dim, xmin, xmax, *buf, *nbuf) != SUCCESS)
goto done; /* FAILURE */
numEval += new_nbuf;
}
done:
free(val1);
free_cachevals(&vc);
return ret;
}
/***************************************************************************/
#define DEFAULT_MAX_NBUF (1U << 20)
int pcubature_v(unsigned fdim, integrand_v f, void *fdata,
unsigned dim, const double *xmin, const double *xmax,
size_t maxEval, double reqAbsError, double reqRelError,
error_norm norm,
double *val, double *err)
{
int ret;
size_t nbuf = 0;
unsigned m[MAXDIM];
double *buf = NULL;
memset(m, 0, sizeof(unsigned) * dim);
ret = pcubature_v_buf(fdim, f, fdata, dim, xmin, xmax,
maxEval, reqAbsError, reqRelError, norm,
m, &buf, &nbuf, DEFAULT_MAX_NBUF, val, err);
free(buf);
return ret;
}
#include "vwrapper.h"
int pcubature(unsigned fdim, integrand f, void *fdata,
unsigned dim, const double *xmin, const double *xmax,
size_t maxEval, double reqAbsError, double reqRelError,
error_norm norm,
double *val, double *err)
{
int ret;
size_t nbuf = 0;
unsigned m[MAXDIM];
double *buf = NULL;
fv_data d;
d.f = f; d.fdata = fdata;
memset(m, 0, sizeof(unsigned) * dim);
ret = pcubature_v_buf(
fdim, fv, &d, dim, xmin, xmax,
maxEval, reqAbsError, reqRelError, norm,
m, &buf, &nbuf, 16 /* max_nbuf > 0 to amortize function overhead */,
val, err);
free(buf);
return ret;
}