1 // -*- C++ -*-
2 /**
3 * @brief
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #include <cmath>
22 #include "openturns/OrthogonalDirection.hxx"
23 #include "openturns/Log.hxx"
24 #include "openturns/PersistentObjectFactory.hxx"
25
26 BEGIN_NAMESPACE_OPENTURNS
27
28 /**
29 * @class OrthogonalDirection
30 */
31
32 CLASSNAMEINIT(OrthogonalDirection)
33
34 static const Factory<OrthogonalDirection> Factory_OrthogonalDirection;
35
36 /* Default constructor */
OrthogonalDirection()37 OrthogonalDirection::OrthogonalDirection()
38 : SamplingStrategyImplementation(0)
39 , size_(1)
40 {
41 // Nothing to do
42 }
43
44 /* Constructor with parameters */
OrthogonalDirection(const UnsignedInteger dimension,const UnsignedInteger size)45 OrthogonalDirection::OrthogonalDirection(const UnsignedInteger dimension,
46 const UnsignedInteger size)
47 : SamplingStrategyImplementation(dimension)
48 , size_(size)
49 {
50 // Nothing to do
51 }
52
53 /* Virtual constructor */
clone() const54 OrthogonalDirection * OrthogonalDirection::clone() const
55 {
56 return new OrthogonalDirection(*this);
57 }
58
59 /* Generate the next permutation of indices in-place in the size_ first elements */
nextCombination(Indices & indices) const60 void OrthogonalDirection::nextCombination(Indices & indices) const
61 {
62 UnsignedInteger i = size_ - 1;
63 while (indices[i] == dimension_ - size_ + i) --i;
64 ++indices[i];
65 for (UnsignedInteger j = i + 1; j < size_; ++j) indices[j] = indices[i] + j - i;
66 }
67
68 /* Generate a random realization of an orientation matrix in SO(dimension) uniformly
69 distributed relatively to the Haar mesure of SO(dimension).
70 The algorithm generate an element of SO(n) with the desired properties from an
71 element of SO(n-1) by the application of a Householder reflexion associated to a
72 uniform random vector on the hypersphere Sn. The starting transformation on
73 SO(1) = {-1, 1} is taken equal to the identity Id or its opposite according to the
74 parity of dimension in order to generate elements of SO(dimension) and not of
75 O(dimension) \ SO(dimension). For an explanation of why it works, please refer to:
76 Francesco Mezzadri, "How to Generate Random Matrices from the Classical Compact Groups",
77 notices of the AMS, Volume 54, Number 5., May 2007, available online at:
78 https://www.ams.org/journals/notices/200705/fea-mezzadri-web.pdf
79 */
getUniformOrientationRealization() const80 Matrix OrthogonalDirection::getUniformOrientationRealization() const
81 {
82 Matrix Q(dimension_, dimension_);
83 // Initialization according to the parity of dimension
84 Q(0, 0) = (dimension_ % 2 == 0 ? -1.0 : 1.0);
85 Matrix column(dimension_, 1);
86 for (UnsignedInteger indexDimension = 1; indexDimension < dimension_; indexDimension++)
87 {
88 Q(indexDimension, indexDimension) = 1.0;
89 Point v(getUniformUnitVectorRealization(indexDimension + 1));
90 for (UnsignedInteger index = 0; index <= indexDimension; ++index) column(index, 0) = v[index];
91 Q = Q - ((2.0 * column) * (column.transpose() * Q));
92 }
93 return Q;
94 }
95
96 /* Add the 2^size linear combinations of columns of Q indicated in the
97 * size first elements of indices by affecting all the choices of sign
98 * to the coefficients of the linear combination */
computePartialSample(const Indices & indices,const Matrix & Q,Sample & result) const99 void OrthogonalDirection::computePartialSample(const Indices & indices,
100 const Matrix & Q,
101 Sample & result) const
102 {
103 // Normalization factor of the linear combination
104 const Scalar factor = 1.0 / sqrt(1.0 * size_);
105 // We have 2^size linear combinations to generate
106 const UnsignedInteger indexLinearCombinationMax = 1 << size_;
107 // For each combination
108 for (UnsignedInteger indexLinearCombination = 0; indexLinearCombination < indexLinearCombinationMax; ++indexLinearCombination)
109 {
110 Point direction(dimension_);
111 // The combination index is used as a mask to select the coefficients equal to 1.0 or to -1.0
112 UnsignedInteger mask = indexLinearCombination;
113 for (UnsignedInteger index = 0; index < size_; ++index)
114 {
115 // Which column of Q corresponds to the index position of indices?
116 const UnsignedInteger column = indices[index];
117 // Sign affected to this column
118 const Scalar sign = 1.0 - 2.0 * (mask % 2);
119 // Summation
120 for (UnsignedInteger row = 0; row < dimension_; ++row) direction[row] += sign * Q(row, column);
121 // Next bit of the mask
122 mask /= 2;
123 }
124 // Normalize the linear combination
125 result.add(factor * direction);
126 }
127 }
128
129 /* Generate a set of directions */
generate() const130 Sample OrthogonalDirection::generate() const
131 {
132 Sample result(0, dimension_);
133 Matrix Q(getUniformOrientationRealization());
134 Indices indices(size_);
135 // Start with the first lexicographic combination
136 indices.fill();
137 computePartialSample(indices, Q, result);
138 while(indices[0] != dimension_ - size_)
139 {
140 nextCombination(indices);
141 computePartialSample(indices, Q, result);
142 }
143 LOGDEBUG(OSS() << "OrthogonalDirection::generate: directions=\n" << result);
144 return result;
145 }
146
147 /* String converter */
__repr__() const148 String OrthogonalDirection::__repr__() const
149 {
150 OSS oss;
151 oss << "class=" << OrthogonalDirection::GetClassName()
152 << " derived from " << SamplingStrategyImplementation::__repr__()
153 << " size=" << size_;
154 return oss;
155 }
156
157 END_NAMESPACE_OPENTURNS
158
159