-
Notifications
You must be signed in to change notification settings - Fork 3
/
multi_channel.dox
71 lines (62 loc) · 4.01 KB
/
multi_channel.dox
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
/**
\defgroup multi_channel_group Multi Channel Integrator
\brief The Multi Channel integration algorithm
In contrast to the VEGAS algorithm, which uses a product of one-dimensional PDFs, the multi channel
integration algorithm uses a sum of \f$ M \f$ user-defined PDFs \f$ p_j ( \vec{y}) \f$ with
automatically adapted weights \f$ \alpha_j \f$ in the form
\f[
p ( \vec{y} ) = \sum_{j=1}^M \alpha_j p_j ( \vec{y} ) \text{,} \quad
\int \mathrm{d}^d y \; p_j \left( \vec{y} \right) = 1 \text{,} \quad
\sum_{j=1}^M \alpha_j = 1 \text{,}
\f]
which is used to sample the function. Multi channel integration performs better compared to VEGAS if
the peak structure of the integrand does not factorize. To avoid explicitly summing over every
*channel* \f$ j \f$ for each randomly chosen point \f$ \vec{x} \in U \equiv [0,1]^d \f$ the MC
integrator also samples over the possible channels. This leads to the final formula
\f[
I = \int \mathrm{d} i \int_U \mathrm{d}^d x \; \left.
\frac{f (\vec{y})}{p (\vec{y})} \right|_{\vec{y} = \vec{y}_i ( \vec{x}
)} \text{,} \qquad
i \sim \left\{ \alpha_1, \alpha_2, \ldots, \alpha_M \right\}
\f]
where the integral over the index \f$ i \f$ is understood to be a Monte Carlo summation in which the
index \f$ i \f$ is randomly chosen according to the specified weights. Since the PDFs \f$ \left\{
p_j \right\}_{j=1}^M \f$ are user-defined, the user also has to specify the CDFs \f$ \left\{
\vec{y}_j ( \vec{x} ) \right\}_{j=1}^M \f$ for each channel.
The multi channel integration uses the following parameters, where `T` denotes the numerical type,
e.g. `double`:
- `dimensions` determines the parameter \f$ d \f$,
- `map_dimensions` determines the size of the vector `coordinates` (see below),
- `channels` determines \f$ M \f$, the number of channels,
- `map` must be the function that calculates *both* the PDFs and the CDFs. Its declaration must
be as follows:
T map(
std::size_t channel,
std::vector<T> const& random_numbers,
std::vector<T>& coordinates,
std::vector<std::size_t> const& enabled_channels,
std::vector<T>& densities,
hep::multi_channel_map action
);
The parameter `action` determines what `map` should do:
- If `action` is \ref hep::multi_channel_map::calculate_coordinates, then the function
must use `random_numbers` to generate a point according to the PDF with the index given
by `channel`. The point itself must be written into `coordinates`, whose size is already
set by the parameter `map_dimensions` (see above). The return value of `map` is ignored.
- If `action` is \ref hep::multi_channel_map::calculate_densities, then the `map` function
must generate the PDFs for the previously generated point, which can be read out of
`coordinates`. The PDFs must be written into `densities`, whose size is determined by
`channels` (see above). The inverse of the return value of the `map` function is
multiplied with every PDF. The vector `enabled_channels` stores the indices of all
channels whose weights are not zero. If a weight is zero, then the PDF for this channel
will be ignored; this can be used to speed up the calculation if there are many disabled
channels.
The reason for the having the parameter `action` is that for some applications the integrand
often evaluates to zero. If this happens, the `densities` are not needed, and then the multi
channels integrator skips the possibly costly call to the `map` function with
`calculate_densities`.
- `function` must be the integrand function that is integrated over. Its declaration must be as
described in \ref integrands. The Monte Carlo point can be captured e.g. using \ref
multi_channel_point or, if `function` needs access to data that has been computed already in
`densities`, it can be captured using \ref multi_channel_point2.
*/