PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
test_cube.c
Go to the documentation of this file.
1
2 /*
3  Usage: ./test_cube <dim> <tol> <integrand> <maxeval>
4
5  where <dim> = # dimensions, <tol> = relative tolerance,
6  <integrand> is either 0/1/2 for the three test integrands (see below),
7  and <maxeval> is the maximum # function evaluations (0 for none).
8 */
9
10 /* author Steven G. Johnson */
11
12 #include "cubature.h"
13 #include <math.h>
14 #include <gsl/gsl_math.h> /* M_PI */
15
16 int count = 0;
18 const double radius = 0.50124145262344534123412; /* random */
19
20 double f_test(unsigned dim, const double *x, void *data)
21 {
22  double val;
23  unsigned i;
24  ++count;
25  switch (which_integrand) {
26  case 0: /* discontinuous objective: volume of hypersphere */
27  val = 0;
28  for (i = 0; i < dim; ++i)
29  val += x[i] * x[i];
31  break;
32  case 1: /* simple smooth (separable) objective: prod. cos(x[i]). */
33  val = 1;
34  for (i = 0; i < dim; ++i)
35  val *= cos(x[i]);
36  break;
37  case 2: { /* integral of exp(-x^2), rescaled to (0,infinity) limits */
38  double scale = 1.0;
39  val = 0;
40  for (i = 0; i < dim; ++i) {
41  double z = (1 - x[i]) / x[i];
42  val += z * z;
43  scale *= M_2_SQRTPI / (x[i] * x[i]);
44  }
45  val = exp(-val) * scale;
46  break;
47  }
48
49  default:
50  fprintf(stderr, "unknown integrand %d\n", which_integrand);
51  exit(EXIT_FAILURE);
52  }
53  /* if (count < 100) printf("%d: f(%g, ...) = %g\n", count, x[0], val); */
54  return val;
55 }
56
57 /* surface area of n-dimensional unit hypersphere */
58 static double S(unsigned n)
59 {
60  double val;
61  int fact = 1;
62  if (n % 2 == 0) { /* n even */
63  val = 2 * pow(M_PI, n * 0.5);
64  n = n / 2;
65  while (n > 1) fact *= (n -= 1);
66  val /= fact;
67  }
68  else { /* n odd */
69  val = (1 << (n/2 + 1)) * pow(M_PI, n/2);
70  while (n > 2) fact *= (n -= 2);
71  val /= fact;
72  }
73  return val;
74 }
75
76 static double exact_integral(unsigned dim, const double *xmax) {
77  unsigned i;
78  double val;
79  switch(which_integrand) {
80  case 0:
81  val = dim == 0 ? 1 : S(dim) * pow(radius * 0.5, dim) / dim;
82  break;
83  case 1:
84  val = 1;
85  for (i = 0; i < dim; ++i)
86  val *= sin(xmax[i]);
87  break;
88  case 2:
89  val = 1;
90  break;
91  default:
92  fprintf(stderr, "unknown integrand %d\n", which_integrand);
93  exit(EXIT_FAILURE);
94  }
95  return val;
96 }
97
98 int main(int argc, char **argv)
99 {
100  double *xmin, *xmax;
101  double tol, val, err;
102  unsigned i, dim, maxEval;
103
104  dim = argc > 1 ? atoi(argv[1]) : 2;
105  tol = argc > 2 ? atof(argv[2]) : 1e-2;
106  which_integrand = argc > 3 ? atoi(argv[3]) : 1;
107  maxEval = argc > 4 ? atoi(argv[4]) : 0;
108
109  xmin = (double *) malloc(dim * sizeof(double));
110  xmax = (double *) malloc(dim * sizeof(double));
111  for (i = 0; i < dim; ++i) {
112  xmin[i] = 0;
113  xmax[i] = 1 + (which_integrand == 2 ? 0 : 0.4 * sin(i));
114  }
115
116  printf("%u-dim integral, tolerance = %g, integrand = %d\n",
117  dim, tol, which_integrand);
118  adapt_integrate(f_test, 0, dim, xmin, xmax, maxEval, 0, tol, &val, &err);
119  printf("integration val = %g, est. err = %g, true err = %g\n",
120  val, err, fabs(val - exact_integral(dim, xmax)));
121  printf("#evals = %d\n", count);
122
123  free(xmax);
124  free(xmin);
125
126  return 0;
127 }
128
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)
Definition: cubature.c:701
#define n
Definition: exactTestM.c:37
std::string printf(const char *format,...)
static double S(unsigned n)
Definition: test_cube.c:58