1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
/*
Grid.c
utility functions for the Vegas grid
this file is part of Vegas
last modified 28 May 10 th
*/
static inline void GetGrid(cThis *t, Grid *grid)
{
count bin, dim;
unsigned const int slot = abs(t->gridno) - 1;
if( t->gridno < 0 && slot < MAXGRIDS ) griddim_[slot] = 0;
if( slot < MAXGRIDS && gridptr_[slot] ) {
if( griddim_[slot] == t->ndim ) {
VecCopy(grid, gridptr_[slot]);
return;
}
free(gridptr_[slot]);
gridptr_[slot] = NULL;
}
for( bin = 0; bin < NBINS; ++bin )
grid[0][bin] = (bin + 1)/(real)NBINS;
for( dim = 1; dim < t->ndim; ++dim )
Copy(&grid[dim], &grid[0], 1);
}
/*********************************************************************/
static inline void PutGrid(cThis *t, Grid *grid)
{
unsigned const int slot = abs(t->gridno) - 1;
if( slot < MAXGRIDS ) {
if( gridptr_[slot] == NULL ) Alloc(gridptr_[slot], t->ndim);
griddim_[slot] = t->ndim;
VecCopy(gridptr_[slot], grid);
}
}
/*********************************************************************/
static void RefineGrid(cThis *t, Grid grid, Grid margsum)
{
real avgperbin, thisbin, newcur, delta;
Grid imp, newgrid;
int bin, newbin;
/* smooth the f^2 value stored for each bin */
real prev = margsum[0];
real cur = margsum[1];
real norm = margsum[0] = .5*(prev + cur);
for( bin = 1; bin < NBINS - 1; ++bin ) {
creal s = prev + cur;
prev = cur;
cur = margsum[bin + 1];
norm += margsum[bin] = (s + cur)/3.;
}
norm += margsum[NBINS - 1] = .5*(prev + cur);
if( norm == 0 ) return;
norm = 1/norm;
/* compute the importance function for each bin */
avgperbin = 0;
for( bin = 0; bin < NBINS; ++bin ) {
real impfun = 0;
if( margsum[bin] > 0 ) {
creal r = margsum[bin]*norm;
avgperbin += impfun = pow((r - 1)/log(r), 1.5);
}
imp[bin] = impfun;
}
avgperbin /= NBINS;
/* redefine the size of each bin */
cur = newcur = 0;
thisbin = 0;
bin = -1;
for( newbin = 0; newbin < NBINS - 1; ++newbin ) {
while( thisbin < avgperbin ) {
thisbin += imp[++bin];
prev = cur;
cur = grid[bin];
}
thisbin -= avgperbin;
delta = (cur - prev)*thisbin;
newgrid[newbin] = SHARPEDGES ?
cur - delta/imp[bin] :
(newcur = Max(newcur,
cur - 2*delta/(imp[bin] + imp[IDim(bin - 1)])));
}
Copy(grid, newgrid, NBINS - 1);
grid[NBINS - 1] = 1;
}
|
