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_multidimquadrature_hpp
21#define quantlib_math_multidimquadrature_hpp
22
23#include <ql/qldefines.hpp>
24
25/* Currently, this doesn't compile under Sun C++ (see
26 https://github.com/lballabio/QuantLib/issues/223). Until that's
27 fixed, we disable it so that the rest of the library can be built.
28*/
29
30#ifndef QL_PATCH_SOLARIS
31
32#include <ql/math/integrals/gaussianquadratures.hpp>
33#include <ql/functional.hpp>
34
35namespace QuantLib {
36
37 /*! \brief Integrates a vector or scalar function of vector domain.
38
39 A template recursion along dimensions avoids calling depth
40 test or virtual functions.
41
42 \todo Add coherence test between the integrand function dimensions (the
43 vector size) and the declared dimension in the constructor.
44
45 \todo Split into integrator classes for functions returning scalar and
46 vector?
47 */
48 class GaussianQuadMultidimIntegrator {
49 private:
50 // Vector integration. Quadrature to functions returning a vector of
51 // real numbers, turns 1D quadratures into ND
52 class VectorIntegrator : public GaussHermiteIntegration {
53 public:
54 explicit VectorIntegrator(Size n, Real mu = 0.0)
55 : GaussHermiteIntegration(n, mu) {}
56
57 template <class F> // todo: fix copies.
58 std::vector<Real> operator()(const F& f) const {
59 //first one, we do not know the size of the vector returned by f
60 Integer i = order()-1;
61 std::vector<Real> term = f(x_[i]);// potential copy! @#$%^!!!
62 std::for_each(term.begin(), term.end(),
63 [&](Real x) -> Real { return x * w_[i]; });
64 std::vector<Real> sum = term;
65
66 for (i--; i >= 0; --i) {
67 term = f(x_[i]);// potential copy! @#$%^!!!
68 // sum[j] += term[j] * w_[i];
69 std::transform(term.begin(), term.end(), sum.begin(),
70 sum.begin(),
71 [&](Real x, Real y) -> Real { return w_[i]*x + y; });
72 }
73 return sum;
74 }
75 };
76
77 public:
78 /*!
79 @param dimension The number of dimensions of the argument of the
80 function we want to integrate.
81 @param quadOrder Quadrature order.
82 @param mu Parameter in the Gauss Hermite weight (i.e. points load).
83 */
84 GaussianQuadMultidimIntegrator(Size dimension, Size quadOrder,
85 Real mu = 0.);
86 //! Integration quadrature order.
87 Size order() const {return integralV_.order();}
88
89 //! Integrates function f over \f$ R^{dim} \f$
90 /* This function is just syntax since the only thing it does is calling
91 to integrate<RetType> which has to exist for the type returned by the
92 function. So theres one redundant call but there should not be any extra
93 cost... up to the compiler. It can not be templated all the way since
94 the integration entries functions can not be templates.
95 Most times integrands will return a scalar or vector but could be a
96 matrix too.
97 */
98 template<class RetType_T>
99 RetType_T operator()(const ext::function<RetType_T (
100 const std::vector<Real>& arg)>& f) const
101 {
102 return integrate<RetType_T>(f);
103 }
104
105
106 //---------------------------------------------------------
107 /* Boost fails on MSVC2008 to recognise the return type when
108 calling op() , its not boost, its me.... FIX ME*/
109
110 // Declare, spezializations follow.
111 template<class RetType_T>
112 RetType_T integrate(const ext::function<RetType_T (
113 const std::vector<Real>& v1)>& f) const;
114
115 private:
116 /* The maximum number of dimensions of the integration variable domain
117 A higher than this number of dimension would presumably be
118 impractical and another integration algorithm (MC) should be
119 considered.
120 \to do Consider moving it to a library configuration variable.
121 */
122 static const Size maxDimensions_ = 15;
123
124 //! \name Integration entry points generation
125 //@{
126 //! Recursive template methods to statically generate (at this
127 // class construction time) handles to the integration entry points
128 template<Size levelSpawn>
129 void spawnFcts() const {
130 integrationEntries_[levelSpawn-1] =
131 [&](ext::function<Real (const std::vector<Real>&)> f, Real x){
132 return scalarIntegrator<levelSpawn>(f, x);
133 };
134 integrationEntriesVR_[levelSpawn-1] =
135 [&](const ext::function<std::vector<Real>(const std::vector<Real>&)>& f, Real x){
136 return vectorIntegratorVR<levelSpawn>(f, x);
137 };
138 spawnFcts<levelSpawn-1>();
139 }
140 //@}
141
142 //---------------------------------------------------------
143
144 template <int intgDepth>
145 Real scalarIntegrator(
146 const ext::function<Real (const std::vector<Real>& arg1)>& f,
147 const Real mFctr) const
148 {
149 varBuffer_[intgDepth-1] = mFctr;
150 return integral_([&](Real x){ return scalarIntegrator<intgDepth-1>(f, x); });
151 }
152
153 template <int intgDepth>
154 std::vector<Real> vectorIntegratorVR(
155 const ext::function<std::vector<Real>(const std::vector<Real>& arg1)>& f,
156 const Real mFctr) const
157 {
158 varBuffer_[intgDepth-1] = mFctr;
159 return integralV_([&](Real x){ return vectorIntegratorVR<intgDepth-1>(f, x); });
160 }
161
162 // Same object for all dimensions poses problems when using the
163 // parallelized integrals version.
164 //! The actual integrators.
165 GaussHermiteIntegration integral_;
166 VectorIntegrator integralV_;
167
168 //! Buffer to allow acces to integrations. We do not know at which
169 // level/dimension we are going to start integration
170 // \todo Declare typedefs for traits
171 mutable std::vector<
172 ext::function<Real (ext::function<Real (
173 const std::vector<Real>& varg2)> f1,
174 const Real r3)> > integrationEntries_;
175 mutable std::vector<
176 ext::function<std::vector<Real> (const ext::function<std::vector<Real>(
177 const std::vector<Real>& vvarg2)>& vf1,
178 const Real vr3)> > integrationEntriesVR_;
179
180 Size dimension_;
181 // integration veriable buffer
182 mutable std::vector<Real> varBuffer_;
183 };
184
185
186 // Template specializations ---------------------------------------------
187
188 template<>
189 inline Real GaussianQuadMultidimIntegrator::operator()(
190 const ext::function<Real (const std::vector<Real>& v1)>& f) const
191 {
192 // integration entry level is selected now
193 return integral_([&](Real x){ return integrationEntries_[dimension_-1](ext::cref(t: f), x); });
194 }
195
196 // Scalar integrand version (merge with vector case?)
197 template<>
198 inline Real GaussianQuadMultidimIntegrator::integrate<Real>(
199 const ext::function<Real (const std::vector<Real>& v1)>& f) const
200 {
201 // integration variables
202 // call vector quadrature integration with the function and start
203 // values, kicks in recursion over the dimensions of the integration
204 // variable.
205 return integral_([&](Real x){ return integrationEntries_[dimension_-1](ext::cref(t: f), x); });
206 }
207
208 // Vector integrand version
209 template<>
210 inline std::vector<Real> GaussianQuadMultidimIntegrator::integrate<std::vector<Real>>(
211 const ext::function<std::vector<Real> (const std::vector<Real>& v1)>& f) const
212 {
213 return integralV_([&](Real x){ return integrationEntriesVR_[dimension_-1](ext::cref(t: f), x); });
214 }
215
216 //! Terminal integrand; scalar function version
217 template<>
218 inline Real GaussianQuadMultidimIntegrator::scalarIntegrator<1>(
219 const ext::function<Real (const std::vector<Real>& arg1)>& f,
220 const Real mFctr) const
221 {
222 varBuffer_[0] = mFctr;
223 return f(varBuffer_);
224 }
225
226 //! Terminal integrand; vector function version
227 template<>
228 inline std::vector<Real>
229 GaussianQuadMultidimIntegrator::vectorIntegratorVR<1>(
230 const ext::function<std::vector<Real> (const std::vector<Real>& arg1)>& f,
231 const Real mFctr) const
232 {
233 varBuffer_[0] = mFctr;
234 return f(varBuffer_);
235 }
236
237 //! Terminal level:
238 template<>
239 inline void GaussianQuadMultidimIntegrator::spawnFcts<1>() const {
240 integrationEntries_[0] = [&](const ext::function<Real(const std::vector<Real>&)>& f,
241 Real x) { return scalarIntegrator<1>(f, mFctr: x); };
242 integrationEntriesVR_[0] =
243 [&](const ext::function<std::vector<Real>(const std::vector<Real>&)>& f, Real x) {
244 return vectorIntegratorVR<1>(f, mFctr: x);
245 };
246 }
247
248}
249
250#endif
251
252#endif
253

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