1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2014 Jose Aparicio
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#ifndef quantlib_math_multidimintegrator_hpp
21#define quantlib_math_multidimintegrator_hpp
22
23#include <ql/types.hpp>
24#include <ql/errors.hpp>
25#include <ql/math/integrals/integral.hpp>
26#include <ql/functional.hpp>
27#include <vector>
28
29namespace QuantLib {
30
31 /*! \brief Integrates a vector or scalar function of vector domain.
32
33 Uses a collection of arbitrary 1D integrators along each of the
34 dimensions. A template recursion along dimensions avoids calling depth
35 test or virtual functions.\par
36 This class generalizes to an arbitrary number of dimensions the
37 functionality in class TwoDimensionalIntegral
38 */
39 class MultidimIntegral {
40 public:
41 explicit MultidimIntegral(
42 const std::vector<ext::shared_ptr<Integrator> >& integrators);
43
44 // scalar variant
45 /*! f is the integrand function; a and b are the lower and
46 upper integration limit domain for each dimension.
47 */
48 Real operator()(
49 const ext::function<Real (const std::vector<Real>&)>& f,
50 const std::vector<Real>& a,
51 const std::vector<Real>& b) const
52 {
53 QL_REQUIRE((a.size()==b.size())&&(b.size()==integrators_.size()),
54 "Incompatible integration problem dimensions");
55 return integrationLevelEntries_[integrators_.size()-1](f, a, b);
56 }
57 // to do: write std::vector<Real> operator()(...) version
58
59 private:
60 static const Size maxDimensions_ = 15;
61
62 /* Here is the tradeoff; this is avoiding the dimension limits checks
63 during integration at the price of these asignments during construction.
64 Explicit template instantiation is of no use, an object is needed
65 (notice 'this' is needed for the asignment.)
66 If not all the dimensions up the maximum number are used the waste goes
67 into storage of the functions (in fact only one is used)
68 */
69 template<Size depth>
70 void spawnFcts() const;
71 // Splits the integration in cross-sections per dimension.
72 template<int T_N>
73 Real vectorBinder (
74 const ext::function<Real (const std::vector<Real>&)>& f,
75 Real z,
76 const std::vector<Real>& a,
77 const std::vector<Real>& b) const ;
78 // actual integration of dimension nT
79 template<int nT>
80 Real integrate(
81 const ext::function<Real (const std::vector<Real>&)>& f,
82 const std::vector<Real>& a,
83 const std::vector<Real>& b) const;
84
85 const std::vector<ext::shared_ptr<Integrator> > integrators_;
86
87 /* typedef (const ext::function<Real
88 (const std::vector<Real>&arg1)>&arg2) integrableFunctType;
89 */
90
91 /* vector of, functions returning reals And taking as argument:
92 1.- a const ref to a function taking vectors
93 2.- a vector, 3. another vector. typedefs eventually...
94 at first sight this might look like mimicking a virtual table, it isnt
95 that. The reason is to be able to select the correct integration
96 dimension at run time, this can not be done before because of the
97 template argument restriction to be constant known at compilation.
98 */
99 mutable std::vector<ext::function<Real (//<- members: integrate<N>
100 // integrable function:
101 const ext::function<Real (const std::vector<Real>&)>&,
102 const std::vector<Real>&, //<- a
103 const std::vector<Real>&) //<- b
104 > >
105 integrationLevelEntries_;
106
107 /* One can avoid the passing around of the ct refs to a and b but the
108 price is to keep a copy of them (they are unknown at construction time)
109 On the other hand the vector integration variable has to be created.*/
110 mutable std::vector<Real> varBuffer_;
111
112 };
113
114 // spez last call/dimension
115 template<>
116 Real inline MultidimIntegral::vectorBinder<0> (
117 const ext::function<Real (const std::vector<Real>&)>& f,
118 Real z,
119 const std::vector<Real>& a,
120 const std::vector<Real>& b) const
121 {
122 varBuffer_[0] = z;
123 return f(varBuffer_);
124 }
125
126 template<>
127 void inline MultidimIntegral::spawnFcts<1>() const {
128 integrationLevelEntries_[0] = [this](const auto& f, const auto& a, const auto& b) {
129 return this->integrate<0>(f, a, b);
130 };
131 }
132
133 template<int nT>
134 inline Real MultidimIntegral::integrate(
135 const ext::function<Real (const std::vector<Real>&)>& f,
136 const std::vector<Real>& a,
137 const std::vector<Real>& b) const
138 {
139 return
140 (*integrators_[nT])([this, &f, &a, &b](auto z) {
141 return this->vectorBinder<nT>(f, z, a, b);
142 }, a[nT], b[nT]);
143 }
144
145 template<int T_N>
146 inline Real MultidimIntegral::vectorBinder (
147 const ext::function<Real (const std::vector<Real>&)>& f,
148 Real z,
149 const std::vector<Real>& a,
150 const std::vector<Real>& b) const
151 {
152 varBuffer_[T_N] = z;
153 return integrate<T_N-1>(f, a, b);
154 }
155
156 template<Size depth>
157 void MultidimIntegral::spawnFcts() const {
158 integrationLevelEntries_[depth-1] = [this](const auto& f, const auto& a, const auto& b) {
159 return this->integrate<depth-1>(f, a, b);
160 };
161 spawnFcts<depth-1>();
162 }
163
164}
165
166#endif
167

source code of quantlib/ql/experimental/math/multidimintegrator.hpp