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 #include <ql/qldefines.hpp>
21 #ifdef BOOST_MSVC
22 # include <ql/auto_link.hpp>
23 #endif
24 #include <ql/experimental/math/multidimintegrator.hpp>
25 #include <ql/experimental/math/multidimquadrature.hpp>
26 #include <ql/math/integrals/trapezoidintegral.hpp>
27 #include <ql/patterns/singleton.hpp>
28 #include <ql/functional.hpp>
29
30
31 #include <iostream>
32 #include <iomanip>
33
34 using namespace QuantLib;
35 using namespace std;
36
37 #if defined(QL_ENABLE_SESSIONS)
38 namespace QuantLib {
39
sessionId()40 ThreadKey sessionId() { return 0; }
41
42 }
43 #endif
44
45 // Correct value is: (e^{-.25} \sqrt{\pi})^{dimension}
46 struct integrand {
operator ()integrand47 Real operator()(const std::vector<Real>& arg) const {
48 Real sum = 1.;
49 for(Size i=0; i<arg.size(); i++)
50 sum *= std::exp(-arg[i]*arg[i]) * std::cos(arg[i]);
51 return sum;
52 }
53 };
54
main()55 int main() {
56
57 try {
58
59 std::cout << std::endl;
60
61 /*
62 Integrates the function above over several dimensions, the size of the
63 vector argument is the dimension one.
64 Both algorithms are not really on the same stand since the quadrature
65 will be incorrect to use if the integrand is not appropiately behaved. Over
66 dimension 3 you might need to modify the points in the integral to retain a
67 sensible computing time.
68 */
69 Size dimension = 3;
70 Real exactSol = std::pow(std::exp(-.25) *
71 std::sqrt(M_PI), static_cast<Real>(dimension));
72
73 ext::function<Real(const std::vector<Real>& arg)> f = integrand();
74
75 #ifndef QL_PATCH_SOLARIS
76 GaussianQuadMultidimIntegrator intg(dimension, 15);
77
78 Real valueQuad = intg(f);
79 #endif
80
81 std::vector<ext::shared_ptr<Integrator> > integrals;
82 for(Size i=0; i<dimension; i++)
83 integrals.push_back(
84 ext::make_shared<TrapezoidIntegral<Default> >(1.e-4, 20));
85 std::vector<Real> a_limits(integrals.size(), -4.);
86 std::vector<Real> b_limits(integrals.size(), 4.);
87 MultidimIntegral testIntg(integrals);
88
89 Real valueGrid = testIntg(f, a_limits, b_limits);
90
91 cout << fixed << setprecision(4);
92 cout << endl << "-------------- " << endl
93 << "Exact: " << exactSol << endl
94 #ifndef QL_PATCH_SOLARIS
95 << "Quad: " << valueQuad << endl
96 #endif
97 << "Grid: " << valueGrid << endl
98 << endl;
99
100 return 0;
101
102 } catch (std::exception& e) {
103 std::cerr << e.what() << std::endl;
104 return 1;
105 } catch (...) {
106 std::cerr << "unknown error" << std::endl;
107 return 1;
108 }
109 }
110