« back to all changes in this revision
Viewing changes to src/integrator.c
- browse files at revision 4
- compare with another revision
- download diff
- download tarball
- view history from revision 4
- Committer: Bazaar Package Importer
- Author(s): Josselin Mouette
- Date: 2006-05-01 20:25:01 UTC
- mfrom: (1.1.2 upstream)
- Revision ID: james.westby@ubuntu.com-20060501202501-lytbmb3oevyoqzxi
Tags: 3.0.1-1
* New upstream release (closes: #361676).
* Major rework of the debian/ directory. Switch to cdbs.
* Migrate Scheme files to a versioned location to allow several
versions to be installed at once.
* Write a Makefile to put with the example.
* Update copyright, the library is now GPL.
* Use gfortran for the F77 wrappers.
* Standards-version is 3.7.0.
* Major rework of the debian/ directory. Switch to cdbs.
* Migrate Scheme files to a versioned location to allow several
versions to be installed at once.
* Write a Makefile to put with the example.
* Update copyright, the library is now GPL.
* Use gfortran for the F77 wrappers.
* Standards-version is 3.7.0.
added
removed
src/integrator.c
1
/*
2
* Copyright (c) 2005 Steven G. Johnson
3
*
4
* Portions (see comments) based on HIntLib (also distributed under
5
* the GNU GPL, v2 or later), copyright (c) 2002-2005 Rudolf Schuerer.
6
* (http://www.cosy.sbg.ac.at/~rschuer/hintlib/)
7
*
8
* Portions (see comments) based on GNU GSL (also distributed under
9
* the GNU GPL, v2 or later), copyright (c) 1996-2000 Brian Gough.
10
* (http://www.gnu.org/software/gsl/)
11
*
12
* This program is free software; you can redistribute it and/or modify
13
* it under the terms of the GNU General Public License as published by
14
* the Free Software Foundation; either version 2 of the License, or
15
* (at your option) any later version.
16
*
17
* This program is distributed in the hope that it will be useful,
18
* but WITHOUT ANY WARRANTY; without even the implied warranty of
19
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20
* GNU General Public License for more details.
21
*
22
* You should have received a copy of the GNU General Public License
23
* along with this program; if not, write to the Free Software
24
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25
*
26
*/
27
28
#include <stdio.h>
29
#include <stdlib.h>
30
#include <math.h>
31
#include <limits.h>
32
#include <float.h>
33
34
/* Adaptive multidimensional integration on hypercubes (or, really,
35
hyper-rectangles) using cubature rules.
36
37
A cubature rule takes a function and a hypercube and evaluates
38
the function at a small number of points, returning an estimate
39
of the integral as well as an estimate of the error, and also
40
a suggested dimension of the hypercube to subdivide.
41
42
Given such a rule, the adaptive integration is simple:
43
44
1) Evaluate the cubature rule on the hypercube(s).
45
Stop if converged.
46
47
2) Pick the hypercube with the largest estimated error,
48
and divide it in two along the suggested dimension.
49
50
3) Goto (1).
51
52
*/
53
54
typedef double (*integrand) (unsigned ndim, const double *x, void *);
55
56
/* Integrate the function f from xmin[dim] to xmax[dim], with at most
57
maxEval function evaluations (0 for no limit), until the given
58
absolute or relative error is achieved. val returns the integral,
59
and err returns the estimate for the absolute error in val. The
60
return value of the function is 0 on success and non-zero if there
61
was an error. */
62
static int adapt_integrate(integrand f, void *fdata,
63
unsigned dim, const double *xmin, const double *xmax,
64
unsigned maxEval,
65
double reqAbsError, double reqRelError,
66
double *val, double *err);
67
68
/***************************************************************************/
69
/* Basic datatypes */
70
71
typedef struct {
72
double val, err;
73
} esterr;
74
75
static double relError(esterr ee)
76
{
77
return (ee.val == 0.0 ? HUGE_VAL : fabs(ee.err / ee.val));
78
}
79
80
typedef struct {
81
unsigned dim;
82
double *data; /* length 2*dim = center followed by half-widths */
83
double vol; /* cache volume = product of widths */
84
} hypercube;
85
86
static double compute_vol(const hypercube *h)
87
{
88
unsigned i;
89
double vol = 1;
90
for (i = 0; i < h->dim; ++i)
91
vol *= 2 * h->data[i + h->dim];
92
return vol;
93
}
94
95
static hypercube make_hypercube(unsigned dim, const double *center, const double *halfwidth)
96
{
97
unsigned i;
98
hypercube h;
99
h.dim = dim;
100
h.data = (double *) malloc(sizeof(double) * dim * 2);
101
for (i = 0; i < dim; ++i) {
102
h.data[i] = center[i];
103
h.data[i + dim] = halfwidth[i];
104
}
105
h.vol = compute_vol(&h);
106
return h;
107
}
108
109
static hypercube make_hypercube_range(unsigned dim, const double *xmin, const double *xmax)
110
{
111
hypercube h = make_hypercube(dim, xmin, xmax);
112
unsigned i;
113
for (i = 0; i < dim; ++i) {
114
h.data[i] = 0.5 * (xmin[i] + xmax[i]);
115
h.data[i + dim] = 0.5 * (xmax[i] - xmin[i]);
116
}
117
h.vol = compute_vol(&h);
118
return h;
119
}
120
121
static void destroy_hypercube(hypercube *h)
122
{
123
free(h->data);
124
h->dim = 0;
125
}
126
127
typedef struct {
128
hypercube h;
129
esterr ee;
130
unsigned splitDim;
131
} region;
132
133
static region make_region(const hypercube *h)
134
{
135
region R;
136
R.h = make_hypercube(h->dim, h->data, h->data + h->dim);
137
R.splitDim = 0;
138
return R;
139
}
140
141
static void destroy_region(region *R)
142
{
143
destroy_hypercube(&R->h);
144
}
145
146
static void cut_region(region *R, region *R2)
147
{
148
unsigned d = R->splitDim, dim = R->h.dim;
149
*R2 = *R;
150
R->h.data[d + dim] *= 0.5;
151
R->h.vol *= 0.5;
152
R2->h = make_hypercube(dim, R->h.data, R->h.data + dim);
153
R->h.data[d] -= R->h.data[d + dim];
154
R2->h.data[d] += R->h.data[d + dim];
155
}
156
157
typedef struct rule_s {
158
unsigned dim; /* the dimensionality */
159
unsigned num_points; /* number of evaluation points */
160
unsigned (*evalError)(struct rule_s *r, integrand f, void *fdata,
161
const hypercube *h, esterr *ee);
162
void (*destroy)(struct rule_s *r);
163
} rule;
164
165
static void destroy_rule(rule *r)
166
{
167
if (r->destroy) r->destroy(r);
168
free(r);
169
}
170
171
static region eval_region(region R, integrand f, void *fdata, rule *r)
172
{
173
R.splitDim = r->evalError(r, f, fdata, &R.h, &R.ee);
174
return R;
175
}
176
177
/***************************************************************************/
178
/* Functions to loop over points in a hypercube. */
179
180
/* Based on orbitrule.cpp in HIntLib-0.0.10 */
181
182
/* ls0 returns the least-significant 0 bit of n (e.g. it returns
183
0 if the LSB is 0, it returns 1 if the 2 LSBs are 01, etcetera). */
184
185
#if (defined(__GNUC__) || defined(__ICC)) && (defined(__i386__) || defined (__x86_64__))
186
/* use x86 bit-scan instruction, based on count_trailing_zeros()
187
macro in GNU GMP's longlong.h. */
188
static unsigned ls0(unsigned n)
189
{
190
unsigned count;
191
n = ~n;
192
__asm__("bsfl %1,%0": "=r"(count):"rm"(n));
193
return count;
194
}
195
#else
196
static unsigned ls0(unsigned n)
197
{
198
const unsigned bits[256] = {
199
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
200
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
201
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
202
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
203
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
204
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
205
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
206
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7,
207
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
208
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
209
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
210
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
211
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
212
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
213
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
214
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8,
215
};
216
unsigned bit = 0;
217
while ((n & 0xff) == 0xff) {
218
n >>= 8;
219
bit += 8;
220
}
221
return bit + bits[n & 0xff];
222
}
223
#endif
224
225
/**
226
* Evaluate the integral on all 2^n points (+/-r,...+/-r)
227
*
228
* A Gray-code ordering is used to minimize the number of coordinate updates
229
* in p.
230
*/
231
static double evalR_Rfs(integrand f, void *fdata, unsigned dim, double *p, const double *c, const double *r)
232
{
233
double sum = 0;
234
unsigned i;
235
unsigned signs = 0; /* 0/1 bit = +/- for corresponding element of r[] */
236
237
/* We start with the point where r is ADDed in every coordinate
238
(this implies signs=0). */
239
for (i = 0; i < dim; ++i)
240
p[i] = c[i] + r[i];
241
242
/* Loop through the points in Gray-code ordering */
243
for (i = 0;; ++i) {
244
unsigned mask, d;
245
246
sum += f(dim, p, fdata);
247
248
d = ls0(i); /* which coordinate to flip */
249
if (d >= dim)
250
break;
251
252
/* flip the d-th bit and add/subtract r[d] */
253
mask = 1U << d;
254
signs ^= mask;
255
p[d] = (signs & mask) ? c[d] - r[d] : c[d] + r[d];
256
}
257
return sum;
258
}
259
260
static double evalRR0_0fs(integrand f, void *fdata, unsigned dim, double *p, const double *c, const double *r)
261
{
262
unsigned i, j;
263
double sum = 0;
264
265
for (i = 0; i < dim - 1; ++i) {
266
p[i] = c[i] - r[i];
267
for (j = i + 1; j < dim; ++j) {
268
p[j] = c[j] - r[j];
269
sum += f(dim, p, fdata);
270
p[i] = c[i] + r[i];
271
sum += f(dim, p, fdata);
272
p[j] = c[j] + r[j];
273
sum += f(dim, p, fdata);
274
p[i] = c[i] - r[i];
275
sum += f(dim, p, fdata);
276
277
p[j] = c[j]; /* Done with j -> Restore p[j] */
278
}
279
p[i] = c[i]; /* Done with i -> Restore p[i] */
280
}
281
return sum;
282
}
283
284
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_)
285
{
286
double maxdiff = 0;
287
unsigned i, dimDiffMax = 0;
288
double sum0, sum1 = 0, sum2 = 0; /* copies for aliasing, performance */
289
290
double ratio = r1[0] / r2[0];
291
292
ratio *= ratio;
293
sum0 = f(dim, p, fdata);
294
295
for (i = 0; i < dim; i++) {
296
double f1a, f1b, f2a, f2b, diff;
297
298
p[i] = c[i] - r1[i];
299
sum1 += (f1a = f(dim, p, fdata));
300
p[i] = c[i] + r1[i];
301
sum1 += (f1b = f(dim, p, fdata));
302
p[i] = c[i] - r2[i];
303
sum2 += (f2a = f(dim, p, fdata));
304
p[i] = c[i] + r2[i];
305
sum2 += (f2b = f(dim, p, fdata));
306
p[i] = c[i];
307
308
diff = fabs(f1a + f1b - 2 * sum0 - ratio * (f2a + f2b - 2 * sum0));
309
310
if (diff > maxdiff) {
311
maxdiff = diff;
312
dimDiffMax = i;
313
}
314
}
315
316
*sum0_ += sum0;
317
*sum1_ += sum1;
318
*sum2_ += sum2;
319
320
return dimDiffMax;
321
}
322
323
#define num0_0(dim) (1U)
324
#define numR0_0fs(dim) (2 * (dim))
325
#define numRR0_0fs(dim) (2 * (dim) * (dim-1))
326
#define numR_Rfs(dim) (1U << (dim))
327
328
/***************************************************************************/
329
/* Based on rule75genzmalik.cpp in HIntLib-0.0.10: An embedded
330
cubature rule of degree 7 (embedded rule degree 5) due to A. C. Genz
331
and A. A. Malik. See:
332
333
A. C. Genz and A. A. Malik, "An imbedded [sic] family of fully
334
symmetric numerical integration rules," SIAM
335
J. Numer. Anal. 20 (3), 580-588 (1983).
336
*/
337
338
typedef struct {
339
rule parent;
340
341
/* temporary arrays of length dim */
342
double *widthLambda, *widthLambda2, *p;
343
344
/* dimension-dependent constants */
345
double weight1, weight3, weight5;
346
double weightE1, weightE3;
347
} rule75genzmalik;
348
349
#define real(x) ((double)(x))
350
#define to_int(n) ((int)(n))
351
352
static int isqr(int x)
353
{
354
return x * x;
355
}
356
357
static void destroy_rule75genzmalik(rule *r_)
358
{
359
rule75genzmalik *r = (rule75genzmalik *) r_;
360
free(r->p);
361
}
362
363
static unsigned rule75genzmalik_evalError(rule *r_, integrand f, void *fdata, const hypercube *h, esterr *ee)
364
{
365
/* lambda2 = sqrt(9/70), lambda4 = sqrt(9/10), lambda5 = sqrt(9/19) */
366
const double lambda2 = 0.3585685828003180919906451539079374954541;
367
const double lambda4 = 0.9486832980505137995996680633298155601160;
368
const double lambda5 = 0.6882472016116852977216287342936235251269;
369
const double weight2 = 980. / 6561.;
370
const double weight4 = 200. / 19683.;
371
const double weightE2 = 245. / 486.;
372
const double weightE4 = 25. / 729.;
373
374
rule75genzmalik *r = (rule75genzmalik *) r_;
375
unsigned i, dimDiffMax, dim = r_->dim;
376
double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, sum4, sum5, result, res5th;
377
const double *center = h->data;
378
const double *halfwidth = h->data + dim;
379
380
for (i = 0; i < dim; ++i)
381
r->p[i] = center[i];
382
383
for (i = 0; i < dim; ++i)
384
r->widthLambda2[i] = halfwidth[i] * lambda2;
385
for (i = 0; i < dim; ++i)
386
r->widthLambda[i] = halfwidth[i] * lambda4;
387
388
/* Evaluate function in the center, in f(lambda2,0,...,0) and
389
f(lambda3=lambda4, 0,...,0). Estimate dimension with largest error */
390
dimDiffMax = evalR0_0fs4d(f, fdata, dim, r->p, center, &sum1, r->widthLambda2, &sum2, r->widthLambda, &sum3);
391
392
/* Calculate sum4 for f(lambda4, lambda4, 0, ...,0) */
393
sum4 = evalRR0_0fs(f, fdata, dim, r->p, center, r->widthLambda);
394
395
/* Calculate sum5 for f(lambda5, lambda5, ..., lambda5) */
396
for (i = 0; i < dim; ++i)
397
r->widthLambda[i] = halfwidth[i] * lambda5;
398
sum5 = evalR_Rfs(f, fdata, dim, r->p, center, r->widthLambda);
399
400
/* Calculate fifth and seventh order results */
401
402
result = h->vol * (r->weight1 * sum1 + weight2 * sum2 + r->weight3 * sum3 + weight4 * sum4 + r->weight5 * sum5);
403
res5th = h->vol * (r->weightE1 * sum1 + weightE2 * sum2 + r->weightE3 * sum3 + weightE4 * sum4);
404
405
ee->val = result;
406
ee->err = fabs(res5th - result);
407
408
return dimDiffMax;
409
}
410
411
static rule *make_rule75genzmalik(unsigned dim)
412
{
413
rule75genzmalik *r;
414
415
if (dim < 2) return 0; /* this rule does not support 1d integrals */
416
417
/* Because of the use of a bit-field in evalR_Rfs, we are limited
418
to be < 32 dimensions (or however many bits are in unsigned).
419
This is not a practical limitation...long before you reach
420
32 dimensions, the Genz-Malik cubature becomes excruciatingly
421
slow and is superseded by other methods (e.g. Monte-Carlo). */
422
if (dim >= sizeof(unsigned) * 8) return 0;
423
424
r = (rule75genzmalik *) malloc(sizeof(rule75genzmalik));
425
r->parent.dim = dim;
426
427
r->weight1 = (real(12824 - 9120 * to_int(dim) + 400 * isqr(to_int(dim)))
428
/ real(19683));
429
r->weight3 = real(1820 - 400 * to_int(dim)) / real(19683);
430
r->weight5 = real(6859) / real(19683) / real(1U << dim);
431
r->weightE1 = (real(729 - 950 * to_int(dim) + 50 * isqr(to_int(dim)))
432
/ real(729));
433
r->weightE3 = real(265 - 100 * to_int(dim)) / real(1458);
434
435
r->p = (double *) malloc(sizeof(double) * dim * 3);
436
r->widthLambda = r->p + dim;
437
r->widthLambda2 = r->p + 2 * dim;
438
439
r->parent.num_points = num0_0(dim) + 2 * numR0_0fs(dim)
440
+ numRR0_0fs(dim) + numR_Rfs(dim);
441
442
r->parent.evalError = rule75genzmalik_evalError;
443
r->parent.destroy = destroy_rule75genzmalik;
444
445
return (rule *) r;
446
}
447
448
/***************************************************************************/
449
/* 1d 15-point Gaussian quadrature rule, based on qk15.c and qk.c in
450
GNU GSL (which in turn is based on QUADPACK). */
451
452
static unsigned rule15gauss_evalError(rule *r, integrand f, void *fdata,
453
const hypercube *h, esterr *ee)
454
{
455
/* Gauss quadrature weights and kronrod quadrature abscissae and
456
weights as evaluated with 80 decimal digit arithmetic by
457
L. W. Fullerton, Bell Labs, Nov. 1981. */
458
const unsigned n = 8;
459
const double xgk[8] = { /* abscissae of the 15-point kronrod rule */
460
0.991455371120812639206854697526329,
461
0.949107912342758524526189684047851,
462
0.864864423359769072789712788640926,
463
0.741531185599394439863864773280788,
464
0.586087235467691130294144838258730,
465
0.405845151377397166906606412076961,
466
0.207784955007898467600689403773245,
467
0.000000000000000000000000000000000
468
/* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
469
xgk[0], xgk[2], ... to optimally extend the 7-point gauss rule */
470
};
471
static const double wg[4] = { /* weights of the 7-point gauss rule */
472
0.129484966168869693270611432679082,
473
0.279705391489276667901467771423780,
474
0.381830050505118944950369775488975,
475
0.417959183673469387755102040816327
476
};
477
static const double wgk[8] = { /* weights of the 15-point kronrod rule */
478
0.022935322010529224963732008058970,
479
0.063092092629978553290700663189204,
480
0.104790010322250183839876322541518,
481
0.140653259715525918745189590510238,
482
0.169004726639267902826583426598550,
483
0.190350578064785409913256402421014,
484
0.204432940075298892414161999234649,
485
0.209482141084727828012999174891714
486
};
487
488
const double center = h->data[0];
489
const double halfwidth = h->data[1];
490
double fv1[7], fv2[7];
491
const double f_center = f(1, ¢er, fdata);
492
double result_gauss = f_center * wg[n/2 - 1];
493
double result_kronrod = f_center * wgk[n - 1];
494
double result_abs = fabs(result_kronrod);
495
double result_asc, mean, err;
496
unsigned j;
497
498
for (j = 0; j < (n - 1) / 2; ++j) {
499
int j2 = 2*j + 1;
500
double x, f1, f2, fsum, w = halfwidth * xgk[j2];
501
x = center - w; fv1[j2] = f1 = f(1, &x, fdata);
502
x = center + w; fv2[j2] = f2 = f(1, &x, fdata);
503
fsum = f1 + f2;
504
result_gauss += wg[j] * fsum;
505
result_kronrod += wgk[j2] * fsum;
506
result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
507
}
508
509
for (j = 0; j < n/2; ++j) {
510
int j2 = 2*j;
511
double x, f1, f2, w = halfwidth * xgk[j2];
512
x = center - w; fv1[j2] = f1 = f(1, &x, fdata);
513
x = center + w; fv2[j2] = f2 = f(1, &x, fdata);
514
result_kronrod += wgk[j2] * (f1 + f2);
515
result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
516
}
517
518
ee->val = result_kronrod * halfwidth;
519
520
/* compute error estimate: */
521
mean = result_kronrod * 0.5;
522
result_asc = wgk[n - 1] * fabs(f_center - mean);
523
for (j = 0; j < n - 1; ++j)
524
result_asc += wgk[j] * (fabs(fv1[j]-mean) + fabs(fv2[j]-mean));
525
err = fabs(result_kronrod - result_gauss) * halfwidth;
526
result_abs *= halfwidth;
527
result_asc *= halfwidth;
528
if (result_asc != 0 && err != 0) {
529
double scale = pow((200 * err / result_asc), 1.5);
530
if (scale < 1)
531
err = result_asc * scale;
532
else
533
err = result_asc;
534
}
535
if (result_abs > DBL_MIN / (50 * DBL_EPSILON)) {
536
double min_err = 50 * DBL_EPSILON * result_abs;
537
if (min_err > err)
538
err = min_err;
539
}
540
ee->err = err;
541
542
return 0; /* no choice but to divide 0th dimension */
543
}
544
545
static rule *make_rule15gauss(unsigned dim)
546
{
547
rule *r;
548
if (dim != 1) return 0; /* this rule is only for 1d integrals */
549
r = (rule *) malloc(sizeof(rule));
550
r->dim = dim;
551
r->num_points = 15;
552
r->evalError = rule15gauss_evalError;
553
r->destroy = 0;
554
return r;
555
}
556
557
/***************************************************************************/
558
/* binary heap implementation (ala _Introduction to Algorithms_ by
559
Cormen, Leiserson, and Rivest), for use as a priority queue of
560
regions to integrate. */
561
562
typedef region heap_item;
563
#define KEY(hi) ((hi).ee.err)
564
565
typedef struct {
566
unsigned n, nalloc;
567
heap_item *items;
568
esterr ee;
569
} heap;
570
571
static void heap_resize(heap *h, unsigned nalloc)
572
{
573
h->nalloc = nalloc;
574
h->items = (heap_item *) realloc(h->items, sizeof(heap_item) * nalloc);
575
}
576
577
static heap heap_alloc(unsigned nalloc)
578
{
579
heap h;
580
h.n = 0;
581
h.nalloc = 0;
582
h.items = 0;
583
h.ee.val = h.ee.err = 0;
584
heap_resize(&h, nalloc);
585
return h;
586
}
587
588
/* note that heap_free does not deallocate anything referenced by the items */
589
static void heap_free(heap *h)
590
{
591
h->n = 0;
592
heap_resize(h, 0);
593
}
594
595
static void heap_push(heap *h, heap_item hi)
596
{
597
int insert;
598
599
h->ee.val += hi.ee.val;
600
h->ee.err += hi.ee.err;
601
insert = h->n;
602
if (++(h->n) > h->nalloc)
603
heap_resize(h, h->n * 2);
604
605
while (insert) {
606
int parent = (insert - 1) / 2;
607
if (KEY(hi) <= KEY(h->items[parent]))
608
break;
609
h->items[insert] = h->items[parent];
610
insert = parent;
611
}
612
h->items[insert] = hi;
613
}
614
615
static heap_item heap_pop(heap *h)
616
{
617
heap_item ret;
618
int i, n, child;
619
620
if (!(h->n)) {
621
fprintf(stderr, "attempted to pop an empty heap\n");
622
exit(EXIT_FAILURE);
623
}
624
625
ret = h->items[0];
626
h->items[i = 0] = h->items[n = --(h->n)];
627
while ((child = i * 2 + 1) < n) {
628
int largest;
629
heap_item swap;
630
631
if (KEY(h->items[child]) <= KEY(h->items[i]))
632
largest = i;
633
else
634
largest = child;
635
if (++child < n && KEY(h->items[largest]) < KEY(h->items[child]))
636
largest = child;
637
if (largest == i)
638
break;
639
swap = h->items[i];
640
h->items[i] = h->items[largest];
641
h->items[i = largest] = swap;
642
}
643
644
h->ee.val -= ret.ee.val;
645
h->ee.err -= ret.ee.err;
646
return ret;
647
}
648
649
/***************************************************************************/
650
651
/* adaptive integration, analogous to adaptintegrator.cpp in HIntLib */
652
653
static int ruleadapt_integrate(rule *r, integrand f, void *fdata, const hypercube *h, unsigned maxEval, double reqAbsError, double reqRelError, esterr *ee)
654
{
655
unsigned maxIter; /* maximum number of adaptive subdivisions */
656
heap regions;
657
unsigned i;
658
int status = -1; /* = ERROR */
659
660
if (maxEval) {
661
if (r->num_points > maxEval)
662
return status; /* ERROR */
663
maxIter = (maxEval - r->num_points) / (2 * r->num_points);
664
}
665
else
666
maxIter = UINT_MAX;
667
668
regions = heap_alloc(1);
669
670
heap_push(®ions, eval_region(make_region(h), f, fdata, r));
671
/* another possibility is to specify some non-adaptive subdivisions:
672
if (initialRegions != 1)
673
partition(h, initialRegions, EQUIDISTANT, ®ions, f,fdata, r); */
674
675
for (i = 0; i < maxIter; ++i) {
676
region R, R2;
677
if (regions.ee.err <= reqAbsError
678
|| relError(regions.ee) <= reqRelError) {
679
status = 0; /* converged! */
680
break;
681
}
682
R = heap_pop(®ions); /* get worst region */
683
cut_region(&R, &R2);
684
heap_push(®ions, eval_region(R, f, fdata, r));
685
heap_push(®ions, eval_region(R2, f, fdata, r));
686
}
687
688
ee->val = ee->err = 0; /* re-sum integral and errors */
689
for (i = 0; i < regions.n; ++i) {
690
ee->val += regions.items[i].ee.val;
691
ee->err += regions.items[i].ee.err;
692
destroy_region(®ions.items[i]);
693
}
694
/* printf("regions.nalloc = %d\n", regions.nalloc); */
695
heap_free(®ions);
696
697
return status;
698
}
699
700
static int adapt_integrate(integrand f, void *fdata,
701
unsigned dim, const double *xmin, const double *xmax,
702
unsigned maxEval, double reqAbsError, double reqRelError,
703
double *val, double *err)
704
{
705
rule *r;
706
hypercube h;
707
esterr ee;
708
int status;
709
710
if (dim == 0) { /* trivial integration */
711
*val = f(0, xmin, fdata);
712
*err = 0;
713
return 0;
714
}
715
r = dim == 1 ? make_rule15gauss(dim) : make_rule75genzmalik(dim);
716
if (!r) { *val = 0; *err = HUGE_VAL; return -2; /* ERROR */ }
717
h = make_hypercube_range(dim, xmin, xmax);
718
status = ruleadapt_integrate(r, f, fdata, &h,
719
maxEval, reqAbsError, reqRelError,
720
&ee);
721
*val = ee.val;
722
*err = ee.err;
723
destroy_hypercube(&h);
724
destroy_rule(r);
725
return status;
726
}
727
728
/***************************************************************************/
729
730
/* Compile with -DTEST_INTEGRATOR for a self-contained test program.
731
732
Usage: ./integrator <dim> <tol> <integrand> <maxeval>
733
734
where <dim> = # dimensions, <tol> = relative tolerance,
735
<integrand> is either 0/1/2 for the three test integrands (see below),
736
and <maxeval> is the maximum # function evaluations (0 for none).
737
*/
738
739
#ifdef TEST_INTEGRATOR
740
741
int count = 0;
742
int which_integrand = 0;
743
const double radius = 0.50124145262344534123412; /* random */
744
745
/* Simple constant function */
746
double
747
fconst (double x[], size_t dim, void *params)
748
{
749
return 1;
750
}
751
752
/*** f0, f1, f2, and f3 are test functions from the Monte-Carlo
753
integration routines in GSL 1.6 (monte/test.c). Copyright (c)
754
1996-2000 Michael Booth, GNU GPL. ****/
755
756
/* Simple product function */
757
double f0 (unsigned dim, const double *x, void *params)
758
{
759
double prod = 1.0;
760
unsigned int i;
761
for (i = 0; i < dim; ++i)
762
prod *= 2.0 * x[i];
763
return prod;
764
}
765
766
/* Gaussian centered at 1/2. */
767
double f1 (unsigned dim, const double *x, void *params)
768
{
769
double a = *(double *)params;
770
double sum = 0.;
771
unsigned int i;
772
for (i = 0; i < dim; i++) {
773
double dx = x[i] - 0.5;
774
sum += dx * dx;
775
}
776
return (pow (M_2_SQRTPI / (2. * a), (double) dim) *
777
exp (-sum / (a * a)));
778
}
779
780
/* double gaussian */
781
double f2 (unsigned dim, const double *x, void *params)
782
{
783
double a = *(double *)params;
784
double sum1 = 0.;
785
double sum2 = 0.;
786
unsigned int i;
787
for (i = 0; i < dim; i++) {
788
double dx1 = x[i] - 1. / 3.;
789
double dx2 = x[i] - 2. / 3.;
790
sum1 += dx1 * dx1;
791
sum2 += dx2 * dx2;
792
}
793
return 0.5 * pow (M_2_SQRTPI / (2. * a), dim)
794
* (exp (-sum1 / (a * a)) + exp (-sum2 / (a * a)));
795
}
796
797
/* Tsuda's example */
798
double f3 (unsigned dim, const double *x, void *params)
799
{
800
double c = *(double *)params;
801
double prod = 1.;
802
unsigned int i;
803
for (i = 0; i < dim; i++)
804
prod *= c / (c + 1) * pow((c + 1) / (c + x[i]), 2.0);
805
return prod;
806
}
807
808
/*** end of GSL test functions ***/
809
810
double f_test(unsigned dim, const double *x, void *data)
811
{
812
double val;
813
unsigned i;
814
++count;
815
switch (which_integrand) {
816
case 0: /* simple smooth (separable) objective: prod. cos(x[i]). */
817
val = 1;
818
for (i = 0; i < dim; ++i)
819
val *= cos(x[i]);
820
break;
821
case 1: { /* integral of exp(-x^2), rescaled to (0,infinity) limits */
822
double scale = 1.0;
823
val = 0;
824
for (i = 0; i < dim; ++i) {
825
double z = (1 - x[i]) / x[i];
826
val += z * z;
827
scale *= M_2_SQRTPI / (x[i] * x[i]);
828
}
829
val = exp(-val) * scale;
830
break;
831
}
832
case 2: /* discontinuous objective: volume of hypersphere */
833
val = 0;
834
for (i = 0; i < dim; ++i)
835
val += x[i] * x[i];
836
val = val < radius * radius;
837
break;
838
case 3:
839
val = f0(dim, x, data);
840
break;
841
case 4:
842
val = f1(dim, x, data);
843
break;
844
case 5:
845
val = f2(dim, x, data);
846
break;
847
case 6:
848
val = f3(dim, x, data);
849
break;
850
default:
851
fprintf(stderr, "unknown integrand %d\n", which_integrand);
852
exit(EXIT_FAILURE);
853
}
854
/* if (count < 100) printf("%d: f(%g, ...) = %g\n", count, x[0], val); */
855
return val;
856
}
857
858
/* surface area of n-dimensional unit hypersphere */
859
static double S(unsigned n)
860
{
861
double val;
862
int fact = 1;
863
if (n % 2 == 0) { /* n even */
864
val = 2 * pow(M_PI, n * 0.5);
865
n = n / 2;
866
while (n > 1) fact *= (n -= 1);
867
val /= fact;
868
}
869
else { /* n odd */
870
val = (1 << (n/2 + 1)) * pow(M_PI, n/2);
871
while (n > 2) fact *= (n -= 2);
872
val /= fact;
873
}
874
return val;
875
}
876
877
static double exact_integral(unsigned dim, const double *xmax) {
878
unsigned i;
879
double val;
880
switch(which_integrand) {
881
case 0:
882
val = 1;
883
for (i = 0; i < dim; ++i)
884
val *= sin(xmax[i]);
885
break;
886
case 2:
887
val = dim == 0 ? 1 : S(dim) * pow(radius * 0.5, dim) / dim;
888
break;
889
default:
890
val = 1.0;
891
}
892
return val;
893
}
894
895
int main(int argc, char **argv)
896
{
897
double *xmin, *xmax;
898
double tol, val, err;
899
unsigned i, dim, maxEval;
900
double fdata;
901
902
dim = argc > 1 ? atoi(argv[1]) : 2;
903
tol = argc > 2 ? atof(argv[2]) : 1e-2;
904
which_integrand = argc > 3 ? atoi(argv[3]) : 0;
905
maxEval = argc > 4 ? atoi(argv[4]) : 0;
906
907
fdata = which_integrand == 6 ? (1.0 + sqrt (10.0)) / 9.0 : 0.1;
908
909
xmin = (double *) malloc(dim * sizeof(double));
910
xmax = (double *) malloc(dim * sizeof(double));
911
for (i = 0; i < dim; ++i) {
912
xmin[i] = 0;
913
xmax[i] = 1 + (which_integrand >= 2 ? 0 : 0.4 * sin(i));
914
}
915
916
printf("%u-dim integral, tolerance = %g, integrand = %d\n",
917
dim, tol, which_integrand);
918
adapt_integrate(f_test, &fdata,
919
dim, xmin, xmax,
920
maxEval, 0, tol, &val, &err);
921
printf("integration val = %g, est. err = %g, true err = %g\n",
922
val, err, fabs(val - exact_integral(dim, xmax)));
923
printf("#evals = %d\n", count);
924
925
free(xmax);
926
free(xmin);
927
928
return 0;
929
}
930
931
#else
932
933
/*************************************************************************/
934
/* libctl interface */
935
936
#include "ctl.h"
937
938
static int adapt_integrate(integrand f, void *fdata,
939
unsigned dim, const double *xmin, const double *xmax,
940
unsigned maxEval,
941
double reqAbsError, double reqRelError,
942
double *val, double *err);
943
944
number adaptive_integration(multivar_func f, number *xmin, number *xmax,
945
integer n, void *fdata,
946
number abstol, number reltol, integer maxnfe,
947
number *esterr, integer *errflag)
948
{
949
double val;
950
*errflag = adapt_integrate((integrand) f, fdata, n, xmin, xmax,
951
maxnfe, abstol, reltol, &val, esterr);
952
return val;
953
}
954
955
/* from subplex.c */
956
extern number f_scm_wrapper(integer n, number *x, void *f_scm_p);
957
958
SCM adaptive_integration_scm(SCM f_scm, SCM xmin_scm, SCM xmax_scm,
959
SCM abstol_scm, SCM reltol_scm, SCM maxnfe_scm)
960
{
961
integer n, maxnfe, errflag, i;
962
number *xmin, *xmax, abstol, reltol, integral;
963
964
n = list_length(xmin_scm);
965
abstol = fabs(gh_scm2double(abstol_scm));
966
reltol = fabs(gh_scm2double(reltol_scm));
967
maxnfe = gh_scm2int(maxnfe_scm);
968
969
if (list_length(xmax_scm) != n) {
970
fprintf(stderr, "adaptive_integration: xmin/xmax dimension mismatch\n");
971
return SCM_UNDEFINED;
972
}
973
974
xmin = (number*) malloc(sizeof(number) * n);
975
xmax = (number*) malloc(sizeof(number) * n);
976
if (!xmin || !xmax) {
977
fprintf(stderr, "adaptive_integration: error, out of memory!\n");
978
exit(EXIT_FAILURE);
979
}
980
981
for (i = 0; i < n; ++i) {
982
xmin[i] = number_list_ref(xmin_scm, i);
983
xmax[i] = number_list_ref(xmax_scm, i);
984
}
985
986
integral = adaptive_integration(f_scm_wrapper, xmin, xmax, n, &f_scm,
987
abstol, reltol, maxnfe,
988
&abstol, &errflag);
989
990
free(xmax);
991
free(xmin);
992
993
switch (errflag) {
994
case 3:
995
fprintf(stderr, "adaptive_integration: invalid inputs\n");
996
return SCM_UNDEFINED;
997
case 1:
998
fprintf(stderr, "adaptive_integration: maxnfe too small\n");
999
break;
1000
case 2:
1001
fprintf(stderr, "adaptive_integration: lenwork too small\n");
1002
break;
1003
}
1004
1005
return gh_cons(gh_double2scm(integral), gh_double2scm(abstol));
1006
}
1007
1008
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1