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