ROL
ROL_MonteCarloGeneratorDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_MONTECARLOGENERATORDEF_HPP
11 #define ROL_MONTECARLOGENERATORDEF_HPP
12 
13 namespace ROL {
14 
15 template<typename Real>
16 Real MonteCarloGenerator<Real>::ierf(Real input) const {
17  std::vector<Real> coeff;
18  const Real pi = ScalarTraits<Real>::pi(), zero(0), one(1), two(2), tol(1e-4);
19  Real c(1);
20  Real tmp = c * (std::sqrt(pi)/two * input);
21  Real val = tmp;
22  coeff.push_back(c);
23  int cnt = 1;
24  while (std::abs(tmp) > tol*std::abs(val)) {
25  c = zero;
26  for ( unsigned i = 0; i < coeff.size(); ++i )
27  c += coeff[i]*coeff[coeff.size()-1-i]/static_cast<Real>((i+1)*(2*i+1));
28  Real ind = static_cast<Real>(cnt);
29  tmp = c/(two*ind+one) * std::pow(std::sqrt(pi)/two*input, two*ind+one);
30  val += tmp;
31  coeff.push_back(c);
32  cnt++;
33  }
34  return val;
35 }
36 
37 template<typename Real>
39  return static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
40 }
41 
42 template<typename Real>
43 std::vector<std::vector<Real>> MonteCarloGenerator<Real>::sample(int nSamp, bool store, bool refine) {
44  if (!refine) srand(seed_);
45  const Real zero(0), one(1), two(2), tol(0.1);
46  int rank = SampleGenerator<Real>::batchID();
47  const int dim = (!useDist_ ? data_.size() : dist_.size());
48  std::vector<Real> pts(nSamp*dim, zero);
49  if (rank == 0) {
50  // Generate samples
51  for (int i = 0; i < nSamp; ++i) {
52  if ( !useDist_ ) {
53  for (int j = 0; j < dim; ++j) {
54  if ( use_normal_ )
55  pts[j + i*dim] = std::sqrt(two*(data_[j])[1])*ierf(two*random()-one) + (data_[j])[0];
56  else
57  pts[j + i*dim] = ((data_[j])[1]-(data_[j])[0])*random()+(data_[j])[0];
58  }
59  }
60  else {
61  for (int j = 0; j < dim; ++j) {
62  pts[j + i*dim] = (dist_[j])->invertCDF(random());
63  while (std::abs(pts[j + i*dim]) > tol*ROL_INF<Real>())
64  pts[j + i*dim] = (dist_[j])->invertCDF(random());
65  }
66  }
67  }
68  }
69  SampleGenerator<Real>::broadcast(&pts[0],nSamp*dim,0);
70  // Separate samples across processes
72  int frac = nSamp / nProc;
73  int rem = nSamp % nProc;
74  int N = frac + ((rank < rem) ? 1 : 0);
75  int offset = 0;
76  for (int i = 0; i < rank; ++i) offset += frac + ((i < rem) ? 1 : 0);
77  std::vector<std::vector<Real>> mypts;
78  std::vector<Real> pt(dim);
79  for (int i = 0; i < N; ++i) {
80  int I = offset+i;
81  for (int j = 0; j < dim; ++j) pt[j] = pts[j + I*dim];
82  mypts.push_back(pt);
83  }
84  if ( store ) {
85  std::vector<Real> mywts(N, one/static_cast<Real>(nSamp));
88  }
89  return mypts;
90 }
91 
92 template<typename Real>
94  const std::vector<Ptr<Distribution<Real>>> &dist,
95  const Ptr<BatchManager<Real>> &bman,
96  bool use_SA, bool adaptive, int numNewSamps, int seed)
97  : SampleGenerator<Real>(bman),
98  dist_(dist),
99  nSamp_(nSamp),
100  use_normal_(false),
101  use_SA_(use_SA),
102  adaptive_(adaptive),
103  numNewSamps_(numNewSamps),
104  useDist_(true),
105  seed_(seed),
106  sum_val_(0),
107  sum_val2_(0),
108  sum_ng_(0),
109  sum_ng2_(0) {
110  int nProc = SampleGenerator<Real>::numBatches();
111  ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
112  ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
113  sample(nSamp_,true,false);
114 }
115 
116 template<typename Real>
118  std::vector<std::vector<Real>> &bounds,
119  const Ptr<BatchManager<Real>> &bman,
120  bool use_SA, bool adaptive, int numNewSamps, int seed)
121  : SampleGenerator<Real>(bman),
122  nSamp_(nSamp),
123  use_normal_(false),
124  use_SA_(use_SA),
125  adaptive_(adaptive),
126  numNewSamps_(numNewSamps),
127  useDist_(false),
128  seed_(seed),
129  sum_val_(0),
130  sum_val2_(0),
131  sum_ng_(0),
132  sum_ng2_(0) {
133  int nProc = SampleGenerator<Real>::numBatches();
134  ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
135  ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
136  unsigned dim = bounds.size();
137  data_.clear();
138  Real tmp(0);
139  for ( unsigned j = 0; j < dim; ++j ) {
140  if ( (bounds[j])[0] > (bounds[j])[1] ) {
141  tmp = (bounds[j])[0];
142  (bounds[j])[0] = (bounds[j])[1];
143  (bounds[j])[1] = tmp;
144  data_.push_back(bounds[j]);
145  }
146  data_.push_back(bounds[j]);
147  }
148  sample(nSamp_,true,false);
149 }
150 
151 template<typename Real>
153  const std::vector<Real> &mean,
154  const std::vector<Real> &std,
155  const Ptr<BatchManager<Real>> &bman,
156  bool use_SA, bool adaptive, int numNewSamps, int seed)
157  : SampleGenerator<Real>(bman),
158  nSamp_(nSamp),
159  use_normal_(true),
160  use_SA_(use_SA),
161  adaptive_(adaptive),
162  numNewSamps_(numNewSamps),
163  useDist_(false),
164  seed_(seed),
165  sum_val_(0),
166  sum_val2_(0),
167  sum_ng_(0),
168  sum_ng2_(0) {
169  int nProc = SampleGenerator<Real>::numBatches();
170  ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
171  ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
172  unsigned dim = mean.size();
173  data_.clear();
174  for ( unsigned j = 0; j < dim; ++j )
175  data_.push_back({mean[j],std[j]});
176  sample(nSamp_,true,false);
177 }
178 
179 template<typename Real>
182  sum_val_ = static_cast<Real>(0);
183  sum_val2_ = static_cast<Real>(0);
184  sum_ng_ = static_cast<Real>(0);
185  sum_ng_ = static_cast<Real>(0);
186  if ( use_SA_ ) sample(nSamp_,true,true);
187 }
188 
189 template<typename Real>
190 Real MonteCarloGenerator<Real>::computeError( std::vector<Real> &vals ) {
191  Real err(0);
192  if ( adaptive_ && !use_SA_ ) {
193  const Real zero(0), one(1), tol(1e-8);
194  // Compute unbiased sample variance
195  int cnt = 0;
196  for ( int i = SampleGenerator<Real>::start(); i < SampleGenerator<Real>::numMySamples(); ++i ) {
197  sum_val_ += vals[cnt];
198  sum_val2_ += vals[cnt]*vals[cnt];
199  cnt++;
200  }
201  Real mymean = sum_val_ / static_cast<Real>(nSamp_);
202  Real mean = zero;
203  SampleGenerator<Real>::sumAll(&mymean,&mean,1);
204 
205  Real myvar = (sum_val2_ - mean*mean)/(static_cast<Real>(nSamp_)-one);
206  Real var = zero;
207  SampleGenerator<Real>::sumAll(&myvar,&var,1);
208  // Return Monte Carlo error
209  vals.clear();
210  err = std::sqrt(var/static_cast<Real>(nSamp_))*tol;
211  }
212  else {
213  vals.clear();
214  }
215  return err;
216 }
217 
218 template<typename Real>
220  const Vector<Real> &x ) {
221  Real err(0);
222  if ( adaptive_ && !use_SA_ ) {
223  const Real zero(0), one(1), tol(1e-4);
224  // Compute unbiased sample variance
225  int cnt = 0;
226  Real ng = zero;
227  for ( int i = SampleGenerator<Real>::start(); i < SampleGenerator<Real>::numMySamples(); ++i ) {
228  ng = (vals[cnt])->norm();
229  sum_ng_ += ng;
230  sum_ng2_ += ng*ng;
231  cnt++;
232  }
233  Real mymean = sum_ng_ / static_cast<Real>(nSamp_);
234  Real mean = zero;
235  SampleGenerator<Real>::sumAll(&mymean,&mean,1);
236 
237  Real myvar = (sum_ng2_ - mean*mean)/(static_cast<Real>(nSamp_)-one);
238  Real var = zero;
239  SampleGenerator<Real>::sumAll(&myvar,&var,1);
240  // Return Monte Carlo error
241  vals.clear();
242  err = std::sqrt(var/static_cast<Real>(nSamp_))*tol;
243  }
244  else {
245  vals.clear();
246  }
247  return err;
248 }
249 
250 template<typename Real>
252  if ( adaptive_ && !use_SA_ ) {
253  std::vector<std::vector<Real>> pts;
254  for ( int i = 0; i < SampleGenerator<Real>::numMySamples(); ++i )
255  pts.push_back(SampleGenerator<Real>::getMyPoint(i));
256  std::vector<std::vector<Real>> pts_new = sample(numNewSamps_,false,true);
257  pts.insert(pts.end(),pts_new.begin(),pts_new.end());
258  nSamp_ += numNewSamps_;
259  std::vector<Real> wts(pts.size(),static_cast<Real>(1)/static_cast<Real>(nSamp_));
263  }
264 }
265 
266 template<typename Real>
268  return nSamp_;
269 }
270 
271 }
272 #endif
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)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
void broadcast(Real *input, int cnt, int root) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
static constexpr Real pi() noexcept
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
Definition: example_05.cpp:15
int numBatches(void) const
void setPoints(std::vector< std::vector< Real > > &p)
constexpr auto dim
std::vector< std::vector< Real > > sample(int nSamp, bool store=true, bool refine=false)
void setWeights(std::vector< Real > &w)
std::vector< std::vector< Real > > data_