ROL
ROL_MonteCarloGenerator.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_MONTECARLOGENERATOR_HPP
45 #define ROL_MONTECARLOGENERATOR_HPP
46 
47 #include "ROL_SampleGenerator.hpp"
48 #include "ROL_Distribution.hpp"
49 
50 namespace ROL {
51 
52 template<class Real>
53 class MonteCarloGenerator : public SampleGenerator<Real> {
54 private:
55  int nSamp_;
56  const bool use_normal_;
57  const bool use_SA_;
58  const bool adaptive_;
59  const int numNewSamps_;
60  std::vector<std::vector<Real> > data_;
61 
62  Real sum_val_;
63  Real sum_val2_;
64  Real sum_ng_;
65  Real sum_ng2_;
66 
67  const bool useDist_;
68  const std::vector<ROL::Ptr<Distribution<Real> > > dist_;
69 
70  Real ierf(Real input) const {
71  std::vector<Real> coeff;
72  Real pi = ROL::ScalarTraits<Real>::pi(), zero(0), one(1), two(2), tol(1e-4);
73  Real c(1);
74  Real tmp = c * (std::sqrt(pi)/two * input);
75  Real val = tmp;
76  coeff.push_back(c);
77  int cnt = 1;
78  while (std::abs(tmp) > tol*std::abs(val)) {
79  c = zero;
80  for ( unsigned i = 0; i < coeff.size(); i++ ) {
81  c += coeff[i]*coeff[coeff.size()-1-i]/((i+1)*(2*i+1));
82  }
83  Real ind = static_cast<Real>(cnt);
84  tmp = c/(two*ind+one) * std::pow(std::sqrt(pi)/two*input, two*ind+one);
85  val += tmp;
86  coeff.push_back(c);
87  cnt++;
88  }
89  return val;
90  }
91 
92  Real random(void) const {
93  return static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
94  }
95 
96  std::vector<std::vector<Real> > sample(int nSamp, bool store = true) {
97  srand(123454321);
98  const Real zero(0), one(1), two(2), tol(0.1);
99  int rank = SampleGenerator<Real>::batchID();
100  const int dim = (!useDist_ ? data_.size() : dist_.size());
101  std::vector<Real> pts(nSamp*dim, zero);
102  if (rank == 0) {
103  // Generate samples
104  for (int i = 0; i < nSamp; ++i) {
105  if ( !useDist_ ) {
106  for (int j = 0; j < dim; ++j) {
107  if ( use_normal_ ) {
108  pts[j + i*dim] = std::sqrt(two*(data_[j])[1])*ierf(two*random()-one) + (data_[j])[0];
109  }
110  else {
111  pts[j + i*dim] = ((data_[j])[1]-(data_[j])[0])*random()+(data_[j])[0];
112  }
113  }
114  }
115  else {
116  for (int j = 0; j < dim; ++j) {
117  pts[j + i*dim] = (dist_[j])->invertCDF(random());
118  while (std::abs(pts[j + i*dim]) > tol*ROL_INF<Real>()) {
119  pts[j + i*dim] = (dist_[j])->invertCDF(random());
120  }
121  }
122  }
123  }
124  }
125  SampleGenerator<Real>::broadcast(&pts[0],nSamp*dim,0);
126  // Separate samples across processes
127  int nProc = SampleGenerator<Real>::numBatches();
128  int frac = nSamp / nProc;
129  int rem = nSamp % nProc;
130  int N = frac + ((rank < rem) ? 1 : 0);
131  int offset = 0;
132  for (int i = 0; i < rank; ++i) {
133  offset += frac + ((i < rem) ? 1 : 0);
134  }
135  std::vector<std::vector<Real> > mypts;
136  std::vector<Real> pt(dim);
137  for (int i = 0; i < N; ++i) {
138  int I = offset+i;
139  for (int j = 0; j < dim; ++j) {
140  pt[j] = pts[j + I*dim];
141  }
142  mypts.push_back(pt);
143  }
144  if ( store ) {
145  std::vector<Real> mywts(N, one/static_cast<Real>(nSamp));
148  }
149  return mypts;
150  }
151 
152  void sample(void) {
153  sample(nSamp_,true);
154  }
155 
156 public:
157  MonteCarloGenerator(const int nSamp,
158  const std::vector<ROL::Ptr<Distribution<Real> > > &dist,
159  const ROL::Ptr<BatchManager<Real> > &bman,
160  const bool use_SA = false,
161  const bool adaptive = false,
162  const int numNewSamps = 0)
163  : SampleGenerator<Real>(bman),
164  nSamp_(nSamp),
165  use_normal_(false),
166  use_SA_(use_SA),
167  adaptive_(adaptive),
168  numNewSamps_(numNewSamps),
169  sum_val_(0),
170  sum_val2_(0),
171  sum_ng_(0),
172  sum_ng2_(0),
173  useDist_(true),
174  dist_(dist) {
175  int nProc = SampleGenerator<Real>::numBatches();
176  ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
177  ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
178  sample();
179  }
180 
181  MonteCarloGenerator(const int nSamp,
182  std::vector<std::vector<Real> > &bounds,
183  const ROL::Ptr<BatchManager<Real> > &bman,
184  const bool use_SA = false,
185  const bool adaptive = false,
186  const int numNewSamps = 0)
187  : SampleGenerator<Real>(bman),
188  nSamp_(nSamp),
189  use_normal_(false),
190  use_SA_(use_SA),
191  adaptive_(adaptive),
192  numNewSamps_(numNewSamps),
193  sum_val_(0),
194  sum_val2_(0),
195  sum_ng_(0),
196  sum_ng2_(0),
197  useDist_(false) {
198  int nProc = SampleGenerator<Real>::numBatches();
199  ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
200  ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
201  unsigned dim = bounds.size();
202  data_.clear();
203  Real tmp(0);
204  for ( unsigned j = 0; j < dim; j++ ) {
205  if ( (bounds[j])[0] > (bounds[j])[1] ) {
206  tmp = (bounds[j])[0];
207  (bounds[j])[0] = (bounds[j])[1];
208  (bounds[j])[1] = tmp;
209  data_.push_back(bounds[j]);
210  }
211  data_.push_back(bounds[j]);
212  }
213  sample();
214  }
215 
216  MonteCarloGenerator(const int nSamp,
217  const std::vector<Real> &mean,
218  const std::vector<Real> &std,
219  const ROL::Ptr<BatchManager<Real> > &bman,
220  const bool use_SA = false,
221  const bool adaptive = false,
222  const int numNewSamps = 0 )
223  : SampleGenerator<Real>(bman),
224  nSamp_(nSamp),
225  use_normal_(true),
226  use_SA_(use_SA),
227  adaptive_(adaptive),
228  numNewSamps_(numNewSamps),
229  sum_val_(0),
230  sum_val2_(0),
231  sum_ng_(0),
232  sum_ng2_(0),
233  useDist_(false) {
234  int nProc = SampleGenerator<Real>::numBatches();
235  ROL_TEST_FOR_EXCEPTION( nSamp_ < nProc, std::invalid_argument,
236  ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
237  unsigned dim = mean.size();
238  data_.clear();
239  std::vector<Real> tmp(2,static_cast<Real>(0));
240  for ( unsigned j = 0; j < dim; j++ ) {
241  tmp[0] = mean[j];
242  tmp[1] = std[j];
243  data_.push_back(tmp);
244  }
245  sample();
246  }
247 
248  void update( const Vector<Real> &x ) {
250  Real zero(0);
251  sum_val_ = zero;
252  sum_val2_ = zero;
253  sum_ng_ = zero;
254  sum_ng_ = zero;
255  if ( use_SA_ ) {
256  sample();
257  }
258  }
259 
260  Real computeError( std::vector<Real> &vals ) {
261  if ( adaptive_ && !use_SA_ ) {
262  Real zero(0), one(1), tol(1e-8);
263  // Compute unbiased sample variance
264  int cnt = 0;
265  for ( int i = SampleGenerator<Real>::start(); i < SampleGenerator<Real>::numMySamples(); i++ ) {
266  sum_val_ += vals[cnt];
267  sum_val2_ += vals[cnt]*vals[cnt];
268  cnt++;
269  }
270  Real mymean = sum_val_ / nSamp_;
271  Real mean = zero;
272  SampleGenerator<Real>::sumAll(&mymean,&mean,1);
273 
274  Real myvar = (sum_val2_ - mean*mean)/(nSamp_-one);
275  Real var = zero;
276  SampleGenerator<Real>::sumAll(&myvar,&var,1);
277  // Return Monte Carlo error
278  vals.clear();
279  return std::sqrt(var/(nSamp_))*tol;
280  }
281  else {
282  vals.clear();
283  return static_cast<Real>(0);
284  }
285  }
286 
287  Real computeError( std::vector<ROL::Ptr<Vector<Real> > > &vals, const Vector<Real> &x ) {
288  if ( adaptive_ && !use_SA_ ) {
289  Real zero(0), one(1), tol(1e-4);
290  // Compute unbiased sample variance
291  int cnt = 0;
292  Real ng = zero;
293  for ( int i = SampleGenerator<Real>::start(); i < SampleGenerator<Real>::numMySamples(); i++ ) {
294  ng = (vals[cnt])->norm();
295  sum_ng_ += ng;
296  sum_ng2_ += ng*ng;
297  cnt++;
298  }
299  Real mymean = sum_ng_ / nSamp_;
300  Real mean = zero;
301  SampleGenerator<Real>::sumAll(&mymean,&mean,1);
302 
303  Real myvar = (sum_ng2_ - mean*mean)/(nSamp_-one);
304  Real var = zero;
305  SampleGenerator<Real>::sumAll(&myvar,&var,1);
306  // Return Monte Carlo error
307  vals.clear();
308  return std::sqrt(var/(nSamp_))*tol;
309  }
310  else {
311  vals.clear();
312  return static_cast<Real>(0);
313  }
314  }
315 
316  void refine(void) {
317  if ( adaptive_ && !use_SA_ ) {
318  Real zero(0), one(1);
319  std::vector<std::vector<Real> > pts;
320  std::vector<Real> pt(data_.size(),zero);
321  for ( int i = 0; i < SampleGenerator<Real>::numMySamples(); i++ ) {
323  pts.push_back(pt);
324  }
325  std::vector<std::vector<Real> > pts_new = sample(numNewSamps_,false);
326  pts.insert(pts.end(),pts_new.begin(),pts_new.end());
327  nSamp_ += numNewSamps_;
328  std::vector<Real> wts(pts.size(),one/((Real)nSamp_));
332  }
333  }
334 
335  int numGlobalSamples(void) const {
336  return nSamp_;
337  }
338 
339 };
340 
341 }
342 
343 #endif
virtual std::vector< Real > getMyPoint(const int i) const
virtual void update(const Vector< Real > &x)
Real computeError(std::vector< ROL::Ptr< Vector< Real > > > &vals, const Vector< Real > &x)
void sumAll(Real *input, Real *output, int dim) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
MonteCarloGenerator(const int nSamp, const std::vector< ROL::Ptr< Distribution< Real > > > &dist, const ROL::Ptr< BatchManager< Real > > &bman, const bool use_SA=false, const bool adaptive=false, const int numNewSamps=0)
Real computeError(std::vector< Real > &vals)
static constexpr Real pi() noexcept
std::vector< std::vector< Real > > sample(int nSamp, bool store=true)
virtual void refine(void)
MonteCarloGenerator(const int nSamp, std::vector< std::vector< Real > > &bounds, const ROL::Ptr< BatchManager< Real > > &bman, const bool use_SA=false, const bool adaptive=false, const int numNewSamps=0)
MonteCarloGenerator(const int nSamp, const std::vector< Real > &mean, const std::vector< Real > &std, const ROL::Ptr< BatchManager< Real > > &bman, const bool use_SA=false, const bool adaptive=false, const int numNewSamps=0)
void setPoints(std::vector< std::vector< Real > > &p)
const std::vector< ROL::Ptr< Distribution< Real > > > dist_
void update(const Vector< Real > &x)
void broadcast(Real *input, int cnt, int root) const
void setWeights(std::vector< Real > &w)
std::vector< std::vector< Real > > data_