ROL
ROL_CDFObjective.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_CDFOBJECTIVE_H
11 #define ROL_CDFOBJECTIVE_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_BatchManager.hpp"
15 #include "ROL_Vector.hpp"
16 #include "ROL_Distribution.hpp"
17 #include "ROL_Ptr.hpp"
18 #include <math.h>
19 
20 namespace ROL {
21 
22 template <class Real>
23 class CDFObjective : public Objective<Real> {
24 private:
25  // Batch manager for parallel computation
26  ROL::Ptr<BatchManager<Real> > bman_;
27 
28  // Distribution information
29  std::vector<ROL::Ptr<Distribution<Real> > > dist_;
30  std::vector<Real> lowerBound_;
31  std::vector<Real> upperBound_;
33 
34  const Real scale_;
35  const Real sqrt2_;
36  const Real sqrtpi_;
37 
38  const bool optProb_;
39  const bool optAtom_;
40 
41  std::vector<Real> pts_;
42  std::vector<Real> wts_;
43 
44  // Number of quadrature points
46 
47  void initializeQuadrature(void) {
48  numPoints_ = 20;
49  pts_.clear(); pts_.resize(numPoints_,0.);
50  wts_.clear(); wts_.resize(numPoints_,0.);
51  wts_[0] = 0.1527533871307258; pts_[0] = -0.0765265211334973;
52  wts_[1] = 0.1527533871307258; pts_[1] = 0.0765265211334973;
53  wts_[2] = 0.1491729864726037; pts_[2] = -0.2277858511416451;
54  wts_[3] = 0.1491729864726037; pts_[3] = 0.2277858511416451;
55  wts_[4] = 0.1420961093183820; pts_[4] = -0.3737060887154195;
56  wts_[5] = 0.1420961093183820; pts_[5] = 0.3737060887154195;
57  wts_[6] = 0.1316886384491766; pts_[6] = -0.5108670019508271;
58  wts_[7] = 0.1316886384491766; pts_[7] = 0.5108670019508271;
59  wts_[8] = 0.1181945319615184; pts_[8] = -0.6360536807265150;
60  wts_[9] = 0.1181945319615184; pts_[9] = 0.6360536807265150;
61  wts_[10] = 0.1019301198172404; pts_[10] = -0.7463319064601508;
62  wts_[11] = 0.1019301198172404; pts_[11] = 0.7463319064601508;
63  wts_[12] = 0.0832767415767048; pts_[12] = -0.8391169718222188;
64  wts_[13] = 0.0832767415767048; pts_[13] = 0.8391169718222188;
65  wts_[14] = 0.0626720483341091; pts_[14] = -0.9122344282513259;
66  wts_[15] = 0.0626720483341091; pts_[15] = 0.9122344282513259;
67  wts_[16] = 0.0406014298003869; pts_[16] = -0.9639719272779138;
68  wts_[17] = 0.0406014298003869; pts_[17] = 0.9639719272779138;
69  wts_[18] = 0.0176140071391521; pts_[18] = -0.9931285991850949;
70  wts_[19] = 0.0176140071391521; pts_[19] = 0.9931285991850949;
71  for (int i = 0; i < numPoints_; i++) {
72  wts_[i] *= 0.5;
73  pts_[i] += 1.; pts_[i] *= 0.5;
74  }
75  }
76 
77  Real valueCDF(const int dim, const Real loc,
78  const ProbabilityVector<Real> &prob,
79  const AtomVector<Real> &atom) const {
80  const int numSamples = prob.getNumMyAtoms();
81  Real val = 0, hs = 0, xpt = 0, xwt = 0, sum = 0, half(0.5), one(1);
82  for (int k = 0; k < numSamples; k++) {
83  xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
84  hs = half * (one + erf((loc-xpt)/(sqrt2_*scale_)));
85  val += xwt * hs;
86  }
87  bman_->sumAll(&val,&sum,1);
88  return sum;
89  }
90 
91  Real gradientCDF(std::vector<Real> &gradx, std::vector<Real> &gradp,
92  const int dim, const Real loc,
93  const ProbabilityVector<Real> &prob,
94  const AtomVector<Real> &atom) const {
95  const int numSamples = prob.getNumMyAtoms();
96  gradx.resize(numSamples,0); gradp.resize(numSamples,0);
97  Real val = 0, hs = 0, xpt = 0, xwt = 0, sum = 0, half(0.5), one(1);
98  for (int k = 0; k < numSamples; k++) {
99  xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
100  hs = half * (one + erf((loc-xpt)/(sqrt2_*scale_)));
101  val += xwt * hs;
102  gradx[k] = -(xwt/(sqrt2_*sqrtpi_*scale_))
103  * std::exp(-std::pow((loc-xpt)/(sqrt2_*scale_),2));
104  gradp[k] = hs;
105  }
106  bman_->sumAll(&val,&sum,1);
107  return sum;
108  }
109 
110  Real hessVecCDF(std::vector<Real> &hvxx, std::vector<Real> &hvxp, std::vector<Real> &hvpx,
111  std::vector<Real> &gradx, std::vector<Real> &gradp,
112  Real &sumx, Real &sump,
113  const int dim, const Real loc,
114  const ProbabilityVector<Real> &prob,
115  const AtomVector<Real> &atom,
116  const ProbabilityVector<Real> &vprob,
117  const AtomVector<Real> &vatom) const {
118  const int numSamples = prob.getNumMyAtoms();
119  hvxx.resize(numSamples,0); hvxp.resize(numSamples,0); hvpx.resize(numSamples,0);
120  gradx.resize(numSamples,0); gradp.resize(numSamples,0);
121  sumx = 0; sump = 0;
122  std::vector<Real> psum(3,0), out(3,0);
123  Real val = 0, hs = 0, dval = 0, scale3 = std::pow(scale_,3);
124  Real xpt = 0, xwt = 0, vpt = 0, vwt = 0, half(0.5), one(1);
125  for (int k = 0; k < numSamples; k++) {
126  xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
127  vpt = (*vatom.getAtom(k))[dim]; vwt = vprob.getProbability(k);
128  hs = half * (one + erf((loc-xpt)/(sqrt2_*scale_)));
129  psum[0] += xwt * hs;
130  dval = std::exp(-std::pow((loc-xpt)/(sqrt2_*scale_),2));
131  gradx[k] = -(xwt/(sqrt2_*sqrtpi_*scale_)) * dval;
132  gradp[k] = hs;
133  hvxx[k] = -(xwt/(sqrt2_*sqrtpi_*scale3)) * dval * (loc-xpt) * vpt;
134  hvxp[k] = -dval/(sqrt2_*sqrtpi_*scale_)*vwt;
135  hvpx[k] = -dval/(sqrt2_*sqrtpi_*scale_)*vpt;
136  psum[1] += vpt*gradx[k];
137  psum[2] += vwt*gradp[k];
138  }
139  bman_->sumAll(&psum[0],&out[0],3);
140  val = out[0]; sumx = out[1]; sump = out[2];
141  return val;
142  }
143 
144 public:
145  CDFObjective(const std::vector<ROL::Ptr<Distribution<Real> > > &dist,
146  const ROL::Ptr<BatchManager<Real> > &bman,
147  const Real scale = 1.e-2,
148  const bool optProb = true, const bool optAtom = true)
149  : Objective<Real>(), bman_(bman), dist_(dist), dimension_(dist.size()),
150  scale_(scale), sqrt2_(std::sqrt(2.)), sqrtpi_(std::sqrt(ROL::ScalarTraits<Real>::pi())),
151  optProb_(optProb), optAtom_(optAtom) {
152  lowerBound_.resize(dimension_,0);
153  upperBound_.resize(dimension_,0);
154  for ( int i = 0; i < dimension_; i++ ) {
155  lowerBound_[i] = dist[i]->lowerBound();
156  upperBound_[i] = dist[i]->upperBound();
157  }
159  }
160 
161  Real value( const Vector<Real> &x, Real &tol ) {
162  const SROMVector<Real> &ex = dynamic_cast<const SROMVector<Real>&>(x);
163  const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
164  const AtomVector<Real> &atom = *(ex.getAtomVector());
165  Real val(0), diff(0), pt(0), wt(0), meas(0), lb(0), one(1);
166  for (int d = 0; d < dimension_; d++) {
167  lb = lowerBound_[d];
168  meas = (upperBound_[d] - lb);
169  meas = ((meas > ROL_EPSILON<Real>()) ? meas : one);
170  for (int k = 0; k < numPoints_; k++) {
171  pt = meas*pts_[k] + lb;
172  wt = wts_[k]/meas;
173  diff = (valueCDF(d,pt,prob,atom)-dist_[d]->evaluateCDF(pt));
174  val += wt*std::pow(diff,2);
175  }
176  }
177  return 0.5*val;
178  }
179 
180  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
181  g.zero();
182  const SROMVector<Real> &ex = dynamic_cast<const SROMVector<Real>&>(x);
183  const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
184  const AtomVector<Real> &atom = *(ex.getAtomVector());
185  const int numSamples = prob.getNumMyAtoms();
186  std::vector<Real> gradx(numSamples,0.), gradp(numSamples,0);
187  Real diff(0), pt(0), wt(0), meas(0), lb(0), val(0), one(1);
188  std::vector<Real> val_wt(numSamples,0), tmp(dimension_,0);
189  std::vector<std::vector<Real> > val_pt(numSamples,tmp);
190  for (int d = 0; d < dimension_; d++) {
191  lb = lowerBound_[d];
192  meas = (upperBound_[d] - lb);
193  meas = ((meas > ROL_EPSILON<Real>()) ? meas : one);
194  for (int k = 0; k < numPoints_; k++) {
195  pt = meas*pts_[k] + lb;
196  wt = wts_[k]/meas;
197  val = gradientCDF(gradx,gradp,d,pt,prob,atom);
198  diff = (val-dist_[d]->evaluateCDF(pt));
199  for (int j = 0; j < numSamples; j++) {
200  (val_pt[j])[d] += wt * diff * gradx[j];
201  val_wt[j] += wt * diff * gradp[j];
202  }
203  }
204  }
205  SROMVector<Real> &eg = dynamic_cast<SROMVector<Real>&>(g);
207  AtomVector<Real> &gatom = *(eg.getAtomVector());
208  for (int k = 0; k < numSamples; k++) {
209  if ( optProb_ ) {
210  gprob.setProbability(k,val_wt[k]);
211  }
212  if ( optAtom_ ) {
213  gatom.setAtom(k,val_pt[k]);
214  }
215  }
216  }
217 
218  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
219  hv.zero();
220  const SROMVector<Real> &ev = dynamic_cast<const SROMVector<Real>&>(v);
221  const ProbabilityVector<Real> &vprob = *(ev.getProbabilityVector());
222  const AtomVector<Real> &vatom = *(ev.getAtomVector());
223  const SROMVector<Real> &ex = dynamic_cast<const SROMVector<Real>&>(x);
224  const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
225  const AtomVector<Real> &atom = *(ex.getAtomVector());
226  const int numSamples = prob.getNumMyAtoms();
227  std::vector<Real> hvxx(numSamples,0), hvxp(numSamples,0), hvpx(numSamples,0);
228  std::vector<Real> gradx(numSamples,0), gradp(numSamples,0);
229  Real diff(0), pt(0), wt(0), meas(0), lb(0), val(0), sumx(0), sump(0), one(1);
230  std::vector<Real> val_wt(numSamples,0), tmp(dimension_,0);
231  std::vector<std::vector<Real> > val_pt(numSamples,tmp);
232  for (int d = 0; d < dimension_; d++) {
233  lb = lowerBound_[d];
234  meas = (upperBound_[d] - lb);
235  meas = ((meas > ROL_EPSILON<Real>()) ? meas : one);
236  for (int k = 0; k < numPoints_; k++) {
237  pt = meas*pts_[k] + lb;
238  wt = wts_[k]/meas;
239  val = hessVecCDF(hvxx,hvxp,hvpx,gradx,gradp,sumx,sump,d,pt,prob,atom,vprob,vatom);
240  diff = (val-dist_[d]->evaluateCDF(pt));
241  for (int j = 0; j < numSamples; j++) {
242  (val_pt[j])[d] += wt * ( (sump + sumx) * gradx[j] + diff * (hvxx[j] + hvxp[j]) );
243  val_wt[j] += wt * ( (sump + sumx) * gradp[j] + diff * hvpx[j] );
244  }
245  }
246  }
247  SROMVector<Real> &ehv = dynamic_cast<SROMVector<Real>&>(hv);
249  AtomVector<Real> &hatom = *(ehv.getAtomVector());
250  for (int k = 0; k < numSamples; k++) {
251  if ( optProb_ ) {
252  hprob.setProbability(k,val_wt[k]);
253  }
254  if ( optAtom_ ) {
255  hatom.setAtom(k,val_pt[k]);
256  }
257  }
258  }
259 }; // class LinearCombinationObjective
260 
261 } // namespace ROL
262 
263 #endif
Provides the interface to evaluate objective functions.
std::vector< Real > lowerBound_
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
const ROL::Ptr< const ProbabilityVector< Real > > getProbabilityVector(void) 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
ROL::Ptr< const std::vector< Real > > getAtom(const int i) const
void setAtom(const int i, const std::vector< Real > &pt)
CDFObjective(const std::vector< ROL::Ptr< Distribution< Real > > > &dist, const ROL::Ptr< BatchManager< Real > > &bman, const Real scale=1.e-2, const bool optProb=true, const bool optAtom=true)
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
std::vector< Real > wts_
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.
std::vector< Real > pts_
std::vector< ROL::Ptr< Distribution< Real > > > dist_
const ROL::Ptr< const AtomVector< Real > > getAtomVector(void) const
ROL::Ptr< BatchManager< Real > > bman_
Real valueCDF(const int dim, const Real loc, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom) const
constexpr auto dim