44 #ifndef ROL_MONTECARLOGENERATOR_HPP
45 #define ROL_MONTECARLOGENERATOR_HPP
60 std::vector<std::vector<Real> >
data_;
68 const std::vector<ROL::Ptr<Distribution<Real> > >
dist_;
72 Real
ierf(Real input)
const {
73 std::vector<Real> coeff;
76 Real tmp = c * (std::sqrt(pi)/two * input);
80 while (std::abs(tmp) > tol*std::abs(val)) {
82 for (
unsigned i = 0; i < coeff.size(); i++ ) {
83 c += coeff[i]*coeff[coeff.size()-1-i]/((i+1)*(2*i+1));
85 Real ind =
static_cast<Real
>(cnt);
86 tmp = c/(two*ind+one) * std::pow(std::sqrt(pi)/two*input, two*ind+one);
95 return static_cast<Real
>(rand())/static_cast<Real>(RAND_MAX);
98 std::vector<std::vector<Real> >
sample(
int nSamp,
bool store =
true) {
100 const Real
zero(0), one(1), two(2), tol(0.1);
103 std::vector<Real> pts(nSamp*dim,
zero);
106 for (
int i = 0; i < nSamp; ++i) {
108 for (
int j = 0; j < dim; ++j) {
118 for (
int j = 0; j < dim; ++j) {
120 while (std::abs(pts[j + i*dim]) > tol*ROL_INF<Real>()) {
130 int frac = nSamp / nProc;
131 int rem = nSamp % nProc;
132 int N = frac + ((rank < rem) ? 1 : 0);
134 for (
int i = 0; i < rank; ++i) {
135 offset += frac + ((i < rem) ? 1 : 0);
137 std::vector<std::vector<Real> > mypts;
138 std::vector<Real> pt(dim);
139 for (
int i = 0; i < N; ++i) {
141 for (
int j = 0; j < dim; ++j) {
142 pt[j] = pts[j + I*dim];
147 std::vector<Real> mywts(N, one/static_cast<Real>(nSamp));
162 const bool use_SA =
false,
163 const bool adaptive =
false,
164 const int numNewSamps = 0,
165 const int seed = 123454321)
180 ROL_TEST_FOR_EXCEPTION(
nSamp_ < nProc, std::invalid_argument,
181 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
186 std::vector<std::vector<Real> > &bounds,
188 const bool use_SA =
false,
189 const bool adaptive =
false,
190 const int numNewSamps = 0,
191 const int seed = 123454321)
205 ROL_TEST_FOR_EXCEPTION(
nSamp_ < nProc, std::invalid_argument,
206 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
207 unsigned dim = bounds.size();
210 for (
unsigned j = 0; j < dim; j++ ) {
211 if ( (bounds[j])[0] > (bounds[j])[1] ) {
212 tmp = (bounds[j])[0];
213 (bounds[j])[0] = (bounds[j])[1];
214 (bounds[j])[1] = tmp;
215 data_.push_back(bounds[j]);
217 data_.push_back(bounds[j]);
223 const std::vector<Real> &mean,
224 const std::vector<Real> &std,
226 const bool use_SA =
false,
227 const bool adaptive =
false,
228 const int numNewSamps = 0,
229 const int seed = 123454321)
243 ROL_TEST_FOR_EXCEPTION(
nSamp_ < nProc, std::invalid_argument,
244 ">>> ERROR (ROL::MonteCarloGenerator): Total number of samples is less than the number of batches!");
245 unsigned dim = mean.size();
247 std::vector<Real> tmp(2,static_cast<Real>(0));
248 for (
unsigned j = 0; j < dim; j++ ) {
251 data_.push_back(tmp);
270 Real
zero(0), one(1), tol(1e-8);
287 return std::sqrt(var/(
nSamp_))*tol;
291 return static_cast<Real
>(0);
297 Real
zero(0), one(1), tol(1e-4);
302 ng = (vals[cnt])->norm();
316 return std::sqrt(var/(
nSamp_))*tol;
320 return static_cast<Real
>(0);
326 Real
zero(0), one(1);
327 std::vector<std::vector<Real> > pts;
329 for (
int i = 0; i < SampleGenerator<Real>::numMySamples(); i++ ) {
334 pts.insert(pts.end(),pts_new.begin(),pts_new.end());
336 std::vector<Real> wts(pts.size(),one/((Real)
nSamp_));
Real ierf(Real input) const
virtual void update(const Vector< Real > &x)
Real computeError(std::vector< ROL::Ptr< Vector< Real > > > &vals, const Vector< Real > &x)
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, const int seed=123454321)
virtual std::vector< Real > getMyPoint(const int i) const
void broadcast(Real *input, int cnt, int root) const
Defines the linear algebra or vector space interface.
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, const int seed=123454321)
void sumAll(Real *input, Real *output, int dim) const
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
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)
int numGlobalSamples(void) const
int numBatches(void) const
void setPoints(std::vector< std::vector< Real > > &p)
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, const int seed=123454321)
const std::vector< ROL::Ptr< Distribution< Real > > > dist_
void update(const Vector< Real > &x)
void setWeights(std::vector< Real > &w)
std::vector< std::vector< Real > > data_