44 #ifndef ROL_CDFOBJECTIVE_H 45 #define ROL_CDFOBJECTIVE_H 51 #include "Teuchos_RCP.hpp" 60 Teuchos::RCP<BatchManager<Real> >
bman_;
63 std::vector<Teuchos::RCP<Distribution<Real> > >
dist_;
83 pts_.clear(); pts_.resize(numPoints_,0.);
84 wts_.clear(); wts_.resize(numPoints_,0.);
85 wts_[0] = 0.1527533871307258; pts_[0] = -0.0765265211334973;
86 wts_[1] = 0.1527533871307258; pts_[1] = 0.0765265211334973;
87 wts_[2] = 0.1491729864726037; pts_[2] = -0.2277858511416451;
88 wts_[3] = 0.1491729864726037; pts_[3] = 0.2277858511416451;
89 wts_[4] = 0.1420961093183820; pts_[4] = -0.3737060887154195;
90 wts_[5] = 0.1420961093183820; pts_[5] = 0.3737060887154195;
91 wts_[6] = 0.1316886384491766; pts_[6] = -0.5108670019508271;
92 wts_[7] = 0.1316886384491766; pts_[7] = 0.5108670019508271;
93 wts_[8] = 0.1181945319615184; pts_[8] = -0.6360536807265150;
94 wts_[9] = 0.1181945319615184; pts_[9] = 0.6360536807265150;
95 wts_[10] = 0.1019301198172404; pts_[10] = -0.7463319064601508;
96 wts_[11] = 0.1019301198172404; pts_[11] = 0.7463319064601508;
97 wts_[12] = 0.0832767415767048; pts_[12] = -0.8391169718222188;
98 wts_[13] = 0.0832767415767048; pts_[13] = 0.8391169718222188;
99 wts_[14] = 0.0626720483341091; pts_[14] = -0.9122344282513259;
100 wts_[15] = 0.0626720483341091; pts_[15] = 0.9122344282513259;
101 wts_[16] = 0.0406014298003869; pts_[16] = -0.9639719272779138;
102 wts_[17] = 0.0406014298003869; pts_[17] = 0.9639719272779138;
103 wts_[18] = 0.0176140071391521; pts_[18] = -0.9931285991850949;
104 wts_[19] = 0.0176140071391521; pts_[19] = 0.9931285991850949;
107 pts_[i] += 1.; pts_[i] *= 0.5;
115 Real val = 0., hs = 0., xpt = 0., xwt = 0., sum = 0.;
116 for (
int k = 0; k < numSamples; k++) {
118 hs = 0.5 * (1. + erf((loc-xpt)/(sqrt2_*scale_)));
121 bman_->sumAll(&val,&sum,1);
125 Real
gradientCDF(std::vector<Real> &gradx, std::vector<Real> &gradp,
126 const int dim,
const Real loc,
130 gradx.resize(numSamples,0.); gradp.resize(numSamples,0.);
131 Real val = 0., hs = 0., xpt = 0., xwt = 0., sum = 0.;
132 for (
int k = 0; k < numSamples; k++) {
134 hs = 0.5 * (1. + erf((loc-xpt)/(sqrt2_*scale_)));
136 gradx[k] = -(xwt/(sqrt2_*sqrtpi_*
scale_))
137 * std::exp(-std::pow((loc-xpt)/(sqrt2_*scale_),2));
140 bman_->sumAll(&val,&sum,1);
144 Real
hessVecCDF(std::vector<Real> &hvxx, std::vector<Real> &hvxp, std::vector<Real> &hvpx,
145 std::vector<Real> &gradx, std::vector<Real> &gradp,
146 Real &sumx, Real &sump,
147 const int dim,
const Real loc,
153 hvxx.resize(numSamples,0.); hvxp.resize(numSamples,0.); hvpx.resize(numSamples,0.);
154 gradx.resize(numSamples,0.); gradp.resize(numSamples,0.);
155 sumx = 0.; sump = 0.;
156 std::vector<Real> psum(3,0.0), out(3,0.0);
157 Real val = 0., hs = 0., dval = 0., scale3 = std::pow(scale_,3);
158 Real xpt = 0., xwt = 0., vpt = 0., vwt = 0.;
159 for (
int k = 0; k < numSamples; k++) {
162 hs = 0.5 * (1. + erf((loc-xpt)/(sqrt2_*scale_)));
164 dval = std::exp(-std::pow((loc-xpt)/(sqrt2_*scale_),2));
165 gradx[k] = -(xwt/(sqrt2_*sqrtpi_*
scale_)) * dval;
167 hvxx[k] = -(xwt/(sqrt2_*sqrtpi_*scale3)) * dval * (loc-xpt) * vpt;
168 hvxp[k] = -dval/(sqrt2_*sqrtpi_*
scale_)*vwt;
169 hvpx[k] = -dval/(sqrt2_*sqrtpi_*
scale_)*vpt;
170 psum[1] += vpt*gradx[k];
171 psum[2] += vwt*gradp[k];
173 bman_->sumAll(&psum[0],&out[0],3);
174 val = out[0]; sumx = out[1]; sump = out[2];
181 const Real scale = 1.e-2,
182 const bool optProb =
true,
const bool optAtom =
true)
183 :
Objective<Real>(), bman_(bman), dist_(dist), dimension_(dist.size()),
184 scale_(scale), sqrt2_(std::sqrt(2.)), sqrtpi_(std::sqrt(M_PI)),
185 optProb_(optProb), optAtom_(optAtom) {
186 lowerBound_.resize(dimension_,0.0);
187 upperBound_.resize(dimension_,0.0);
189 lowerBound_[i] = dist[i]->lowerBound();
190 upperBound_[i] = dist[i]->upperBound();
199 Real val = 0., diff = 0., pt = 0., wt = 0., meas = 0., lb = 0.;
202 meas = (upperBound_[d] - lb);
204 pt = meas*pts_[k] + lb;
206 diff = (
valueCDF(d,pt,prob,atom)-dist_[d]->evaluateCDF(pt));
207 val += wt*std::pow(diff,2);
219 std::vector<Real> gradx(numSamples,0.), gradp(numSamples,0.);
220 Real diff = 0., pt = 0., wt = 0., meas = 0., lb = 0., val = 0.;
221 std::vector<Real> val_wt(numSamples,0.), tmp(dimension_,0.);
222 std::vector<std::vector<Real> > val_pt(numSamples,tmp);
225 meas = (upperBound_[d] - lb);
227 pt = meas*pts_[k] + lb;
230 diff = (val-dist_[d]->evaluateCDF(pt));
231 for (
int j = 0; j < numSamples; j++) {
232 (val_pt[j])[d] += wt * diff * gradx[j];
233 val_wt[j] += wt * diff * gradp[j];
240 for (
int k = 0; k < numSamples; k++) {
242 gprob.setProbability(k,val_wt[k]);
259 std::vector<Real> hvxx(numSamples,0.), hvxp(numSamples,0.), hvpx(numSamples,0.);
260 std::vector<Real> gradx(numSamples,0.), gradp(numSamples,0.);
261 Real diff = 0., pt = 0., wt = 0., meas = 0., lb = 0., val = 0., sumx = 0., sump = 0.;
262 std::vector<Real> val_wt(numSamples,0.), tmp(dimension_,0.);
263 std::vector<std::vector<Real> > val_pt(numSamples,tmp);
266 meas = (upperBound_[d] - lb);
268 pt = meas*pts_[k] + lb;
270 val =
hessVecCDF(hvxx,hvxp,hvpx,gradx,gradp,sumx,sump,d,pt,prob,atom,vprob,vatom);
271 diff = (val-dist_[d]->evaluateCDF(pt));
272 for (
int j = 0; j < numSamples; j++) {
273 (val_pt[j])[d] += wt * ( (sump + sumx) * gradx[j] + diff * (hvxx[j] + hvxp[j]) );
274 val_wt[j] += wt * ( (sump + sumx) * gradp[j] + diff * hvpx[j] );
281 for (
int k = 0; k < numSamples; k++) {
283 hprob.setProbability(k,val_wt[k]);
Provides the interface to evaluate objective functions.
std::vector< Real > lowerBound_
std::vector< Teuchos::RCP< Distribution< Real > > > dist_
Real hessVecCDF(std::vector< Real > &hvxx, std::vector< Real > &hvxp, std::vector< Real > &hvpx, std::vector< Real > &gradx, std::vector< Real > &gradp, Real &sumx, Real &sump, const int dim, const Real loc, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom, const ProbabilityVector< Real > &vprob, const AtomVector< Real > &vatom) const
Provides the std::vector implementation of the ROL::Vector interface.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
const Real getProbability(const int i) const
void setAtom(const int i, const std::vector< Real > &pt)
const Teuchos::RCP< const ProbabilityVector< Real > > getProbabilityVector(void) const
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::RCP< BatchManager< Real > > bman_
Teuchos::RCP< const std::vector< Real > > getAtom(const int i) const
const Teuchos::RCP< const AtomVector< Real > > getAtomVector(void) const
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void initializeQuadrature(void)
Provides the std::vector implementation of the ROL::Vector interface.
Real gradientCDF(std::vector< Real > &gradx, std::vector< Real > &gradp, const int dim, const Real loc, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom) const
std::vector< Real > upperBound_
Provides the std::vector implementation of the ROL::Vector interface.
int getNumMyAtoms(void) const
Real valueCDF(const int dim, const Real loc, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom) const
CDFObjective(const std::vector< Teuchos::RCP< Distribution< Real > > > &dist, const Teuchos::RCP< BatchManager< Real > > &bman, const Real scale=1.e-2, const bool optProb=true, const bool optAtom=true)
int getNumMyAtoms(void) const