LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

upgrade of cuba to version 3.2. Merge in from BMW.

Browse Source
This commit is contained in:
suter_a 2013-12-20 09:58:23 +00:00
parent 1e6f1c45be
commit 2818739ec6
51 changed files with 2454 additions and 2216 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ChangeLog
View File

@ -60,6 +60,7 @@ FIXED 2012-09-23 fixed wrong chisq output in musrview if expected chisq is
present.
FIXED 2012-05-30 fixed RRF bug in single histo plotting.
FIXED 2012-05-18 fixed wrong forward/backward tag for ROOT-PPC (MUSR-215)
CHANGED 2013-12-20 upgrade of cuba to version 3.2. Merge in from BMW
CHANGED 2013-12-16 prettyfied the Noakes-Kalvius formulae. Furthermore added
a sub-folder with cross checks for these formulae.
CHANGED 2013-11-12 changed normalization of log max likelihood according to S.

26
configure.ac
View File

@ -69,8 +69,8 @@ PLUGIN_MINOR_VERSION=0
PLUGIN_MICRO_VERSION=0
#release versioning
CUBA_MAJOR_VERSION=2
CUBA_MINOR_VERSION=1
CUBA_MAJOR_VERSION=3
CUBA_MINOR_VERSION=2
CUBA_MICRO_VERSION=0
#API version
@ -200,9 +200,9 @@ dnl -----------------------------------------------
dnl Automake initialization and program checks
dnl -----------------------------------------------
AM_INIT_AUTOMAKE(@PACKAGE_NAME@, @PACKAGE_VERSION@)
m4_ifdef([AM_SILENT_RULES],
[AM_SILENT_RULES([yes])])
AM_INIT_AUTOMAKE
# m4_ifdef([AM_SILENT_RULES],
# [AM_SILENT_RULES([yes])])
AC_CONFIG_HEADER([config.h])
AC_LANG([C++])
AC_PROG_LN_S
@ -657,6 +657,17 @@ AC_ARG_ENABLE([BMWlibs], [AC_HELP_STRING([--enable-BMWlibs],[build optional BMW
AC_CHECK_FUNCS([powl])
AC_CHECK_FUNCS([erf])
# AC_FUNC_FORK
AC_FUNC_ALLOCA
AC_DEFUN([chk_shmget], [dnl
AC_REQUIRE([AC_CANONICAL_HOST])
AS_CASE([$host_os],
[*cygwin*], [],
[AC_CHECK_FUNCS([shmget])])
])
chk_shmget
AC_CHECK_FUNCS([getloadavg])
MAXDIM=${MAXDIM:-16}
AC_ARG_WITH(maxdim,
@ -1112,6 +1123,11 @@ AC_CONFIG_FILES([Makefile \
src/external/libCuba/Makefile \
src/external/libCuba/src/Makefile \
src/external/libCuba/src/cuba.pc \
src/external/libCuba/src/cuhre/Makefile \
src/external/libCuba/src/divonne/Makefile \
src/external/libCuba/src/suave/Makefile \
src/external/libCuba/src/vegas/Makefile \
src/external/libCuba/src/common/Makefile \
src/external/BMWtools/Makefile \
src/external/libFitPofB/Makefile \
src/external/libFitPofB/classes/Makefile \

17
src/external/BMWtools/BMWIntegrator.cpp vendored
View File

@ -34,6 +34,7 @@
#define USERDATA NULL
#define SEED 0
#define STATEFILE NULL
std::vector<double> TDWaveGapIntegralCuhre::fPar;
@ -60,7 +61,7 @@ double TDWaveGapIntegralCuhre::IntegrateFunc()
Cuhre(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
KEY, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -111,7 +112,7 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
Cuhre(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
KEY, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -162,7 +163,7 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
Cuhre(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
KEY, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -213,7 +214,7 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
Cuhre(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
KEY, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -272,7 +273,7 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
Divonne(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE, SEED, MINEVAL, MAXEVAL,
KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -324,7 +325,7 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc()
Suave(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, SEED, MINEVAL, MAXEVAL,
NNEW, FLATNESS,
NNEW, FLATNESS, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -375,7 +376,7 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
Cuhre(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
KEY, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
@ -426,7 +427,7 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
Cuhre(fNDim, NCOMP, Integrand, USERDATA,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
KEY, STATEFILE,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];

15
src/external/libCuba/src/Makefile.am vendored
View File

@ -1,25 +1,18 @@
## Process this file with automake to create Makefile.in
## $Id$
h_sources = cuba.h
c_sources = cuhre/Cuhre.c \
divonne/Divonne.c \
suave/Suave.c \
vegas/Vegas.c
SUBDIRS = cuhre divonne suave vegas common
include_HEADERS = cuba.h
AM_CPPFLAGS = -I. -Icommon
AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
CLEANFILES = common/*~ cuhre/*~ divonne/*~ suave/*~ vegas/*~ *~ core
lib_LTLIBRARIES = libcuba.la
libcuba_la_SOURCES = $(h_sources) $(c_sources)
libcuba_la_LIBADD = -lm
libcuba_la_SOURCES =
libcuba_la_LIBADD = cuhre/libcuhre.la divonne/libdivonne.la suave/libsuave.la vegas/libvegas.la -lm
libcuba_la_LDFLAGS = -version-info $(CUBA_LIBRARY_VERSION) -release $(CUBA_RELEASE) $(AM_LDFLAGS)
pkgconfigdir = $(libdir)/pkgconfig

67
src/external/libCuba/src/common/CSample.c vendored Normal file
View File

@ -0,0 +1,67 @@
/*
CSample.c
the serial sampling routine
for the C versions of the Cuba routines
by Thomas Hahn
last modified 19 Dec 11 th
*/
static inline number SampleRaw(cThis *t, number n, creal *x, real *f
VES_ONLY(, creal *w, ccount iter))
{
for( ; n; --n ) {
if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata
VES_ONLY(, w++, &iter)
DIV_ONLY(, &t->phase)) == ABORT ) return -1;
x += t->ndim;
f += t->ncomp;
}
return 0;
}
/*********************************************************************/
static inline void DoSampleSerial(This *t, cnumber n, creal *x, real *f
VES_ONLY(, creal *w, ccount iter))
{
t->neval += n;
if( SampleRaw(t, n, x, f VES_ONLY(, w, iter)) )
longjmp(t->abort, -99);
}
/*********************************************************************/
#ifdef HAVE_FORK
static void DoSample(This *t, number n, creal *x, real *f
VES_ONLY(, creal *w, ccount iter));
DIV_ONLY(static int Explore(This *t, cint iregion);)
#else
#define DoSample DoSampleSerial
#define Explore ExploreSerial
#define ForkCores(t)
#define WaitCores(t)
#endif
#ifdef DIVONNE
static inline count SampleExtra(This *t, cBounds *b)
{
number n = t->nextra;
t->peakfinder(&t->ndim, b, &n, t->xextra);
DoSample(t, n, t->xextra, t->fextra);
return n;
}
#endif
#include "common.c"
#ifdef HAVE_FORK
#include "Fork.c"
#endif
#include "Integrate.c"

20
src/external/libCuba/src/common/ChiSquare.c vendored
View File

@ -6,26 +6,6 @@
last modified 9 Feb 05 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#ifdef HAVE_ERF
#define Erf erf
#else

20
src/external/libCuba/src/common/Erf.c vendored
View File

@ -7,26 +7,6 @@
last modified 8 Feb 05 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
static real Erfc(creal x)
{

422
src/external/libCuba/src/common/Fork.c vendored Normal file
View File

@ -0,0 +1,422 @@
/*
Fork.c
the parallel sampling routine
for the C versions of the Cuba routines
by Thomas Hahn
last modified 25 Sep 13 th
*/
#define MINSLICE 10
#define MINCORES 1
/*#define MINCORES 2*/
typedef struct {
number n, m, i;
VES_ONLY(count iter;)
DIV_ONLY(int phase SHM_ONLY(, shmid);)
} Slice;
workerini cubaini;
#if defined HAVE_SHMGET && (defined SUAVE || defined DIVONNE)
#define FRAMECOPY
#endif
#ifdef DEBUG
#define TERM_RED "\e[31m"
#define TERM_BLUE "\e[34m"
#define TERM_RESET "\e[0m\n"
#define MASTER(s, ...) \
fprintf(stderr, TERM_RED ROUTINE " master %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__)
#define WORKER(s, ...) \
fprintf(stderr, TERM_BLUE ROUTINE " worker %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__)
#else
#define MASTER(s, ...)
#define WORKER(s, ...)
#endif
/*********************************************************************/
#ifndef MSG_WAITALL
/* Windows */
#define MSG_WAITALL 0
#endif
static inline int readsock(int fd, void *data, size_t n)
{
ssize_t got;
size_t remain = n;
do got = recv(fd, data, remain, MSG_WAITALL);
while( got > 0 && (data += got, remain -= got) > 0 );
return got;
}
static inline int writesock(int fd, const void *data, size_t n)
{
ssize_t got;
size_t remain = n;
do got = send(fd, data, remain, MSG_WAITALL);
while( got > 0 && (data += got, remain -= got) > 0 );
return got;
}
/*********************************************************************/
static void DoSample(This *t, number n, creal *x, real *f
VES_ONLY(, creal *w, ccount iter))
{
cint ncores = IMin(t->ncores, n/MINSLICE);
if( ncores < MINCORES ) DoSampleSerial(t, n, x, f VES_ONLY(, w, iter));
else {
Slice slice;
int core, abort;
number nx;
char s[128];
t->neval += n;
nx = n % ncores;
slice.m = slice.n = n/ncores + 1;
if( VERBOSE > 2 ) {
sprintf(s, "sampling " NUMBER " points each on %d cores",
slice.n, ncores);
Print(s);
}
slice.i = 0;
VES_ONLY(slice.iter = iter;)
DIV_ONLY(slice.phase = t->phase;)
#ifdef DIVONNE
if( n > t->nframe ) {
FrameFree(t, ShmRm(t));
t->nframe = n;
FrameAlloc(t);
}
SHM_ONLY(slice.shmid = t->shmid;)
#endif
SHM_ONLY(if( t->shmid != -1 ) {
slice.m = n;
#ifdef FRAMECOPY
VES_ONLY(Copy(t->frame, w, n);)
Copy(t->frame + n*NW, x, n*t->ndim);
#endif
})
for( core = 0; core < ncores; ++core ) {
cint fd = t->child[core];
slice.n -= (core == nx);
MASTER("sending samples (sli:%lu[+" VES_ONLY(NUMBER "w:%lu+")
NUMBER "x:%lu]) to fd %d",
sizeof slice, VES_ONLY(slice.n, sizeof *w,)
slice.n, t->ndim*sizeof *x, fd);
writesock(fd, &slice, sizeof slice);
SHM_ONLY(if( t->shmid == -1 )) {
VES_ONLY(writesock(fd, w, slice.n*sizeof *w);
w += slice.n;)
writesock(fd, x, slice.n*t->ndim*sizeof *x);
x += slice.n*t->ndim;
}
slice.i += slice.n;
}
abort = 0;
for( core = 0; core < ncores; ++core ) {
cint fd = t->child[core];
readsock(fd, &slice, sizeof slice);
MASTER("reading samples (sli:%lu[+" NUMBER "f:%lu]) from fd %d",
sizeof slice, slice.n, t->ncomp*sizeof *f, fd);
if( slice.n == -1 ) abort = 1;
else SHM_ONLY(if( t->shmid == -1 )) readsock(fd,
f + slice.i*t->ncomp, slice.n*t->ncomp*sizeof *f);
}
if( abort ) longjmp(t->abort, -99);
#ifdef FRAMECOPY
if( t->shmid != -1 )
Copy(f, t->frame + slice.m*(NW + t->ndim), slice.m*t->ncomp);
#endif
}
}
/*********************************************************************/
#ifdef DIVONNE
static inline int ReadyCore(cThis *t)
{
fd_set ready;
int core, n = 0;
FD_ZERO(&ready);
for( core = 0; core < t->ncores; ++core ) {
FD_SET(t->child[core], &ready);
n = IMax(n, t->child[core]);
}
select(n + 1, &ready, NULL, NULL, NULL);
for( core = 0; core < t->ncores; ++core )
if( FD_ISSET(t->child[core], &ready) ) break;
return core;
}
/*********************************************************************/
typedef struct {
number neval, neval_opt, neval_cut;
count nregions, iregion, retval;
} ExploreResult;
static int Explore(This *t, cint iregion)
{
csize_t regionsize = RegionSize;
Region *region;
int ireg = iregion, core = t->running;
Vector(Totals, totals, NCOMP);
if( t->ncores < MINCORES ) return ExploreSerial(t, iregion);
if( t->running >= ((iregion < 0) ? 1 : t->ncores) ) {
cint fd = t->child[core = ReadyCore(t)];
ExploreResult res;
count comp, succ;
--t->running;
MASTER("reading res + region (res:%lu+reg:%lu) from fd %d",
sizeof res, regionsize, fd);
readsock(fd, &res, sizeof res);
ireg = res.iregion;
region = RegionPtr(ireg);
succ = ireg + region->next;
readsock(fd, region, regionsize);
if( --res.nregions > 0 ) {
region->next = t->nregions - ireg;
EnlargeRegions(t, res.nregions);
MASTER("reading regions (%dreg:%lu) from fd %d",
res.nregions, regionsize, fd);
readsock(fd, RegionPtr(t->nregions), res.nregions*regionsize);
t->nregions += res.nregions;
RegionPtr(t->nregions-1)->next = succ - t->nregions + 1;
}
MASTER("reading totals (tot:%lu) from fd %d",
t->ncomp*sizeof(Totals), fd);
readsock(fd, totals, t->ncomp*sizeof(Totals));
for( comp = 0; comp < t->ncomp; ++comp )
t->totals[comp].secondspread =
Max(t->totals[comp].secondspread, totals[comp].secondspread);
t->neval += res.neval;
t->neval_opt += res.neval_opt;
t->neval_cut += res.neval_cut;
if( res.retval == -1 ) return -1;
}
if( iregion >= 0 ) {
Slice slice;
cint fd = t->child[core];
slice.n = 0;
slice.i = iregion;
slice.phase = t->phase;
region = RegionPtr(iregion);
MASTER("writing region (sli:%lu+sam:%lu+reg:%lu+tot:%lu) to fd %d",
sizeof slice, sizeof(Samples), regionsize,
t->ncomp*sizeof(Totals), fd);
writesock(fd, &slice, sizeof slice);
writesock(fd, &t->samples[region->isamples], sizeof(Samples));
writesock(fd, region, regionsize);
writesock(fd, t->totals, t->ncomp*sizeof(Totals));
region->depth = 0;
++t->running;
}
return ireg;
}
#endif
/*********************************************************************/
static void DoChild(This *t, cint core, cint fd)
{
Slice slice;
#ifdef DIVONNE
csize_t regionsize = RegionSize;
Vector(Totals, totals, NCOMP);
ExploreResult res;
t->totals = totals;
t->ncores = 0; /* no recursive forks */
t->size = 2*t->ndim + 2;
AllocRegions(t);
#endif
#ifdef SUAVE
SHM_ONLY(if( t->shmid == -1 ))
MemAlloc(t->frame, t->nframe*SAMPLESIZE);
#endif
if( cubaini.initfun ) cubaini.initfun(cubaini.initarg);
while( readsock(fd, &slice, sizeof slice) ) {
number n = slice.n;
DIV_ONLY(t->phase = slice.phase;)
if( n > 0 ) {
real VES_ONLY(*w,) *x, *f;
WORKER("reading samples (sli:%lu[+" VES_ONLY(NUMBER "w:%lu+")
NUMBER "x:%lu]) from fd %d",
sizeof slice, VES_ONLY(n, sizeof *w,) n, t->ndim*sizeof *x, fd);
#ifdef DIVONNE
if( slice.m > t->nframe ) {
FrameFree(t);
t->nframe = slice.m;
SHM_ONLY(t->shmid = slice.shmid; ShmMap(t) else)
MemAlloc(t->frame, t->nframe*SAMPLESIZE);
}
#endif
VES_ONLY(w = t->frame;)
x = t->frame + slice.m*NW;
f = x + slice.m*t->ndim;
SHM_ONLY(if( t->shmid != -1 ) {
VES_ONLY(w += slice.i;)
x += slice.i*t->ndim;
f += slice.i*t->ncomp;
}
else) {
VES_ONLY(readsock(fd, w, n*sizeof *w);)
readsock(fd, x, n*t->ndim*sizeof *x);
}
slice.n |= SampleRaw(t, n, x, f VES_ONLY(, w, slice.iter));
WORKER("writing samples (sli:%lu[+" NUMBER "f:%lu]) to fd %d",
sizeof slice, slice.n, t->ncomp*sizeof *f, fd);
writesock(fd, &slice, sizeof slice);
if( SHM_ONLY(t->shmid == -1 &&) slice.n != -1 )
writesock(fd, f, slice.n*t->ncomp*sizeof *f);
}
#ifdef DIVONNE
else {
Samples *samples, psamples;
WORKER("reading region (sli:%lu+sam:%lu+reg:%lu+tot:%lu) from fd %d",
sizeof slice, sizeof(Samples), regionsize,
t->ncomp*sizeof(Totals), fd);
readsock(fd, &psamples, sizeof(Samples));
readsock(fd, t->region, regionsize);
readsock(fd, totals, t->ncomp*sizeof(Totals));
t->nregions = 1;
t->neval = t->neval_opt = t->neval_cut = 0;
samples = &t->samples[RegionPtr(0)->isamples];
if( psamples.n != samples->n ) {
SamplesFree(samples);
*samples = psamples;
SamplesAlloc(t, samples);
}
res.retval = ExploreSerial(t, 0);
res.neval = t->neval;
res.neval_opt = t->neval_opt;
res.neval_cut = t->neval_cut;
res.nregions = t->nregions;
res.iregion = slice.i;
WORKER("writing regions (res:%lu+%dreg:%lu+tot:%lu) to fd %d",
sizeof res, t->nregions, regionsize,
t->ncomp*sizeof(Totals), fd);
writesock(fd, &res, sizeof res);
writesock(fd, t->region, t->nregions*regionsize);
writesock(fd, totals, t->ncomp*sizeof(Totals));
}
#endif
}
WORKER("wrapping up");
if( cubaini.exitfun ) cubaini.exitfun(cubaini.exitarg);
exit(0);
}
/*********************************************************************/
#ifdef HAVE_GETLOADAVG
double cubaloadavg_;
#endif
static inline void ForkCores(This *t)
{
int core;
char s[128];
cchar *env = getenv("CUBACORES");
t->ncores = env ? atoi(env) : sysconf(_SC_NPROCESSORS_ONLN);
#ifdef HAVE_GETLOADAVG
if( env == NULL || t->ncores < 0 ) {
if( cubaloadavg_ < 0 ) getloadavg(&cubaloadavg_, 1);
t->ncores = abs(t->ncores) - floor(cubaloadavg_);
}
#endif
DIV_ONLY(t->running = 0;)
if( t->ncores < MINCORES ) return;
if( VERBOSE ) {
sprintf(s, "using %d cores via "
#ifdef HAVE_SHMGET
"shared memory",
#else
"pipes",
#endif
t->ncores);
Print(s);
}
fflush(NULL); /* make sure all buffers are flushed,
or else buffered content will be written
out multiply, at each child's exit(0) */
Alloc(t->child, t->ncores);
for( core = 0; core < t->ncores; ++core ) {
int fd[2];
pid_t pid;
assert(
socketpair(AF_LOCAL, SOCK_STREAM, 0, fd) != -1 &&
(pid = fork()) != -1 );
if( pid == 0 ) {
close(fd[0]);
DoChild(t, core, fd[1]);
}
MASTER("forked pid %d pipe %d(master) -> %d(worker)",
pid, fd[0], fd[1]);
close(fd[1]);
t->child[core] = fd[0];
}
}
/*********************************************************************/
static inline void WaitCores(cThis *t)
{
if( t->ncores >= MINCORES ) {
int core;
pid_t pid;
for( core = 0; core < t->ncores; ++core ) {
MASTER("closing fd %d", t->child[core]);
close(t->child[core]);
}
free(t->child);
for( core = 0; core < t->ncores; ++core ) {
MASTER("waiting for child");
wait(&pid);
MASTER("pid %d terminated", pid);
}
}
}

90
src/external/libCuba/src/common/MSample.c vendored Normal file
View File

@ -0,0 +1,90 @@
/*
MSample.c
the sampling routine for the
Mathematica versions of the Cuba routines
by Thomas Hahn
last modified 19 Mar 12 th
*/
static void DoSample(This *t, cnumber n, real *x, real *f
VES_ONLY(, real *w, ccount iter))
{
real *mma_f;
long mma_n;
if( MLAbort ) longjmp(t->abort, -99);
MLPutFunction(stdlink, "EvaluatePacket", 1);
MLPutFunction(stdlink, "Cuba`" ROUTINE "`sample", 1 VES_ONLY(+2) DIV_ONLY(+1));
MLPutRealList(stdlink, x, n*t->ndim);
VES_ONLY(MLPutRealList(stdlink, w, n);
MLPutInteger(stdlink, iter);)
DIV_ONLY(MLPutInteger(stdlink, t->phase);)
MLEndPacket(stdlink);
MLNextPacket(stdlink);
if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
MLClearError(stdlink);
MLNewPacket(stdlink);
longjmp(t->abort, -99);
}
t->neval += mma_n;
if( mma_n != n*t->ncomp ) {
MLDisownRealList(stdlink, mma_f, mma_n);
longjmp(t->abort, -3);
}
Copy(f, mma_f, n*t->ncomp);
MLDisownRealList(stdlink, mma_f, mma_n);
}
/*********************************************************************/
#ifdef DIVONNE
#define Explore ExploreSerial
static count SampleExtra(This *t, cBounds *b)
{
count n, nget;
real *mma_f;
long mma_n;
MLPutFunction(stdlink, "EvaluatePacket", 1);
MLPutFunction(stdlink, "Cuba`Divonne`findpeak", 2);
MLPutRealList(stdlink, (real *)b, 2*t->ndim);
MLPutInteger(stdlink, t->phase);
MLEndPacket(stdlink);
MLNextPacket(stdlink);
if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
MLClearError(stdlink);
MLNewPacket(stdlink);
longjmp(t->abort, -99);
}
t->neval += nget = mma_n/(t->ndim + t->ncomp);
n = IMin(nget, t->nextra);
if( n ) {
Copy(t->xextra, mma_f, n*t->ndim);
Copy(t->fextra, mma_f + nget*t->ndim, n*t->ncomp);
}
MLDisownRealList(stdlink, mma_f, mma_n);
return n;
}
#endif
/*********************************************************************/
#include "common.c"
#define ForkCores(t)
#define WaitCores(t)
#include "Integrate.c"

13
src/external/libCuba/src/common/Makefile.am vendored Normal file
View File

@ -0,0 +1,13 @@
## Process this file with automake to create Makefile.in
## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $
c_sources = WorkerIni.c
AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
noinst_LTLIBRARIES = libworkerini.la
libworkerini_la_SOURCES = $(c_sources)
libworkerini_la_LDFLAGS = $(AM_LDFLAGS)

34
src/external/libCuba/src/common/Random.c vendored
View File

@ -1,29 +1,9 @@
/*
Random.c
quasi- and pseudo-random-number generation
last modified 8 Jun 10 th
last modified 7 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
/*
PART 1: Sobol quasi-random-number generator
@ -142,7 +122,7 @@ static inline void SobolIni(This *t)
}
t->rng.sobol.seq = 0;
VecClear(t->rng.sobol.prev);
XClear(t->rng.sobol.prev);
t->rng.getrandom = SobolGet;
t->rng.skiprandom = SobolSkip;
@ -300,10 +280,10 @@ static void RanluxGet(This *t, real *x)
cint nskip = (--t->rng.ranlux.n24 >= 0) ? 0 :
(t->rng.ranlux.n24 = 24, t->rng.ranlux.nskip);
cint s = RanluxInt(t, 1 + nskip);
x[dim] = s*0x1p-24;
x[dim] = ldexp(s, -24);
/* small numbers (with less than 12 significant bits) are "padded" */
if( s < (1 << 12) )
x[dim] += t->rng.ranlux.state[t->rng.ranlux.j24]*0x1p-48;
x[dim] += ldexp(t->rng.ranlux.state[t->rng.ranlux.j24], -48);
}
}
@ -320,11 +300,11 @@ static inline void RanluxIni(This *t)
cint skip[] = {24, 48, 97, 223, 389,
223, 223, 223, 223, 223, 223, 223, 223, 223, 223,
223, 223, 223, 223, 223, 223, 223, 223, 223, 223};
state_t seed = t->seed;
state_t level = RNG;
int seed = t->seed;
int level = RNG;
count i;
if( level < sizeof skip ) level = skip[level];
if( level < Elements(skip) ) level = skip[level];
t->rng.ranlux.nskip = level - 24;
t->rng.ranlux.i24 = 23;

37
src/external/libCuba/src/common/WorkerIni.c vendored Normal file
View File

@ -0,0 +1,37 @@
/*
WorkerIni.c
set/run the init/exit functions for worker processes
by Thomas Hahn
last modified 6 Sep 12 th
*/
#include "stddecl.h"
extern workerini cubaini;
Extern void SUFFIX(cubasetinit)(subroutine f, void *arg)
{
cubaini.initfun = f;
cubaini.initarg = arg;
}
Extern void SUFFIX(cubasetexit)(subroutine f, void *arg)
{
cubaini.exitfun = f;
cubaini.exitarg = arg;
}
Extern void SUFFIX(cubaruninit)()
{
if( cubaini.initfun ) cubaini.initfun(cubaini.initarg);
}
Extern void SUFFIX(cubarunexit)()
{
if( cubaini.exitfun ) cubaini.exitfun(cubaini.exitarg);
}

105
src/external/libCuba/src/common/debug.c vendored
View File

@ -1,105 +0,0 @@
/***************************************************************************
* Copyright (C) 2004-2009 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define DEBFile stdout
#define DEBF "%9.5f"
//#define DEBF "%21.17f"
//#define DEBF "%a"
#ifdef DEBFile
#define DEBOpen()
#define DEBClose()
#else
FILE *DEBFile;
static inline void DEBOpen()
{
#ifdef MLVERSION
DEBFile = fopen("log-mma", "w");
#else
DEBFile = fopen("log-c", "w");
#endif
}
static inline void DEBClose()
{
fclose(DEBFile);
}
#endif
#define DEB(...) fprintf(DEBFile, __VA_ARGS__); fflush(DEBFile)
static inline void DEBVec(const char *s, creal *d)
{
char space[strlen(s) + 2];
count dim;
memset(space, ' ', sizeof(space));
space[sizeof(space) - 1] = 0;
DEB("%s=" DEBF "\n", s, d[0]);
for( dim = 1; dim < ndim_; ++dim )
DEB("%s" DEBF "\n", space, d[dim]);
}
/*
static inline void DEBRegion(const char *s, cBounds *b)
{
char space[strlen(s) + 3];
count dim;
memset(space, ' ', sizeof(space));
space[sizeof(space) - 1] = 0;
DEB("%s: " DEBF " - " DEBF "\n", s, b[0].lower, b[0].upper);
for( dim = 1; dim < ndim_; ++dim )
DEB("%s" DEBF " - " DEBF "\n", space, b[dim].lower, b[dim].upper);
}
*/
static inline void DEBMem(const char *s)
{
int kbytes = -1;
FILE *f;
char procfile[128];
sprintf(procfile, "/proc/%d/status", getpid());
f = fopen(procfile, "r");
while( !feof(f) ) {
char s[128];
*s = 0;
fgets(s, sizeof(s), f);
if( sscanf(s, "VmSize: %d", &kbytes) == 1 ) break;
}
fclose(f);
DEB("MEM %s: %d kbytes\n", s, kbytes);
}

319
src/external/libCuba/src/common/stddecl.h vendored
View File

@ -1,37 +1,20 @@
/*
stddecl.h
Type declarations common to all Cuba routines
last modified 16 Jun 10 th
declarations common to all Cuba routines
last modified 17 Sep 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#ifndef __stddecl_h__
#define __stddecl_h__
#ifndef _stddecl_h_
#define _stddecl_h_
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#define _BSD_SOURCE
#define _XOPEN_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
@ -39,28 +22,83 @@
#include <float.h>
#include <limits.h>
#include <unistd.h>
#include <assert.h>
#include <fcntl.h>
#include <setjmp.h>
#include <sys/stat.h>
#include <sys/types.h>
#ifdef HAVE_FORK
#include <sys/wait.h>
#include <sys/socket.h>
#ifdef HAVE_SHMGET
#include <sys/ipc.h>
#include <sys/shm.h>
#endif
#endif
#ifdef HAVE_ALLOCA_H
#include <alloca.h>
#elif defined __GNUC__
#define alloca __builtin_alloca
#elif defined _AIX
#define alloca __alloca
#elif defined _MSC_VER
#include <malloc.h>
#define alloca _alloca
#else
#include <stddef.h>
#ifdef __cplusplus
extern "C"
#endif
void *alloca (size_t);
#endif
#ifndef NDIM
#define NDIM t->ndim
#endif
#ifndef NCOMP
#define NCOMP t->ncomp
#define MAXDIM 1024
#else
#define MAXDIM NDIM
#endif
#ifndef NCOMP
#define NCOMP t->ncomp
#define MAXCOMP 1024
#else
#define MAXCOMP NCOMP
#endif
#if defined(VEGAS) || defined(SUAVE)
#define VES_ONLY(...) __VA_ARGS__
#define NW 1
#else
#define VES_ONLY(...)
#define NW 0
#endif
#ifdef DIVONNE
#define DIV_ONLY(...) __VA_ARGS__
#else
#define DIV_ONLY(...)
#endif
#define SAMPLESIZE (NW + t->ndim + t->ncomp)*sizeof(real)
#define VERBOSE (t->flags & 3)
#define LAST (t->flags & 4)
#define SHARPEDGES (t->flags & 8)
#define KEEPFILE (t->flags & 16)
#define REGIONS (t->flags & 128)
#define RNG (t->flags >> 8)
#define INFTY DBL_MAX
#define NOTZERO 0x1p-104
#if __STDC_VERSION__ >= 199901L
#define POW2(n) 0x1p-##n
#else
#define POW2(n) ldexp(1., -n)
#endif
#define NOTZERO POW2(104)
#define ABORT -999
@ -68,15 +106,17 @@
#define Copy(d, s, n) memcpy(d, s, (n)*sizeof(*(d)))
#define VecCopy(d, s) Copy(d, s, t->ndim)
#define Move(d, s, n) memmove(d, s, (n)*sizeof(*(d)))
#define ResCopy(d, s) Copy(d, s, t->ncomp)
#define XCopy(d, s) Copy(d, s, t->ndim)
#define FCopy(d, s) Copy(d, s, t->ncomp)
#define Clear(d, n) memset(d, 0, (n)*sizeof(*(d)))
#define VecClear(d) Clear(d, t->ndim)
#define XClear(d) Clear(d, t->ndim)
#define ResClear(d) Clear(d, t->ncomp)
#define FClear(d) Clear(d, t->ncomp)
#define Zap(d) memset(d, 0, sizeof(d))
@ -90,14 +130,155 @@
#define reallocset(p, n) (p = realloc(p, n))
#endif
#define ChkAlloc(r) if( r == NULL ) { \
fprintf(stderr, "Out of memory in " __FILE__ " line %d.\n", __LINE__); \
exit(1); \
#define Abort(s) abort1(s, __LINE__)
#define abort1(s, line) abort2(s, line)
#define abort2(s, line) { perror(s " " __FILE__ "(" #line ")"); exit(1); }
#define Die(p) if( (p) == NULL ) Abort("malloc")
#define MemAlloc(p, n) Die(mallocset(p, n))
#define ReAlloc(p, n) Die(reallocset(p, n))
#define Alloc(p, n) MemAlloc(p, (n)*sizeof(*p))
#if __STDC_VERSION__ >= 199901L
#define Sized(type, var, size) char var##_[size]; type *var = (type *)var##_
#define Vector(type, var, n1) type var[n1]
#define Array(type, var, n1, n2) type var[n1][n2]
#else
#define Sized(type, var, size) type *var = alloca(size)
#define Vector(type, var, n1) type *var = alloca((n1)*sizeof(type))
#define Array(type, var, n1, n2) type (*var)[n2] = alloca((n1)*(n2)*sizeof(type))
#endif
#define FORK_ONLY(...)
#define SHM_ONLY(...)
#define ShmAlloc(...)
#define ShmFree(...)
#ifdef MLVERSION
#define ML_ONLY(...) __VA_ARGS__
#else
#define ML_ONLY(...)
#ifdef HAVE_FORK
#undef FORK_ONLY
#define FORK_ONLY(...) __VA_ARGS__
#ifdef HAVE_SHMGET
#undef SHM_ONLY
#define SHM_ONLY(...) __VA_ARGS__
#define ShmMap(t, ...) if( t->shmid != -1 ) { \
t->frame = shmat(t->shmid, NULL, 0); \
if( t->frame == (void *)-1 ) Abort("shmat"); \
__VA_ARGS__ \
}
#define Alloc(p, n) MemAlloc(p, (n)*sizeof(*p))
#define MemAlloc(p, n) ChkAlloc(mallocset(p, n))
#define ReAlloc(p, n) ChkAlloc(reallocset(p, n))
#define ShmRm(t) shmctl(t->shmid, IPC_RMID, NULL);
#undef ShmAlloc
#define ShmAlloc(t, ...) \
t->shmid = shmget(IPC_PRIVATE, t->nframe*SAMPLESIZE, IPC_CREAT | 0600); \
ShmMap(t, __VA_ARGS__)
#undef ShmFree
#define ShmFree(t, ...) if( t->shmid != -1 ) { \
shmdt(t->frame); \
__VA_ARGS__ \
}
#endif
#endif
#endif
#define FrameAlloc(t, ...) \
SHM_ONLY(ShmAlloc(t, __VA_ARGS__) else) \
MemAlloc(t->frame, t->nframe*SAMPLESIZE);
#define FrameFree(t, ...) DIV_ONLY(if( t->nframe )) { \
SHM_ONLY(ShmFree(t, __VA_ARGS__) else) \
free(t->frame); \
}
#define StateDecl \
char *statefile_tmp = NULL, *statefile_XXXXXX = NULL; \
int statemsg = VERBOSE; \
ssize_t ini = 1; \
struct stat st
#define StateSetup(t) if( (t)->statefile ) { \
if( *(t)->statefile == 0 ) (t)->statefile = NULL; \
else { \
ccount len = strlen((t)->statefile); \
statefile_tmp = alloca(len + 8); \
strcpy(statefile_tmp, (t)->statefile); \
statefile_XXXXXX = statefile_tmp + len; \
} \
}
typedef long long int signature_t;
#define StateSignature(t, i) (0x41425543 + \
((signature_t)(i) << 60) + \
((signature_t)(t)->ncomp << 48) + \
((signature_t)(t)->ndim << 32))
#define StateReadTest(t) (t)->statefile && \
stat((t)->statefile, &st) == 0 && (st.st_mode & 0400)
#define StateReadOpen(t, fd) do { \
int fd; \
if( (fd = open((t)->statefile, O_RDONLY)) != -1 ) { \
do
#define StateRead(fd, buf, size) \
ini += size - read(fd, buf, size)
#define StateReadClose(t, fd) \
while( (--ini, 0) ); \
close(fd); \
} \
if( ini | statemsg ) { \
char s[512]; \
sprintf(s, ini ? \
"\nError restoring state from %s, starting from scratch." : \
"\nRestored state from %s.", (t)->statefile); \
Print(s); \
} \
} while( 0 )
#define StateWriteTest(t) ((t)->statefile)
#define StateWriteOpen(t, fd) do { \
ssize_t fail = 1; \
int fd; \
strcpy(statefile_XXXXXX, "-XXXXXX"); \
if( (fd = mkstemp(statefile_tmp)) != -1 ) { \
do
#define StateWrite(fd, buf, size) \
fail += size - write(fd, buf, size)
#define StateWriteClose(t, fd) \
while( (--fail, 0) ); \
close(fd); \
if( fail == 0 ) fail |= rename(statefile_tmp, (t)->statefile); \
} \
if( fail | statemsg ) { \
char s[512]; \
sprintf(s, fail ? \
"\nError saving state to %s." : \
"\nSaved state to %s.", (t)->statefile); \
Print(s); \
statemsg &= fail & -2; \
} \
} while( 0 )
#define StateRemove(t) \
if( fail == 0 && (t)->statefile && KEEPFILE == 0 ) unlink((t)->statefile)
#ifdef __cplusplus
@ -115,6 +296,8 @@ typedef const int cint;
typedef const long clong;
typedef const size_t csize_t;
#define COUNT "%d"
typedef /*unsigned*/ int count;
typedef const count ccount;
@ -144,6 +327,15 @@ typedef /*long*/ double real;
typedef const real creal;
typedef void (*subroutine)();
typedef struct {
subroutine initfun;
void *initarg;
subroutine exitfun;
void *exitarg;
} workerini;
struct _this;
@ -179,33 +371,42 @@ typedef struct {
} RNGState;
#ifdef UNDERSCORE
#define SUFFIX(s) s##_
#else
#if NOUNDERSCORE
#define SUFFIX(s) s
#else
#define SUFFIX(s) s##_
#endif
#define EXPORT(s) EXPORT_(PREFIX(s))
#define EXPORT_(s) SUFFIX(s)
static inline real Sq(creal x)
{
#define CString(cs, fs, len) { \
char *_s = NULL; \
if( fs ) { \
int _l = len; \
while( _l > 0 && fs[_l - 1] == ' ' ) --_l; \
if( _l > 0 && (_s = alloca(_l + 1)) ) { \
memcpy(_s, fs, _l); \
_s[_l] = 0; \
} \
} \
cs = _s; \
}
static inline real Sq(creal x) {
return x*x;
}
static inline real Min(creal a, creal b)
{
static inline real Min(creal a, creal b) {
return (a < b) ? a : b;
}
static inline real Max(creal a, creal b)
{
static inline real Max(creal a, creal b) {
return (a > b) ? a : b;
}
static inline real Weight(creal sum, creal sqsum, cnumber n)
{
static inline real Weight(creal sum, creal sqsum, cnumber n) {
creal w = sqrt(sqsum*n);
return (n - 1)/Max((w + sum)*(w - sum), NOTZERO);
}
@ -235,5 +436,25 @@ static inline real Weight(creal sum, creal sqsum, cnumber n)
/* abs(a) + (a == 0) */
#define Abs1(a) (((a) ^ NegQ(a)) - NegQ((a) - 1))
#ifdef MLVERSION
static inline void Print(MLCONST char *s)
{
MLPutFunction(stdlink, "EvaluatePacket", 1);
MLPutFunction(stdlink, "Print", 1);
MLPutString(stdlink, s);
MLEndPacket(stdlink);
MLNextPacket(stdlink);
MLNewPacket(stdlink);
}
#else
#define Print(s) puts(s); fflush(stdout)
#endif
#endif

33
src/external/libCuba/src/cuba.h vendored
View File

@ -2,29 +2,9 @@
cuba.h
Prototypes for the Cuba library
this file is part of Cuba
last modified 16 Jun 10 th
last modified 30 Apr 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#ifdef __cplusplus
extern "C" {
#endif
@ -68,6 +48,7 @@ void Suave(const int ndim, const int ncomp,
const int flags, const int seed,
const int mineval, const int maxeval,
const int nnew, const double flatness,
const char *statefile,
int *nregions, int *neval, int *fail,
double integral[], double error[], double prob[]);
@ -77,6 +58,7 @@ void llSuave(const int ndim, const int ncomp,
const int flags, const int seed,
const long long int mineval, const long long int maxeval,
const long long int nnew, const double flatness,
const char *statefile,
int *nregions, long long int *neval, int *fail,
double integral[], double error[], double prob[]);
@ -89,6 +71,7 @@ void Divonne(const int ndim, const int ncomp,
const double border, const double maxchisq, const double mindeviation,
const int ngiven, const int ldxgiven, double xgiven[],
const int nextra, peakfinder_t peakfinder,
const char *statefile,
int *nregions, int *neval, int *fail,
double integral[], double error[], double prob[]);
@ -102,6 +85,7 @@ void llDivonne(const int ndim, const int ncomp,
const long long int ngiven, const int ldxgiven, double xgiven[],
const long long int nextra,
void (*peakfinder)(const int *, const double [], int *, double []),
const char *statefile,
int *nregions, long long int *neval, int *fail,
double integral[], double error[], double prob[]);
@ -110,6 +94,7 @@ void Cuhre(const int ndim, const int ncomp,
const double epsrel, const double epsabs,
const int flags, const int mineval, const int maxeval,
const int key,
const char *statefile,
int *nregions, int *neval, int *fail,
double integral[], double error[], double prob[]);
@ -119,9 +104,15 @@ void llCuhre(const int ndim, const int ncomp,
const int flags,
const long long int mineval, const long long int maxeval,
const int key,
const char *statefile,
int *nregions, long long int *neval, int *fail,
double integral[], double error[], double prob[]);
void cubasetinit(void (*)(), void *);
void cubasetexit(void (*)(), void *);
void cubaruninit(void);
void cubaruninit(void);
#ifdef __cplusplus
}
#endif

54
src/external/libCuba/src/cuhre/Cuhre.c vendored
View File

@ -2,56 +2,23 @@
Cuhre.c
Adaptive integration using cubature rules
by Thomas Hahn
last modified 16 Jun 10 th
last modified 17 Sep 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define CUHRE
#define ROUTINE "Cuhre"
#include "decl.h"
#define Print(s) puts(s); fflush(stdout)
#include "CSample.c"
/*********************************************************************/
static inline void DoSample(This *t, count n, creal *x, real *f)
{
t->neval += n;
while( n-- ) {
if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata) == ABORT )
longjmp(t->abort, 1);
x += t->ndim;
f += t->ncomp;
}
}
/*********************************************************************/
#include "common.c"
Extern void EXPORT(Cuhre)(ccount ndim, ccount ncomp,
Integrand integrand, void *userdata,
creal epsrel, creal epsabs,
cint flags, cnumber mineval, cnumber maxeval,
ccount key,
ccount key, cchar *statefile,
count *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
{
@ -66,8 +33,7 @@ Extern void EXPORT(Cuhre)(ccount ndim, ccount ncomp,
t.mineval = mineval;
t.maxeval = maxeval;
t.key = key;
t.nregions = 0;
t.neval = 0;
t.statefile = statefile;
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
@ -80,9 +46,9 @@ Extern void EXPORT(cuhre)(ccount *pndim, ccount *pncomp,
Integrand integrand, void *userdata,
creal *pepsrel, creal *pepsabs,
cint *pflags, cnumber *pmineval, cnumber *pmaxeval,
ccount *pkey,
ccount *pkey, cchar *statefile,
count *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
real *integral, real *error, real *prob, cint statefilelen)
{
This t;
t.ndim = *pndim;
@ -95,11 +61,9 @@ Extern void EXPORT(cuhre)(ccount *pndim, ccount *pncomp,
t.mineval = *pmineval;
t.maxeval = *pmaxeval;
t.key = *pkey;
t.nregions = 0;
t.neval = 0;
CString(t.statefile, statefile, statefilelen);
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
*pneval = t.neval;
}

251
src/external/libCuba/src/cuhre/Integrate.c vendored
View File

@ -2,158 +2,174 @@
Integrate.c
integrate over the unit hypercube
this file is part of Cuhre
last modified 7 Jun 10 th
checkpointing by B. Chokoufe
last modified 17 Sep 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define POOLSIZE 1024
typedef struct pool {
struct pool *next;
char region[];
} Pool;
typedef struct {
signature_t signature;
count nregions, ncur;
number neval;
Totals totals[];
} State;
static int Integrate(This *t, real *integral, real *error, real *prob)
{
TYPEDEFREGION;
typedef struct pool {
struct pool *next;
Region region[POOLSIZE];
} Pool;
StateDecl;
csize_t statesize = sizeof(State) + NCOMP*sizeof(Totals);
Sized(State, state, statesize);
csize_t regionsize = RegionSize;
csize_t poolsize = sizeof(Pool) + POOLSIZE*regionsize;
Vector(Result, result, NCOMP);
Vector(char, out, 128*NCOMP + 256);
count dim, comp, ncur, ipool, npool;
int fail = -99;
Totals totals[NCOMP];
Totals *tot, *Tot = state->totals + t->ncomp;
Result *res, *resL, *resR;
Bounds *b, *B;
Pool *cur = NULL, *pool;
Region *region;
count comp, ipool, npool;
int fail;
if( VERBOSE > 1 ) {
char s[256];
sprintf(s, "Cuhre input parameters:\n"
sprintf(out, "Cuhre input parameters:\n"
" ndim " COUNT "\n ncomp " COUNT "\n"
" epsrel " REAL "\n epsabs " REAL "\n"
" flags %d\n mineval " NUMBER "\n maxeval " NUMBER "\n"
" key " COUNT,
" key " COUNT "\n"
" statefile \"%s\"",
t->ndim, t->ncomp,
t->epsrel, t->epsabs,
t->flags, t->mineval, t->maxeval,
t->key);
Print(s);
t->key,
t->statefile);
Print(out);
}
if( BadComponent(t) ) return -2;
if( BadDimension(t) ) return -1;
RuleAlloc(t);
t->epsabs = Max(t->epsabs, NOTZERO);
RuleAlloc(t);
t->mineval = IMax(t->mineval, t->rule.n + 1);
FrameAlloc(t, ShmRm(t));
ForkCores(t);
if( setjmp(t->abort) ) goto abort;
if( (fail = setjmp(t->abort)) ) goto abort;
Alloc(cur, 1);
cur->next = NULL;
ncur = 1;
StateSetup(t);
region = cur->region;
region->div = 0;
for( dim = 0; dim < t->ndim; ++dim ) {
Bounds *b = &region->bounds[dim];
b->lower = 0;
b->upper = 1;
if( StateReadTest(t) ) {
StateReadOpen(t, fd) {
Pool *prev = NULL;
int size;
if( read(fd, state, statesize) != statesize ||
state->signature != StateSignature(t, 4) ) break;
t->neval = state->neval;
t->nregions = state->nregions;
do {
MemAlloc(cur, poolsize);
cur->next = prev;
prev = cur;
size = read(fd, cur, poolsize);
} while( size == poolsize );
if( size != state->ncur*regionsize ) break;
} StateReadClose(t, fd);
}
Sample(t, region);
if( ini ) {
MemAlloc(cur, poolsize);
cur->next = NULL;
state->ncur = t->nregions = 1;
for( comp = 0; comp < t->ncomp; ++comp ) {
Totals *tot = &totals[comp];
Result *r = &region->result[comp];
tot->avg = tot->lastavg = tot->guess = r->avg;
tot->err = tot->lasterr = r->err;
tot->weightsum = 1/Max(Sq(r->err), NOTZERO);
tot->avgsum = tot->weightsum*r->avg;
tot->chisq = tot->chisqsum = tot->chisum = 0;
region = (Region *)cur->region;
region->div = 0;
for( B = (b = region->bounds) + t->ndim; b < B; ++b ) {
b->lower = 0;
b->upper = 1;
}
t->neval = 0;
Sample(t, region);
for( res = RegionResult(region), tot = state->totals;
tot < Tot; ++res, ++tot ) {
tot->avg = tot->lastavg = tot->guess = res->avg;
tot->err = tot->lasterr = res->err;
tot->weightsum = 1/Max(Sq(res->err), NOTZERO);
tot->avgsum = tot->weightsum*res->avg;
tot->chisq = tot->chisqsum = tot->chisum = 0;
}
}
for( t->nregions = 1; ; ++t->nregions ) {
/* main iteration loop */
for( ; ; ) {
count maxcomp, bisectdim;
real maxratio, maxerr;
Result result[NCOMP];
Region *regionL, *regionR;
Bounds *bL, *bR;
if( VERBOSE ) {
char s[128 + 128*NCOMP], *p = s;
p += sprintf(p, "\n"
char *oe = out + sprintf(out, "\n"
"Iteration " COUNT ": " NUMBER " integrand evaluations so far",
t->nregions, t->neval);
for( comp = 0; comp < t->ncomp; ++comp ) {
cTotals *tot = &totals[comp];
p += sprintf(p, "\n[" COUNT "] "
for( tot = state->totals, comp = 0; tot < Tot; ++tot )
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)",
comp + 1, tot->avg, tot->err, tot->chisq, t->nregions - 1);
}
Print(s);
++comp, tot->avg, tot->err, tot->chisq, t->nregions - 1);
Print(out);
}
maxratio = -INFTY;
maxcomp = 0;
for( comp = 0; comp < t->ncomp; ++comp ) {
creal ratio = totals[comp].err/MaxErr(totals[comp].avg);
for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
creal ratio = tot->err/MaxErr(tot->avg);
if( ratio > maxratio ) {
maxratio = ratio;
maxcomp = comp;
}
}
if( maxratio <= 1 && t->neval >= t->mineval ) {
fail = 0;
if( maxratio <= 1 && t->neval >= t->mineval ) break;
if( t->neval >= t->maxeval ) {
fail = 1;
break;
}
if( t->neval >= t->maxeval ) break;
maxerr = -INFTY;
regionL = cur->region;
npool = ncur;
regionL = (Region *)cur->region;
npool = state->ncur;
for( pool = cur; pool; npool = POOLSIZE, pool = pool->next )
for( ipool = 0; ipool < npool; ++ipool ) {
Region *region = &pool->region[ipool];
creal err = region->result[maxcomp].err;
Region *region = RegionPtr(pool, ipool);
creal err = RegionResult(region)[maxcomp].err;
if( err > maxerr ) {
maxerr = err;
regionL = region;
}
}
if( ncur == POOLSIZE ) {
if( state->ncur == POOLSIZE ) {
Pool *prev = cur;
Alloc(cur, 1);
MemAlloc(cur, poolsize);
cur->next = prev;
ncur = 0;
state->ncur = 0;
}
regionR = &cur->region[ncur++];
regionR = RegionPtr(cur, state->ncur++);
regionR->div = ++regionL->div;
ResCopy(result, regionL->result);
VecCopy(regionR->bounds, regionL->bounds);
FCopy(result, RegionResult(regionL));
XCopy(regionR->bounds, regionL->bounds);
bisectdim = result[maxcomp].bisectdim;
bL = &regionL->bounds[bisectdim];
@ -163,25 +179,25 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
Sample(t, regionL);
Sample(t, regionR);
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &result[comp];
Result *rL = &regionL->result[comp];
Result *rR = &regionR->result[comp];
Totals *tot = &totals[comp];
for( res = result,
resL = RegionResult(regionL),
resR = RegionResult(regionR),
tot = state->totals;
tot < Tot; ++res, ++resL, ++resR, ++tot ) {
real diff, err, w, avg, sigsq;
tot->lastavg += diff = rL->avg + rR->avg - r->avg;
tot->lastavg += diff = resL->avg + resR->avg - res->avg;
diff = fabs(.25*diff);
err = rL->err + rR->err;
err = resL->err + resR->err;
if( err > 0 ) {
creal c = 1 + 2*diff/err;
rL->err *= c;
rR->err *= c;
resL->err *= c;
resR->err *= c;
}
rL->err += diff;
rR->err += diff;
tot->lasterr += rL->err + rR->err - r->err;
resL->err += diff;
resR->err += diff;
tot->lasterr += resL->err + resR->err - res->err;
tot->weightsum += w = 1/Max(Sq(tot->lasterr), NOTZERO);
sigsq = 1/tot->weightsum;
@ -200,10 +216,22 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
tot->err = sqrt(sigsq);
}
}
++t->nregions;
if( StateWriteTest(t) ) {
StateWriteOpen(t, fd) {
Pool *prev = cur;
state->signature = StateSignature(t, 4);
state->nregions = t->nregions;
state->neval = t->neval;
StateWrite(fd, state, statesize);
while( (prev = prev->next) ) StateWrite(fd, prev, poolsize);
StateWrite(fd, cur, state->ncur*regionsize);
} StateWriteClose(t, fd);
}
}
for( comp = 0; comp < t->ncomp; ++comp ) {
cTotals *tot = &totals[comp];
for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
integral[comp] = tot->avg;
error[comp] = tot->err;
prob[comp] = ChiSquare(tot->chisq, t->nregions - 1);
@ -214,27 +242,20 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
MLPutFunction(stdlink, "List", 2);
MLPutFunction(stdlink, "List", t->nregions);
npool = ncur;
npool = state->ncur;
for( pool = cur; pool; npool = POOLSIZE, pool = pool->next )
for( ipool = 0; ipool < npool; ++ipool ) {
Region const *region = &pool->region[ipool];
real lower[NDIM], upper[NDIM];
Region const *region = RegionPtr(pool, ipool);
Result *Res;
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
lower[dim] = b->lower;
upper[dim] = b->upper;
}
MLPutFunction(stdlink, "Cuba`Cuhre`region", 3);
MLPutRealList(stdlink, lower, t->ndim);
MLPutRealList(stdlink, upper, t->ndim);
MLPutFunction(stdlink, "Cuba`Cuhre`region", 2);
MLPutRealList(stdlink, (real *)region->bounds, 2*t->ndim);
MLPutFunction(stdlink, "List", t->ncomp);
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
real res[] = {r->avg, r->err};
MLPutRealList(stdlink, res, Elements(res));
for( Res = (res = RegionResult(region)) + t->ncomp;
res < Res; ++res ) {
real r[] = {res->avg, res->err};
MLPutRealList(stdlink, r, Elements(r));
}
}
}
@ -246,8 +267,12 @@ abort:
free(pool);
}
WaitCores(t);
FrameFree(t);
RuleFree(t);
StateRemove(t);
return fail;
}

13
src/external/libCuba/src/cuhre/Makefile.am vendored Normal file
View File

@ -0,0 +1,13 @@
## Process this file with automake to create Makefile.in
## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $
c_sources = Cuhre.c
AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
noinst_LTLIBRARIES = libcuhre.la
libcuhre_la_SOURCES = $(c_sources)
libcuhre_la_LDFLAGS = $(AM_LDFLAGS)

210
src/external/libCuba/src/cuhre/Rule.c vendored
View File

@ -3,38 +3,12 @@
integration with cubature rules
code lifted with minor modifications from DCUHRE
by J. Berntsen, T. Espelid, and A. Genz
this file is part of Divonne
last modified 8 Jun 10 th
this file is part of Cuhre
last modified 5 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
enum { nrules = 5 };
#define TYPEDEFSET \
typedef struct { \
count n; \
real weight[nrules], scale[nrules], norm[nrules]; \
real gen[NDIM]; \
} Set
#define NextSet(p) p = (Set *)((char *)p + setsize)
/*********************************************************************/
@ -83,86 +57,84 @@ static void Rule13Alloc(This *t)
enum { nsets = 14, ndim = 2 };
TYPEDEFSET;
count n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
Copy(last->weight, w[0], nrules);
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[1], nrules);
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[2], nrules);
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[3], nrules);
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[4], nrules);
last->gen[0] = g[3];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[5], nrules);
last->gen[0] = g[4];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[6], nrules);
last->gen[0] = g[5];
last->gen[1] = g[5];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[7], nrules);
last->gen[0] = g[6];
last->gen[1] = g[6];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[8], nrules);
last->gen[0] = g[7];
last->gen[1] = g[7];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[9], nrules);
last->gen[0] = g[8];
last->gen[1] = g[8];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[10], nrules);
last->gen[0] = g[9];
last->gen[1] = g[9];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
Copy(last->weight, w[11], nrules);
last->gen[0] = g[10];
last->gen[1] = g[11];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
Copy(last->weight, w[12], nrules);
last->gen[0] = g[12];
last->gen[1] = g[13];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
Copy(last->weight, w[13], nrules);
last->gen[0] = g[14];
@ -175,12 +147,12 @@ static void Rule13Alloc(This *t)
t->rule.errcoeff[2] = 5;
t->rule.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -231,84 +203,82 @@ static void Rule11Alloc(This *t)
enum { nsets = 13, ndim = 3 };
TYPEDEFSET;
count n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
Copy(last->weight, w[0], nrules);
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[1], nrules);
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[2], nrules);
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[3], nrules);
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[4], nrules);
last->gen[0] = g[3];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[5], nrules);
last->gen[0] = g[4];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[6], nrules);
last->gen[0] = g[5];
last->gen[1] = g[5];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[7], nrules);
last->gen[0] = g[6];
last->gen[1] = g[6];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
Copy(last->weight, w[8], nrules);
last->gen[0] = g[7];
last->gen[1] = g[7];
last->gen[2] = g[7];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
Copy(last->weight, w[9], nrules);
last->gen[0] = g[8];
last->gen[1] = g[8];
last->gen[2] = g[8];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
Copy(last->weight, w[10], nrules);
last->gen[0] = g[9];
last->gen[1] = g[9];
last->gen[2] = g[9];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2);
Copy(last->weight, w[11], nrules);
last->gen[0] = g[10];
last->gen[1] = g[11];
last->gen[2] = g[11];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2);
Copy(last->weight, w[12], nrules);
last->gen[0] = g[12];
@ -322,12 +292,12 @@ static void Rule11Alloc(This *t)
t->rule.errcoeff[2] = 3;
t->rule.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -368,15 +338,13 @@ static void Rule9Alloc(This *t)
enum { nsets = 9 };
TYPEDEFSET;
ccount ndim = t->ndim;
ccount twondim = 1 << ndim;
count dim, n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
@ -386,7 +354,7 @@ static void Rule9Alloc(This *t)
last->weight[3] = ndim*(ndim*w[9] + w[10]) - 1 + last->weight[0];
last->weight[4] = ndim*w[11] + 1 - last->weight[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = ndim*(ndim*w[12] + w[13]) + w[14];
last->weight[1] = ndim*(ndim*w[15] + w[16]) + w[17];
@ -395,7 +363,7 @@ static void Rule9Alloc(This *t)
last->weight[4] = w[21] - last->weight[0];
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = ndim*w[22] + w[23];
last->weight[1] = ndim*w[24] + w[25];
@ -404,7 +372,7 @@ static void Rule9Alloc(This *t)
last->weight[4] = -last->weight[0];
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = w[29];
last->weight[1] = w[30];
@ -413,12 +381,12 @@ static void Rule9Alloc(This *t)
last->weight[4] = -w[29];
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[2] = w[32];
last->gen[0] = g[3];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
last->weight[0] = w[33] - ndim*w[12];
last->weight[1] = w[34] - ndim*w[15];
@ -428,7 +396,7 @@ static void Rule9Alloc(This *t)
last->gen[0] = g[0];
last->gen[1] = g[0];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
last->weight[0] = w[36];
last->weight[1] = w[37];
@ -438,7 +406,7 @@ static void Rule9Alloc(This *t)
last->gen[0] = g[0];
last->gen[1] = g[1];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
last->weight[0] = w[39];
last->weight[1] = w[40];
@ -449,7 +417,7 @@ static void Rule9Alloc(This *t)
last->gen[1] = g[0];
last->gen[2] = g[0];
++last;
NextSet(last);
n += last->n = twondim;
last->weight[0] = w[41]/twondim;
last->weight[1] = w[7]/twondim;
@ -466,12 +434,12 @@ static void Rule9Alloc(This *t)
t->rule.errcoeff[2] = 5;
t->rule.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -501,15 +469,13 @@ static void Rule7Alloc(This *t)
enum { nsets = 6 };
TYPEDEFSET;
ccount ndim = t->ndim;
ccount twondim = 1 << ndim;
count dim, n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
@ -519,7 +485,7 @@ static void Rule7Alloc(This *t)
last->weight[3] = ndim*(ndim*w[7] + w[8]) - w[9];
last->weight[4] = 1 - last->weight[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = w[10];
last->weight[1] = w[11];
@ -528,7 +494,7 @@ static void Rule7Alloc(This *t)
last->weight[4] = -w[10];
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = w[13] - ndim*w[0];
last->weight[1] = w[14] - ndim*w[3];
@ -537,12 +503,12 @@ static void Rule7Alloc(This *t)
last->weight[4] = -last->weight[0];
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[2] = w[17];
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
last->weight[0] = -w[7];
last->weight[1] = w[18];
@ -552,7 +518,7 @@ static void Rule7Alloc(This *t)
last->gen[0] = g[0];
last->gen[1] = g[0];
++last;
NextSet(last);
n += last->n = twondim;
last->weight[0] = w[20]/twondim;
last->weight[1] = w[5]/twondim;
@ -569,12 +535,12 @@ static void Rule7Alloc(This *t)
t->rule.errcoeff[2] = 5;
t->rule.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -594,15 +560,12 @@ static inline void RuleAlloc(This *t)
else if( t->ndim == 3 ) Rule11Alloc(t);
else Rule9Alloc(t);
}
Alloc(t->rule.x, t->rule.n*(t->ndim + t->ncomp));
t->rule.f = t->rule.x + t->rule.n*t->ndim;
}
/*********************************************************************/
static inline void RuleFree(cThis *t)
{
free(t->rule.x);
free(t->rule.first);
}
@ -657,25 +620,23 @@ next:
/*********************************************************************/
static void Sample(This *t, void *voidregion)
static void Sample(This *t, Region *region)
{
TYPEDEFREGION;
TYPEDEFSET;
Region *const region = (Region *)voidregion;
csize_t setsize = SetSize;
creal vol = ldexp(1., -region->div);
real *x = t->rule.x, *f = t->rule.f;
Set *first = (Set *)t->rule.first, *last = (Set *)t->rule.last, *s;
real *x = t->frame, *f = x + t->rule.n*t->ndim;
Set *first = t->rule.first, *last = t->rule.last, *s;
Bounds *b, *B = region->bounds + t->ndim;
Result *result = RegionResult(region), *res, *Res = result + t->ncomp;
creal *errcoeff = t->rule.errcoeff;
creal ratio = Sq(first[2].gen[0]/first[1].gen[0]);
ccount offset = 2*t->ndim*t->ncomp;
count dim, comp, rul, n, maxdim = 0;
count dim, rul, n, maxdim = 0;
real maxrange = 0;
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
for( b = region->bounds, dim = 0; b < B; ++b, ++dim ) {
creal range = b->upper - b->lower;
if( range > maxrange ) {
maxrange = range;
@ -683,13 +644,12 @@ static void Sample(This *t, void *voidregion)
}
}
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
if( s->n ) x = ExpandFS(t, region->bounds, s->gen, x);
DoSample(t, t->rule.n, t->rule.x, f);
DoSample(t, t->rule.n, t->frame, f);
for( comp = 0; comp < t->ncomp; ++comp ) {
Result *r = &region->result[comp];
for( res = result; res < Res; ++res ) {
real sum[nrules];
creal *f1 = f;
creal base = *f1*2*(1 - ratio);
@ -707,11 +667,11 @@ static void Sample(This *t, void *voidregion)
bisectdim = dim;
}
}
r->bisectdim = bisectdim;
res->bisectdim = bisectdim;
f1 = f++;
Zap(sum);
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( n = s->n; n; --n ) {
creal fun = *f1;
f1 += t->ncomp;
@ -726,37 +686,35 @@ static void Sample(This *t, void *voidregion)
for( rul = 1; rul < nrules - 1; ++rul ) {
real maxerr = 0;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
maxerr = Max(maxerr,
fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
sum[rul] = maxerr;
}
r->avg = vol*sum[0];
r->err = vol*(
res->avg = vol*sum[0];
res->err = vol*(
(errcoeff[0]*sum[1] <= sum[2] && errcoeff[0]*sum[2] <= sum[3]) ?
errcoeff[1]*sum[1] :
errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3]) );
}
if( VERBOSE > 2 ) {
char s[64*NDIM + 128*NCOMP], *p = s;
Vector(char, out, 64*NDIM + 128*NCOMP);
char *oe = out;
count comp;
cchar *msg = "\nRegion (" REALF ") - (" REALF ")";
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
p += sprintf(p,
(dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" :
"\n (" REALF ") - (" REALF ")",
b->lower, b->upper);
for( b = region->bounds; b < B; ++b ) {
oe += sprintf(oe, msg, b->lower, b->upper);
msg = "\n (" REALF ") - (" REALF ")";
}
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
p += sprintf(p, "\n[" COUNT "] "
REAL " +- " REAL, comp + 1, r->avg, r->err);
}
for( res = result, comp = 0; res < Res; ++res )
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL, ++comp, res->avg, res->err);
Print(s);
Print(out);
}
}

28
src/external/libCuba/src/cuhre/common.c vendored
View File

@ -2,44 +2,22 @@
common.c
includes most of the modules
this file is part of Cuhre
last modified 7 Jun 10 th
last modified 2 Aug 13 11 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "ChiSquare.c"
#include "Rule.c"
static inline bool BadDimension(cThis *t)
{
if( t->ndim > NDIM ) return true;
if( t->ndim > MAXDIM ) return true;
return t->ndim < 2;
}
static inline bool BadComponent(cThis *t)
{
if( t->ncomp > NCOMP ) return true;
if( t->ncomp > MAXCOMP ) return true;
return t->ncomp < 1;
}
#include "Integrate.c"

59
src/external/libCuba/src/cuhre/decl.h vendored
View File

@ -2,27 +2,9 @@
decl.h
Type declarations
this file is part of Cuhre
last modified 7 Jun 10 th */
last modified 26 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "stddecl.h"
@ -47,9 +29,18 @@ typedef struct {
typedef const Bounds cBounds;
enum { nrules = 5 };
typedef struct {
real *x, *f;
void *first, *last;
count n;
real weight[nrules], scale[nrules], norm[nrules];
real gen[];
} Set;
#define SetSize (sizeof(Set) + t->ndim*sizeof(real))
typedef struct {
Set *first, *last;
real errcoeff[3];
count n;
} Rule;
@ -63,22 +54,34 @@ typedef struct _this {
#ifndef MLVERSION
Integrand integrand;
void *userdata;
#ifdef HAVE_FORK
int ncores, *child;
SHM_ONLY(int shmid;)
#endif
#endif
real *frame;
real epsrel, epsabs;
int flags;
number mineval, maxeval;
count key, nregions;
cchar *statefile;
number neval;
Rule rule;
jmp_buf abort;
} This;
#define nframe rule.n
typedef const This cThis;
#define TYPEDEFREGION \
typedef struct region { \
count div; \
Result result[NCOMP]; \
Bounds bounds[NDIM]; \
} Region
typedef struct region {
count div;
Bounds bounds[];
} Region;
#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*sizeof(Result))
#define RegionResult(r) ((Result *)(r->bounds + t->ndim))
#define RegionPtr(p, n) ((Region *)((char *)p->region + (n)*regionsize))

37
src/external/libCuba/src/cuhre/util.c vendored
View File

@ -1,37 +0,0 @@
/*
util.c
Utility functions
this file is part of Cuhre
last modified 9 Jan 05 th
*/
/***************************************************************************
* Copyright (C) 2004-2009 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "decl.h"
static count ndim_, ncomp_, nregions_;
static number neval_;
#ifdef DEBUG
#include "debug.c"
#endif

100
src/external/libCuba/src/divonne/Divonne.c vendored
View File

@ -4,83 +4,17 @@
originally by J.H. Friedman and M.H. Wright
(CERNLIB subroutine D151)
this version by Thomas Hahn
last modified 16 Jun 10 th
last modified 17 Sep 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define DIVONNE
#define ROUTINE "Divonne"
#include "decl.h"
#define Print(s) puts(s); fflush(stdout)
#include "CSample.c"
/*********************************************************************/
static inline void DoSample(This *t, number n, ccount ldx, creal *x, real *f)
{
t->neval += n;
while( n-- ) {
if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata, &t->phase) == ABORT )
longjmp(t->abort, 1);
x += ldx;
f += t->ncomp;
}
}
/*********************************************************************/
static inline count SampleExtra(This *t, cBounds *b)
{
number n = t->nextra;
t->peakfinder(&t->ndim, b, &n, t->xextra);
DoSample(t, n, t->ldxgiven, t->xextra, t->fextra);
return n;
}
/*********************************************************************/
static inline void AllocGiven(This *t, creal *xgiven)
{
if( t->ngiven | t->nextra ) {
cnumber nxgiven = t->ngiven*(t->ldxgiven = IMax(t->ldxgiven, t->ndim));
cnumber nxextra = t->nextra*t->ldxgiven;
cnumber nfgiven = t->ngiven*t->ncomp;
cnumber nfextra = t->nextra*t->ncomp;
Alloc(t->xgiven, nxgiven + nxextra + nfgiven + nfextra);
t->xextra = t->xgiven + nxgiven;
t->fgiven = t->xextra + nxextra;
t->fextra = t->fgiven + nfgiven;
if( nxgiven ) {
t->phase = 0;
Copy(t->xgiven, xgiven, nxgiven);
DoSample(t, t->ngiven, t->ldxgiven, t->xgiven, t->fgiven);
}
}
}
/*********************************************************************/
#include "common.c"
Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
Integrand integrand, void *userdata,
creal epsrel, creal epsabs,
@ -88,8 +22,9 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
cnumber mineval, cnumber maxeval,
cint key1, cint key2, cint key3, ccount maxpass,
creal border, creal maxchisq, creal mindeviation,
cnumber ngiven, ccount ldxgiven, creal *xgiven,
cnumber ngiven, ccount ldxgiven, real *xgiven,
cnumber nextra, PeakFinder peakfinder,
cchar *statefile,
int *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
{
@ -112,20 +47,15 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
t.maxchisq = maxchisq;
t.mindeviation = mindeviation;
t.ngiven = ngiven;
t.xgiven = NULL;
t.xgiven = xgiven;
t.ldxgiven = ldxgiven;
t.nextra = nextra;
t.peakfinder = peakfinder;
t.nregions = 0;
t.neval = 0;
AllocGiven(&t, xgiven);
t.statefile = statefile;
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
*pneval = t.neval;
free(t.xgiven);
}
/*********************************************************************/
@ -137,10 +67,11 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
cnumber *pmineval, cnumber *pmaxeval,
cint *pkey1, cint *pkey2, cint *pkey3, ccount *pmaxpass,
creal *pborder, creal *pmaxchisq, creal *pmindeviation,
cnumber *pngiven, ccount *pldxgiven, creal *xgiven,
cnumber *pngiven, ccount *pldxgiven, real *xgiven,
cnumber *pnextra, PeakFinder peakfinder,
cchar *statefile,
int *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
real *integral, real *error, real *prob, cint statefilelen)
{
This t;
t.ndim = *pndim;
@ -161,19 +92,14 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
t.maxchisq = *pmaxchisq;
t.mindeviation = *pmindeviation;
t.ngiven = *pngiven;
t.xgiven = NULL;
t.xgiven = xgiven;
t.ldxgiven = *pldxgiven;
t.nextra = *pnextra;
t.peakfinder = peakfinder;
t.nregions = 0;
t.neval = 0;
AllocGiven(&t, xgiven);
CString(t.statefile, statefile, statefilelen);
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
*pneval = t.neval;
free(t.xgiven);
}

96
src/external/libCuba/src/divonne/Explore.c vendored
View File

@ -2,29 +2,9 @@
Explore.c
sample region, determine min and max, split if necessary
this file is part of Divonne
last modified 8 Jun 10 th
last modified 2 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
typedef struct {
real fmin, fmax;
@ -33,38 +13,21 @@ typedef struct {
/*********************************************************************/
static bool Explore(This *t, count iregion, cSamples *samples,
cint depth, cint flags)
static int ExploreSerial(This *t, ccount iregion)
{
#define SPLICE (flags & 1)
#define HAVESAMPLES (flags & 2)
csize_t regionsize = RegionSize;
Region *region = RegionPtr(iregion);
cBounds *bounds = region->bounds;
Result *result = RegionResult(region);
TYPEDEFREGION;
count n, dim, comp, maxcomp;
Extrema extrema[NCOMP];
Result *r;
Vector(Extrema, extrema, NCOMP);
Vector(real, xtmp, NDIM);
Result *r, *r0;
creal *x;
real *f;
real halfvol, maxerr;
Region *region;
Bounds *bounds;
Result *result;
/* needed as of gcc 3.3 to make gcc correctly address region #@$&! */
sizeof(*region);
if( SPLICE ) {
if( t->nregions == t->size ) {
t->size += CHUNKSIZE;
ReAlloc(t->voidregion, t->size*sizeof(Region));
}
VecCopy(RegionPtr(t->nregions)->bounds, RegionPtr(iregion)->bounds);
iregion = t->nregions++;
}
region = RegionPtr(iregion);
bounds = region->bounds;
result = region->result;
count n, dim, comp, maxcomp;
cSamples *samples = &t->samples[region->isamples];
for( comp = 0; comp < t->ncomp; ++comp ) {
Extrema *e = &extrema[comp];
@ -73,7 +36,7 @@ static bool Explore(This *t, count iregion, cSamples *samples,
e->xmin = e->xmax = NULL;
}
if( !HAVESAMPLES ) {
if( region->isamples == 0 ) { /* others already sampled */
real vol = 1;
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &bounds[dim];
@ -108,7 +71,7 @@ skip:
f += t->ncomp;
}
samples->sampler(t, samples, bounds, vol);
samples->sampler(t, iregion);
}
x = samples->x;
@ -131,30 +94,27 @@ skip:
for( comp = 0; comp < t->ncomp; ++comp ) {
Extrema *e = &extrema[comp];
Result *r = &result[comp];
real xtmp[NDIM], ftmp, err;
real ftmp, err;
if( e->xmin ) { /* not all NaNs */
t->selectedcomp = comp;
VecCopy(xtmp, e->xmin);
XCopy(xtmp, e->xmin);
ftmp = FindMinimum(t, bounds, xtmp, e->fmin);
if( ftmp < r->fmin ) {
r->fmin = ftmp;
VecCopy(r->xmin, xtmp);
XCopy(&r->xminmax[0], xtmp);
}
t->selectedcomp = Tag(comp);
VecCopy(xtmp, e->xmax);
XCopy(xtmp, e->xmax);
ftmp = -FindMinimum(t, bounds, xtmp, -e->fmax);
if( ftmp > r->fmax ) {
r->fmax = ftmp;
VecCopy(r->xmax, xtmp);
XCopy(&r->xminmax[t->ndim], xtmp);
}
}
r->avg = samples->avg[comp];
r->err = samples->err[comp];
r->spread = halfvol*(r->fmax - r->fmin);
err = r->spread/Max(fabs(r->avg), NOTZERO);
if( err > maxerr ) {
maxerr = err;
@ -166,26 +126,25 @@ skip:
if( maxcomp == -1 ) { /* all NaNs */
region->depth = 0;
return false;
return -1;
}
region->cutcomp = maxcomp;
r = &region->result[maxcomp];
r0 = RegionResult(region);
r = r0 + maxcomp;
if( halfvol*(r->fmin + r->fmax) > r->avg ) {
region->fminor = r->fmin;
region->fmajor = r->fmax;
region->xmajor = r->xmax - (real *)region->result;
region->xmajor = &r->xminmax[t->ndim] - (real *)r0;
}
else {
region->fminor = r->fmax;
region->fmajor = r->fmin;
region->xmajor = r->xmin - (real *)region->result;
region->xmajor = &r->xminmax[0] - (real *)r0;
}
region->depth = IDim(depth);
if( !HAVESAMPLES ) {
if( samples->weight*r->spread < r->err ||
if( region->isamples == 0 ) {
if( (region->depth < INIDEPTH && r->spread < samples->neff*r->err) ||
r->spread < t->totals[maxcomp].secondspread ) region->depth = 0;
if( region->depth == 0 )
for( comp = 0; comp < t->ncomp; ++comp )
@ -193,7 +152,8 @@ skip:
Max(t->totals[comp].secondspread, result[comp].spread);
}
if( region->depth ) Split(t, iregion, region->depth);
return true;
if( region->depth ) Split(t, iregion);
return iregion;
}

77
src/external/libCuba/src/divonne/FindMinimum.c vendored
View File

@ -2,37 +2,17 @@
FindMinimum.c
find minimum (maximum) of hyperrectangular region
this file is part of Divonne
last modified 8 Jun 10 th
last modified 7 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define EPS POW2(52)
#define RTEPS POW2(26)
#define QEPS POW2(13)
#define EPS 0x1p-52
#define RTEPS 0x1p-26
#define QEPS 0x1p-13
#define DELTA 0x1p-16
#define RTDELTA 0x1p-8
#define QDELTA 0x1p-4
#define DELTA POW2(16)
#define RTDELTA POW2(8)
#define QDELTA POW2(4)
/*
#define DELTA 1e-5
@ -162,7 +142,8 @@ static void UpdateCholesky(cThis *t, ccount n, real *hessian,
static inline void BFGS(cThis *t, ccount n, real *hessian,
creal *gnew, creal *g, real *p, creal dx)
{
real y[NDIM], c;
Vector(real, y, NDIM);
real c;
count i, j;
for( i = 0; i < n; ++i )
@ -210,7 +191,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
real tol = ftol, tol2 = tol + tol;
Point cur = {0, fini};
VecCopy(x, xini);
XCopy(x, xini);
/* don't even try if
a) we'd walk backwards,
@ -360,9 +341,10 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
static real LocalSearch(This *t, ccount nfree, ccount *ifree,
cBounds *b, creal *x, creal fx, real *z)
{
Vector(real, y, NDIM);
Vector(real, p, NDIM);
real delta, smax, sopp, spmax, snmax;
real y[NDIM], fy, fz, ftest;
real p[NDIM];
real fy, fz, ftest;
int sign;
count i;
@ -389,7 +371,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
/* Move along p until the integrand changes appreciably
or we come close to a border. */
VecCopy(y, x);
XCopy(y, x);
ftest = SUFTOL*(1 + fabs(fx));
delta = RTDELTA/5;
do {
@ -442,7 +424,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
/* Move along p' until the integrand changes appreciably
or we come close to a border. */
VecCopy(z, y);
XCopy(z, y);
ftest = SUFTOL*(1 + fabs(fy));
delta = RTDELTA/5;
do {
@ -468,7 +450,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
grad = (fy - fz)/delta;
range = sopp/.9 + delta;
step = Min(delta + delta, sopp);
VecCopy(y, z);
XCopy(y, z);
fy = fz;
for( i = 0; i < nfree; ++i )
p[i] = -p[i];
@ -535,15 +517,18 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
{
real hessian[NDIM*NDIM];
real gfree[NDIM], p[NDIM];
real tmp[NDIM], ftmp, fini = fmin;
Vector(real, hessian, NDIM*NDIM);
Vector(real, gfree, NDIM);
Vector(real, p, NDIM);
Vector(real, tmp, NDIM);
Vector(count, ifree, NDIM);
Vector(count, ifix, NDIM);
real ftmp, fini = fmin;
ccount maxeval = t->neval + 50*t->ndim;
count nfree, nfix;
count ifree[NDIM], ifix[NDIM];
count dim, local;
Zap(hessian);
Clear(hessian, t->ndim*t->ndim);
for( dim = 0; dim < t->ndim; ++dim )
Hessian(dim, dim) = 1;
@ -580,7 +565,7 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
if( ftmp > fmin - (1 + fabs(fmin))*RTEPS )
goto releasebounds;
fmin = ftmp;
VecCopy(xmin, tmp);
XCopy(xmin, tmp);
}
while( t->neval <= maxeval ) {
@ -632,13 +617,13 @@ fixbound:
tmp[i] = Hessian(i + 1, mini);
for( i = mini; i < nfree; ++i ) {
Copy(&Hessian(i, 0), &Hessian(i + 1, 0), i);
Move(&Hessian(i, 0), &Hessian(i + 1, 0), i);
Hessian(i, i) = Hessian(i + 1, i + 1);
}
RenormalizeCholesky(t, nfree, hessian, tmp, diag);
Copy(&ifree[mini], &ifree[mini + 1], nfree - mini);
Copy(&gfree[mini], &gfree[mini + 1], nfree - mini);
Move(&ifree[mini], &ifree[mini + 1], nfree - mini);
Move(&gfree[mini], &gfree[mini + 1], nfree - mini);
}
continue;
}
@ -651,11 +636,11 @@ fixbound:
real fdiff;
fmin = low.f;
VecCopy(xmin, tmp);
XCopy(xmin, tmp);
Gradient(t, nfree, ifree, b, xmin, fmin, tmp);
BFGS(t, nfree, hessian, tmp, gfree, p, low.dx);
VecCopy(gfree, tmp);
XCopy(gfree, tmp);
if( fabs(low.dx - minstep) < QEPS*minstep ) goto fixbound;
@ -696,7 +681,7 @@ releasebounds:
++nfree;
--nfix;
Copy(&ifix[mini], &ifix[mini + 1], nfix - mini);
Move(&ifix[mini], &ifix[mini + 1], nfix - mini);
continue;
}
}

541
src/external/libCuba/src/divonne/Integrate.c vendored
View File

@ -4,64 +4,56 @@
has approximately equal spread = 1/2 vol (max - min),
then do a main integration over all regions
this file is part of Divonne
last modified 8 Jun 10 th
checkpointing by B. Chokoufe
last modified 5 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define INIDEPTH 3
#define DEPTH 5
#define POSTDEPTH 15
/*********************************************************************/
typedef struct {
signature_t signature;
number neval, neval_opt, neval_cut;
number nmin, nrand;
count iregion, nregions;
count phase, iter, pass, size;
Totals totals[];
} State;
static int Integrate(This *t, real *integral, real *error, real *prob)
{
TYPEDEFREGION;
StateDecl;
csize_t statesize = sizeof(State) + NCOMP*sizeof(Totals);
Sized(State, state, statesize);
csize_t regionsize = RegionSize;
Vector(char, out, 64*NDIM + 256*NCOMP);
Totals totals[NCOMP];
real nneed, weight;
count dim, comp, iter, pass = 0, err, iregion;
number nwant, nmin = INT_MAX;
int fail = -99;
Totals *tot, *Tot = state->totals + t->ncomp;
Bounds *b, *B;
Result *res;
count comp, err, iregion;
number nwant;
real nneed;
ML_ONLY(number neff;)
int fail;
if( VERBOSE > 1 ) {
char s[512];
sprintf(s, "Divonne input parameters:\n"
sprintf(out, "Divonne input parameters:\n"
" ndim " COUNT "\n ncomp " COUNT "\n"
" epsrel " REAL "\n epsabs " REAL "\n"
" flags %d\n seed %d\n"
" mineval " NUMBER "\n maxeval " NUMBER "\n"
" key1 %d\n key2 %d\n key3 %d\n maxpass " COUNT "\n"
" border " REAL "\n maxchisq " REAL "\n mindeviation " REAL "\n"
" ngiven " NUMBER "\n nextra " NUMBER "\n",
" ngiven " NUMBER "\n nextra " NUMBER "\n"
" statefile \"%s\"",
t->ndim, t->ncomp,
t->epsrel, t->epsabs,
t->flags, t->seed,
t->mineval, t->maxeval,
t->key1, t->key2, t->key3, t->maxpass,
t->border.lower, t->maxchisq, t->mindeviation,
t->ngiven, t->nextra);
Print(s);
t->ngiven, t->nextra,
t->statefile);
Print(out);
}
if( BadComponent(t) ) return -2;
@ -69,124 +61,169 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
BadDimension(t, t->key2) ||
((t->key3 & -2) && BadDimension(t, t->key3)) ) return -1;
t->neval_opt = t->neval_cut = 0;
FORK_ONLY(t->nframe = 0;)
t->size = CHUNKSIZE;
MemAlloc(t->voidregion, t->size*sizeof(Region));
for( dim = 0; dim < t->ndim; ++dim ) {
Bounds *b = &RegionPtr(0)->bounds[dim];
b->lower = 0;
b->upper = 1;
}
RuleIni(&t->rule7);
RuleIni(&t->rule9);
RuleIni(&t->rule11);
RuleIni(&t->rule13);
AllocGiven(t);
SamplesIni(&t->samples[0]);
SamplesIni(&t->samples[1]);
SamplesIni(&t->samples[2]);
if( setjmp(t->abort) ) goto abort;
t->epsabs = Max(t->epsabs, NOTZERO);
/* Step 1: partition the integration region */
if( VERBOSE ) Print("Partitioning phase:");
if( IsSobol(t->key1) || IsSobol(t->key2) || IsSobol(t->key3) )
RuleAlloc(t);
if( IsSobol(t->key1) | IsSobol(t->key2) | IsSobol(t->key3) )
IniRandom(t);
t->epsabs = Max(t->epsabs, NOTZERO);
t->totals = state->totals;
ForkCores(t);
if( (fail = setjmp(t->abort)) ) goto abort;
SamplesLookup(t, &t->samples[0], t->key1,
(number)47, (number)INT_MAX, (number)0);
SamplesAlloc(t, &t->samples[0]);
t->totals = totals;
Zap(totals);
t->phase = 1;
StateSetup(t);
Explore(t, 0, &t->samples[0], INIDEPTH, 1);
if( StateReadTest(t) ) {
StateReadOpen(t, fd) {
if( read(fd, state, statesize) != statesize ||
state->signature != StateSignature(t, 3) ) break;
t->nregions = state->nregions;
if( st.st_size != statesize + t->nregions*regionsize ) break;
t->neval = state->neval;
t->neval_opt = state->neval_opt;
t->neval_cut = state->neval_cut;
t->nrand = state->nrand;
t->phase = state->phase;
t->size = state->size;
AllocRegions(t);
StateRead(fd, t->region, t->nregions*regionsize);
} StateReadClose(t, fd);
for( iter = 1; ; ++iter ) {
Totals *maxtot;
count valid;
if( IsSobol(t->key1) | IsSobol(t->key2) | IsSobol(t->key3) )
t->rng.skiprandom(t, t->nrand);
}
for( comp = 0; comp < t->ncomp; ++comp ) {
Totals *tot = &totals[comp];
tot->avg = tot->spreadsq = 0;
tot->spread = tot->secondspread = -INFTY;
if( ini ) {
#if MLVERSION
t->neval = t->ngiven;
#else
t->neval = 0;
#endif
t->neval_opt = 0;
t->neval_cut = 0;
t->nrand = 0;
t->size = CHUNKSIZE;
AllocRegions(t);
for( B = (b = t->region->bounds) + t->ndim; b < B; ++b ) {
b->lower = 0;
b->upper = 1;
}
t->nregions = 1;
for( iregion = 0; iregion < t->nregions; ++iregion ) {
Region *region = RegionPtr(iregion);
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
Totals *tot = &totals[comp];
tot->avg += r->avg;
tot->spreadsq += Sq(r->spread);
if( r->spread > tot->spread ) {
tot->secondspread = tot->spread;
tot->spread = r->spread;
tot->iregion = iregion;
t->phase = 1;
state->iter = 1;
state->pass = 0;
state->nmin = INT_MAX;
state->iregion = 0;
FClear(state->totals);
}
/* Step 1: partition the integration region */
if( t->phase == 1 ) {
if( VERBOSE ) Print("Partitioning phase:");
if( ini ) Iterate(t, 0, INIDEPTH, 0, NULL);
for( ; ; ++state->iter ) {
Totals *maxtot;
count valid;
for( tot = state->totals; tot < Tot; ++tot ) {
tot->avg = tot->spreadsq = 0;
tot->spread = tot->secondspread = -INFTY;
}
for( iregion = 0; iregion < t->nregions; ++iregion )
for( res = RegionResult(RegionPtr(iregion)), tot = state->totals;
tot < Tot; ++res, ++tot ) {
tot->avg += res->avg;
tot->spreadsq += Sq(res->spread);
if( res->spread > tot->spread ) {
tot->secondspread = tot->spread;
tot->spread = res->spread;
tot->iregion = iregion;
}
else if( res->spread > tot->secondspread )
tot->secondspread = res->spread;
}
else if( r->spread > tot->secondspread )
tot->secondspread = r->spread;
valid = 0;
for( maxtot = tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
integral[comp] = tot->avg;
valid += tot->avg == tot->avg;
if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot;
tot->spread = sqrt(tot->spreadsq);
error[comp] = tot->spread/t->samples[0].neff;
}
}
maxtot = totals;
valid = 0;
for( comp = 0; comp < t->ncomp; ++comp ) {
Totals *tot = &totals[comp];
integral[comp] = tot->avg;
valid += tot->avg == tot->avg;
if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot;
tot->spread = sqrt(tot->spreadsq);
error[comp] = tot->spread*t->samples[0].weight;
}
#define WriteState(t) \
if( StateWriteTest(t) ) { \
StateWriteOpen(t, fd) { \
state->signature = StateSignature(t, 3); \
state->neval = t->neval; \
state->neval_opt = t->neval_opt; \
state->neval_cut = t->neval_cut; \
state->nrand = t->nrand; \
state->nregions = t->nregions; \
state->phase = t->phase; \
state->size = t->size; \
StateWrite(fd, state, statesize); \
StateWrite(fd, t->region, t->nregions*regionsize); \
} StateWriteClose(t, fd); \
}
if( VERBOSE ) {
char s[128 + 64*NCOMP], *p = s;
WriteState(t);
p += sprintf(p, "\n"
"Iteration " COUNT " (pass " COUNT "): " COUNT " regions\n"
NUMBER7 " integrand evaluations so far,\n"
NUMBER7 " in optimizing regions,\n"
NUMBER7 " in finding cuts",
iter, pass, t->nregions, t->neval, t->neval_opt, t->neval_cut);
for( comp = 0; comp < t->ncomp; ++comp )
p += sprintf(p, "\n[" COUNT "] "
REAL " +- " REAL,
comp + 1, integral[comp], error[comp]);
Print(s);
}
if( valid == 0 ) goto abort; /* all NaNs */
if( t->neval > t->maxeval ) break;
nneed = maxtot->spread/MaxErr(maxtot->avg);
if( nneed < MAXPRIME ) {
cnumber n = t->neval + t->nregions*(number)ceil(nneed);
if( n < nmin ) {
nmin = n;
pass = 0;
if( VERBOSE ) {
char *oe = out + sprintf(out, "\n"
"Iteration " COUNT " (pass " COUNT "): " COUNT " regions\n"
NUMBER7 " integrand evaluations so far,\n"
NUMBER7 " in optimizing regions,\n"
NUMBER7 " in finding cuts",
state->iter, state->pass, t->nregions,
t->neval, t->neval_opt, t->neval_cut);
for( comp = 0; comp < t->ncomp; ++comp )
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL,
comp + 1, integral[comp], error[comp]);
Print(out);
}
else if( ++pass > t->maxpass && n >= t->mineval ) break;
}
Split(t, maxtot->iregion, DEPTH);
if( valid == 0 ) goto abort; /* all NaNs */
if( t->neval > t->maxeval ) break;
nneed = maxtot->spread/MaxErr(maxtot->avg);
if( nneed < MAXPRIME ) {
cnumber n = t->neval + t->nregions*(number)ceil(nneed);
if( n < state->nmin ) {
state->nmin = n;
state->pass = 0;
}
else if( ++state->pass > t->maxpass && n >= t->mineval ) break;
}
Iterate(t, maxtot->iregion, DEPTH, -1, NULL);
}
}
/* Step 2: do a "full" integration on each region */
/* nneed = t->samples[0].neff + 1; */
nneed = 2*t->samples[0].neff;
for( comp = 0; comp < t->ncomp; ++comp ) {
Totals *tot = &totals[comp];
for( tot = state->totals; tot < Tot; ++tot ) {
creal maxerr = MaxErr(tot->avg);
tot->nneed = tot->spread/maxerr;
nneed = Max(nneed, tot->nneed);
@ -205,85 +242,73 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
if( VERBOSE ) Print("\nNot enough samples left for main integration.");
for( comp = 0; comp < t->ncomp; ++comp )
prob[comp] = -999;
weight = t->samples[0].weight;
ML_ONLY(neff = t->samples[0].neff;)
}
else {
bool can_adjust = (t->key3 == 1 && t->samples[1].sampler != SampleRule &&
bool can_adjust = (t->key3 == 1 &&
t->samples[1].sampler != SampleRule &&
(t->key2 < 0 || t->samples[1].neff < MAXPRIME));
count df, nlimit;
SamplesAlloc(t, &t->samples[1]);
if( VERBOSE ) {
char s[128];
sprintf(s, "\nMain integration on " COUNT
sprintf(out, "\nMain integration on " COUNT
" regions with " NUMBER " samples per region.",
t->nregions, t->samples[1].neff);
Print(s);
Print(out);
}
ResClear(integral);
ResClear(error);
ResClear(prob);
nlimit = t->maxeval - t->nregions*t->samples[1].n;
df = 0;
for( iregion = 0; iregion < t->nregions; ++iregion ) {
Region *region = RegionPtr(iregion);
char s[64*NDIM + 256*NCOMP], *p = s;
#define CopyPhaseResults(f) \
for( res = RegionResult(region), tot = state->totals; tot < Tot; ++res, ++tot ) { \
PhaseResult *p = &tot->phase[f]; \
p->avg = res->avg; \
p->err = res->err; \
}
#define Var2(f, r) Sq((r)->err ? (r)->err : res->spread/t->samples[f].neff)
#define Var(f) Var2(f, &tot->phase[f])
while( state->iregion < t->nregions ) {
Region *region;
char *oe = out;
int todo;
refine:
region = RegionPtr(state->iregion);
CopyPhaseResults(0);
t->phase = 2;
t->samples[1].sampler(t, &t->samples[1], region->bounds, region->vol);
region->isamples = 1;
t->samples[1].sampler(t, state->iregion);
CopyPhaseResults(1);
if( can_adjust )
for( comp = 0; comp < t->ncomp; ++comp )
totals[comp].spreadsq -= Sq(region->result[comp].spread);
for( res = RegionResult(region), tot = state->totals;
tot < Tot; ++res, ++tot )
tot->spreadsq -= Sq(res->spread);
nlimit += t->samples[1].n;
todo = 0;
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
Totals *tot = &totals[comp];
t->samples[0].avg[comp] = r->avg;
t->samples[0].err[comp] = r->err;
for( res = RegionResult(region), tot = state->totals;
tot < Tot; ++tot ) {
if( t->neval < nlimit ) {
creal avg2 = t->samples[1].avg[comp];
creal err2 = t->samples[1].err[comp];
creal diffsq = Sq(avg2 - r->avg);
creal avg2 = tot->phase[1].avg;
creal diffsq = Sq(avg2 - tot->phase[0].avg);
#define Var(s) Sq((s.err[comp] == 0) ? r->spread*s.weight : s.err[comp])
if( err2*tot->nneed > r->spread ||
diffsq > Max(t->maxchisq*(Var(t->samples[0]) + Var(t->samples[1])),
EPS*Sq(avg2)) ) {
if( res->err*tot->nneed > res->spread ||
diffsq > Max(t->maxchisq*(Var(0) + Var(1)), EPS*Sq(avg2)) ) {
if( t->key3 && diffsq > tot->mindevsq ) {
if( t->key3 == 1 ) {
ccount xregion = t->nregions;
if( VERBOSE > 2 ) Print("\nSplit");
t->phase = 1;
Explore(t, iregion, &t->samples[1], POSTDEPTH, 2);
Iterate(t, state->iregion, POSTDEPTH, 1, state->totals);
if( can_adjust ) {
number nnew;
count ireg, xreg;
for( ireg = iregion, xreg = xregion;
ireg < t->nregions; ireg = xreg++ ) {
cResult *result = RegionPtr(ireg)->result;
count c;
for( c = 0; c < t->ncomp; ++c )
totals[c].spreadsq += Sq(result[c].spread);
}
nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ?
cnumber nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ?
MARKMASK :
(number)ceil(sqrt(tot->spreadsq/tot->maxerrsq));
if( nnew > nwant + nwant/64 ) {
@ -298,15 +323,13 @@ refine:
can_adjust = false;
if( VERBOSE > 2 ) {
char s[128];
sprintf(s, "Sampling remaining " COUNT
sprintf(out, "Sampling remaining " COUNT
" regions with " NUMBER " points per region.",
t->nregions, t->samples[1].neff);
Print(s);
Print(out);
}
}
}
goto refine;
}
todo |= 3;
@ -316,15 +339,15 @@ refine:
}
}
if( can_adjust ) {
for( comp = 0; comp < t->ncomp; ++comp )
totals[comp].maxerrsq -=
Sq(region->result[comp].spread*t->samples[1].weight);
}
if( can_adjust )
for( res = RegionResult(region), tot = state->totals;
tot < Tot; ++res, ++tot )
tot->maxerrsq -= Sq(res->spread/t->samples[1].neff);
switch( todo ) {
case 1: /* get spread right */
Explore(t, iregion, &t->samples[1], 0, 2);
region->isamples = 1;
ExploreSerial(t, state->iregion);
break;
case 3: /* sample region again with more points */
@ -334,144 +357,142 @@ refine:
SamplesAlloc(t, &t->samples[2]);
}
t->phase = 3;
t->samples[2].sampler(t, &t->samples[2], region->bounds, region->vol);
Explore(t, iregion, &t->samples[2], 0, 2);
region->isamples = 2;
t->samples[2].sampler(t, state->iregion);
ExploreSerial(t, state->iregion);
++region->depth; /* misused for df here */
++df;
}
++region->depth; /* misused for df here */
if( VERBOSE > 2 ) {
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
p += sprintf(p,
(dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" :
"\n (" REALF ") - (" REALF ")",
b->lower, b->upper);
cchar *msg = "\nRegion (" REALF ") - (" REALF ")";
for( B = (b = region->bounds) + t->ndim; b < B; ++b ) {
oe += sprintf(oe, msg, b->lower, b->upper);
msg = "\n (" REALF ") - (" REALF ")";
}
}
for( comp = 0; comp < t->ncomp; ++comp ) {
Result *r = &region->result[comp];
creal x1 = t->samples[0].avg[comp];
creal s1 = Var(t->samples[0]);
creal x2 = t->samples[1].avg[comp];
creal s2 = Var(t->samples[1]);
creal r2 = (s1 == 0) ? Sq(t->samples[1].neff*t->samples[0].weight) : s2/s1;
for( tot = state->totals, res = RegionResult(region), comp = 0;
tot < Tot; ++tot, ++res ) {
creal x1 = tot->phase[0].avg;
creal v1 = Var(0);
creal x2 = tot->phase[1].avg;
creal v2 = Var(1);
creal r2 = v1 ? v2/v1 :
Sq(t->samples[1].neff/(real)t->samples[0].neff);
real norm = 1 + r2;
real avg = x2 + r2*x1;
real sigsq = s2;
real sigsq = v2;
real chisq = Sq(x2 - x1);
real chiden = s1 + s2;
real chiden = v1 + v2;
if( todo == 3 ) {
creal x3 = t->samples[2].avg[comp];
creal s3 = Var(t->samples[2]);
creal r3 = (s2 == 0) ? Sq(t->samples[2].neff*t->samples[1].weight) : s3/s2;
creal x3 = res->avg;
creal v3 = Var2(2, res);
creal r3 = v2 ? v3/v2 :
Sq(t->samples[2].neff/(real)t->samples[1].neff);
norm = 1 + r3*norm;
avg = x3 + r3*avg;
sigsq = s3;
chisq = s1*Sq(x3 - x2) + s2*Sq(x3 - x1) + s3*chisq;
chiden = s1*s2 + s3*chiden;
sigsq = v3;
chisq = v1*Sq(x3 - x2) + v2*Sq(x3 - x1) + v3*chisq;
chiden = v1*v2 + v3*chiden;
}
avg = LAST ? r->avg : (sigsq *= norm = 1/norm, avg*norm);
avg = LAST ? res->avg : (sigsq *= norm = 1/norm, avg*norm);
if( chisq > EPS ) chisq /= Max(chiden, NOTZERO);
#define Out(s) s.avg[comp], r->spread*s.weight, s.err[comp]
if( VERBOSE > 2 ) {
p += sprintf(p, "\n[" COUNT "] "
#define Out2(f, r) (r)->avg, res->spread/t->samples[f].neff, (r)->err
#define Out(f) Out2(f, &tot->phase[f])
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL "(" REAL ")\n "
REAL " +- " REAL "(" REAL ")",
comp + 1, Out(t->samples[0]), Out(t->samples[1]));
if( todo == 3 ) p += sprintf(p, "\n "
REAL " +- " REAL "(" REAL ")",
Out(t->samples[2]));
p += sprintf(p, " \tchisq " REAL, chisq);
REAL " +- " REAL "(" REAL ")", ++comp, Out(0), Out(1));
if( todo == 3 ) oe += sprintf(oe, "\n "
REAL " +- " REAL "(" REAL ")", Out2(2, res));
oe += sprintf(oe, " \tchisq " REAL, chisq);
}
integral[comp] += avg;
error[comp] += sigsq;
prob[comp] += chisq;
tot->integral += avg;
tot->sigsq += sigsq;
tot->chisq += chisq;
r->avg = avg;
r->spread = sqrt(sigsq);
r->chisq = chisq;
res->avg = avg;
res->spread = sqrt(sigsq);
res->chisq = chisq;
}
if( VERBOSE > 2 ) Print(s);
}
if( VERBOSE > 2 ) Print(out);
++state->iregion;
for( comp = 0; comp < t->ncomp; ++comp )
error[comp] = sqrt(error[comp]);
WriteState(t);
}
df += t->nregions;
if( VERBOSE > 2 ) {
char s[16 + 128*NCOMP], *p = s;
p += sprintf(p, "\nTotals:");
for( comp = 0; comp < t->ncomp; ++comp )
p += sprintf(p, "\n[" COUNT "] "
REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)",
comp + 1, integral[comp], error[comp], prob[comp], df);
Print(s);
for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
integral[comp] = tot->integral;
error[comp] = sqrt(tot->sigsq);
prob[comp] = ChiSquare(tot->chisq, df);
}
for( comp = 0; comp < t->ncomp; ++comp )
prob[comp] = ChiSquare(prob[comp], df);
if( VERBOSE > 2 ) {
char *oe = out + sprintf(out, "\nTotals:");
for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp )
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)",
comp + 1, integral[comp], error[comp], tot->chisq, df);
Print(out);
}
weight = 1;
ML_ONLY(neff = 1;)
}
#ifdef MLVERSION
if( REGIONS ) {
Vector(real, bounds, t->ndim*2);
MLPutFunction(stdlink, "List", 2);
MLPutFunction(stdlink, "List", t->nregions);
for( iregion = 0; iregion < t->nregions; ++iregion ) {
Region *region = RegionPtr(iregion);
cBounds *b = region->bounds;
real lower[NDIM], upper[NDIM];
cResult *Res;
real *d = bounds;
for( dim = 0; dim < t->ndim; ++dim ) {
lower[dim] = b[dim].lower;
upper[dim] = b[dim].upper;
for( B = (b = region->bounds) + t->ndim; b < B; ++b ) {
*d++ = b->lower;
*d++ = b->upper;
}
MLPutFunction(stdlink, "Cuba`Divonne`region", 4);
MLPutRealList(stdlink, lower, t->ndim);
MLPutRealList(stdlink, upper, t->ndim);
MLPutRealList(stdlink, bounds, 2*t->ndim);
MLPutFunction(stdlink, "List", t->ncomp);
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
real res[] = {r->avg, r->spread*weight, r->chisq};
MLPutRealList(stdlink, res, Elements(res));
for( Res = (res = RegionResult(region)) + t->ncomp; res < Res; ++res ) {
real r[] = {res->avg, res->spread/neff, res->chisq};
MLPutRealList(stdlink, r, Elements(r));
}
MLPutInteger(stdlink, region->depth); /* misused for df */
MLPutInteger(stdlink, region->depth + 1); /* misused for df */
}
}
#endif
abort:
WaitCores(t);
FORK_ONLY(FrameFree(t, ShmRm(t));)
RuleFree(t);
SamplesFree(&t->samples[2]);
SamplesFree(&t->samples[1]);
SamplesFree(&t->samples[0]);
RuleFree(&t->rule13);
RuleFree(&t->rule11);
RuleFree(&t->rule9);
RuleFree(&t->rule7);
free(t->region);
free(t->xgiven);
free(t->voidregion);
StateRemove(t);
return fail;
}

132
src/external/libCuba/src/divonne/Iterate.c vendored Normal file
View File

@ -0,0 +1,132 @@
/*
Iterate.c
recursion over regions
this file is part of Divonne
last modified 2 Aug 13 th
*/
static void Iterate(This *t, count iregion, cint depth, cint isamples,
Totals *totals)
{
csize_t regionsize = RegionSize;
Region *parent, *region;
typedef struct {
real avg, err, spread, spreadsq;
} Corr;
Vector(Corr, corr, NCOMP);
Corr *c, *C = corr + t->ncomp;
Result *res;
count ireg, mreg = iregion;
count comp, maxsplit;
int last, idest, isrc;
region = RegionPtr(iregion);
region->depth = depth;
region->next = -iregion - 1;
if( isamples < 0 ) Split(t, iregion);
else {
region->isamples = isamples;
ExploreSerial(t, iregion);
}
ireg = iregion + RegionPtr(iregion)->next;
do {
region = RegionPtr(ireg);
if( region->depth > 0 ) {
--region->depth;
FORK_ONLY(more:)
ireg = Explore(t, ireg);
if( ireg == -1 ) return;
region = RegionPtr(ireg);
}
if( region->depth < 0 ) mreg = IMax(mreg, ireg);
ireg += region->next;
} while( ireg > 0 );
FORK_ONLY(if( t->running ) goto more;)
maxsplit = 1;
for( ireg = mreg; ireg >= iregion; --ireg ) {
parent = RegionPtr(ireg);
maxsplit -= NegQ(parent->depth);
if( parent->depth < 0 ) {
count xreg;
struct {
count from, to;
} todo[maxsplit], *tdmax = todo, *td;
count nsplit = 0;
real norm;
FClear(corr);
tdmax->from = ireg + parent->next;
tdmax->to = tdmax->from - parent->depth;
++tdmax;
for( td = todo; td < tdmax; ++td ) {
for( xreg = td->from; xreg < td->to; ++xreg ) {
Region *region = RegionPtr(xreg);
if( region->depth < 0 ) {
tdmax->from = xreg + region->next;
tdmax->to = tdmax->from - region->depth;
++tdmax;
}
else {
++nsplit;
for( res = RegionResult(region), c = corr; c < C; ++res, ++c ) {
c->avg += res->avg;
c->err += res->err;
c->spread += Sq(res->spread);
}
}
}
}
norm = 1./nsplit--;
for( res = RegionResult(parent), c = corr; c < C; ++res, ++c ) {
creal diff = fabs(res->avg - c->avg)*norm;
c->avg = diff*norm*nsplit;
c->err = (c->err == 0) ? 1 : 1 + diff/c->err;
c->spread = (c->spread == 0) ? 1 : 1 + diff/sqrt(c->spread);
}
for( td = todo; td < tdmax; ++td )
for( xreg = td->from; xreg < td->to; ++xreg ) {
Region *region = RegionPtr(xreg);
if( region->depth >= 0 ) {
cnumber neff = t->samples[region->isamples].neff;
for( res = RegionResult(region), c = corr; c < C; ++res, ++c ) {
if( res->err > 0 ) res->err = res->err*c->err + c->avg;
res->spread = res->spread*c->spread + c->avg*neff;
c->spreadsq += Sq(res->spread);
}
}
}
}
}
if( totals )
for( comp = 0; comp < t->ncomp; ++comp )
totals[comp].spreadsq += corr[comp].spreadsq;
for( last = -1, idest = isrc = iregion; iregion <= mreg; ++iregion ) {
Region *region = RegionPtr(iregion);
cint cur = NegQ(region->depth);
switch( cur - last ) {
case -1:
memmove(RegionPtr(idest), RegionPtr(isrc),
(iregion - isrc)*regionsize);
idest += iregion - isrc;
break;
case 1:
isrc = iregion;
}
last = cur;
}
memmove(RegionPtr(idest), RegionPtr(iregion),
(t->nregions - iregion)*regionsize);
t->nregions += idest - iregion;
}

20
src/external/libCuba/src/divonne/KorobovCoeff.c vendored
View File

@ -1,23 +1,3 @@
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define KOROBOV_MINDIM 2
#define KOROBOV_MAXDIM 33
#define MAXPRIME 9689

13
src/external/libCuba/src/divonne/Makefile.am vendored Normal file
View File

@ -0,0 +1,13 @@
## Process this file with automake to create Makefile.in
## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $
c_sources = Divonne.c
AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
noinst_LTLIBRARIES = libdivonne.la
libdivonne_la_SOURCES = $(c_sources)
libdivonne_la_LDFLAGS = $(AM_LDFLAGS)

207
src/external/libCuba/src/divonne/Rule.c vendored
View File

@ -4,37 +4,11 @@
code lifted with minor modifications from DCUHRE
by J. Berntsen, T. Espelid, and A. Genz
this file is part of Divonne
last modified 16 Jun 10 th
last modified 26 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
enum { nrules = 5 };
#define TYPEDEFSET \
typedef struct { \
count n; \
real weight[nrules], scale[nrules], norm[nrules]; \
real gen[NDIM]; \
} Set
#define NextSet(p) p = (Set *)((char *)p + setsize)
/*********************************************************************/
@ -83,86 +57,84 @@ static void Rule13Alloc(This *t)
enum { nsets = 14, ndim = 2 };
TYPEDEFSET;
count n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
Copy(last->weight, w[0], nrules);
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[1], nrules);
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[2], nrules);
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[3], nrules);
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[4], nrules);
last->gen[0] = g[3];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[5], nrules);
last->gen[0] = g[4];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[6], nrules);
last->gen[0] = g[5];
last->gen[1] = g[5];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[7], nrules);
last->gen[0] = g[6];
last->gen[1] = g[6];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[8], nrules);
last->gen[0] = g[7];
last->gen[1] = g[7];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[9], nrules);
last->gen[0] = g[8];
last->gen[1] = g[8];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[10], nrules);
last->gen[0] = g[9];
last->gen[1] = g[9];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
Copy(last->weight, w[11], nrules);
last->gen[0] = g[10];
last->gen[1] = g[11];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
Copy(last->weight, w[12], nrules);
last->gen[0] = g[12];
last->gen[1] = g[13];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
Copy(last->weight, w[13], nrules);
last->gen[0] = g[14];
@ -175,12 +147,12 @@ static void Rule13Alloc(This *t)
t->rule13.errcoeff[2] = 5;
t->rule13.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -231,84 +203,82 @@ static void Rule11Alloc(This *t)
enum { nsets = 13, ndim = 3 };
TYPEDEFSET;
count n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
Copy(last->weight, w[0], nrules);
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[1], nrules);
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[2], nrules);
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[3], nrules);
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[4], nrules);
last->gen[0] = g[3];
++last;
NextSet(last);
n += last->n = 2*ndim;
Copy(last->weight, w[5], nrules);
last->gen[0] = g[4];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[6], nrules);
last->gen[0] = g[5];
last->gen[1] = g[5];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
Copy(last->weight, w[7], nrules);
last->gen[0] = g[6];
last->gen[1] = g[6];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
Copy(last->weight, w[8], nrules);
last->gen[0] = g[7];
last->gen[1] = g[7];
last->gen[2] = g[7];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
Copy(last->weight, w[9], nrules);
last->gen[0] = g[8];
last->gen[1] = g[8];
last->gen[2] = g[8];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
Copy(last->weight, w[10], nrules);
last->gen[0] = g[9];
last->gen[1] = g[9];
last->gen[2] = g[9];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2);
Copy(last->weight, w[11], nrules);
last->gen[0] = g[10];
last->gen[1] = g[11];
last->gen[2] = g[11];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2);
Copy(last->weight, w[12], nrules);
last->gen[0] = g[12];
@ -322,12 +292,12 @@ static void Rule11Alloc(This *t)
t->rule11.errcoeff[2] = 3;
t->rule11.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -368,15 +338,13 @@ static void Rule9Alloc(This *t)
enum { nsets = 9 };
TYPEDEFSET;
ccount ndim = t->ndim;
ccount twondim = 1 << ndim;
count dim, n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
@ -386,7 +354,7 @@ static void Rule9Alloc(This *t)
last->weight[3] = ndim*(ndim*w[9] + w[10]) - 1 + last->weight[0];
last->weight[4] = ndim*w[11] + 1 - last->weight[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = ndim*(ndim*w[12] + w[13]) + w[14];
last->weight[1] = ndim*(ndim*w[15] + w[16]) + w[17];
@ -395,7 +363,7 @@ static void Rule9Alloc(This *t)
last->weight[4] = w[21] - last->weight[0];
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = ndim*w[22] + w[23];
last->weight[1] = ndim*w[24] + w[25];
@ -404,7 +372,7 @@ static void Rule9Alloc(This *t)
last->weight[4] = -last->weight[0];
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = w[29];
last->weight[1] = w[30];
@ -413,12 +381,12 @@ static void Rule9Alloc(This *t)
last->weight[4] = -w[29];
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[2] = w[32];
last->gen[0] = g[3];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
last->weight[0] = w[33] - ndim*w[12];
last->weight[1] = w[34] - ndim*w[15];
@ -428,7 +396,7 @@ static void Rule9Alloc(This *t)
last->gen[0] = g[0];
last->gen[1] = g[0];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1);
last->weight[0] = w[36];
last->weight[1] = w[37];
@ -438,7 +406,7 @@ static void Rule9Alloc(This *t)
last->gen[0] = g[0];
last->gen[1] = g[1];
++last;
NextSet(last);
n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
last->weight[0] = w[39];
last->weight[1] = w[40];
@ -449,7 +417,7 @@ static void Rule9Alloc(This *t)
last->gen[1] = g[0];
last->gen[2] = g[0];
++last;
NextSet(last);
n += last->n = twondim;
last->weight[0] = w[41]/twondim;
last->weight[1] = w[7]/twondim;
@ -466,12 +434,12 @@ static void Rule9Alloc(This *t)
t->rule9.errcoeff[2] = 5;
t->rule9.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -501,15 +469,13 @@ static void Rule7Alloc(This *t)
enum { nsets = 6 };
TYPEDEFSET;
ccount ndim = t->ndim;
ccount twondim = 1 << ndim;
count dim, n, r;
Set *first, *last, *s, *x;
csize_t setsize = SetSize;
Alloc(first, nsets);
Clear(first, nsets);
Die(first = calloc(nsets, setsize));
last = first;
n = last->n = 1;
@ -519,7 +485,7 @@ static void Rule7Alloc(This *t)
last->weight[3] = ndim*(ndim*w[7] + w[8]) - w[9];
last->weight[4] = 1 - last->weight[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = w[10];
last->weight[1] = w[11];
@ -528,7 +494,7 @@ static void Rule7Alloc(This *t)
last->weight[4] = -w[10];
last->gen[0] = g[1];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[0] = w[13] - ndim*w[0];
last->weight[1] = w[14] - ndim*w[3];
@ -537,12 +503,12 @@ static void Rule7Alloc(This *t)
last->weight[4] = -last->weight[0];
last->gen[0] = g[0];
++last;
NextSet(last);
n += last->n = 2*ndim;
last->weight[2] = w[17];
last->gen[0] = g[2];
++last;
NextSet(last);
n += last->n = 2*ndim*(ndim - 1);
last->weight[0] = -w[7];
last->weight[1] = w[18];
@ -552,7 +518,7 @@ static void Rule7Alloc(This *t)
last->gen[0] = g[0];
last->gen[1] = g[0];
++last;
NextSet(last);
n += last->n = twondim;
last->weight[0] = w[20]/twondim;
last->weight[1] = w[5]/twondim;
@ -569,12 +535,12 @@ static void Rule7Alloc(This *t)
t->rule7.errcoeff[2] = 5;
t->rule7.n = n;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
for( r = 1; r < nrules - 1; ++r ) {
creal scale = (s->weight[r] == 0) ? 100 :
-s->weight[r + 1]/s->weight[r];
real sum = 0;
for( x = first; x <= last; ++x )
for( x = first; x <= last; NextSet(x) )
sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
s->scale[r] = scale;
s->norm[r] = 1/sum;
@ -583,23 +549,33 @@ static void Rule7Alloc(This *t)
/*********************************************************************/
static inline void RuleIni(Rule *rule)
static inline void RuleAlloc(This *t)
{
rule->first = NULL;
if( (t->ndim - 2) | (t->key1 - 13)*(t->key2 - 13)*(t->key3 - 13) )
t->rule13.first = NULL;
else Rule13Alloc(t);
if( (t->ndim - 3) | (t->key1 - 11)*(t->key2 - 11)*(t->key3 - 11) )
t->rule11.first = NULL;
else Rule11Alloc(t);
if( (t->key1 - 9)*(t->key2 - 9)*(t->key3 - 9) )
t->rule9.first = NULL;
else Rule9Alloc(t);
if( (t->key1 - 7)*(t->key2 - 7)*(t->key3 - 7) )
t->rule7.first = NULL;
else Rule7Alloc(t);
}
/*********************************************************************/
static inline bool RuleIniQ(Rule *rule)
static inline void RuleFree(cThis *t)
{
return rule->first == NULL;
}
/*********************************************************************/
static inline void RuleFree(Rule *rule)
{
free(rule->first);
free(t->rule7.first);
free(t->rule9.first);
free(t->rule11.first);
free(t->rule13.first);
}
/*********************************************************************/
@ -653,29 +629,28 @@ next:
/*********************************************************************/
static void SampleRule(This *t, cSamples *samples, cBounds *b, creal vol)
static void SampleRule(This *t, ccount iregion)
{
TYPEDEFSET;
real *x = samples->x, *f = samples->f;
Set *first = (Set *)samples->rule->first;
Set *last = (Set *)samples->rule->last;
SAMPLERDEFS;
Set *first = samples->rule->first;
Set *last = samples->rule->last;
Set *s;
creal *errcoeff = samples->rule->errcoeff;
count comp, rul, n;
count comp, rul, sn;
csize_t setsize = SetSize;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
if( s->n ) x = ExpandFS(t, b, s->gen, x);
DoSample(t, samples->n, t->ndim, samples->x, f);
DoSample(t, n, samples->x, f);
for( comp = 0; comp < t->ncomp; ++comp ) {
real sum[nrules];
creal *f1 = f++;
Zap(sum);
for( s = first; s <= last; ++s )
for( n = s->n; n; --n ) {
for( s = first; s <= last; NextSet(s) )
for( sn = s->n; sn; --sn ) {
creal fun = *f1;
f1 += t->ncomp;
for( rul = 0; rul < nrules; ++rul )
@ -689,17 +664,17 @@ static void SampleRule(This *t, cSamples *samples, cBounds *b, creal vol)
for( rul = 1; rul < nrules - 1; ++rul ) {
real maxerr = 0;
for( s = first; s <= last; ++s )
for( s = first; s <= last; NextSet(s) )
maxerr = Max(maxerr,
fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
sum[rul] = maxerr;
}
samples->avg[comp] = vol*sum[0];
samples->err[comp] = vol*(
res[comp].avg = region->vol*sum[0];
res[comp].err = region->vol*(
(errcoeff[0]*sum[1] <= sum[2] && errcoeff[0]*sum[2] <= sum[3]) ?
errcoeff[1]*sum[1] :
errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3]));
errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3]) );
}
}

102
src/external/libCuba/src/divonne/Sample.c vendored
View File

@ -2,29 +2,9 @@
Sample.c
most of what is related to sampling
this file is part of Divonne
last modified 16 Jun 10 th
last modified 30 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define MARKMASK 0xfffffff
#define Marked(x) ((x) & ~MARKMASK)
@ -38,6 +18,7 @@
static inline void SamplesIni(Samples *samples)
{
samples->x = NULL;
samples->n = 0;
}
/*********************************************************************/
@ -56,11 +37,11 @@ static inline void SamplesFree(cSamples *samples)
/*********************************************************************/
static void SampleSobol(This *t, cSamples *samples, cBounds *b, creal vol)
static void SampleSobol(This *t, ccount iregion)
{
creal norm = vol*samples->weight;
real *x = samples->x, *f = samples->f, *avg = samples->avg;
cnumber n = samples->n;
SAMPLERDEFS;
Vector(real, avg, NCOMP);
real norm;
number i;
count dim, comp;
@ -70,35 +51,40 @@ static void SampleSobol(This *t, cSamples *samples, cBounds *b, creal vol)
*x = b[dim].lower + *x*(b[dim].upper - b[dim].lower);
}
DoSample(t, n, t->ndim, samples->x, f);
t->nrand += n;
ResCopy(avg, f);
DoSample(t, n, samples->x, f);
FCopy(avg, f);
f += t->ncomp;
for( i = 1; i < n; ++i )
for( i = 2; i < n; ++i )
for( comp = 0; comp < t->ncomp; ++comp )
avg[comp] += *f++;
for( comp = 0; comp < t->ncomp; ++comp )
avg[comp] *= norm;
norm = region->vol/samples->neff;
for( comp = 0; comp < t->ncomp; ++comp ) {
res[comp].avg = norm*avg[comp];
res[comp].err = 0;
}
}
/*********************************************************************/
static void SampleKorobov(This *t, cSamples *samples, cBounds *b, creal vol)
static void SampleKorobov(This *t, ccount iregion)
{
creal norm = vol*samples->weight;
real *x = samples->x, *xlast = x + t->ndim;
real *f = samples->f, *flast = f + t->ncomp;
real *avg = samples->avg;
cnumber n = samples->n, neff = samples->neff;
number nextra = n, i;
SAMPLERDEFS;
real *xlast = x + t->ndim, *flast = f + t->ncomp;
Vector(real, avg, NCOMP);
real norm;
cnumber neff = samples->neff;
number nextra = 0, i;
real dist = 0;
count dim, comp;
for( i = 1; i < n; ++i ) {
number c = i;
for( dim = 0; dim < t->ndim; ++dim ) {
creal dx = abs(2*c - neff)*samples->weight;
creal dx = abs(2*c - neff)/(real)neff;
*xlast++ = b[dim].lower + dx*(b[dim].upper - b[dim].lower);
c = c*samples->coeff % neff;
}
@ -119,26 +105,29 @@ static void SampleKorobov(This *t, cSamples *samples, cBounds *b, creal vol)
}
xlast[dim] = x2;
}
++nextra;
nextra = 1;
}
DoSample(t, nextra, t->ndim, x, f);
DoSample(t, n + nextra, x, f);
ResCopy(avg, flast);
FCopy(avg, flast);
flast += t->ncomp;
for( i = 2; i < n; ++i )
for( comp = 0; comp < t->ncomp; ++comp )
avg[comp] += *flast++;
if( nextra > n ) {
if( nextra ) {
for( comp = 0; comp < t->ncomp; ++comp )
f[comp] += dist*(f[comp] - flast[comp]);
for( dim = 0; dim < t->ndim; ++dim )
x[dim] = b[dim].upper;
}
for( comp = 0; comp < t->ncomp; ++comp )
avg[comp] = (avg[comp] + avg[comp] + f[comp])*norm;
norm = region->vol/samples->neff;
for( comp = 0; comp < t->ncomp; ++comp ) {
res[comp].avg = norm*(avg[comp] + avg[comp] + f[comp]);
res[comp].err = 0;
}
}
/*********************************************************************/
@ -165,25 +154,21 @@ static count SamplesLookup(This *t, Samples *samples, cint key,
number n;
if( key == 13 && t->ndim == 2 ) {
if( RuleIniQ(&t->rule13) ) Rule13Alloc(t);
samples->rule = &t->rule13;
samples->n = n = nmin = t->rule13.n;
samples->sampler = SampleRule;
}
else if( key == 11 && t->ndim == 3 ) {
if( RuleIniQ(&t->rule11) ) Rule11Alloc(t);
samples->rule = &t->rule11;
samples->n = n = nmin = t->rule11.n;
samples->sampler = SampleRule;
}
else if( key == 9 ) {
if( RuleIniQ(&t->rule9) ) Rule9Alloc(t);
samples->rule = &t->rule9;
samples->n = n = nmin = t->rule9.n;
samples->sampler = SampleRule;
}
else if( key == 7 ) {
if( RuleIniQ(&t->rule7) ) Rule7Alloc(t);
samples->rule = &t->rule7;
samples->n = n = nmin = t->rule7.n;
samples->sampler = SampleRule;
@ -234,21 +219,18 @@ static void SamplesAlloc(cThis *t, Samples *samples)
Alloc(samples->x, nx + nf + t->ncomp + t->ncomp);
samples->f = samples->x + nx;
samples->avg = samples->f + nf;
samples->err = samples->avg + t->ncomp;
ResClear(samples->err);
samples->weight = 1./samples->neff;
}
/*********************************************************************/
static real Sample(This *t, creal *x0)
{
real xtmp[2*NDIM], ftmp[2*NCOMP], *xlast = xtmp, f;
Vector(real, xtmp, 2*NDIM);
Vector(real, ftmp, 2*NCOMP);
real *xlast = xtmp, f;
real dist = 0;
count dim, comp;
number nextra = 1;
number n = 1;
for( dim = 0; dim < t->ndim; ++dim ) {
creal x1 = *xlast++ = Min(Max(*x0++, 0.), 1.);
@ -268,14 +250,16 @@ static real Sample(This *t, creal *x0)
}
*xlast++ = x2;
}
nextra = 2;
n = 2;
}
DoSample(t, nextra, t->ndim, xtmp, ftmp);
DoSample(t, n, xtmp, ftmp);
#define fin(x) Min(Max(x, -.5*INFTY), .5*INFTY)
comp = Untag(t->selectedcomp);
f = ftmp[comp];
if( nextra > 1 ) f += dist*(f - ftmp[comp + t->ncomp]);
f = fin(ftmp[comp]);
if( n > 1 ) f += dist*(f - fin(ftmp[comp + t->ncomp]));
return Sign(t->selectedcomp)*f;
}

179
src/external/libCuba/src/divonne/Split.c vendored
View File

@ -2,33 +2,13 @@
Split.c
determine optimal cuts for splitting a region
this file is part of Divonne
last modified 20 Jul 10 th
last modified 31 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define BNDTOL .05
#define FRACT .5
#define BIG 1e10
#define SMALL 1e-10
#define SINGTOL 1e-4
#define LHSTOL .1
@ -47,26 +27,21 @@ typedef struct {
real lhs, row, sol;
} Cut;
typedef struct {
real diff, err, spread;
} Errors;
typedef const Errors cErrors;
/*********************************************************************/
static inline real Div(creal a, creal b)
{
return (b != 0 && fabs(b) < BIG*fabs(a)) ? a/b : a;
return (b != 0 /*&& fabs(a) > SMALL*fabs(b)*/) ? a/b : a;
}
/*********************************************************************/
static void SomeCut(This *t, Cut *cut, Bounds *b)
{
Vector(real, xmid, NDIM);
real ymid, maxdev;
count dim, maxdim;
static count nextdim = 0;
real xmid[NDIM], ymid, maxdev;
for( dim = 0; dim < t->ndim; ++dim )
xmid[dim] = .5*(b[dim].upper + b[dim].lower);
@ -111,10 +86,10 @@ static inline real Volume(cThis *t, creal *delta)
/*********************************************************************/
static inline real SetupEqs(Cut *cut, ccount ncut, real f)
static inline real SetupEqs(Cut *cut, ccount ncuts, real f)
{
real sqsum = 0;
Cut *c = &cut[ncut];
Cut *c = &cut[ncuts];
while( --c >= cut ) {
sqsum += Sq(c->lhs = f - c->f);
f = c->f;
@ -124,7 +99,7 @@ static inline real SetupEqs(Cut *cut, ccount ncut, real f)
/*********************************************************************/
static inline void SolveEqs(Cut *cut, count ncut,
static inline void SolveEqs(Cut *cut, count ncuts,
creal *delta, creal diff)
{
real last = 0;
@ -135,7 +110,7 @@ static inline void SolveEqs(Cut *cut, count ncut,
ccount dim = Dim(c->i);
c->row = r -=
Div(diff, (delta[Lower(dim)] + delta[Upper(dim)])*c->df);
if( --ncut == 0 ) break;
if( --ncuts == 0 ) break;
last += r*c->lhs;
}
@ -157,8 +132,8 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
{
cint sign = (fdiff < 0) ? -1 : 1;
count ncut = 0, icut;
real delta[2*NDIM];
count ncuts = 0, icut;
Vector(real, delta, 2*NDIM);
real gamma, fgamma, lhssq;
count dim, div;
@ -167,7 +142,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
creal xsave = xmajor[dim];
real dist = b->upper - xsave;
if( dist >= BNDTOL*(b->upper - b->lower) ) {
Cut *c = &cut[ncut++];
Cut *c = &cut[ncuts++];
c->i = Upper(dim);
c->save = dist;
xmajor[dim] += dist *= FRACT;
@ -182,7 +157,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
creal xsave = xmajor[dim];
real dist = xsave - b->lower;
if( dist >= BNDTOL*(b->upper - b->lower) ) {
Cut *c = &cut[ncut++];
Cut *c = &cut[ncuts++];
c->i = Lower(dim);
c->save = dist;
xmajor[dim] -= dist *= FRACT;
@ -192,7 +167,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
delta[Lower(dim)] = dist;
}
if( ncut == 0 ) {
if( ncuts == 0 ) {
SomeCut(t, cut, bounds);
return 1;
}
@ -201,7 +176,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
real mindiff = INFTY;
Cut *mincut = cut;
for( icut = 0; icut < ncut; ++icut ) {
for( icut = 0; icut < ncuts; ++icut ) {
Cut *c = &cut[icut];
creal diff = fabs(fmajor - c->f);
if( diff <= mindiff ) {
@ -215,30 +190,30 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
if( sign*(mincut->f - fgamma) < 0 ) break;
if( --ncut == 0 ) {
if( --ncuts == 0 ) {
SomeCut(t, cut, bounds);
return 1;
}
delta[mincut->i] = mincut->save;
memcpy(mincut, mincut + 1, (char *)&cut[ncut] - (char *)mincut);
memmove(mincut, mincut + 1, (char *)&cut[ncuts] - (char *)mincut);
}
for( icut = 0; icut < ncut; ++icut ) {
for( icut = 0; icut < ncuts; ++icut ) {
Cut *c = &cut[icut];
c->fold = c->f;
c->df = (c->f - fmajor)/delta[c->i];
}
lhssq = SetupEqs(cut, ncut, fgamma);
lhssq = SetupEqs(cut, ncuts, fgamma);
repeat:
SolveEqs(cut, ncut, delta, gamma*fdiff);
SolveEqs(cut, ncuts, delta, gamma*fdiff);
for( div = 1; div <= 16; div *= 4 ) {
real gammanew, lhssqnew;
for( icut = 0; icut < ncut; ++icut ) {
for( icut = 0; icut < ncuts; ++icut ) {
Cut *c = &cut[icut];
real *x = &xmajor[Dim(c->i)];
creal xsave = *x;
@ -250,7 +225,7 @@ repeat:
gammanew = Volume(t, delta)/vol;
fgamma = fmajor + (gammanew - 1)*fdiff;
lhssqnew = SetupEqs(cut, ncut, fgamma);
lhssqnew = SetupEqs(cut, ncuts, fgamma);
if( lhssqnew <= lhssq ) {
real fmax;
@ -259,12 +234,12 @@ repeat:
gamma = gammanew;
fmax = fabs(fgamma);
for( icut = 0; icut < ncut; ++icut ) {
for( icut = 0; icut < ncuts; ++icut ) {
Cut *c = &cut[icut];
creal dfmin = SINGTOL*c->df;
creal sol = c->sol/div;
real df = c->f - c->fold;
df = (fabs(sol) < BIG*fabs(df)) ? df/sol : 1;
df = (fabs(df) > SMALL*fabs(sol)) ? df/sol : 1;
c->df = (fabs(df) < fabs(dfmin)) ? dfmin : df;
fmax = Max(fmax, fabs(c->f));
c->fold = c->f;
@ -276,90 +251,62 @@ repeat:
}
}
for( icut = 0; icut < ncut; ++icut ) {
for( icut = 0; icut < ncuts; ++icut ) {
Cut *c = &cut[icut];
real *b = (real *)bounds + c->i;
c->save = *b;
*b = xmajor[Dim(c->i)] + SignedDelta(c->i);
}
return ncut;
return ncuts;
}
/*********************************************************************/
static void Split(This *t, count iregion, int depth)
static void Split(This *t, ccount iregion)
{
TYPEDEFREGION;
csize_t regionsize = RegionSize;
Region *region = RegionPtr(iregion);
Vector(Cut, cut, 2*NDIM);
Cut *c;
count ncuts, succ;
int depth;
real *b;
Cut cut[2*NDIM];
Errors errors[NCOMP];
count comp, ncut, nsplit, xregion, ireg, xreg;
real tmp;
{
Region *const region = RegionPtr(iregion);
t->selectedcomp = region->cutcomp;
t->neval_cut -= t->neval;
ncut = FindCuts(t, cut, region->bounds, region->vol,
(real *)region->result + region->xmajor, region->fmajor,
ncuts = FindCuts(t, cut, region->bounds, region->vol,
(real *)RegionResult(region) + region->xmajor, region->fmajor,
region->fmajor - region->fminor);
t->neval_cut += t->neval;
for( comp = 0; comp < t->ncomp; ++comp ) {
Errors *e = &errors[comp];
e->diff = region->result[comp].avg;
e->spread = e->err = 0;
}
}
xregion = t->nregions;
depth -= ncut;
if( Explore(t, iregion, &t->samples[0], depth, 1) ) {
Cut *c;
for( c = cut; ncut--; ++c ) {
real *b = (real *)RegionPtr(iregion)->bounds;
ccount c0 = c->i, c1 = c0 ^ 1;
creal tmp = b[c1];
b[c1] = b[c0];
b[c0] = c->save;
if( !Explore(t, iregion, &t->samples[0], depth++, ncut != 0) ) break;
if( ncut ) ((real *)RegionPtr(iregion)->bounds)[c1] = tmp;
}
}
nsplit = t->nregions - xregion + 1;
for( ireg = iregion, xreg = xregion; ireg < t->nregions; ireg = xreg++ ) {
cResult *result = RegionPtr(ireg)->result;
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &result[comp];
Errors *e = &errors[comp];
e->diff -= r->avg;
e->err += r->err;
e->spread += Sq(r->spread);
}
}
tmp = 1./nsplit;
for( comp = 0; comp < t->ncomp; ++comp ) {
Errors *e = &errors[comp];
e->diff = tmp*fabs(e->diff);
e->err = (e->err == 0) ? 1 : 1 + e->diff/e->err;
e->spread = (e->spread == 0) ? 1 : 1 + e->diff/sqrt(e->spread);
}
tmp = 1 - tmp;
for( ireg = iregion, xreg = xregion; ireg < t->nregions; ireg = xreg++ ) {
Result *result = RegionPtr(ireg)->result;
for( comp = 0; comp < t->ncomp; ++comp ) {
Result *r = &result[comp];
cErrors *e = &errors[comp];
creal c = tmp*e->diff;
if( r->err > 0 ) r->err = r->err*e->err + c;
r->spread = r->spread*e->spread + c*t->samples[0].neff;
}
depth = region->depth - ncuts;
EnlargeRegions(t, ++ncuts);
region = RegionPtr(iregion);
region->depth = -ncuts;
succ = iregion + region->next;
region->next = t->nregions - iregion;
b = (real *)region->bounds;
region = RegionPtr(t->nregions);
XCopy(region->bounds, b);
region->depth = IDim(depth) + 1;
region->next = 1;
region->isamples = 0;
for( c = cut; --ncuts; ++c ) {
ccount ci = c->i;
creal tmp = b[ci ^ 1];
b[ci ^ 1] = b[ci];
b[ci] = c->save;
region = RegionPtr(++t->nregions);
XCopy(region->bounds, b);
region->depth = IDim(depth) + 1;
region->next = 1;
region->isamples = 0;
++depth;
b[ci ^ 1] = tmp;
}
region->next = succ - t->nregions++;
}

76
src/external/libCuba/src/divonne/common.c vendored
View File

@ -2,47 +2,22 @@
common.c
includes most of the modules
this file is part of Divonne
last modified 8 Jun 10 th
last modified 29 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "Random.c"
#include "ChiSquare.c"
#include "Rule.c"
#include "Sample.c"
#include "FindMinimum.c"
static bool Explore(This *t, count iregion, cSamples *samples,
cint depth, cint flags);
static void Split(This *t, count iregion, int depth);
#include "Explore.c"
#include "Split.c"
#include "Explore.c"
#include "Iterate.c"
static inline bool BadDimension(cThis *t, ccount key)
{
if( t->ndim > NDIM ) return true;
if( t->ndim > MAXDIM ) return true;
if( IsSobol(key) ) return
t->ndim < SOBOL_MINDIM || (t->seed == 0 && t->ndim > SOBOL_MAXDIM);
if( IsRule(key, t->ndim) ) return t->ndim < 1;
@ -51,9 +26,48 @@ static inline bool BadDimension(cThis *t, ccount key)
static inline bool BadComponent(cThis *t)
{
if( t->ncomp > NCOMP ) return true;
if( t->ncomp > MAXCOMP ) return true;
return t->ncomp < 1;
}
#include "Integrate.c"
static inline void AllocGiven(This *t)
{
real *xgiven = NULL, *fgiven = NULL;
if( t->ngiven | t->nextra ) {
cnumber nxgiven = t->ngiven*t->ndim;
cnumber nxextra = t->nextra*t->ndim;
cnumber nfgiven = t->ngiven*t->ncomp;
cnumber nfextra = t->nextra*t->ncomp;
Alloc(xgiven, nxgiven + nxextra + nfgiven + nfextra);
t->xextra = xgiven + nxgiven;
fgiven = t->xextra + nxextra;
t->fextra = fgiven + nfgiven;
if( nxgiven ) {
#ifdef MLVERSION
Copy(xgiven, t->xgiven, nxgiven);
Copy(fgiven, t->fgiven, nfgiven);
#else
if( t->ldxgiven == t->ndim )
Copy(xgiven, t->xgiven, nxgiven);
else {
number i;
real *sgiven = t->xgiven, *dgiven = xgiven;
for( i = 0; i < t->ngiven; ++i ) {
Copy(dgiven, sgiven, t->ndim);
sgiven += t->ldxgiven;
dgiven += t->ndim;
}
}
t->phase = 0;
DoSample(t, t->ngiven, xgiven, fgiven);
#endif
}
}
t->xgiven = xgiven;
t->fgiven = fgiven;
}

125
src/external/libCuba/src/divonne/decl.h vendored
View File

@ -2,32 +2,16 @@
decl.h
Type declarations
this file is part of Divonne
last modified 8 Jun 10 th
last modified 26 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "stddecl.h"
#define INIDEPTH 3
#define DEPTH 5
#define POSTDEPTH 15
#define Tag(x) ((x) | INT_MIN)
#define Untag(x) ((x) & INT_MAX)
@ -37,15 +21,31 @@ typedef struct {
typedef const Bounds cBounds;
typedef struct {
real avg, err;
} PhaseResult;
typedef struct {
real avg, spreadsq;
real spread, secondspread;
real nneed, maxerrsq, mindevsq;
real integral, sigsq, chisq;
PhaseResult phase[2];
int iregion;
} Totals;
enum { nrules = 5 };
typedef struct {
void *first, *last;
count n;
real weight[nrules], scale[nrules], norm[nrules];
real gen[];
} Set;
#define SetSize (sizeof(Set) + t->ndim*sizeof(real))
typedef struct {
Set *first, *last;
real errcoeff[3];
count n;
} Rule;
@ -53,16 +53,45 @@ typedef struct {
typedef const Rule cRule;
typedef struct samples {
real weight;
real *x, *f, *avg, *err;
void (*sampler)(struct _this *t, const struct samples *, cBounds *, creal);
real *x, *f;
void (*sampler)(struct _this *t, ccount);
cRule *rule;
count coeff;
number n, neff;
count coeff;
} Samples;
typedef const Samples cSamples;
typedef struct {
real diff, err, spread;
} Errors;
typedef const Errors cErrors;
typedef struct {
real avg, err, spread, chisq;
real fmin, fmax;
real xminmax[];
} Result;
typedef const Result cResult;
#define ResultSize (sizeof(Result) + t->ndim*2*sizeof(real))
typedef struct region {
int depth, next;
count isamples, cutcomp, xmajor;
real fmajor, fminor, vol;
Bounds bounds[];
} Region;
#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*ResultSize)
#define RegionResult(r) ((Result *)(r->bounds + t->ndim))
#define RegionPtr(n) ((Region *)((char *)t->region + (n)*regionsize))
typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *, cint *);
typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *);
@ -73,6 +102,12 @@ typedef struct _this {
Integrand integrand;
void *userdata;
PeakFinder peakfinder;
#ifdef HAVE_FORK
int ncores, running, *child;
real *frame;
number nframe;
SHM_ONLY(int shmid;)
#endif
#endif
real epsrel, epsabs;
int flags, seed;
@ -82,37 +117,39 @@ typedef struct _this {
Bounds border;
real maxchisq, mindeviation;
number ngiven, nextra;
real *xgiven, *xextra, *fgiven, *fextra;
real *xgiven, *fgiven;
real *xextra, *fextra;
count ldxgiven;
count nregions;
number neval, neval_opt, neval_cut;
int phase;
cchar *statefile;
number neval, neval_opt, neval_cut, nrand;
count phase;
count selectedcomp, size;
Samples samples[3];
Totals *totals;
Rule rule7, rule9, rule11, rule13;
RNGState rng;
void *voidregion;
Region *region;
jmp_buf abort;
} This;
typedef const This cThis;
#define CHUNKSIZE 4096
#define TYPEDEFREGION \
typedef struct { \
real avg, err, spread, chisq; \
real fmin, fmax; \
real xmin[NDIM], xmax[NDIM]; \
} Result; \
typedef const Result cResult; \
typedef struct region { \
count cutcomp, depth, xmajor; \
real fmajor, fminor, vol; \
Bounds bounds[NDIM]; \
Result result[NCOMP]; \
} Region
#define AllocRegions(t) \
MemAlloc((t)->region, (t)->size*regionsize)
#define RegionPtr(n) (&((Region *)t->voidregion)[n])
#define EnlargeRegions(t, n) if( (t)->nregions + n > (t)->size ) \
ReAlloc((t)->region, ((t)->size += CHUNKSIZE)*regionsize)
#define SAMPLERDEFS \
csize_t regionsize = RegionSize; \
Region *region = RegionPtr(iregion); \
cBounds *b = region->bounds; \
Result *res = RegionResult(region); \
cSamples *samples = &t->samples[region->isamples]; \
real *x = samples->x, *f = samples->f; \
cnumber n = samples->n

51
src/external/libCuba/src/divonne/util.c vendored
View File

@ -1,51 +0,0 @@
/*
util.c
Utility functions
this file is part of Divonne
last modified 9 Apr 09 th
*/
/***************************************************************************
* Copyright (C) 2004-2009 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "decl.h"
static count ndim_, ncomp_, selectedcomp_, nregions_;
static number neval_, neval_opt_, neval_cut_;
static int sign_, phase_;
static Bounds border_;
static Samples samples_[3];
static Rule rule7_, rule9_, rule11_, rule13_;
static real *xgiven_, *fgiven_, *xextra_, *fextra_;
static count ldxgiven_;
static number ngiven_, nextra_;
static Totals *totals_;
static void *voidregion_;
#define region_ ((Region *)voidregion_)
static count size_;
#ifdef DEBUG
#include "debug.c"
#endif

46
src/external/libCuba/src/suave/Fluct.c vendored
View File

@ -2,48 +2,32 @@
Fluct.c
compute the fluctuation in the left and right half
this file is part of Suave
last modified 9 Feb 05 th
last modified 29 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL)
typedef long double xdouble;
typedef long double realx;
#define XDBL_MAX_EXP LDBL_MAX_EXP
#define XDBL_MAX LDBL_MAX
#define powx powl
#define ldexpx ldexpl
#else
typedef double xdouble;
typedef double realx;
#define XDBL_MAX_EXP DBL_MAX_EXP
#define XDBL_MAX DBL_MAX
#define powx pow
#define ldexpx ldexp
#endif
typedef const realx crealx;
typedef struct {
xdouble fluct;
realx fluct;
number n;
} Var;
@ -54,23 +38,23 @@ static void Fluct(cThis *t, Var *var,
{
creal *x = w + n;
creal *f = x + n*t->ndim + comp;
creal flat = 2/3./t->flatness;
creal max = ldexp(1., (int)((XDBL_MAX_EXP - 2)/t->flatness));
creal norm = 1/(err*Max(fabs(avg), err));
count nvar = 2*t->ndim;
creal norm = 1/(err*Max(fabs(avg), err));
creal flat = 2/3./t->flatness;
crealx max = ldexpx(1., (int)((XDBL_MAX_EXP - 2)/t->flatness));
Clear(var, nvar);
while( n-- ) {
count dim;
const xdouble ft =
powx(Min(1 + fabs(*w++)*Sq(*f - avg)*norm, max), t->flatness);
crealx arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
crealx ft = powx(arg < max ? arg : max, t->flatness);
f += t->ncomp;
for( dim = 0; dim < t->ndim; ++dim ) {
Var *v = &var[2*dim + (*x++ >= b[dim].mid)];
const xdouble f = v->fluct + ft;
Var *v = &var[2*dim + (*x++ >= .5*(b[dim].lower + b[dim].upper))];
crealx f = v->fluct + ft;
v->fluct = (f > XDBL_MAX/2) ? XDBL_MAX/2 : f;
++v->n;
}

29
src/external/libCuba/src/suave/Grid.c vendored
View File

@ -2,28 +2,9 @@
Grid.c
utility functions for the Vegas grid
this file is part of Suave
last modified 1 Jun 10 th
last modified 7 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
static void RefineGrid(cThis *t, Grid grid, Grid margsum)
{
@ -72,7 +53,7 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum)
delta = (cur - prev)*thisbin;
newgrid[newbin] = SHARPEDGES ?
cur - delta/imp[bin] :
(newcur = Max(newcur + 0x1p-48,
(newcur = Max(newcur + 16*DBL_EPSILON,
cur - 2*delta/(imp[bin] + imp[IDim(bin - 1)])));
}
Copy(grid, newgrid, NBINS - 1);
@ -84,8 +65,8 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum)
static void Reweight(cThis *t, Bounds *b,
creal *w, creal *f, creal *lastf, cResult *total)
{
Grid margsum[NDIM];
real scale[NCOMP];
Vector(Grid, margsum, NDIM);
Vector(real, scale, NCOMP);
cbin_t *bin = (cbin_t *)lastf;
count dim, comp;
@ -95,7 +76,7 @@ static void Reweight(cThis *t, Bounds *b,
scale[comp] = (total[comp].avg == 0) ? 0 : 1/total[comp].avg;
}
Zap(margsum);
XClear(margsum);
while( f < lastf ) {
real fsq = 0;

306
src/external/libCuba/src/suave/Integrate.c vendored
View File

@ -2,90 +2,122 @@
Integrate.c
integrate over the unit hypercube
this file is part of Suave
last modified 16 Jun 10 th
checkpointing by B. Chokoufe
last modified 5 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
typedef struct {
signature_t signature;
count nregions, df;
number neval;
Result totals[];
} State;
static int Integrate(This *t, real *integral, real *error, real *prob)
{
TYPEDEFREGION;
StateDecl;
csize_t statesize = sizeof(State) + NCOMP*sizeof(Result);
Sized(State, state, statesize);
Array(Var, var, NDIM, 2);
Vector(char, out, 128*NCOMP + 256);
count dim, comp, df;
int fail = -99;
Result totals[NCOMP];
Region *anchor = NULL, *region = NULL;
Result *tot, *Tot = state->totals + t->ncomp;
Result *res, *resL, *resR;
Bounds *b, *B;
count dim, comp;
int fail;
if( VERBOSE > 1 ) {
char s[256];
sprintf(s, "Suave input parameters:\n"
sprintf(out, "Suave input parameters:\n"
" ndim " COUNT "\n ncomp " COUNT "\n"
" epsrel " REAL "\n epsabs " REAL "\n"
" flags %d\n seed %d\n"
" mineval " NUMBER "\n maxeval " NUMBER "\n"
" nnew " NUMBER "\n flatness " REAL,
" nnew " NUMBER "\n flatness " REAL "\n"
" statefile \"%s\"\n",
t->ndim, t->ncomp,
t->epsrel, t->epsabs,
t->flags, t->seed,
t->mineval, t->maxeval,
t->nnew, t->flatness);
Print(s);
t->nnew, t->flatness,
t->statefile);
Print(out);
}
if( BadComponent(t) ) return -2;
if( BadDimension(t) ) return -1;
if( setjmp(t->abort) ) goto abort;
ShmAlloc(t, ShmRm(t));
ForkCores(t);
if( (fail = setjmp(t->abort)) ) goto abort;
t->epsabs = Max(t->epsabs, NOTZERO);
IniRandom(t);
RegionAlloc(anchor, t->nnew, t->nnew);
anchor->next = NULL;
anchor->div = 0;
StateSetup(t);
for( dim = 0; dim < t->ndim; ++dim ) {
Bounds *b = &anchor->bounds[dim];
b->lower = 0;
b->upper = 1;
b->mid = .5;
if( dim == 0 ) {
count bin;
/* define the initial distribution of bins */
for( bin = 0; bin < NBINS; ++bin )
b->grid[bin] = (bin + 1)/(real)NBINS;
}
else Copy(b->grid, anchor->bounds[0].grid, NBINS);
if( StateReadTest(t) ) {
size_t *regsize = NULL;
StateReadOpen(t, fd) {
count iregion;
size_t totsize;
Region **last = &anchor;
if( read(fd, state, statesize) != statesize ||
state->signature != StateSignature(t, 2) ) break;
t->nregions = state->nregions;
totsize = t->nregions*sizeof *regsize;
MemAlloc(regsize, totsize);
StateRead(fd, regsize, totsize);
for( iregion = 0; iregion < t->nregions; ++iregion )
totsize += regsize[iregion];
if( st.st_size != statesize + totsize ) break;
for( iregion = 0; iregion < t->nregions; ++iregion ) {
Region *reg;
MemAlloc(reg, regsize[iregion]);
StateRead(fd, reg, regsize[iregion]);
*last = reg;
last = &reg->next;
}
*last = NULL;
} StateReadClose(t, fd);
free(regsize);
t->neval = state->neval;
t->rng.skiprandom(t, t->neval);
}
Sample(t, t->nnew, anchor, anchor->w,
anchor->w + t->nnew,
anchor->w + t->nnew + t->ndim*t->nnew);
df = anchor->df;
ResCopy(totals, anchor->result);
if( ini ) {
count bin;
Bounds *b0;
for( t->nregions = 1; ; ++t->nregions ) {
Var var[NDIM][2], *vLR;
anchor = RegionAlloc(t, t->nnew, t->nnew);
anchor->next = NULL;
anchor->div = 0;
t->nregions = 1;
b0 = RegionBounds(anchor);
b0->lower = 0;
b0->upper = 1;
/* define the initial distribution of bins */
for( bin = 0; bin < NBINS; ++bin )
b0->grid[bin] = (bin + 1)/(real)NBINS;
for( b = b0 + 1, B = b0 + t->ndim; b < B; ++b )
Copy(b, b0, 1);
t->neval = 0;
Sample(t, t->nnew, anchor, RegionW(anchor),
RegionW(anchor) + t->nnew,
RegionW(anchor) + t->nnew + t->ndim*t->nnew);
state->df = anchor->df;
FCopy(state->totals, anchor->result);
}
/* main iteration loop */
for( ; ; ) {
Var *vLR;
real maxratio, maxerr, minfluct, bias, mid;
Region *regionL, *regionR, *reg, **parent, **par;
Bounds *bounds, *boundsL, *boundsR;
@ -94,39 +126,33 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
real *w, *wL, *wR, *x, *xL, *xR, *f, *fL, *fR, *wlast, *flast;
if( VERBOSE ) {
char s[128 + 128*NCOMP], *p = s;
p += sprintf(p, "\n"
char *oe = out + sprintf(out, "\n"
"Iteration " COUNT ": " NUMBER " integrand evaluations so far",
t->nregions, t->neval);
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *tot = &totals[comp];
p += sprintf(p, "\n[" COUNT "] "
for( tot = state->totals, comp = 0; tot < Tot; ++tot )
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)",
comp + 1, tot->avg, tot->err, tot->chisq, df);
}
Print(s);
++comp, tot->avg, tot->err, tot->chisq, state->df);
Print(out);
}
maxratio = -INFTY;
maxcomp = 0;
for( comp = 0; comp < t->ncomp; ++comp ) {
creal ratio = totals[comp].err/MaxErr(totals[comp].avg);
for( tot = state->totals, comp = 0; tot < Tot; ++tot ) {
creal ratio = tot->err/MaxErr(tot->avg);
if( ratio > maxratio ) {
maxratio = ratio;
maxcomp = comp;
}
}
if( maxratio <= 1 && t->neval >= t->mineval ) {
fail = 0;
if( maxratio <= 1 && t->neval >= t->mineval ) break;
if( t->neval >= t->maxeval ) {
fail = 1;
break;
}
if( t->neval >= t->maxeval ) break;
maxerr = -INFTY;
parent = &anchor;
region = anchor;
@ -139,8 +165,11 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
}
}
Fluct(t, var[0],
region->bounds, region->w, region->n, maxcomp,
bounds = RegionBounds(region);
w = RegionW(region);
/* find the bisectdim which minimizes the fluctuations */
Fluct(t, var[0], bounds, w, region->n, maxcomp,
region->result[maxcomp].avg, Max(maxerr, t->epsabs));
bias = (t->epsrel < 1e-50) ? 2 :
@ -148,15 +177,15 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
minfluct = INFTY;
bisectdim = 0;
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
creal fluct = (var[dim][0].fluct + var[dim][1].fluct)*
(bias - b->upper + b->lower);
(bias - bounds[dim].upper + bounds[dim].lower);
if( fluct < minfluct ) {
minfluct = fluct;
bisectdim = dim;
}
}
/* apply stratified sampling to distribute points in bisected region */
vLR = var[bisectdim];
minfluct = vLR[0].fluct + vLR[1].fluct;
nnewL = IMax(
@ -166,37 +195,35 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
nnewR = IMax(t->nnew - nnewL, MINSAMPLES);
nR = vLR[1].n + nnewR;
RegionAlloc(regionL, nL, nnewL);
RegionAlloc(regionR, nR, nnewR);
regionL = RegionAlloc(t, nL, nnewL);
regionR = RegionAlloc(t, nR, nnewR);
*parent = regionL;
regionL->next = regionR;
regionR->next = region->next;
regionL->div = regionR->div = region->div + 1;
bounds = &region->bounds[bisectdim];
mid = bounds->mid;
mid = .5*(bounds[bisectdim].lower + bounds[bisectdim].upper);
n = region->n;
w = wlast = region->w; x = w + n; f = flast = x + n*t->ndim;
wL = regionL->w; xL = wL + nL; fL = xL + nL*t->ndim;
wR = regionR->w; xR = wR + nR; fR = xR + nR*t->ndim;
wlast = w; x = w + n; f = flast = x + n*t->ndim;
wL = RegionW(regionL); xL = wL + nL; fL = xL + nL*t->ndim;
wR = RegionW(regionR); xR = wR + nR; fR = xR + nR*t->ndim;
while( n-- ) {
cbool final = (*w < 0);
if( x[bisectdim] < mid ) {
if( final && wR > regionR->w ) *(wR - 1) = -fabs(*(wR - 1));
if( final && wR > RegionW(regionR) ) wR[-1] = -fabs(wR[-1]);
*wL++ = *w++;
VecCopy(xL, x);
XCopy(xL, x);
xL += t->ndim;
ResCopy(fL, f);
FCopy(fL, f);
fL += t->ncomp;
}
else {
if( final && wL > regionL->w ) *(wL - 1) = -fabs(*(wL - 1));
if( final && wL > RegionW(regionL) ) wL[-1] = -fabs(wL[-1]);
*wR++ = *w++;
VecCopy(xR, x);
XCopy(xR, x);
xR += t->ndim;
ResCopy(fR, f);
FCopy(fR, f);
fR += t->ncomp;
}
x += t->ndim;
@ -204,80 +231,100 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
if( n && final ) wlast = w, flast = f;
}
Reweight(t, region->bounds, wlast, flast, f, totals);
VecCopy(regionL->bounds, region->bounds);
VecCopy(regionR->bounds, region->bounds);
Reweight(t, bounds, wlast, flast, f, state->totals);
boundsL = RegionBounds(regionL);
XCopy(boundsL, bounds);
boundsR = RegionBounds(regionR);
XCopy(boundsR, bounds);
boundsL = &regionL->bounds[bisectdim];
boundsR = &regionR->bounds[bisectdim];
boundsL->mid = .5*(boundsL->lower + (boundsL->upper = mid));
boundsR->mid = .5*((boundsR->lower = mid) + boundsR->upper);
StretchGrid(bounds->grid, boundsL->grid, boundsR->grid);
boundsL[bisectdim].upper = mid;
boundsR[bisectdim].lower = mid;
StretchGrid(bounds[bisectdim].grid,
boundsL[bisectdim].grid, boundsR[bisectdim].grid);
Sample(t, nnewL, regionL, wL, xL, fL);
Sample(t, nnewR, regionR, wR, xR, fR);
df += regionL->df + regionR->df - region->df;
state->df += regionL->df + regionR->df - region->df;
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
Result *rL = &regionL->result[comp];
Result *rR = &regionR->result[comp];
Result *tot = &totals[comp];
for( res = region->result,
resL = regionL->result,
resR = regionR->result,
tot = state->totals;
tot < Tot; ++res, ++resL, ++resR, ++tot ) {
real diff, sigsq;
tot->avg += diff = rL->avg + rR->avg - r->avg;
tot->avg += diff = resL->avg + resR->avg - res->avg;
diff = Sq(.25*diff);
sigsq = rL->sigsq + rR->sigsq;
sigsq = resL->sigsq + resR->sigsq;
if( sigsq > 0 ) {
creal c = Sq(1 + sqrt(diff/sigsq));
rL->sigsq *= c;
rR->sigsq *= c;
resL->sigsq *= c;
resR->sigsq *= c;
}
rL->err = sqrt(rL->sigsq += diff);
rR->err = sqrt(rR->sigsq += diff);
resL->err = sqrt(resL->sigsq += diff);
resR->err = sqrt(resR->sigsq += diff);
tot->sigsq += rL->sigsq + rR->sigsq - r->sigsq;
tot->sigsq += resL->sigsq + resR->sigsq - res->sigsq;
tot->err = sqrt(tot->sigsq);
tot->chisq += rL->chisq + rR->chisq - r->chisq;
tot->chisq += resL->chisq + resR->chisq - res->chisq;
}
free(region);
region = NULL;
++t->nregions;
if( StateWriteTest(t) ) {
StateWriteOpen(t, fd) {
Region *reg;
size_t totsize, *regsize, *s;
state->signature = StateSignature(t, 2);
state->neval = t->neval;
state->nregions = t->nregions;
StateWrite(fd, state, statesize);
MemAlloc(regsize, totsize = t->nregions*sizeof(size_t));
s = regsize;
for( reg = anchor; reg; reg = reg->next )
*s++ = reg->size;
StateWrite(fd, regsize, totsize);
free(regsize);
for( reg = anchor; reg; reg = reg->next )
StateWrite(fd, reg, reg->size);
} StateWriteClose(t, fd);
}
}
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *tot = &totals[comp];
for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
integral[comp] = tot->avg;
error[comp] = tot->err;
prob[comp] = ChiSquare(tot->chisq, df);
prob[comp] = ChiSquare(tot->chisq, state->df);
}
#ifdef MLVERSION
if( REGIONS ) {
Vector(real, bounds, NDIM*2);
MLPutFunction(stdlink, "List", 2);
MLPutFunction(stdlink, "List", t->nregions);
for( region = anchor; region; region = region->next ) {
real lower[NDIM], upper[NDIM];
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
lower[dim] = b->lower;
upper[dim] = b->upper;
cResult *Res;
real *d = bounds;
for( B = (b = RegionBounds(region)) + t->ndim; b < B; ++b ) {
*d++ = b->lower;
*d++ = b->upper;
}
MLPutFunction(stdlink, "Cuba`Suave`region", 4);
MLPutRealList(stdlink, lower, t->ndim);
MLPutRealList(stdlink, upper, t->ndim);
MLPutFunction(stdlink, "Cuba`Suave`region", 3);
MLPutRealList(stdlink, bounds, 2*t->ndim);
MLPutFunction(stdlink, "List", t->ncomp);
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
real res[] = {r->avg, r->err, r->chisq};
MLPutRealList(stdlink, res, Elements(res));
for( Res = (res = region->result) + t->ncomp; res < Res; ++res ) {
real r[] = {res->avg, res->err, res->chisq};
MLPutRealList(stdlink, r, Elements(r));
}
MLPutInteger(stdlink, region->df);
@ -287,11 +334,14 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
abort:
free(region);
while( (region = anchor) ) {
anchor = anchor->next;
free(region);
}
WaitCores(t);
ShmFree(t);
StateRemove(t);
return fail;
}

13
src/external/libCuba/src/suave/Makefile.am vendored Normal file
View File

@ -0,0 +1,13 @@
## Process this file with automake to create Makefile.in
## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $
c_sources = Suave.c
AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
noinst_LTLIBRARIES = libsuave.la
libsuave_la_SOURCES = $(c_sources)
libsuave_la_LDFLAGS = $(AM_LDFLAGS)

81
src/external/libCuba/src/suave/Sample.c vendored
View File

@ -2,29 +2,9 @@
Sample.c
the sampling step of Suave
this file is part of Suave
last modified 13 Sep 10 th
last modified 30 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
typedef struct {
real sum, sqsum;
@ -34,16 +14,16 @@ typedef struct {
/*********************************************************************/
static void Sample(This *t, cnumber nnew, void *voidregion,
static void Sample(This *t, cnumber nnew, Region *region,
real *lastw, real *lastx, real *lastf)
{
TYPEDEFREGION;
Region *const region = (Region *)voidregion;
count comp, dim, df;
count comp, df;
number n;
Cumulants cumul[NCOMP];
char **ss, *s;
Vector(Cumulants, cumul, NCOMP);
Cumulants *C = cumul + t->ncomp, *c;
Bounds *bounds = RegionBounds(region), *B = bounds + t->ndim, *b;
Result *res;
char **ss = NULL, *s = NULL;
ccount chars = 128*(region->div + 1);
creal jacobian = 1/ldexp((real)nnew, region->div);
@ -55,8 +35,7 @@ static void Sample(This *t, cnumber nnew, void *voidregion,
t->rng.getrandom(t, f);
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
for( b = bounds; b < B; ++b ) {
creal pos = *f*NBINS;
ccount ipos = (count)pos;
creal prev = (ipos == 0) ? 0 : b->grid[ipos - 1];
@ -69,11 +48,11 @@ static void Sample(This *t, cnumber nnew, void *voidregion,
*w++ = weight;
}
DoSample(t, nnew, lastw, lastx, lastf, region->div + 1);
DoSample(t, nnew, lastx, lastf, lastw, region->div + 1);
w[-1] = -w[-1];
lastw = w;
w = region->w;
w = RegionW(region);
region->n = lastw - w;
if( VERBOSE > 2 ) {
@ -87,7 +66,7 @@ static void Sample(This *t, cnumber nnew, void *voidregion,
}
}
Zap(cumul);
FClear(cumul);
df = n = 0;
while( w < lastw ) {
@ -95,9 +74,7 @@ static void Sample(This *t, cnumber nnew, void *voidregion,
creal weight = fabs(*w++);
++n;
for( comp = 0; comp < t->ncomp; ++comp ) {
Cumulants *c = &cumul[comp];
for( c = cumul, comp = 0; c < C; ++c ) {
creal wfun = weight*(*f++);
c->sum += wfun;
c->sqsum += Sq(wfun);
@ -135,23 +112,21 @@ static void Sample(This *t, cnumber nnew, void *voidregion,
region->df = --df;
for( comp = 0; comp < t->ncomp; ++comp ) {
Result *r = &region->result[comp];
Cumulants *c = &cumul[comp];
for( c = cumul, res = region->result; c < C; ++c, ++res ) {
creal sigsq = 1/c->weightsum;
creal avg = sigsq*c->avgsum;
if( LAST ) {
r->sigsq = 1/c->weight;
r->avg = r->sigsq*c->avg;
res->sigsq = 1/c->weight;
res->avg = res->sigsq*c->avg;
}
else {
r->sigsq = sigsq;
r->avg = avg;
res->sigsq = sigsq;
res->avg = avg;
}
r->err = sqrt(r->sigsq);
res->err = sqrt(res->sigsq);
r->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum;
res->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum;
/* This catches the special case where the integrand is constant
over the entire region. Unless that constant is zero, only the
first set of samples will have zero variance, and hence weight
@ -168,19 +143,17 @@ static void Sample(This *t, cnumber nnew, void *voidregion,
if( VERBOSE > 2 ) {
char *p = s;
char *p0 = p + t->ndim*64;
char *msg = "\nRegion (" REALF ") - (" REALF ")";
for( dim = 0; dim < t->ndim; ++dim ) {
cBounds *b = &region->bounds[dim];
p += sprintf(p,
(dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" :
"\n (" REALF ") - (" REALF ")",
b->lower, b->upper);
for( b = bounds; b < B; ++b ) {
p += sprintf(p, msg, b->lower, b->upper);
msg = "\n (" REALF ") - (" REALF ")";
}
for( comp = 0; comp < t->ncomp; ++comp ) {
cResult *r = &region->result[comp];
for( comp = 0, res = region->result;
comp < t->ncomp; ++comp, ++res ) {
p += sprintf(p, "%s \tchisq " REAL " (" COUNT " df)",
p0, r->chisq, df);
p0, res->chisq, df);
p0 += chars;
}

53
src/external/libCuba/src/suave/Suave.c vendored
View File

@ -2,59 +2,25 @@
Suave.c
Subregion-adaptive Vegas Monte-Carlo integration
by Thomas Hahn
last modified 13 Sep 10 th
last modified 17 Sep 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define SUAVE
#define ROUTINE "Suave"
#include "decl.h"
#define Print(s) puts(s); fflush(stdout)
#include "CSample.c"
/*********************************************************************/
static inline void DoSample(This *t, number n,
creal *w, creal *x, real *f, cint iter)
{
t->neval += n;
while( n-- ) {
if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata,
w++, &iter) == ABORT )
longjmp(t->abort, 1);
x += t->ndim;
f += t->ncomp;
}
}
/*********************************************************************/
#include "common.c"
Extern void EXPORT(Suave)(ccount ndim, ccount ncomp,
Integrand integrand, void *userdata,
creal epsrel, creal epsabs,
cint flags, cint seed,
cnumber mineval, cnumber maxeval,
cnumber nnew, creal flatness,
cchar *statefile,
count *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
{
@ -71,8 +37,7 @@ Extern void EXPORT(Suave)(ccount ndim, ccount ncomp,
t.maxeval = maxeval;
t.nnew = nnew;
t.flatness = flatness;
t.nregions = 0;
t.neval = 0;
t.statefile = statefile;
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
@ -87,8 +52,9 @@ Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp,
cint *pflags, cint *pseed,
cnumber *pmineval, cnumber *pmaxeval,
cnumber *pnnew, creal *pflatness,
cchar *statefile,
count *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
real *integral, real *error, real *prob, cint statefilelen)
{
This t;
t.ndim = *pndim;
@ -103,8 +69,7 @@ Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp,
t.maxeval = *pmaxeval;
t.nnew = *pnnew;
t.flatness = *pflatness;
t.nregions = 0;
t.neval = 0;
CString(t.statefile, statefile, statefilelen);
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;

44
src/external/libCuba/src/suave/common.c vendored
View File

@ -2,45 +2,34 @@
common.c
includes most of the modules
this file is part of Suave
last modified 2 Jun 10 th
last modified 29 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define RegionAlloc(p, n, nnew) MemAlloc(p, \
sizeof(Region) + \
(n)*(t->ndim + t->ncomp + 1)*sizeof(real) + \
(nnew)*t->ndim*sizeof(bin_t))
static inline Region *RegionAlloc(cThis *t, cnumber n, cnumber nnew)
{
csize_t size = sizeof(Region) +
t->ncomp*sizeof(Result) +
t->ndim*sizeof(Bounds) +
t->ncomp*t->ndim*2*sizeof(real) +
n*SAMPLESIZE +
nnew*t->ndim*sizeof(bin_t);
Region *p;
MemAlloc(p, size);
p->size = size;
return p;
}
static inline bool BadDimension(cThis *t)
{
if( t->ndim > NDIM ) return true;
if( t->ndim > MAXDIM ) return true;
return t->ndim < SOBOL_MINDIM ||
(t->seed == 0 && t->ndim > SOBOL_MAXDIM);
}
static inline bool BadComponent(cThis *t)
{
if( t->ncomp > NCOMP ) return true;
if( t->ncomp > MAXCOMP ) return true;
return t->ncomp < 1;
}
@ -49,5 +38,4 @@ static inline bool BadComponent(cThis *t)
#include "Grid.c"
#include "Sample.c"
#include "Fluct.c"
#include "Integrate.c"

52
src/external/libCuba/src/suave/decl.h vendored
View File

@ -2,29 +2,9 @@
decl.h
Type declarations
this file is part of Suave
last modified 13 Sep 10 th
last modified 29 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "stddecl.h"
@ -48,7 +28,7 @@ typedef struct {
typedef const Result cResult;
typedef struct {
real lower, upper, mid;
real lower, upper;
Grid grid;
} Bounds;
@ -62,28 +42,36 @@ typedef struct _this {
#ifndef MLVERSION
Integrand integrand;
void *userdata;
#ifdef HAVE_FORK
int ncores, *child;
real *frame;
SHM_ONLY(int shmid;)
#endif
#endif
real epsrel, epsabs;
int flags, seed;
number mineval, maxeval;
number nnew;
real flatness;
cchar *statefile;
count nregions;
number neval;
RNGState rng;
jmp_buf abort;
} This;
#define nframe nnew
typedef const This cThis;
#define TYPEDEFREGION \
typedef struct region { \
struct region *next; \
count div, df; \
number n; \
Result result[NCOMP]; \
Bounds bounds[NDIM]; \
real fluct[NCOMP][NDIM][2]; \
real w[]; \
} Region
typedef struct region {
struct region *next;
size_t size;
count div, df;
number n;
Result result[];
} Region;
#define RegionBounds(r) ((Bounds *)(r->result + t->ncomp))
#define RegionW(r) ((real *)(RegionBounds(r) + t->ndim))

43
src/external/libCuba/src/suave/util.c vendored
View File

@ -1,43 +0,0 @@
/*
util.c
Utility functions
this file is part of Suave
last modified 9 Feb 05 th
*/
/***************************************************************************
* Copyright (C) 2004-2009 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "decl.h"
static count ndim_, ncomp_, nregions_;
static number neval_;
#define RegionAlloc(p, n, nnew) \
MemAlloc(p, sizeof(Region) + \
(n)*(ndim_ + ncomp_ + 1)*sizeof(real) + \
(nnew)*ndim_*sizeof(bin_t))
#ifdef DEBUG
#include "debug.c"
#endif

26
src/external/libCuba/src/vegas/Grid.c vendored
View File

@ -2,29 +2,9 @@
Grid.c
utility functions for the Vegas grid
this file is part of Vegas
last modified 28 May 10 th
last modified 13 Dec 11 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
static inline void GetGrid(cThis *t, Grid *grid)
{
@ -35,7 +15,7 @@ static inline void GetGrid(cThis *t, Grid *grid)
if( slot < MAXGRIDS && gridptr_[slot] ) {
if( griddim_[slot] == t->ndim ) {
VecCopy(grid, gridptr_[slot]);
XCopy(grid, gridptr_[slot]);
return;
}
free(gridptr_[slot]);
@ -57,7 +37,7 @@ static inline void PutGrid(cThis *t, Grid *grid)
if( slot < MAXGRIDS ) {
if( gridptr_[slot] == NULL ) Alloc(gridptr_[slot], t->ndim);
griddim_[slot] = t->ndim;
VecCopy(gridptr_[slot], grid);
XCopy(gridptr_[slot], grid);
}
}

201
src/external/libCuba/src/vegas/Integrate.c vendored
View File

@ -2,110 +2,94 @@
Integrate.c
integrate over the unit hypercube
this file is part of Vegas
last modified 13 Sep 10 th
last modified 8 Aug 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
typedef struct {
signature_t signature;
count niter;
number nsamples, neval;
Cumulants cumul[];
} State;
static int Integrate(This *t, real *integral, real *error, real *prob)
{
real *sample;
bin_t *bins;
count dim, comp;
int fail = -99;
struct {
count niter;
number nsamples, neval;
Cumulants cumul[NCOMP];
Grid grid[NDIM];
} state;
int statemsg = VERBOSE;
struct stat st;
int fail;
StateDecl;
csize_t statesize = sizeof(State) +
NCOMP*sizeof(Cumulants) + NDIM*sizeof(Grid);
Sized(State, state, statesize);
Cumulants *c, *C = state->cumul + t->ncomp;
Grid *state_grid = (Grid *)C;
Array(Grid, margsum, NCOMP, NDIM);
Vector(char, out, 128*NCOMP + 256);
if( VERBOSE > 1 ) {
char s[512];
sprintf(s, "Vegas input parameters:\n"
sprintf(out, "Vegas input parameters:\n"
" ndim " COUNT "\n ncomp " COUNT "\n"
" epsrel " REAL "\n epsabs " REAL "\n"
" flags %d\n seed %d\n"
" mineval " NUMBER "\n maxeval " NUMBER "\n"
" nstart " NUMBER "\n nincrease " NUMBER "\n"
" nbatch " NUMBER "\n gridno %d\n"
" statefile \"%s\"\n",
" statefile \"%s\"",
t->ndim, t->ncomp,
t->epsrel, t->epsabs,
t->flags, t->seed,
t->mineval, t->maxeval,
t->nstart, t->nincrease, t->nbatch,
t->gridno, t->statefile);
Print(s);
Print(out);
}
if( BadComponent(t) ) return -2;
if( BadDimension(t) ) return -1;
SamplesAlloc(sample);
FrameAlloc(t, ShmRm(t));
ForkCores(t);
Alloc(bins, t->nbatch*t->ndim);
if( setjmp(t->abort) ) goto abort;
if( (fail = setjmp(t->abort)) ) goto abort;
IniRandom(t);
if( t->statefile && *t->statefile &&
stat(t->statefile, &st) == 0 &&
st.st_size == sizeof state && (st.st_mode & 0400) ) {
cint h = open(t->statefile, O_RDONLY);
read(h, &state, sizeof state);
close(h);
t->rng.skiprandom(t, t->neval = state.neval);
StateSetup(t);
if( VERBOSE ) {
char s[256];
sprintf(s, "\nRestoring state from %s.", t->statefile);
Print(s);
}
if( StateReadTest(t) ) {
StateReadOpen(t, fd) {
if( read(fd, state, statesize) != statesize ||
state->signature != StateSignature(t, 1) ) break;
} StateReadClose(t, fd);
t->neval = state->neval;
t->rng.skiprandom(t, t->neval);
}
else {
t->statefile = NULL;
state.niter = 0;
state.nsamples = t->nstart;
Zap(state.cumul);
GetGrid(t, state.grid);
if( ini ) {
state->niter = 0;
state->nsamples = t->nstart;
FClear(state->cumul);
GetGrid(t, state_grid);
t->neval = 0;
}
/* main iteration loop */
for( ; ; ) {
number nsamples = state.nsamples;
number nsamples = state->nsamples;
creal jacobian = 1./nsamples;
Grid margsum[NCOMP][NDIM];
Zap(margsum);
FClear(margsum);
for( ; nsamples > 0; nsamples -= t->nbatch ) {
cnumber n = IMin(t->nbatch, nsamples);
real *w = sample;
real *w = t->frame;
real *x = w + n;
real *f = x + n*t->ndim;
real *lastf = f + n*t->ncomp;
bin_t *bin = (bin_t *)lastf;
bin_t *bin = bins;
while( x < f ) {
real weight = jacobian;
@ -115,8 +99,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
for( dim = 0; dim < t->ndim; ++dim ) {
creal pos = *x*NBINS;
ccount ipos = (count)pos;
creal prev = (ipos == 0) ? 0 : state.grid[dim][ipos - 1];
creal diff = state.grid[dim][ipos] - prev;
creal prev = (ipos == 0) ? 0 : state_grid[dim][ipos - 1];
creal diff = state_grid[dim][ipos] - prev;
*x++ = prev + (pos - ipos)*diff;
*bin++ = ipos;
weight *= diff*NBINS;
@ -125,25 +109,24 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
*w++ = weight;
}
DoSample(t, n, sample, w, f, state.niter + 1);
DoSample(t, n, w, f, t->frame, state->niter + 1);
w = sample;
bin = (bin_t *)lastf;
bin = bins;
w = t->frame;
while( f < lastf ) {
creal weight = *w++;
Grid *m = &margsum[0][0];
for( comp = 0; comp < t->ncomp; ++comp ) {
for( c = state->cumul; c < C; ++c ) {
real wfun = weight*(*f++);
if( wfun ) {
Cumulants *c = &state.cumul[comp];
Grid *m = margsum[comp];
c->sum += wfun;
c->sqsum += wfun *= wfun;
for( dim = 0; dim < t->ndim; ++dim )
m[dim][bin[dim]] += wfun;
}
m += t->ndim;
}
bin += t->ndim;
@ -154,19 +137,16 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
/* compute the integral and error values */
for( comp = 0; comp < t->ncomp; ++comp ) {
Cumulants *c = &state.cumul[comp];
real avg, sigsq;
real w = Weight(c->sum, c->sqsum, state.nsamples);
sigsq = 1/(c->weightsum += w);
avg = sigsq*(c->avgsum += w*c->sum);
for( c = state->cumul; c < C; ++c ) {
real w = Weight(c->sum, c->sqsum, state->nsamples);
real sigsq = 1/(c->weightsum += w);
real avg = sigsq*(c->avgsum += w*c->sum);
c->avg = LAST ? (sigsq = 1/w, c->sum) : avg;
c->err = sqrt(sigsq);
fail |= (c->err > MaxErr(c->avg));
if( state.niter == 0 ) c->guess = c->sum;
if( state->niter == 0 ) c->guess = c->sum;
else {
c->chisum += w *= c->sum - c->guess;
c->chisqsum += w*c->sum;
@ -177,38 +157,29 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
}
if( VERBOSE ) {
char s[128 + 128*NCOMP], *p = s;
p += sprintf(p, "\n"
char *oe = out + sprintf(out, "\n"
"Iteration " COUNT ": " NUMBER " integrand evaluations so far",
state.niter + 1, t->neval);
for( comp = 0; comp < t->ncomp; ++comp ) {
cCumulants *c = &state.cumul[comp];
p += sprintf(p, "\n[" COUNT "] "
state->niter + 1, t->neval);
for( c = state->cumul, comp = 0; c < C; ++c )
oe += sprintf(oe, "\n[" COUNT "] "
REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)",
comp + 1, c->avg, c->err, c->chisq, state.niter);
}
Print(s);
++comp, c->avg, c->err, c->chisq, state->niter);
Print(out);
}
if( fail == 0 && t->neval >= t->mineval ) {
if( t->statefile ) unlink(t->statefile);
break;
}
if( fail == 0 && t->neval >= t->mineval ) break;
if( t->neval >= t->maxeval && t->statefile == NULL ) break;
if( t->neval >= t->maxeval && !StateWriteTest(t) ) break;
if( t->ncomp == 1 )
for( dim = 0; dim < t->ndim; ++dim )
RefineGrid(t, state.grid[dim], margsum[0][dim]);
RefineGrid(t, state_grid[dim], margsum[0][dim]);
else {
for( dim = 0; dim < t->ndim; ++dim ) {
Grid wmargsum;
Zap(wmargsum);
for( comp = 0; comp < t->ncomp; ++comp ) {
real w = state.cumul[comp].avg;
real w = state->cumul[comp].avg;
if( w != 0 ) {
creal *m = margsum[comp][dim];
count bin;
@ -217,41 +188,37 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
wmargsum[bin] += w*m[bin];
}
}
RefineGrid(t, state.grid[dim], wmargsum);
RefineGrid(t, state_grid[dim], wmargsum);
}
}
++state.niter;
state.nsamples += t->nincrease;
++state->niter;
state->nsamples += t->nincrease;
if( t->statefile ) {
cint h = creat(t->statefile, 0666);
if( h != -1 ) {
state.neval = t->neval;
write(h, &state, sizeof state);
close(h);
if( statemsg ) {
char s[256];
sprintf(s, "\nSaving state to %s.", t->statefile);
Print(s);
statemsg = false;
}
}
if( StateWriteTest(t) ) {
state->signature = StateSignature(t, 1);
state->neval = t->neval;
StateWriteOpen(t, fd) {
StateWrite(fd, state, statesize);
} StateWriteClose(t, fd);
if( t->neval >= t->maxeval ) break;
}
}
for( comp = 0; comp < t->ncomp; ++comp ) {
cCumulants *c = &state.cumul[comp];
cCumulants *c = &state->cumul[comp];
integral[comp] = c->avg;
error[comp] = c->err;
prob[comp] = ChiSquare(c->chisq, state.niter);
prob[comp] = ChiSquare(c->chisq, state->niter);
}
abort:
free(sample);
PutGrid(t, state.grid);
PutGrid(t, state_grid);
free(bins);
WaitCores(t);
FrameFree(t);
StateRemove(t);
return fail;
}

13
src/external/libCuba/src/vegas/Makefile.am vendored Normal file
View File

@ -0,0 +1,13 @@
## Process this file with automake to create Makefile.in
## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $
c_sources = Vegas.c
AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
noinst_LTLIBRARIES = libvegas.la
libvegas_la_SOURCES = $(c_sources)
libvegas_la_LDFLAGS = $(AM_LDFLAGS)

62
src/external/libCuba/src/vegas/Vegas.c vendored
View File

@ -2,53 +2,18 @@
Vegas.c
Vegas Monte-Carlo integration
by Thomas Hahn
last modified 13 Sep 10 th
last modified 17 Sep 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#define VEGAS
#define ROUTINE "Vegas"
#include "decl.h"
#define Print(s) puts(s); fflush(stdout)
#include "CSample.c"
/*********************************************************************/
static inline void DoSample(This *t, number n,
creal *w, creal *x, real *f, cint iter)
{
t->neval += n;
while( n-- ) {
if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata,
w++, &iter) == ABORT )
longjmp(t->abort, 1);
x += t->ndim;
f += t->ncomp;
}
}
/*********************************************************************/
#include "common.c"
Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
Integrand integrand, void *userdata,
creal epsrel, creal epsabs, cint flags, cint seed,
@ -74,7 +39,6 @@ Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
t.nbatch = nbatch;
t.gridno = gridno;
t.statefile = statefile;
t.neval = 0;
*pfail = Integrate(&t, integral, error, prob);
*pneval = t.neval;
@ -89,9 +53,8 @@ Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
cnumber *pnstart, cnumber *pnincrease,
cnumber *pnbatch, cint *pgridno, cchar *statefile,
number *pneval, int *pfail,
real *integral, real *error, real *prob, int statefilelen)
real *integral, real *error, real *prob, cint statefilelen)
{
char *s = NULL;
This t;
t.ndim = *pndim;
t.ncomp = *pncomp;
@ -107,22 +70,9 @@ Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
t.nincrease = *pnincrease;
t.nbatch = *pnbatch;
t.gridno = *pgridno;
t.neval = 0;
if( statefile ) {
/* strip trailing spaces */
while( statefilelen > 0 && statefile[statefilelen - 1] == ' ' )
--statefilelen;
if( statefilelen > 0 && (s = malloc(statefilelen + 1)) ) {
memcpy(s, statefile, statefilelen);
s[statefilelen] = 0;
}
}
t.statefile = s;
CString(t.statefile, statefile, statefilelen);
*pfail = Integrate(&t, integral, error, prob);
*pneval = t.neval;
free(s);
}

32
src/external/libCuba/src/vegas/common.c vendored
View File

@ -2,50 +2,24 @@
common.c
Code common to Vegas.c and Vegas.tm
this file is part of Vegas
last modified 6 Jun 10 th
last modified 29 Jul 13 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "Random.c"
#include "ChiSquare.c"
#include "Grid.c"
#define SamplesAlloc(p) MemAlloc(p, \
t->nbatch*((t->ndim + t->ncomp + 1)*sizeof(real) + \
t->ndim*sizeof(bin_t)))
static inline bool BadDimension(cThis *t)
{
if( t->ndim > NDIM ) return true;
if( t->ndim > MAXDIM ) return true;
return t->ndim < SOBOL_MINDIM ||
(t->seed == 0 && t->ndim > SOBOL_MAXDIM);
}
static inline bool BadComponent(cThis *t)
{
if( t->ncomp > NCOMP ) return true;
if( t->ncomp > MAXCOMP ) return true;
return t->ncomp < 1;
}
#include "Integrate.c"

29
src/external/libCuba/src/vegas/decl.h vendored
View File

@ -2,29 +2,9 @@
decl.h
Type declarations
this file is part of Vegas
last modified 13 Sep 10 th
last modified 21 Dec 11 th
*/
/***************************************************************************
* Copyright (C) 2004-2010 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "stddecl.h"
@ -56,7 +36,12 @@ typedef struct _this {
#ifndef MLVERSION
Integrand integrand;
void *userdata;
#ifdef HAVE_FORK
int ncores, *child;
SHM_ONLY(int shmid;)
#endif
#endif
real *frame;
real epsrel, epsabs;
int flags, seed;
number mineval, maxeval;
@ -68,6 +53,8 @@ typedef struct _this {
jmp_buf abort;
} This;
#define nframe nbatch
typedef const This cThis;
static Grid *gridptr_[MAXGRIDS];

46
src/external/libCuba/src/vegas/util.c vendored
View File

@ -1,46 +0,0 @@
/*
util.c
Utility functions
this file is part of Vegas
last modified 2 Mar 06 th
*/
/***************************************************************************
* Copyright (C) 2004-2009 by Thomas Hahn *
* hahn@feynarts.de *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Lesser General Public *
* License as published by the Free Software Foundation; either *
* version 2.1 of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public *
* License along with this library; if not, write to the *
* Free Software Foundation, Inc., *
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
***************************************************************************/
#include "decl.h"
static count ndim_, ncomp_;
static number neval_;
static Grid *gridptr_[MAXGRIDS];
static count griddim_[MAXGRIDS];
int EXPORT(vegasnbatch) = 1000;
int EXPORT(vegasgridno) = 0;
char EXPORT(vegasstate)[MAXSTATESIZE] = "";
#define SamplesAlloc(p, n) \
MemAlloc(p, (n)*((ndim_ + ncomp_ + 1)*sizeof(real) + ndim_*sizeof(bin_t)))
#ifdef DEBUG
#include "debug.c"
#endif