Imported Upstream version 3.0+20111124upstream/3.0+20111124upstream

6 files changed, 609 insertions, 0 deletions

diff --git a/demo/cuba.F b/demo/cuba.F
new file mode 100644
index 0000000..a227d4f
--- /dev/null
+++ b/demo/cuba.F
@@ -0,0 +1,104 @@
+* cuba.F
+* Fortran chooser for the Cuba routines
+* last modified 3 Feb 05 th
+
+#define VEGAS 1
+#define SUAVE 2
+#define DIVONNE 3
+#define CUHRE 4
+
+
+	subroutine Cuba(method, ndim, ncomp, integrand,
+     &    integral, error, prob)
+	implicit none
+	integer method, ndim, ncomp
+	external integrand
+	double precision integral(*), error(*), prob(*)
+
+	character*7 name(4)
+	data name /"Vegas", "Suave", "Divonne", "Cuhre"/
+
+	integer mineval, maxeval, verbose, last
+	double precision epsrel, epsabs
+	parameter (epsrel = 1D-3)
+	parameter (epsabs = 1D-12)
+	parameter (verbose = 2)
+	parameter (last = 4)
+	parameter (mineval = 0)
+	parameter (maxeval = 50000)
+
+	integer nstart, nincrease
+	parameter (nstart = 1000)
+	parameter (nincrease = 500)
+
+	integer nnew
+	double precision flatness
+	parameter (nnew = 1000)
+	parameter (flatness = 25D0)
+
+	integer key1, key2, key3, maxpass
+	double precision border, maxchisq, mindeviation
+	integer ngiven, ldxgiven, nextra
+	parameter (key1 = 47)
+	parameter (key2 = 1)
+	parameter (key3 = 1)
+	parameter (maxpass = 5)
+	parameter (border = 0D0)
+	parameter (maxchisq = 10D0)
+	parameter (mindeviation = .25D0)
+	parameter (ngiven = 0)
+	parameter (ldxgiven = ndim)
+	parameter (nextra = 0)
+
+	integer key
+	parameter (key = 0)
+
+	integer nregions, neval, fail
+
+
+	if( method .eq. VEGAS ) then
+
+	  call vegas(ndim, ncomp, integrand,
+     &      epsrel, epsabs, verbose, mineval, maxeval,
+     &      nstart, nincrease,
+     &      neval, fail, integral, error, prob)
+	  nregions = 1
+
+	else if( method .eq. SUAVE ) then
+
+	  call suave(ndim, ncomp, integrand,
+     &      epsrel, epsabs, verbose + last, mineval, maxeval,
+     &      nnew, flatness,
+     &      nregions, neval, fail, integral, error, prob)
+
+	else if( method .eq. DIVONNE ) then
+
+	  call divonne(ndim, ncomp, integrand,
+     &      epsrel, epsabs, verbose, mineval, maxeval,
+     &      key1, key2, key3, maxpass,
+     &      border, maxchisq, mindeviation,
+     &      ngiven, ldxgiven, 0, nextra, 0,
+     &      nregions, neval, fail, integral, error, prob)
+
+	else if( method .eq. CUHRE ) then
+
+	  call cuhre(ndim, ncomp, integrand,
+     &      epsrel, epsabs, verbose + last, mineval, maxeval,
+     &      key,
+     &      nregions, neval, fail, integral, error, prob)
+
+	else
+
+	  print *, "invalid method ", method
+	  return
+
+	endif
+
+	print *, "method   =", name(method)
+	print *, "nregions =", nregions
+	print *, "neval    =", neval
+	print *, "fail     =", fail
+	print '(G20.12," +- ",G20.12,"   p = ",F8.3)',
+     &    (integral(c), error(c), prob(c), c = 1, ncomp)
+	end
+
diff --git a/demo/demo-c.c b/demo/demo-c.c
new file mode 100644
index 0000000..300db4d
--- /dev/null
+++ b/demo/demo-c.c
@@ -0,0 +1,152 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "cuba.h"
+
+
+static inline double Sq(double x)
+{
+  return x*x;
+}
+
+
+static int Integrand(const int *ndim, const double xx[],
+  const int *ncomp, double ff[], void *userdata)
+{
+#define x xx[0]
+#define y xx[1]
+#define z xx[2]
+#define f ff[0]
+
+#define FUN 1
+
+  const double rsq = Sq(x) + Sq(y) + Sq(z);
+
+#if FUN == 1
+  f = sin(x)*cos(y)*exp(z);
+#elif FUN == 2
+  f = 1/(Sq(x + y) + .003)*cos(y)*exp(z);
+#elif FUN == 3
+  f = 1/(3.75 - cos(M_PI*x) - cos(M_PI*y) - cos(M_PI*z));
+#elif FUN == 4
+  f = fabs(rsq - .125);
+#elif FUN == 5
+  f = exp(-rsq);
+#elif FUN == 6
+  f = 1/(1 - x*y*z + 1e-10);
+#elif FUN == 7
+  f = sqrt(fabs(x - y - z));
+#elif FUN == 8
+  f = exp(-x*y*z);
+#elif FUN == 9
+  f = Sq(x)/(cos(x + y + z + 1) + 5);
+#elif FUN == 10
+  f = (x > .5) ? 1/sqrt(x*y*z + 1e-5) : sqrt(x*y*z);
+#else
+  f = (rsq < 1) ? 1 : 0;
+#endif
+
+  return 0;
+}
+
+/*********************************************************************/
+
+#define NDIM 3
+#define NCOMP 1
+#define USERDATA NULL
+#define EPSREL 1e-3
+#define EPSABS 1e-12
+#define LAST 4
+#define SEED 0
+#define MINEVAL 0
+#define MAXEVAL 50000
+
+#define NSTART 1000
+#define NINCREASE 500
+#define NBATCH 1000
+#define GRIDNO 0
+#define STATEFILE NULL
+
+#define NNEW 1000
+#define FLATNESS 25.
+
+#define KEY1 47
+#define KEY2 1
+#define KEY3 1
+#define MAXPASS 5
+#define BORDER 0.
+#define MAXCHISQ 10.
+#define MINDEVIATION .25
+#define NGIVEN 0
+#define LDXGIVEN NDIM
+#define NEXTRA 0
+
+#define KEY 0
+
+int main()
+{
+  int verbose, comp, nregions, neval, fail;
+  double integral[NCOMP], error[NCOMP], prob[NCOMP];
+
+  const char *env = getenv("CUBAVERBOSE");
+  verbose = 2;
+  if( env ) verbose = atoi(env);
+
+  printf("-------------------- Vegas test --------------------\n");
+
+  Vegas(NDIM, NCOMP, Integrand, USERDATA,
+    EPSREL, EPSABS, verbose, SEED,
+    MINEVAL, MAXEVAL, NSTART, NINCREASE, NBATCH,
+    GRIDNO, STATEFILE,
+    &neval, &fail, integral, error, prob);
+
+  printf("VEGAS RESULT:\tneval %d\tfail %d\n",
+    neval, fail);
+  for( comp = 0; comp < NCOMP; ++comp )
+    printf("VEGAS RESULT:\t%.8f +- %.8f\tp = %.3f\n",
+      integral[comp], error[comp], prob[comp]);
+
+  printf("\n-------------------- Suave test --------------------\n");
+
+  Suave(NDIM, NCOMP, Integrand, USERDATA,
+    EPSREL, EPSABS, verbose | LAST, SEED,
+    MINEVAL, MAXEVAL, NNEW, FLATNESS,
+    &nregions, &neval, &fail, integral, error, prob);
+
+  printf("SUAVE RESULT:\tnregions %d\tneval %d\tfail %d\n",
+    nregions, neval, fail);
+  for( comp = 0; comp < NCOMP; ++comp )
+    printf("SUAVE RESULT:\t%.8f +- %.8f\tp = %.3f\n",
+      integral[comp], error[comp], prob[comp]);
+
+  printf("\n------------------- Divonne test -------------------\n");
+
+  Divonne(NDIM, NCOMP, Integrand, USERDATA,
+    EPSREL, EPSABS, verbose, SEED,
+    MINEVAL, MAXEVAL, KEY1, KEY2, KEY3, MAXPASS,
+    BORDER, MAXCHISQ, MINDEVIATION,
+    NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
+    &nregions, &neval, &fail, integral, error, prob);
+
+  printf("DIVONNE RESULT:\tnregions %d\tneval %d\tfail %d\n",
+    nregions, neval, fail);
+  for( comp = 0; comp < NCOMP; ++comp )
+    printf("DIVONNE RESULT:\t%.8f +- %.8f\tp = %.3f\n",
+      integral[comp], error[comp], prob[comp]);
+
+  printf("\n-------------------- Cuhre test --------------------\n");
+
+  Cuhre(NDIM, NCOMP, Integrand, USERDATA,
+    EPSREL, EPSABS, verbose | LAST,
+    MINEVAL, MAXEVAL, KEY,
+    &nregions, &neval, &fail, integral, error, prob);
+
+  printf("CUHRE RESULT:\tnregions %d\tneval %d\tfail %d\n",
+    nregions, neval, fail);
+  for( comp = 0; comp < NCOMP; ++comp )
+    printf("CUHRE RESULT:\t%.8f +- %.8f\tp = %.3f\n",
+      integral[comp], error[comp], prob[comp]);
+
+  return 0;
+}
+
diff --git a/demo/demo-c.out b/demo/demo-c.out
new file mode 100644
index 0000000..531971d
--- /dev/null
+++ b/demo/demo-c.out
@@ -0,0 +1,8 @@
+VEGAS RESULT:	neval 10000	fail 0
+VEGAS RESULT:	0.66481073 +- 0.00049218	p = 0.089
+SUAVE RESULT:	nregions 7	neval 7000	fail 0
+SUAVE RESULT:	0.66444529 +- 0.00056577	p = 0.210
+DIVONNE RESULT:	nregions 14	neval 3052	fail 0
+DIVONNE RESULT:	0.66461951 +- 0.00063503	p = 0.000
+CUHRE RESULT:	nregions 2	neval 381	fail 0
+CUHRE RESULT:	0.66466968 +- 0.00000000	p = 0.000
diff --git a/demo/demo-fortran.F b/demo/demo-fortran.F
new file mode 100644
index 0000000..4900a86
--- /dev/null
+++ b/demo/demo-fortran.F
@@ -0,0 +1,176 @@
+* demo-fortran.F
+* test program for the Cuba library
+* last modified 12 Aug 11 th
+
+
+	program CubaTest
+	implicit none
+
+	integer ndim, ncomp, last, seed, mineval, maxeval
+	double precision epsrel, epsabs, userdata
+	parameter (ndim = 3)
+	parameter (ncomp = 1)
+	parameter (userdata = 0)
+	parameter (epsrel = 1D-3)
+	parameter (epsabs = 1D-12)
+	parameter (last = 4)
+	parameter (seed = 0)
+	parameter (mineval = 0)
+	parameter (maxeval = 50000)
+
+	integer nstart, nincrease, nbatch, gridno
+	character*(*) statefile	
+	parameter (nstart = 1000)
+	parameter (nincrease = 500)
+	parameter (nbatch = 1000)
+	parameter (gridno = 0)
+	parameter (statefile = "")
+
+	integer nnew
+	double precision flatness
+	parameter (nnew = 1000)
+	parameter (flatness = 25D0)
+
+	integer key1, key2, key3, maxpass
+	double precision border, maxchisq, mindeviation
+	integer ngiven, ldxgiven, nextra
+	parameter (key1 = 47)
+	parameter (key2 = 1)
+	parameter (key3 = 1)
+	parameter (maxpass = 5)
+	parameter (border = 0D0)
+	parameter (maxchisq = 10D0)
+	parameter (mindeviation = .25D0)
+	parameter (ngiven = 0)
+	parameter (ldxgiven = ndim)
+	parameter (nextra = 0)
+
+	integer key
+	parameter (key = 0)
+
+	external integrand
+
+	double precision integral(ncomp), error(ncomp), prob(ncomp)
+	integer verbose, nregions, neval, fail
+	character*16 env
+
+	integer c
+
+	call getenv("CUBAVERBOSE", env)
+	verbose = 2
+	read(env, *, iostat=fail, end=999, err=999) verbose
+999	continue
+
+	print *, "-------------------- Vegas test --------------------"
+
+	call vegas(ndim, ncomp, integrand, userdata,
+     &    epsrel, epsabs, verbose, seed,
+     &    mineval, maxeval, nstart, nincrease, nbatch,
+     &    gridno, statefile,
+     &    neval, fail, integral, error, prob)
+
+	print *, "neval    =", neval
+	print *, "fail     =", fail
+	print '(F20.12," +- ",F20.12,"   p = ",F8.3)',
+     &    (integral(c), error(c), prob(c), c = 1, ncomp)
+
+	print *, " "
+	print *, "-------------------- Suave test --------------------"
+
+	call suave(ndim, ncomp, integrand, userdata,
+     &    epsrel, epsabs, verbose + last, seed,
+     &    mineval, maxeval, nnew, flatness,
+     &    nregions, neval, fail, integral, error, prob)
+
+	print *, "nregions =", nregions
+	print *, "neval    =", neval
+	print *, "fail     =", fail
+	print '(F20.12," +- ",F20.12,"   p = ",F8.3)',
+     &    (integral(c), error(c), prob(c), c = 1, ncomp)
+
+	print *, " "
+	print *, "------------------- Divonne test -------------------"
+
+	call divonne(ndim, ncomp, integrand, userdata,
+     &    epsrel, epsabs, verbose, seed,
+     &    mineval, maxeval, key1, key2, key3, maxpass,
+     &    border, maxchisq, mindeviation,
+     &    ngiven, ldxgiven, 0, nextra, 0,
+     &    nregions, neval, fail, integral, error, prob)
+
+	print *, "nregions =", nregions
+	print *, "neval    =", neval
+	print *, "fail     =", fail
+	print '(F20.12," +- ",F20.12,"   p = ",F8.3)',
+     &    (integral(c), error(c), prob(c), c = 1, ncomp)
+
+	print *, " "
+	print *, "-------------------- Cuhre test --------------------"
+
+	call cuhre(ndim, ncomp, integrand, userdata,
+     &    epsrel, epsabs, verbose + last,
+     &    mineval, maxeval, key,
+     &    nregions, neval, fail, integral, error, prob)
+
+	print *, "nregions =", nregions
+	print *, "neval    =", neval
+	print *, "fail     =", fail
+	print '(F20.12," +- ",F20.12,"   p = ",F8.3)',
+     &    (integral(c), error(c), prob(c), c = 1, ncomp)
+	end
+
+
+************************************************************************
+
+	integer function integrand(ndim, xx, ncomp, ff)
+	implicit none
+	integer ndim, ncomp
+	double precision xx(*), ff(*)
+
+#define x xx(1)
+#define y xx(2)
+#define z xx(3)
+#define f ff(1)
+
+#define FUN 1
+
+	double precision pi, rsq
+	parameter (pi = 3.14159265358979323846D0)
+
+	rsq = x**2 + y**2 + z**2
+
+#if FUN == 1
+	f = sin(x)*cos(y)*exp(z)
+#elif FUN == 2
+	f = 1/((x + y)**2 + .003D0)*cos(y)*exp(z)
+#elif FUN == 3
+	f = 1/(3.75D0 - cos(pi*x) - cos(pi*y) - cos(pi*z))
+#elif FUN == 4
+	f = abs(rsq - .125D0)
+#elif FUN == 5
+	f = exp(-rsq)
+#elif FUN == 6
+	f = 1/(1 - x*y*z + 1D-10)
+#elif FUN == 7
+	f = sqrt(abs(x - y - z))
+#elif FUN == 8
+	f = exp(-x*y*z)
+#elif FUN == 9
+	f = x**2/(cos(x + y + z + 1) + 5)
+#elif FUN == 10
+	if( x .gt. .5D0 ) then
+	  f = 1/sqrt(x*y*z + 1D-5)
+	else
+	  f = sqrt(x*y*z)
+	endif
+#else
+	if( rsq .lt. 1 ) then
+	  f = 1
+	else
+	  f = 0
+	endif
+#endif
+
+	integrand = 0
+	end
+
diff --git a/demo/demo-math.m b/demo/demo-math.m
new file mode 100644
index 0000000..2fe5a3d
--- /dev/null
+++ b/demo/demo-math.m
@@ -0,0 +1,36 @@
+Install["Vegas"]
+
+Install["Suave"]
+
+Install["Divonne"]
+
+Install["Cuhre"]
+
+
+test[n_] := {t[n, Vegas], t[n, Suave], t[n, Divonne], t[n, Cuhre]}
+
+t[n_, int_] := int[f[n][x, y, z], {x,0,1}, {y,0,1}, {z,0,1}]
+
+
+f[1][x_, y_, z_] := Sin[x] Cos[y] Exp[z]
+
+f[2][x_, y_, z_] := 1/((x + y)^2 + .003) Cos[y] Exp[z]
+
+f[3][x_, y_, z_] := 1/(3.75 - Cos[Pi x] - Cos[Pi y] - Cos[Pi z])
+
+f[4][x_, y_, z_] := Abs[x^2 + y^2 + z^2 - .125]
+
+f[5][x_, y_, z_] := Exp[-x^2 - y^2 - z^2]
+
+f[6][x_, y_, z_] := 1/(1 - x y z + 10^-10)
+
+f[7][x_, y_, z_] := Sqrt[Abs[x - y - z]]
+
+f[8][x_, y_, z_] := Exp[-x y z]
+
+f[9][x_, y_, z_] := x^2/(Cos[x + y + z + 1] + 5)
+
+f[10][x_, y_, z_] := If[ x > .5, 1/Sqrt[x y z + 10^-5], Sqrt[x y z] ]
+
+f[11][x_, y_, z_] := If[ x^2 + y^2 + z^2 < 1, 1, 0 ]
+
diff --git a/demo/testsuite.m b/demo/testsuite.m
new file mode 100644
index 0000000..9804b03
--- /dev/null
+++ b/demo/testsuite.m
@@ -0,0 +1,133 @@
+(* Test suite of Genz, used also by Sloan and Joe, and Novak and Ritter *)
+
+seed = 4711
+
+maxpoints = 150000
+
+repeat = 20
+
+
+(* Family 1: Oscillatory *)
+
+f[1][x_, c_, w_] := Cos[2 Pi w[[1]] + c.x]
+
+
+(* Family 2: Product peak *)
+
+f[2][x_, c_, w_] := Times@@ MapThread[f2a, {x, c, w}]
+
+f2a[xi_, ci_, wi_] := 1/(ci^-2 + (xi - wi)^2)
+
+
+(* Family 3: Corner peak *)
+
+f[3][x_, c_, w_] := (1 + c.x)^(-(Length[x] + 1))
+
+
+(* Family 4: Gaussian *)
+
+f[4][x_, c_, w_] := Exp[Plus@@ MapThread[f4a, {x, c, w}]]
+
+f4a[xi_, ci_, wi_] := -ci^2 (xi - wi)^2
+
+
+(* Family 5: Exponential *)
+
+f[5][x_, c_, w_] := Exp[Plus@@ MapThread[f5a, {x, c, w}]]
+
+f5a[xi_, ci_, wi_] := -ci Abs[xi - wi]
+
+
+(* Family 6: Discontinuous *)
+
+f[6][x_, c_, w_] := 0 /; x[[1]] > w[[1]] || x[[2]] > w[[2]]
+
+f[6][x_, c_, w_] := Exp[c.x]
+
+
+(* Novak & Ritter use
+difficulty[fam_] := {9.00, 7.25, 1.85, 7.03, 2.04, 4.30}[[fam]]
+*)
+
+(* Sloan & Joe use
+scale[dim_] := dim^Min[Max[.2 dim, 1], 2]
+
+SetOptions[Interpolation, InterpolationOrder -> 2]
+
+ifun[1] = Interpolation[{{5, 145.7}, {8, 354.0}, {10,  900.0}}];
+ifun[2] = Interpolation[{{5, 261.0}, {8, 545.0}, {10, 1760.0}}];
+ifun[3] = Interpolation[{{5, 433.0}, {8, 193.0}, {10,  185.0}}];
+ifun[4] = Interpolation[{{5, 155.0}, {8, 382.0}, {10, 1230.0}}];
+ifun[5] = Interpolation[{{5, 217.0}, {8, 674.0}, {10, 2040.0}}];
+ifun[6] = Interpolation[{{5,  90.0}, {8, 240.0}, {10, 1470.0}}];
+
+difficulty[fam_] := ifun[fam][ndim]/scale[ndim]
+*)
+
+difficulty[fam_] := {6.0, 18.0, 2.2, 15.2, 16.1, 16.4}[[fam]]
+
+c[fam_] := Block[{r = w}, r difficulty[fam]/Plus@@ r]
+
+
+Install["Vegas"]
+
+Install["Suave"]
+
+Install["Divonne"]
+
+Install["Cuhre"]
+
+
+
+SetAll[opt__] := (
+  SetOptions[Vegas, opt];
+  SetOptions[Suave, opt];
+  SetOptions[Divonne, opt];
+  SetOptions[Cuhre, opt];
+  SetOptions[NIntegrate, opt];
+)
+
+SetAll[PrecisionGoal -> 3, MaxPoints -> maxpoints]
+
+SetOptions[Divonne, Key1 -> -200]
+
+
+def[f_][{x__}][{r__}] := (
+  Attributes[idef] = {HoldAll};
+  idef[int_, NIntegrate] := (
+    int := Module[{count = 0, res},
+             res = NIntegrate[f, r, EvaluationMonitor :> (++count)];
+             {count, res}]
+  ) /; $VersionNumber >= 5;
+  idef[int_, Int_] := (
+    int := Module[{count = 0, res},
+             res = Int[(++count; f), r];
+             {count, res}]
+  );
+  idef[vegas, Vegas];
+  idef[suave, Suave];
+  idef[divonne, Divonne];
+  idef[cuhre, Cuhre];
+  idef[nint, NIntegrate];
+)
+
+vars = Table[Unique["x"], {20}]
+
+test[ndim_, fam_] :=
+Block[ {w, xs = Take[vars, ndim]},
+  w := Table[Random[], {ndim}];
+  def[f[fam][xs, c[fam], w]][xs][{#, 0, 1}&/@ xs];
+  {vegas, suave, divonne, cuhre, nint}
+]
+
+
+dotest[ndim_, from_:1, to_:6] :=
+Block[ {dir = ToString[ndim]},
+  If[ FileType[dir] =!= Directory, CreateDirectory[dir] ];
+  Do[
+    SeedRandom[seed];
+    Put[ Table[test[ndim, fam], {repeat}],
+         ToFileName[dir, "fam" <> ToString[fam]] ],
+  {fam, from, to}]
+]
+