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);
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);
129 while (std::abs(p[j]) > 0.1*ROL::ROL_OVERFLOW<Real>()) {
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);
180 while (std::abs(p[j]) > 0.1*ROL::ROL_OVERFLOW<Real>()) {
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)
215 TEUCHOS_TEST_FOR_EXCEPTION(
nSamp_ < nProc, std::invalid_argument,
216 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
221 std::vector<std::vector<Real> > &bounds,
223 const bool use_SA =
false,
224 const bool adaptive =
false,
225 const int numNewSamps = 0)
238 TEUCHOS_TEST_FOR_EXCEPTION(
nSamp_ < nProc, std::invalid_argument,
239 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
240 unsigned dim = bounds.size();
243 for (
unsigned j = 0; j < dim; j++ ) {
244 if ( (bounds[j])[0] > (bounds[j])[1] ) {
245 tmp = (bounds[j])[0];
246 (bounds[j])[0] = (bounds[j])[1];
247 (bounds[j])[1] = tmp;
248 data_.push_back(bounds[j]);
250 data_.push_back(bounds[j]);
256 const std::vector<Real> &mean,
257 const std::vector<Real> &std,
259 const bool use_SA =
false,
260 const bool adaptive =
false,
261 const int numNewSamps = 0 )
274 TEUCHOS_TEST_FOR_EXCEPTION(
nSamp_ < nProc, std::invalid_argument,
275 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
276 unsigned dim = mean.size();
278 std::vector<Real> tmp(2,0.0);
279 for (
unsigned j = 0; j < dim; j++ ) {
282 data_.push_back(tmp);
316 return std::sqrt(var/(
nSamp_))*1.e-8;
330 ng = (vals[cnt])->norm();
344 return std::sqrt(var/(
nSamp_))*1.e-4;
354 std::vector<std::vector<Real> > pts;
355 std::vector<Real> pt(
data_.size(),0.0);
356 for (
int i = 0; i < SampleGenerator<Real>::numMySamples(); i++ ) {
361 pts.insert(pts.end(),pts_new.begin(),pts_new.end());
363 std::vector<Real> wts(pts.size(),1.0/((Real)
nSamp_));
virtual std::vector< Real > getMyPoint(const int i) const
virtual void update(const Vector< Real > &x)
Real ierf(Real input) const
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)
void sumAll(Real *input, Real *output, int dim) const
Real computeError(std::vector< Teuchos::RCP< Vector< Real > > > &vals, const Vector< Real > &x)
Defines the linear algebra or vector space interface.
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_
int numBatches(void) const
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)
void setPoints(std::vector< std::vector< Real > > &p)
void update(const Vector< Real > &x)
void setWeights(std::vector< Real > &w)
std::vector< std::vector< Real > > data_