ROL
ROL_CDFObjective.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_CDFOBJECTIVE_H
45#define ROL_CDFOBJECTIVE_H
46
47#include "ROL_Objective.hpp"
48#include "ROL_BatchManager.hpp"
49#include "ROL_Vector.hpp"
50#include "ROL_Distribution.hpp"
51#include "ROL_Ptr.hpp"
52#include <math.h>
53
54namespace ROL {
55
56template <class Real>
57class CDFObjective : public Objective<Real> {
58private:
59 // Batch manager for parallel computation
60 ROL::Ptr<BatchManager<Real> > bman_;
61
62 // Distribution information
63 std::vector<ROL::Ptr<Distribution<Real> > > dist_;
64 std::vector<Real> lowerBound_;
65 std::vector<Real> upperBound_;
67
68 const Real scale_;
69 const Real sqrt2_;
70 const Real sqrtpi_;
71
72 const bool optProb_;
73 const bool optAtom_;
74
75 std::vector<Real> pts_;
76 std::vector<Real> wts_;
77
78 // Number of quadrature points
80
82 numPoints_ = 20;
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;
105 for (int i = 0; i < numPoints_; i++) {
106 wts_[i] *= 0.5;
107 pts_[i] += 1.; pts_[i] *= 0.5;
108 }
109 }
110
111 Real valueCDF(const int dim, const Real loc,
112 const ProbabilityVector<Real> &prob,
113 const AtomVector<Real> &atom) const {
114 const int numSamples = prob.getNumMyAtoms();
115 Real val = 0, hs = 0, xpt = 0, xwt = 0, sum = 0, half(0.5), one(1);
116 for (int k = 0; k < numSamples; k++) {
117 xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
118 hs = half * (one + erf((loc-xpt)/(sqrt2_*scale_)));
119 val += xwt * hs;
120 }
121 bman_->sumAll(&val,&sum,1);
122 return sum;
123 }
124
125 Real gradientCDF(std::vector<Real> &gradx, std::vector<Real> &gradp,
126 const int dim, const Real loc,
127 const ProbabilityVector<Real> &prob,
128 const AtomVector<Real> &atom) const {
129 const int numSamples = prob.getNumMyAtoms();
130 gradx.resize(numSamples,0); gradp.resize(numSamples,0);
131 Real val = 0, hs = 0, xpt = 0, xwt = 0, sum = 0, half(0.5), one(1);
132 for (int k = 0; k < numSamples; k++) {
133 xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
134 hs = half * (one + erf((loc-xpt)/(sqrt2_*scale_)));
135 val += xwt * hs;
136 gradx[k] = -(xwt/(sqrt2_*sqrtpi_*scale_))
137 * std::exp(-std::pow((loc-xpt)/(sqrt2_*scale_),2));
138 gradp[k] = hs;
139 }
140 bman_->sumAll(&val,&sum,1);
141 return sum;
142 }
143
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,
148 const ProbabilityVector<Real> &prob,
149 const AtomVector<Real> &atom,
150 const ProbabilityVector<Real> &vprob,
151 const AtomVector<Real> &vatom) const {
152 const int numSamples = prob.getNumMyAtoms();
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), out(3,0);
157 Real val = 0, hs = 0, dval = 0, scale3 = std::pow(scale_,3);
158 Real xpt = 0, xwt = 0, vpt = 0, vwt = 0, half(0.5), one(1);
159 for (int k = 0; k < numSamples; k++) {
160 xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
161 vpt = (*vatom.getAtom(k))[dim]; vwt = vprob.getProbability(k);
162 hs = half * (one + erf((loc-xpt)/(sqrt2_*scale_)));
163 psum[0] += xwt * hs;
164 dval = std::exp(-std::pow((loc-xpt)/(sqrt2_*scale_),2));
165 gradx[k] = -(xwt/(sqrt2_*sqrtpi_*scale_)) * dval;
166 gradp[k] = hs;
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];
172 }
173 bman_->sumAll(&psum[0],&out[0],3);
174 val = out[0]; sumx = out[1]; sump = out[2];
175 return val;
176 }
177
178public:
179 CDFObjective(const std::vector<ROL::Ptr<Distribution<Real> > > &dist,
180 const ROL::Ptr<BatchManager<Real> > &bman,
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(ROL::ScalarTraits<Real>::pi())),
185 optProb_(optProb), optAtom_(optAtom) {
186 lowerBound_.resize(dimension_,0);
187 upperBound_.resize(dimension_,0);
188 for ( int i = 0; i < dimension_; i++ ) {
189 lowerBound_[i] = dist[i]->lowerBound();
190 upperBound_[i] = dist[i]->upperBound();
191 }
193 }
194
195 Real value( const Vector<Real> &x, Real &tol ) {
196 const SROMVector<Real> &ex = dynamic_cast<const SROMVector<Real>&>(x);
197 const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
198 const AtomVector<Real> &atom = *(ex.getAtomVector());
199 Real val(0), diff(0), pt(0), wt(0), meas(0), lb(0), one(1);
200 for (int d = 0; d < dimension_; d++) {
201 lb = lowerBound_[d];
202 meas = (upperBound_[d] - lb);
203 meas = ((meas > ROL_EPSILON<Real>()) ? meas : one);
204 for (int k = 0; k < numPoints_; k++) {
205 pt = meas*pts_[k] + lb;
206 wt = wts_[k]/meas;
207 diff = (valueCDF(d,pt,prob,atom)-dist_[d]->evaluateCDF(pt));
208 val += wt*std::pow(diff,2);
209 }
210 }
211 return 0.5*val;
212 }
213
214 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
215 g.zero();
216 const SROMVector<Real> &ex = dynamic_cast<const SROMVector<Real>&>(x);
217 const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
218 const AtomVector<Real> &atom = *(ex.getAtomVector());
219 const int numSamples = prob.getNumMyAtoms();
220 std::vector<Real> gradx(numSamples,0.), gradp(numSamples,0);
221 Real diff(0), pt(0), wt(0), meas(0), lb(0), val(0), one(1);
222 std::vector<Real> val_wt(numSamples,0), tmp(dimension_,0);
223 std::vector<std::vector<Real> > val_pt(numSamples,tmp);
224 for (int d = 0; d < dimension_; d++) {
225 lb = lowerBound_[d];
226 meas = (upperBound_[d] - lb);
227 meas = ((meas > ROL_EPSILON<Real>()) ? meas : one);
228 for (int k = 0; k < numPoints_; k++) {
229 pt = meas*pts_[k] + lb;
230 wt = wts_[k]/meas;
231 val = gradientCDF(gradx,gradp,d,pt,prob,atom);
232 diff = (val-dist_[d]->evaluateCDF(pt));
233 for (int j = 0; j < numSamples; j++) {
234 (val_pt[j])[d] += wt * diff * gradx[j];
235 val_wt[j] += wt * diff * gradp[j];
236 }
237 }
238 }
239 SROMVector<Real> &eg = dynamic_cast<SROMVector<Real>&>(g);
241 AtomVector<Real> &gatom = *(eg.getAtomVector());
242 for (int k = 0; k < numSamples; k++) {
243 if ( optProb_ ) {
244 gprob.setProbability(k,val_wt[k]);
245 }
246 if ( optAtom_ ) {
247 gatom.setAtom(k,val_pt[k]);
248 }
249 }
250 }
251
252 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
253 hv.zero();
254 const SROMVector<Real> &ev = dynamic_cast<const SROMVector<Real>&>(v);
255 const ProbabilityVector<Real> &vprob = *(ev.getProbabilityVector());
256 const AtomVector<Real> &vatom = *(ev.getAtomVector());
257 const SROMVector<Real> &ex = dynamic_cast<const SROMVector<Real>&>(x);
258 const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
259 const AtomVector<Real> &atom = *(ex.getAtomVector());
260 const int numSamples = prob.getNumMyAtoms();
261 std::vector<Real> hvxx(numSamples,0), hvxp(numSamples,0), hvpx(numSamples,0);
262 std::vector<Real> gradx(numSamples,0), gradp(numSamples,0);
263 Real diff(0), pt(0), wt(0), meas(0), lb(0), val(0), sumx(0), sump(0), one(1);
264 std::vector<Real> val_wt(numSamples,0), tmp(dimension_,0);
265 std::vector<std::vector<Real> > val_pt(numSamples,tmp);
266 for (int d = 0; d < dimension_; d++) {
267 lb = lowerBound_[d];
268 meas = (upperBound_[d] - lb);
269 meas = ((meas > ROL_EPSILON<Real>()) ? meas : one);
270 for (int k = 0; k < numPoints_; k++) {
271 pt = meas*pts_[k] + lb;
272 wt = wts_[k]/meas;
273 val = hessVecCDF(hvxx,hvxp,hvpx,gradx,gradp,sumx,sump,d,pt,prob,atom,vprob,vatom);
274 diff = (val-dist_[d]->evaluateCDF(pt));
275 for (int j = 0; j < numSamples; j++) {
276 (val_pt[j])[d] += wt * ( (sump + sumx) * gradx[j] + diff * (hvxx[j] + hvxp[j]) );
277 val_wt[j] += wt * ( (sump + sumx) * gradp[j] + diff * hvpx[j] );
278 }
279 }
280 }
281 SROMVector<Real> &ehv = dynamic_cast<SROMVector<Real>&>(hv);
283 AtomVector<Real> &hatom = *(ehv.getAtomVector());
284 for (int k = 0; k < numSamples; k++) {
285 if ( optProb_ ) {
286 hprob.setProbability(k,val_wt[k]);
287 }
288 if ( optAtom_ ) {
289 hatom.setAtom(k,val_pt[k]);
290 }
291 }
292 }
293}; // class LinearCombinationObjective
294
295} // namespace ROL
296
297#endif
Provides the std::vector implementation of the ROL::Vector interface.
ROL::Ptr< const std::vector< Real > > getAtom(const int i) const
void setAtom(const int i, const std::vector< Real > &pt)
void initializeQuadrature(void)
std::vector< Real > wts_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
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)
std::vector< Real > upperBound_
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
ROL::Ptr< BatchManager< Real > > bman_
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
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
std::vector< ROL::Ptr< Distribution< Real > > > dist_
std::vector< Real > lowerBound_
Real valueCDF(const int dim, const Real loc, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom) const
std::vector< Real > pts_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Provides the interface to evaluate objective functions.
Provides the std::vector implementation of the ROL::Vector interface.
const Real getProbability(const int i) const
Provides the std::vector implementation of the ROL::Vector interface.
const ROL::Ptr< const ProbabilityVector< Real > > getProbabilityVector(void) const
const ROL::Ptr< const AtomVector< Real > > getAtomVector(void) const
Defines the linear algebra or vector space interface.
virtual void zero()
Set to zero vector.
constexpr auto dim