QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
multidimquadrature.hpp
Go to the documentation of this file.
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
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 */
49 private:
50 // Vector integration. Quadrature to functions returning a vector of
51 // real numbers, turns 1D quadratures into ND
53 public:
54 explicit VectorIntegrator(Size n, Real mu = 0.0)
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>
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.
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,
175 mutable std::vector<
176 ext::function<std::vector<Real> (const ext::function<std::vector<Real>(
177 const std::vector<Real>& vvarg2)>& vf1,
179
181 // integration veriable buffer
182 mutable std::vector<Real> varBuffer_;
183 };
184
185
186 // Template specializations ---------------------------------------------
187
188 template<>
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(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(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(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, x); };
243 [&](const ext::function<std::vector<Real>(const std::vector<Real>&)>& f, Real x) {
244 return vectorIntegratorVR<1>(f, x);
245 };
246 }
247
248}
249
250#endif
251
252#endif
generalized Gauss-Hermite integration
Integrates a vector or scalar function of vector domain.
std::vector< ext::function< Real(ext::function< Real(const std::vector< Real > &varg2)> f1, const Real r3)> > integrationEntries_
Buffer to allow acces to integrations. We do not know at which.
RetType_T operator()(const ext::function< RetType_T(const std::vector< Real > &arg)> &f) const
Integrates function f over .
Real scalarIntegrator(const ext::function< Real(const std::vector< Real > &arg1)> &f, const Real mFctr) const
RetType_T integrate(const ext::function< RetType_T(const std::vector< Real > &v1)> &f) const
void spawnFcts() const
Recursive template methods to statically generate (at this.
Size order() const
Integration quadrature order.
std::vector< ext::function< std::vector< Real >(const ext::function< std::vector< Real >(const std::vector< Real > &vvarg2)> &vf1, const Real vr3)> integrationEntriesVR_
GaussHermiteIntegration integral_
The actual integrators.
std::vector< Real > vectorIntegratorVR(const ext::function< std::vector< Real >(const std::vector< Real > &arg1)> &f, const Real mFctr) const
Maps function, bind and cref to either the boost or std implementation.
Integral of a 1-dimensional function using the Gauss quadratures.
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
Global definitions and compiler switches.
Real F
Definition: sabr.cpp:200