ttest_cube.c - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
 (HTM) git clone git://src.adamsgaard.dk/pism
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       ttest_cube.c (3597B)
       ---
            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;
           17 int which_integrand = 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];
           30           val = val < radius * radius;
           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