1 // -*- C++ -*-
2 /**
3 * @brief Implementation of the Sobol' sequence
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 <cstdlib>
22
23 #include "openturns/SobolSequence.hxx"
24 #include "openturns/ResourceMap.hxx"
25 #include "openturns/Exception.hxx"
26 #include "openturns/PersistentObjectFactory.hxx"
27
28 BEGIN_NAMESPACE_OPENTURNS
29
30 CLASSNAMEINIT(SobolSequence)
31 static const Factory<SobolSequence> Factory_SobolSequence;
32
33
34 const UnsignedInteger SobolSequence::MaximumBase2Logarithm = 62;
35 const Scalar SobolSequence::Epsilon = 1.0 / power2(MaximumBase2Logarithm);
36
37 // Not in the openturns subdir as it has not to be visible
38 #include "SobolSequenceDirections.hxx"
39
40 /* Constructor with parameters */
SobolSequence(const UnsignedInteger dimension)41 SobolSequence::SobolSequence(const UnsignedInteger dimension)
42 : LowDiscrepancySequenceImplementation(dimension)
43 {
44 initialize(dimension);
45 }
46
47
48 /* Virtual constructor */
clone() const49 SobolSequence * SobolSequence::clone() const
50 {
51 return new SobolSequence(*this);
52 }
53
54
55 /* Initialize the sequence */
initialize(const UnsignedInteger dimension)56 void SobolSequence::initialize(const UnsignedInteger dimension)
57 {
58 LowDiscrepancySequenceImplementation::initialize(dimension);
59 if(!(dimension <= MaximumDimension))
60 throw InvalidDimensionException(HERE) << "Dimension must be in range [0-" << MaximumDimension << "], here dimension=" << dimension << ".";
61 // copy initial direction numbers
62 base_ = Unsigned64BitsIntegerCollection(dimension_ * MaximumBase2Logarithm, 0);
63
64 for ( UnsignedInteger i = 0; i < dimension_; ++ i )
65 for ( UnsignedInteger j = 0; j < MaximumInitialDegree; ++ j )
66 base_[i * MaximumBase2Logarithm + j] = InitialBase[i * MaximumInitialDegree + j];
67
68 // initialize row 0 (first dimension)
69 for ( UnsignedInteger j = 0; j < MaximumBase2Logarithm; ++ j )
70 base_[j] = 1;
71 // initialize remaining direction numbers, for each dimension <-> rows of base_[][]
72 for ( UnsignedInteger i = 1; i < dimension_; ++ i )
73 {
74 // number of bits of PrimitivePolynomial[i]
75 UnsignedInteger polynomialCoefficientDegree = 0;
76 Unsigned64BitsInteger polynomialCoefficient = PrimitivePolynomial[i];
77
78 // get number of bits of the PrimitivePolynomial[i] coefficient (could have used log2)
79 while(polynomialCoefficient > 1)
80 {
81 polynomialCoefficient /= 2;
82 ++ polynomialCoefficientDegree;
83 }
84
85 // generate remaining direction numbers
86 for ( UnsignedInteger j = polynomialCoefficientDegree; j < MaximumBase2Logarithm; ++ j )
87 {
88 base_[i * MaximumBase2Logarithm + j] = base_[i * MaximumBase2Logarithm + j - polynomialCoefficientDegree];
89 for ( UnsignedInteger k = 1; k <= polynomialCoefficientDegree; ++ k )
90 {
91 if((PrimitivePolynomial[i] & power2(polynomialCoefficientDegree - k)) > 0)
92 {
93 base_[i * MaximumBase2Logarithm + j] ^= base_[i * MaximumBase2Logarithm + j - k] * power2(k);
94 }
95 } // k
96 } // j
97 } // i
98
99 // multiply columns of directionNumber[][] by appropriate power of 2
100 for ( UnsignedInteger j = 0; j < MaximumBase2Logarithm - 1; ++ j )
101 for ( UnsignedInteger i = 0; i < dimension_; ++ i )
102 base_[i * MaximumBase2Logarithm + j] *= power2(MaximumBase2Logarithm - j - 1);
103
104 // initialize integer coefficients of the sequence : first column of directionNumber[][]
105 coefficients_ = Unsigned64BitsIntegerCollection(dimension_);
106 for ( UnsignedInteger i = 0; i < dimension_; ++ i )
107 coefficients_[i] = base_[i * MaximumBase2Logarithm];
108
109 // set the seed
110 seed_ = ResourceMap::GetAsUnsignedInteger( "SobolSequence-InitialSeed" );
111 }
112
113
114 /* Generate a pseudo-random vector of independant numbers uniformly distributed over [0, 1[ */
generate() const115 Point SobolSequence::generate() const
116 {
117 if (seed_ == 0)
118 {
119 ++seed_;
120 return Point(dimension_);
121 }
122 // initialize a point with values 2^-MaximumBase2Logarithm
123 Point sequencePoint(dimension_, Epsilon);
124
125 // compute the position of the lowest 0 bit in the binary representation of seed_
126 const UnsignedInteger positionOfLowest0BitOfSeed = computePositionOfLowest0Bit(seed_);
127
128 // for each dimension
129 for(UnsignedInteger i = 0; i < dimension_; ++ i )
130 {
131 // compute sequence from integer coefficients
132 sequencePoint[i] *= coefficients_[i];
133
134 // update integer coefficients for next generation
135 coefficients_[i] ^= base_[i * MaximumBase2Logarithm + positionOfLowest0BitOfSeed - 1];
136 }
137 // increment seed for next generation
138 ++seed_;
139 return sequencePoint;
140 }
141
142
143 /* String converter */
__repr__() const144 String SobolSequence::__repr__() const
145 {
146 OSS oss;
147 oss << "class=" << SobolSequence::GetClassName()
148 << " coefficients=" << coefficients_
149 << " seed=" << seed_;
150 return oss;
151 }
152
153
154 /* return 2^n */
power2(const UnsignedInteger n)155 Unsigned64BitsInteger SobolSequence::power2(const UnsignedInteger n)
156 {
157 return (Unsigned64BitsInteger)1 << n;
158 }
159
160
161 /* returns the position of the lowest '0' in the binary representation of an integer */
computePositionOfLowest0Bit(const Unsigned64BitsInteger number)162 UnsignedInteger SobolSequence::computePositionOfLowest0Bit(const Unsigned64BitsInteger number)
163 {
164 if (number == 0) return 0;
165 UnsignedInteger base2Logarithm = 0;
166 while((number & power2(base2Logarithm)) && (base2Logarithm <= MaximumBase2Logarithm))
167 {
168 ++ base2Logarithm;
169 }
170 return base2Logarithm + 1;// 1 <=> first bit
171 }
172
173
save(Advocate & adv) const174 void SobolSequence::save(Advocate & adv) const
175 {
176 LowDiscrepancySequenceImplementation::save(adv);
177 adv.saveAttribute( "base_", base_);
178 adv.saveAttribute( "seed_", seed_);
179 adv.saveAttribute( "coefficients_", coefficients_);
180 }
181
182
183 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)184 void SobolSequence::load(Advocate & adv)
185 {
186 LowDiscrepancySequenceImplementation::load(adv);
187 adv.loadAttribute( "base_", base_);
188 adv.loadAttribute( "seed_", seed_);
189 adv.loadAttribute( "coefficients_", coefficients_);
190 }
191
192
193 END_NAMESPACE_OPENTURNS
194