1 // -*- C++ -*-
2 /**
3 * @brief Abstract top-level view of an monteCarloExperiment plane
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 "openturns/LowDiscrepancyExperiment.hxx"
22 #include "openturns/SobolSequence.hxx"
23 #include "openturns/Exception.hxx"
24 #include "openturns/IndependentCopula.hxx"
25 #include "openturns/DistributionTransformation.hxx"
26 #include "openturns/PersistentObjectFactory.hxx"
27 #include "openturns/RandomGenerator.hxx"
28
29 BEGIN_NAMESPACE_OPENTURNS
30
31 CLASSNAMEINIT(LowDiscrepancyExperiment)
32
33 static const Factory<LowDiscrepancyExperiment> Factory_LowDiscrepancyExperiment;
34
35 /* Default constructor */
LowDiscrepancyExperiment()36 LowDiscrepancyExperiment::LowDiscrepancyExperiment()
37 : WeightedExperimentImplementation()
38 , sequence_(SobolSequence())
39 , restart_(true)
40 , randomize_(false)
41 {
42 // Build the iso-probabilistic transformation
43 setDistribution(distribution_);
44 }
45
46 /* Constructor with parameters */
LowDiscrepancyExperiment(const UnsignedInteger size,const Bool restart)47 LowDiscrepancyExperiment::LowDiscrepancyExperiment(const UnsignedInteger size,
48 const Bool restart)
49 : WeightedExperimentImplementation(size)
50 , sequence_(SobolSequence())
51 , restart_(restart)
52 , randomize_(false)
53 {
54 // Build the iso-probabilistic transformation
55 setDistribution(distribution_);
56 }
57
58 /* Constructor with parameters */
LowDiscrepancyExperiment(const LowDiscrepancySequence & sequence,const Distribution & distribution,const UnsignedInteger size,const Bool restart)59 LowDiscrepancyExperiment::LowDiscrepancyExperiment(const LowDiscrepancySequence & sequence,
60 const Distribution & distribution,
61 const UnsignedInteger size,
62 const Bool restart)
63 : WeightedExperimentImplementation(size)
64 , sequence_(sequence)
65 , restart_(restart)
66 , randomize_(false)
67 {
68 // Warning! The distribution must not be given to the upper class directly
69 // because the correct initialization of the sequence depends on a test on
70 // its dimension
71 setDistribution(distribution);
72 }
73
74 /* Constructor with parameters */
LowDiscrepancyExperiment(const LowDiscrepancySequence & sequence,const UnsignedInteger size,const Bool restart)75 LowDiscrepancyExperiment::LowDiscrepancyExperiment(const LowDiscrepancySequence & sequence,
76 const UnsignedInteger size,
77 const Bool restart)
78 : WeightedExperimentImplementation(size)
79 , sequence_(sequence)
80 , restart_(restart)
81 , randomize_(false)
82 {
83 // Warning! The distribution must not be given to the upper class directly
84 // because the correct initialization of the sequence depends on a test on
85 // its dimension
86 setDistribution(IndependentCopula(sequence.getDimension()));
87 }
88
89 /* Virtual constructor */
clone() const90 LowDiscrepancyExperiment * LowDiscrepancyExperiment::clone() const
91 {
92 return new LowDiscrepancyExperiment(*this);
93 }
94
95 /* String converter */
__repr__() const96 String LowDiscrepancyExperiment::__repr__() const
97 {
98 OSS oss;
99 oss << "class=" << GetClassName()
100 << " name=" << getName()
101 << " sequence=" << sequence_
102 << " distribution=" << distribution_
103 << " size=" << size_
104 << " restart=" << restart_
105 << " randomize=" << randomize_;
106 return oss;
107 }
108
__str__(const String &) const109 String LowDiscrepancyExperiment::__str__(const String & ) const
110 {
111 OSS oss;
112 oss << GetClassName()
113 << "(sequence=" << sequence_
114 << ", distribution=" << distribution_
115 << ", size" << size_
116 << ", restart=" << restart_
117 << ", randomize=" << randomize_;
118 return oss;
119 }
120
121 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const122 void LowDiscrepancyExperiment::save(Advocate & adv) const
123 {
124 WeightedExperimentImplementation::save(adv);
125 adv.saveAttribute("sequence_", sequence_);
126 adv.saveAttribute("restart_", restart_);
127 adv.saveAttribute("randomize_", randomize_);
128 }
129
130 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)131 void LowDiscrepancyExperiment::load(Advocate & adv)
132 {
133 WeightedExperimentImplementation::load(adv);
134 adv.loadAttribute("sequence_", sequence_);
135 adv.loadAttribute("restart_", restart_);
136 adv.loadAttribute("randomize_", randomize_);
137 setDistribution(distribution_);
138 }
139
140 /* Distribution accessor */
setDistribution(const Distribution & distribution)141 void LowDiscrepancyExperiment::setDistribution(const Distribution & distribution)
142 {
143 const UnsignedInteger dimension = distribution.getDimension();
144 // restart the low-discrepancy sequence if asked for or mandatory (dimension changed)
145 if (restart_ || (dimension != getDistribution().getDimension()))
146 sequence_.initialize(dimension);
147 // Build the iso-probabilistic transformation
148 // For distributions with non-indepedent copula, it resorts to using the method
149 // described in:
150 // Mathieu Cambou, Marius Hofert, Christiane Lemieux, "Quasi-Random numbers for copula models", Statistics and Computing, September 2017, Volume 27, Issue 5, pp 1307–1329
151 // preprint here: https://arxiv.org/pdf/1508.03483.pdf
152 transformation_ = DistributionTransformation(IndependentCopula(dimension), distribution);
153 WeightedExperimentImplementation::setDistribution(distribution);
154 }
155
156 /* Sequence accessor */
getSequence() const157 LowDiscrepancySequence LowDiscrepancyExperiment::getSequence() const
158 {
159 return sequence_;
160 }
161
162 /* Restart accessor */
getRestart() const163 Bool LowDiscrepancyExperiment::getRestart() const
164 {
165 return restart_;
166 }
167
setRestart(const Bool restart)168 void LowDiscrepancyExperiment::setRestart(const Bool restart)
169 {
170 restart_ = restart;
171 }
172
173
174 /* Randomization accessor */
getRandomize() const175 Bool LowDiscrepancyExperiment::getRandomize() const
176 {
177 return randomize_;
178 }
179
setRandomize(const Bool randomize)180 void LowDiscrepancyExperiment::setRandomize(const Bool randomize)
181 {
182 randomize_ = randomize;
183 }
184
185
186 /* Sample generation */
generateWithWeights(Point & weights) const187 Sample LowDiscrepancyExperiment::generateWithWeights(Point & weights) const
188 {
189 Sample sample(sequence_.generate(size_));
190 sample.setDescription(distribution_.getDescription());
191 const UnsignedInteger dimension = distribution_.getDimension();
192 Scalar tmp = -1.0;
193 if (randomize_)
194 {
195 const Point shift(RandomGenerator::Generate(dimension));
196 for (UnsignedInteger i = 0; i < size_; ++ i)
197 {
198 for (UnsignedInteger j = 0; j < dimension; ++ j)
199 {
200 // with a cyclic scrambling of the low discrepancy point as in
201 // L’Ecuyer P., Lemieux C. (2005) Recent Advances in Randomized Quasi-Monte Carlo Methods. In: Dror M., L’Ecuyer P., Szidarovszky F. (eds) Modeling Uncertainty. International Series in Operations Research & Management Science, vol 46. Springer, Boston, MA
202 sample(i, j) = std::modf(sample(i, j) + shift[j], &tmp);
203 } // j
204 } // i
205 } // randomize
206 // In-place transformation to reduce memory consumption
207 sample = transformation_(sample);
208 weights = Point(size_, 1.0 / size_);
209 return sample;
210 }
211
212 END_NAMESPACE_OPENTURNS
213