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
20#include <ql/math/matrixutilities/basisincompleteordered.hpp>
21#include <algorithm>
22
23namespace QuantLib {
24
25 BasisIncompleteOrdered::BasisIncompleteOrdered(Size euclideanDimension)
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
36 if (currentBasis_.size()==euclideanDimension_)
37 return false;
38
39 for (auto& currentBasi : currentBasis_) {
40 Real innerProd =
41 std::inner_product(first1: newVector_.begin(), last1: newVector_.end(), first2: currentBasi.begin(), init: 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(x: std::inner_product(first1: newVector_.begin(),
48 last1: newVector_.end(),
49 first2: newVector_.begin(), init: 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(x: newVector_);
58
59 return true;
60 }
61
62 Size BasisIncompleteOrdered::basisSize() const {
63 return currentBasis_.size();
64 }
65
66 Size BasisIncompleteOrdered::euclideanDimension() const {
67 return euclideanDimension_;
68 }
69
70
71 Matrix BasisIncompleteOrdered::getBasisAsRowsInMatrix() const {
72 Matrix basis(currentBasis_.size(), euclideanDimension_);
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(x: 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
111 OrthogonalProjections::OrthogonalProjections(const Matrix& originalVectors,
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)
131 orthoNormalizedVectors_[k][m] = originalVectors_[k][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(v: orthoNormalizedVectors_, row1: k, w: orthoNormalizedVectors_,row2: l);
142 for (Size n=0; n < dimension_; ++n)
143 orthoNormalizedVectors_[k][n] -= dotProduct*orthoNormalizedVectors_[l][n];
144 }
145
146 }
147
148 Real normBeforeScaling= norm(v: orthoNormalizedVectors_,row: 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(v: originalVectors_, row: j);
169
170
171 for (Size r=0; r < numberVectors_; ++r)
172 if (validVectors_[r] && r != j)
173 {
174 Real dotProduct = innerProduct(v: orthoNormalizedVectors_, row1: j, w: orthoNormalizedVectors_,row2: r);
175
176 for (Size s=0; s < dimension_; ++s)
177 orthoNormalizedVectors_[j][s] -= dotProduct*orthoNormalizedVectors_[r][s];
178
179 }
180
181 Real projectionOnOriginalDirection = innerProduct(v: originalVectors_,row1: j,w: orthoNormalizedVectors_,row2: j);
182 Real sizeMultiplier = prevNormSquared/projectionOnOriginalDirection;
183
184 if (std::fabs(x: 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(x: currentVector);
197
198
199 } //end of j loop
200
201 numberValidVectors_ =0;
202 for (Size i=0; i < numberVectors_; ++i)
203 numberValidVectors_ += validVectors_[i] ? 1 : 0;
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
220 Size OrthogonalProjections::numberValidVectors() const
221 {
222 return numberValidVectors_;
223 }
224
225
226
227
228}
229

source code of quantlib/ql/math/matrixutilities/basisincompleteordered.cpp