ROL
ROL_MonteCarloGeneratorDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_MONTECARLOGENERATORDEF_HPP
45#define ROL_MONTECARLOGENERATORDEF_HPP
46
47namespace ROL {
48
49template<typename Real>
50Real MonteCarloGenerator<Real>::ierf(Real input) const {
51 std::vector<Real> coeff;
52 const Real pi = ScalarTraits<Real>::pi(), zero(0), one(1), two(2), tol(1e-4);
53 Real c(1);
54 Real tmp = c * (std::sqrt(pi)/two * input);
55 Real val = tmp;
56 coeff.push_back(c);
57 int cnt = 1;
58 while (std::abs(tmp) > tol*std::abs(val)) {
59 c = zero;
60 for ( unsigned i = 0; i < coeff.size(); ++i )
61 c += coeff[i]*coeff[coeff.size()-1-i]/static_cast<Real>((i+1)*(2*i+1));
62 Real ind = static_cast<Real>(cnt);
63 tmp = c/(two*ind+one) * std::pow(std::sqrt(pi)/two*input, two*ind+one);
64 val += tmp;
65 coeff.push_back(c);
66 cnt++;
67 }
68 return val;
69}
70
71template<typename Real>
73 return static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
74}
75
76template<typename Real>
77std::vector<std::vector<Real>> MonteCarloGenerator<Real>::sample(int nSamp, bool store, bool refine) {
78 if (!refine) srand(seed_);
79 const Real zero(0), one(1), two(2), tol(0.1);
81 const int dim = (!useDist_ ? data_.size() : dist_.size());
82 std::vector<Real> pts(nSamp*dim, zero);
83 if (rank == 0) {
84 // Generate samples
85 for (int i = 0; i < nSamp; ++i) {
86 if ( !useDist_ ) {
87 for (int j = 0; j < dim; ++j) {
88 if ( use_normal_ )
89 pts[j + i*dim] = std::sqrt(two*(data_[j])[1])*ierf(two*random()-one) + (data_[j])[0];
90 else
91 pts[j + i*dim] = ((data_[j])[1]-(data_[j])[0])*random()+(data_[j])[0];
92 }
93 }
94 else {
95 for (int j = 0; j < dim; ++j) {
96 pts[j + i*dim] = (dist_[j])->invertCDF(random());
97 while (std::abs(pts[j + i*dim]) > tol*ROL_INF<Real>())
98 pts[j + i*dim] = (dist_[j])->invertCDF(random());
99 }
100 }
101 }
102 }
103 SampleGenerator<Real>::broadcast(&pts[0],nSamp*dim,0);
104 // Separate samples across processes
106 int frac = nSamp / nProc;
107 int rem = nSamp % nProc;
108 int N = frac + ((rank < rem) ? 1 : 0);
109 int offset = 0;
110 for (int i = 0; i < rank; ++i) offset += frac + ((i < rem) ? 1 : 0);
111 std::vector<std::vector<Real>> mypts;
112 std::vector<Real> pt(dim);
113 for (int i = 0; i < N; ++i) {
114 int I = offset+i;
115 for (int j = 0; j < dim; ++j) pt[j] = pts[j + I*dim];
116 mypts.push_back(pt);
117 }
118 if ( store ) {
119 std::vector<Real> mywts(N, one/static_cast<Real>(nSamp));
122 }
123 return mypts;
124}
125
126template<typename Real>
128 const std::vector<Ptr<Distribution<Real>>> &dist,
129 const Ptr<BatchManager<Real>> &bman,
130 bool use_SA, bool adaptive, int numNewSamps, int seed)
131 : SampleGenerator<Real>(bman),
132 dist_(dist),
133 nSamp_(nSamp),
134 use_normal_(false),
135 use_SA_(use_SA),
136 adaptive_(adaptive),
137 numNewSamps_(numNewSamps),
138 useDist_(true),
139 seed_(seed),
140 sum_val_(0),
141 sum_val2_(0),
142 sum_ng_(0),
143 sum_ng2_(0) {
145 ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
146 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
147 sample(nSamp_,true,false);
148}
149
150template<typename Real>
152 std::vector<std::vector<Real>> &bounds,
153 const Ptr<BatchManager<Real>> &bman,
154 bool use_SA, bool adaptive, int numNewSamps, int seed)
155 : SampleGenerator<Real>(bman),
156 nSamp_(nSamp),
157 use_normal_(false),
158 use_SA_(use_SA),
159 adaptive_(adaptive),
160 numNewSamps_(numNewSamps),
161 useDist_(false),
162 seed_(seed),
163 sum_val_(0),
164 sum_val2_(0),
165 sum_ng_(0),
166 sum_ng2_(0) {
168 ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
169 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
170 unsigned dim = bounds.size();
171 data_.clear();
172 Real tmp(0);
173 for ( unsigned j = 0; j < dim; ++j ) {
174 if ( (bounds[j])[0] > (bounds[j])[1] ) {
175 tmp = (bounds[j])[0];
176 (bounds[j])[0] = (bounds[j])[1];
177 (bounds[j])[1] = tmp;
178 data_.push_back(bounds[j]);
179 }
180 data_.push_back(bounds[j]);
181 }
182 sample(nSamp_,true,false);
183}
184
185template<typename Real>
187 const std::vector<Real> &mean,
188 const std::vector<Real> &std,
189 const Ptr<BatchManager<Real>> &bman,
190 bool use_SA, bool adaptive, int numNewSamps, int seed)
191 : SampleGenerator<Real>(bman),
192 nSamp_(nSamp),
193 use_normal_(true),
194 use_SA_(use_SA),
195 adaptive_(adaptive),
196 numNewSamps_(numNewSamps),
197 useDist_(false),
198 seed_(seed),
199 sum_val_(0),
200 sum_val2_(0),
201 sum_ng_(0),
202 sum_ng2_(0) {
204 ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
205 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
206 unsigned dim = mean.size();
207 data_.clear();
208 for ( unsigned j = 0; j < dim; ++j )
209 data_.push_back({mean[j],std[j]});
210 sample(nSamp_,true,false);
211}
212
213template<typename Real>
216 sum_val_ = static_cast<Real>(0);
217 sum_val2_ = static_cast<Real>(0);
218 sum_ng_ = static_cast<Real>(0);
219 sum_ng_ = static_cast<Real>(0);
220 if ( use_SA_ ) sample(nSamp_,true,true);
221}
222
223template<typename Real>
224Real MonteCarloGenerator<Real>::computeError( std::vector<Real> &vals ) {
225 Real err(0);
226 if ( adaptive_ && !use_SA_ ) {
227 const Real zero(0), one(1), tol(1e-8);
228 // Compute unbiased sample variance
229 int cnt = 0;
230 for ( int i = SampleGenerator<Real>::start(); i < SampleGenerator<Real>::numMySamples(); ++i ) {
231 sum_val_ += vals[cnt];
232 sum_val2_ += vals[cnt]*vals[cnt];
233 cnt++;
234 }
235 Real mymean = sum_val_ / static_cast<Real>(nSamp_);
236 Real mean = zero;
237 SampleGenerator<Real>::sumAll(&mymean,&mean,1);
238
239 Real myvar = (sum_val2_ - mean*mean)/(static_cast<Real>(nSamp_)-one);
240 Real var = zero;
241 SampleGenerator<Real>::sumAll(&myvar,&var,1);
242 // Return Monte Carlo error
243 vals.clear();
244 err = std::sqrt(var/static_cast<Real>(nSamp_))*tol;
245 }
246 else {
247 vals.clear();
248 }
249 return err;
250}
251
252template<typename Real>
254 const Vector<Real> &x ) {
255 Real err(0);
256 if ( adaptive_ && !use_SA_ ) {
257 const Real zero(0), one(1), tol(1e-4);
258 // Compute unbiased sample variance
259 int cnt = 0;
260 Real ng = zero;
261 for ( int i = SampleGenerator<Real>::start(); i < SampleGenerator<Real>::numMySamples(); ++i ) {
262 ng = (vals[cnt])->norm();
263 sum_ng_ += ng;
264 sum_ng2_ += ng*ng;
265 cnt++;
266 }
267 Real mymean = sum_ng_ / static_cast<Real>(nSamp_);
268 Real mean = zero;
269 SampleGenerator<Real>::sumAll(&mymean,&mean,1);
270
271 Real myvar = (sum_ng2_ - mean*mean)/(static_cast<Real>(nSamp_)-one);
272 Real var = zero;
273 SampleGenerator<Real>::sumAll(&myvar,&var,1);
274 // Return Monte Carlo error
275 vals.clear();
276 err = std::sqrt(var/static_cast<Real>(nSamp_))*tol;
277 }
278 else {
279 vals.clear();
280 }
281 return err;
282}
283
284template<typename Real>
286 if ( adaptive_ && !use_SA_ ) {
287 std::vector<std::vector<Real>> pts;
288 for ( int i = 0; i < SampleGenerator<Real>::numMySamples(); ++i )
289 pts.push_back(SampleGenerator<Real>::getMyPoint(i));
290 std::vector<std::vector<Real>> pts_new = sample(numNewSamps_,false,true);
291 pts.insert(pts.end(),pts_new.begin(),pts_new.end());
292 nSamp_ += numNewSamps_;
293 std::vector<Real> wts(pts.size(),static_cast<Real>(1)/static_cast<Real>(nSamp_));
297 }
298}
299
300template<typename Real>
302 return nSamp_;
303}
304
305}
306#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Real computeError(std::vector< Real > &vals) override
std::vector< std::vector< Real > > data_
MonteCarloGenerator(int nSamp, const std::vector< Ptr< Distribution< Real > > > &dist, const Ptr< BatchManager< Real > > &bman, bool use_SA=false, bool adaptive=false, int numNewSamps=0, int seed=123454321)
void update(const Vector< Real > &x) override
int numGlobalSamples(void) const override
std::vector< std::vector< Real > > sample(int nSamp, bool store=true, bool refine=false)
void setPoints(std::vector< std::vector< Real > > &p)
void setWeights(std::vector< Real > &w)
void broadcast(Real *input, int cnt, int root) const
virtual void refine(void)
void sumAll(Real *input, Real *output, int dim) const
virtual void update(const Vector< Real > &x)
Defines the linear algebra or vector space interface.
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
static constexpr Real pi() noexcept
constexpr auto dim