10 #ifndef ROL_MONTECARLOGENERATORDEF_HPP
11 #define ROL_MONTECARLOGENERATORDEF_HPP
15 template<
typename Real>
17 std::vector<Real> coeff;
20 Real tmp = c * (std::sqrt(pi)/two * input);
24 while (std::abs(tmp) > tol*std::abs(val)) {
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);
37 template<
typename Real>
39 return static_cast<Real
>(rand())/static_cast<Real>(RAND_MAX);
42 template<
typename Real>
44 if (!refine) srand(seed_);
45 const Real
zero(0), one(1), two(2), tol(0.1);
47 const int dim = (!useDist_ ? data_.size() : dist_.size());
48 std::vector<Real> pts(nSamp*dim,
zero);
51 for (
int i = 0; i < nSamp; ++i) {
53 for (
int j = 0; j <
dim; ++j) {
55 pts[j + i*
dim] = std::sqrt(two*(data_[j])[1])*ierf(two*
random()-one) + (data_[j])[0];
57 pts[j + i*
dim] = ((data_[j])[1]-(data_[j])[0])*
random()+(data_[j])[0];
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());
72 int frac = nSamp / nProc;
73 int rem = nSamp % nProc;
74 int N = frac + ((rank < rem) ? 1 : 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) {
81 for (
int j = 0; j <
dim; ++j) pt[j] = pts[j + I*dim];
85 std::vector<Real> mywts(N, one/static_cast<Real>(nSamp));
92 template<
typename Real>
96 bool use_SA,
bool adaptive,
int numNewSamps,
int seed)
103 numNewSamps_(numNewSamps),
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!");
116 template<
typename Real>
118 std::vector<std::vector<Real>> &bounds,
120 bool use_SA,
bool adaptive,
int numNewSamps,
int seed)
126 numNewSamps_(numNewSamps),
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();
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]);
146 data_.push_back(bounds[j]);
151 template<
typename Real>
153 const std::vector<Real> &mean,
154 const std::vector<Real> &std,
156 bool use_SA,
bool adaptive,
int numNewSamps,
int seed)
162 numNewSamps_(numNewSamps),
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();
174 for (
unsigned j = 0; j <
dim; ++j )
175 data_.push_back({mean[j],std[j]});
176 sample(nSamp_,
true,
false);
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);
189 template<
typename Real>
192 if ( adaptive_ && !use_SA_ ) {
193 const Real
zero(0), one(1), tol(1e-8);
197 sum_val_ += vals[cnt];
198 sum_val2_ += vals[cnt]*vals[cnt];
201 Real mymean = sum_val_ /
static_cast<Real
>(nSamp_);
205 Real myvar = (sum_val2_ - mean*mean)/(static_cast<Real>(nSamp_)-one);
210 err = std::sqrt(var/static_cast<Real>(nSamp_))*tol;
218 template<
typename Real>
222 if ( adaptive_ && !use_SA_ ) {
223 const Real
zero(0), one(1), tol(1e-4);
228 ng = (vals[cnt])->norm();
233 Real mymean = sum_ng_ /
static_cast<Real
>(nSamp_);
237 Real myvar = (sum_ng2_ - mean*mean)/(static_cast<Real>(nSamp_)-one);
242 err = std::sqrt(var/static_cast<Real>(nSamp_))*tol;
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 )
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_));
266 template<
typename Real>
Real ierf(Real input) const
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.
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)
int numBatches(void) const
void setPoints(std::vector< std::vector< Real > > &p)
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_