QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
basisincompleteordered.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2007, 2008 Mark Joshi
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy 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
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
21#include <algorithm>
22
23namespace QuantLib {
24
26 : euclideanDimension_(euclideanDimension) {}
27
28 bool BasisIncompleteOrdered::addVector(const Array& newVector1) {
29
30 QL_REQUIRE(newVector1.size() == euclideanDimension_,
31 "missized vector passed to "
32 "BasisIncompleteOrdered::addVector");
33
34 newVector_ = newVector1;
35
37 return false;
38
39 for (auto& currentBasi : currentBasis_) {
40 Real innerProd =
41 std::inner_product(newVector_.begin(), newVector_.end(), currentBasi.begin(), Real(0.0));
42
43 for (Size k=0; k<euclideanDimension_; ++k)
44 newVector_[k] -= innerProd * currentBasi[k];
45 }
46
47 Real norm = std::sqrt(std::inner_product(newVector_.begin(),
49 newVector_.begin(), Real(0.0)));
50
51 if (norm<1e-12) // maybe this should be a tolerance
52 return false;
53
54 for (Size l=0; l<euclideanDimension_; ++l)
55 newVector_[l]/=norm;
56
57 currentBasis_.push_back(newVector_);
58
59 return true;
60 }
61
63 return currentBasis_.size();
64 }
65
68 }
69
70
73 for (Size i=0; i<basis.rows(); ++i)
74 for (Size j=0; j<basis.columns(); ++j)
75 basis[i][j] = currentBasis_[i][j];
76
77 return basis;
78 }
79
80 namespace
81 {
82 Real normSquared(const Matrix& v, Size row)
83 {
84 Real x=0.0;
85 for (Size i=0; i < v.columns(); ++i)
86 x += v[row][i]*v[row][i];
87
88 return x;
89 }
90
91
92 Real norm(const Matrix& v, Size row)
93 {
94 return std::sqrt(normSquared( v, row));
95 }
96
97 Real innerProduct(const Matrix& v, Size row1, const Matrix& w, Size row2)
98 {
99
100 Real x=0.0;
101 for (Size i=0; i < v.columns(); ++i)
102 x += v[row1][i]*w[row2][i];
103
104 return x;
105 }
106
107 }
108
109
110
112 Real multiplierCutoff,
113 Real tolerance)
114 : originalVectors_(originalVectors),
115 multiplierCutoff_(multiplierCutoff),
116 numberVectors_(originalVectors.rows()),
117 dimension_(originalVectors.columns()),
118 validVectors_(true,originalVectors.rows()), // opposite way round from vector constructor
119 orthoNormalizedVectors_(originalVectors.rows(),
120 originalVectors.columns())
121 {
122 std::vector<Real> currentVector(dimension_);
123 for (Size j=0; j < numberVectors_; ++j)
124 {
125
126 if (validVectors_[j])
127 {
128 for (Size k=0; k< numberVectors_; ++k) // create an orthormal basis not containing j
129 {
130 for (Size m=0; m < dimension_; ++m)
132
133 if ( k !=j && validVectors_[k])
134 {
135
136
137 for (Size l=0; l < k; ++l)
138 {
139 if (validVectors_[l] && l !=j)
140 {
141 Real dotProduct = innerProduct(orthoNormalizedVectors_, k, orthoNormalizedVectors_,l);
142 for (Size n=0; n < dimension_; ++n)
144 }
145
146 }
147
148 Real normBeforeScaling= norm(orthoNormalizedVectors_,k);
149
150 if (normBeforeScaling < tolerance)
151 {
152 validVectors_[k] = false;
153 }
154 else
155 {
156 Real normBeforeScalingRecip = 1.0/normBeforeScaling;
157 for (Size m=0; m < dimension_; ++m)
158 orthoNormalizedVectors_[k][m] *= normBeforeScalingRecip;
159
160 } // end of else (norm < tolerance)
161
162 } // end of if k !=j && validVectors_[k])
163
164 }// end of for (Size k=0; k< numberVectors_; ++k)
165
166 // we now have an o.n. basis for everything except j
167
168 Real prevNormSquared = normSquared(originalVectors_, j);
169
170
171 for (Size r=0; r < numberVectors_; ++r)
172 if (validVectors_[r] && r != j)
173 {
174 Real dotProduct = innerProduct(orthoNormalizedVectors_, j, orthoNormalizedVectors_,r);
175
176 for (Size s=0; s < dimension_; ++s)
178
179 }
180
181 Real projectionOnOriginalDirection = innerProduct(originalVectors_,j,orthoNormalizedVectors_,j);
182 Real sizeMultiplier = prevNormSquared/projectionOnOriginalDirection;
183
184 if (std::fabs(sizeMultiplier) < multiplierCutoff_)
185 {
186 for (Size t=0; t < dimension_; ++t)
187 currentVector[t] = orthoNormalizedVectors_[j][t]*sizeMultiplier;
188
189 }
190 else
191 validVectors_[j] = false;
192
193
194 } // end of if (validVectors_[j])
195
196 projectedVectors_.push_back(currentVector);
197
198
199 } //end of j loop
200
202 for (Size i=0; i < numberVectors_; ++i)
204
205
206 } // end of constructor
207
208 const std::valarray<bool>& OrthogonalProjections::validVectors() const
209 {
210 return validVectors_;
211
212 }
213
214 const std::vector<Real>& OrthogonalProjections::GetVector(Size index) const
215 {
216 return projectedVectors_[index];
217 }
218
219
221 {
222 return numberValidVectors_;
223 }
224
225
226
227
228}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
BasisIncompleteOrdered(Size euclideanDimension)
bool addVector(const Array &newVector)
return value indicates if the vector was linearly independent
Matrix used in linear algebra.
Definition: matrix.hpp:41
Size rows() const
Definition: matrix.hpp:504
Size columns() const
Definition: matrix.hpp:508
OrthogonalProjections(const Matrix &originalVectors, Real multiplierCutOff, Real tolerance)
const std::valarray< bool > & validVectors() const
const std::vector< Real > & GetVector(Size index) const
std::valarray< bool > validVectors_
outputs
std::vector< std::vector< Real > > projectedVectors_
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v