44 #ifndef ROL_MONTECARLOGENERATOR_HPP 45 #define ROL_MONTECARLOGENERATOR_HPP 60 std::vector<std::vector<Real> >
data_;
68 const std::vector<Teuchos::RCP<ROL::Distribution<Real> > >
dist_;
70 Real
ierf(Real input)
const {
71 std::vector<Real> coeff;
73 Real tmp = c * (std::sqrt(M_PI)/2.0 * input);
77 while (std::abs(tmp) > 1.e-4*std::abs(val)) {
79 for (
unsigned i = 0; i < coeff.size(); i++ ) {
80 c += coeff[i]*coeff[coeff.size()-1-i]/((i+1)*(2*i+1));
82 tmp = c/(2.0*(Real)cnt+1.0) * std::pow(std::sqrt(M_PI)/2.0 * input,2.0*(Real)cnt+1.0);
95 int frac = nSamp_ / nProc;
96 int rem = nSamp_ % nProc;
97 unsigned N = (unsigned)frac;
98 unsigned sumN = N*(unsigned)rank;
99 for (
int i = 0; i < rank; i++) {
108 std::vector<std::vector<Real> > pts;
111 for (
unsigned i = 0; i < N; i++ ) {
112 srand(123456*(sumN + i + 1));
114 p.resize(data_.size(),0.0);
115 for (
unsigned j = 0; j < data_.size(); j++ ) {
117 p[j] = std::sqrt(2.0*(data_[j])[1])*
ierf(2.0*((Real)rand())/((Real)RAND_MAX)-1.0) +
121 p[j] = ((data_[j])[1]-(data_[j])[0])*((Real)rand())/((Real)RAND_MAX)+(data_[j])[0];
126 p.resize(dist_.size(),0.0);
127 for (
unsigned j = 0; j < dist_.size(); j++ ) {
128 p[j] = (dist_[j])->invertCDF((Real)rand()/(Real)RAND_MAX);
130 p[j] = (dist_[j])->invertCDF((Real)rand()/(Real)RAND_MAX);
136 std::vector<Real> wts(N,1.0/((Real)nSamp_));
141 std::vector<std::vector<Real> >
sample(
int nSamp,
bool store =
true) {
146 int frac = nSamp / nProc;
147 int rem = nSamp % nProc;
148 unsigned N = (unsigned)frac;
149 unsigned sumN = N*(unsigned)rank;
150 for (
int i = 0; i < rank; i++) {
159 std::vector<std::vector<Real> > pts;
162 for (
unsigned i = 0; i < N; i++ ) {
163 srand(123456*(sumN + i + 1));
165 p.resize(data_.size(),0.0);
166 for (
unsigned j = 0; j < data_.size(); j++ ) {
168 p[j] = std::sqrt(2.0*(data_[j])[1])*
ierf(2.0*((Real)rand())/((Real)RAND_MAX)-1.0) +
172 p[j] = ((data_[j])[1]-(data_[j])[0])*((Real)rand())/((Real)RAND_MAX)+(data_[j])[0];
177 p.resize(dist_.size(),0.0);
178 for (
unsigned j = 0; j < dist_.size(); j++ ) {
179 p[j] = (dist_[j])->invertCDF((Real)rand()/(Real)RAND_MAX);
181 p[j] = (dist_[j])->invertCDF((Real)rand()/(Real)RAND_MAX);
188 std::vector<Real> wts(N,1.0/((Real)nSamp));
199 const bool use_SA =
false,
200 const bool adaptive =
false,
201 const int numNewSamps = 0)
207 numNewSamps_(numNewSamps),
218 std::vector<std::vector<Real> > &bounds,
220 const bool use_SA =
false,
221 const bool adaptive =
false,
222 const int numNewSamps = 0)
228 numNewSamps_(numNewSamps),
234 unsigned dim = bounds.size();
237 for (
unsigned j = 0; j < dim; j++ ) {
238 if ( (bounds[j])[0] > (bounds[j])[1] ) {
239 tmp = (bounds[j])[0];
240 (bounds[j])[0] = (bounds[j])[1];
241 (bounds[j])[1] = tmp;
242 data_.push_back(bounds[j]);
244 data_.push_back(bounds[j]);
250 const std::vector<Real> &mean,
251 const std::vector<Real> &std,
253 const bool use_SA =
false,
254 const bool adaptive =
false,
255 const int numNewSamps = 0 )
261 numNewSamps_(numNewSamps),
267 unsigned dim = mean.size();
269 std::vector<Real> tmp(2,0.0);
270 for (
unsigned j = 0; j < dim; j++ ) {
273 data_.push_back(tmp);
290 if ( adaptive_ && !use_SA_ ) {
294 sum_val_ += vals[cnt];
295 sum_val2_ += vals[cnt]*vals[cnt];
298 Real mymean = sum_val_ /
nSamp_;
302 Real myvar = (sum_val2_ - mean*mean)/(nSamp_-1.0);
307 return std::sqrt(var/(nSamp_))*1.e-8;
316 if ( adaptive_ && !use_SA_ ) {
321 ng = (vals[cnt])->norm();
326 Real mymean = sum_ng_ /
nSamp_;
330 Real myvar = (sum_ng2_ - mean*mean)/(nSamp_-1.0);
335 return std::sqrt(var/(nSamp_))*1.e-4;
344 if ( adaptive_ && !use_SA_ ) {
345 std::vector<std::vector<Real> > pts;
346 std::vector<Real> pt(data_.size(),0.0);
347 for (
int i = 0; i < SampleGenerator<Real>::numMySamples(); i++ ) {
351 std::vector<std::vector<Real> > pts_new =
sample(numNewSamps_,
false);
352 pts.insert(pts.end(),pts_new.begin(),pts_new.end());
354 std::vector<Real> wts(pts.size(),1.0/((Real)nSamp_));
Real ierf(Real input) const
virtual void update(const Vector< Real > &x)
MonteCarloGenerator(const int nSamp, const std::vector< Teuchos::RCP< Distribution< Real > > > &dist, const Teuchos::RCP< BatchManager< Real > > &bman, const bool use_SA=false, const bool adaptive=false, const int numNewSamps=0)
virtual std::vector< Real > getMyPoint(const int i) const
Real computeError(std::vector< Teuchos::RCP< Vector< Real > > > &vals, const Vector< Real > &x)
Defines the linear algebra or vector space interface.
void sumAll(Real *input, Real *output, int dim) const
Real computeError(std::vector< Real > &vals)
std::vector< std::vector< Real > > sample(int nSamp, bool store=true)
virtual void refine(void)
const std::vector< Teuchos::RCP< ROL::Distribution< Real > > > dist_
MonteCarloGenerator(const int nSamp, std::vector< std::vector< Real > > &bounds, const Teuchos::RCP< 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 Teuchos::RCP< BatchManager< Real > > &bman, const bool use_SA=false, const bool adaptive=false, const int numNewSamps=0)
int numBatches(void) const
void setPoints(std::vector< std::vector< Real > > &p)
static const double ROL_OVERFLOW
Platform-dependent maximum double.
void update(const Vector< Real > &x)
void setWeights(std::vector< Real > &w)
std::vector< std::vector< Real > > data_