ROL
ROL_BPOE.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_BPOE_HPP
45#define ROL_BPOE_HPP
46
48
63namespace ROL {
64
65template<class Real>
66class BPOE : public RandVarFunctional<Real> {
67private:
69 Real order_;
70
71 std::vector<Real> hvec_;
72 ROL::Ptr<Vector<Real> > dualVec1_, dualVec2_;
73
75
76 using RandVarFunctional<Real>::val_;
77 using RandVarFunctional<Real>::gv_;
78 using RandVarFunctional<Real>::g_;
79 using RandVarFunctional<Real>::hv_;
81
82 using RandVarFunctional<Real>::point_;
84
89
90public:
91 BPOE(const Real threshold, const Real order=1)
92 : RandVarFunctional<Real>(), threshold_(threshold), order_(order), firstResetBPOE_(true) {
93 hvec_.resize(5);
94 }
95
96 BPOE(ROL::ParameterList &parlist) : RandVarFunctional<Real>(), firstResetBPOE_(true) {
97 ROL::ParameterList &list = parlist.sublist("SOL").sublist("Probability").sublist("bPOE");
98 threshold_ = list.get<Real>("Threshold");
99 order_ = list.get<Real>("Moment Order");
100 hvec_.resize(5);
101 }
102
103 void initialize(const Vector<Real> &x) {
105 if ( firstResetBPOE_ ) {
106 dualVec1_ = x.dual().clone();
107 dualVec2_ = x.dual().clone();
108 firstResetBPOE_ = false;
109 }
110 dualVec1_->zero();
111 dualVec2_->zero();
112 hvec_.assign(5,0);
113 }
114
116 const Vector<Real> &x,
117 const std::vector<Real> &xstat,
118 Real &tol) {
119 const Real zero(0), one(1);
120 Real val = computeValue(obj,x,tol);
121 Real bp = xstat[0]*(val-threshold_)+one;
122 if ( bp > zero ) {
123 val_ += weight_*((order_==one) ? bp : std::pow(bp,order_));
124 }
125 }
126
127 Real getValue(const Vector<Real> &x,
128 const std::vector<Real> &xstat,
129 SampleGenerator<Real> &sampler) {
130 const Real one(1);
131 Real bpoe(0);
132 sampler.sumAll(&val_,&bpoe,1);
133 return ((order_==one) ? bpoe : std::pow(bpoe,one/order_));
134 }
135
137 const Vector<Real> &x,
138 const std::vector<Real> &xstat,
139 Real &tol) {
140 const Real zero(0), one(1), two(2);
141 Real val = computeValue(obj,x,tol);
142 Real bp = xstat[0]*(val-threshold_)+one;
143 if ( bp > zero ) {
144 computeGradient(*dualVector_,obj,x,tol);
145 Real pvalp0 = ((order_==one) ? bp : std::pow(bp,order_));
146 Real pvalp1 = ((order_==one) ? one : ((order_==two) ? bp : std::pow(bp,order_-one)));
147 val_ += weight_ * pvalp0;
148 gv_ += weight_ * pvalp1 * (val - threshold_);
149 g_->axpy(weight_ * pvalp1, *dualVector_);
150 }
151 }
152
154 std::vector<Real> &gstat,
155 const Vector<Real> &x,
156 const std::vector<Real> &xstat,
157 SampleGenerator<Real> &sampler) {
158 const Real zero(0), one(1);
159 std::vector<Real> myvals(2), gvals(2);
160 myvals[0] = val_; myvals[1] = gv_;
161 sampler.sumAll(&myvals[0],&gvals[0],2);
162 if ( gvals[0] > zero) {
163 sampler.sumAll(*g_,g);
164 Real norm = std::pow(gvals[0],(order_-one)/order_);
165 g.scale(xstat[0]/norm);
166 gstat[0] = gvals[1]/norm;
167 }
168 else {
169 g.zero();
170 gstat[0] = zero;
171 }
172 }
173
175 const Vector<Real> &v,
176 const std::vector<Real> &vstat,
177 const Vector<Real> &x,
178 const std::vector<Real> &xstat,
179 Real &tol) {
180 const Real zero(0), one(1), two(2), three(3);
181 Real val = computeValue(obj,x,tol);
182 Real bp = xstat[0]*(val-threshold_)+one;
183 if ( bp > zero ) {
184 // Gradient only
185 Real gv = computeGradVec(*dualVector_,obj,v,x,tol);
186 Real pvalp0 = ((order_==one) ? bp : std::pow(bp,order_));
187 Real pvalp1 = ((order_==one) ? one
188 : ((order_==two) ? bp : std::pow(bp,order_-one)));
189 Real pvalp2 = ((order_==one) ? zero
190 : ((order_==two) ? one
191 : ((order_==three) ? bp : std::pow(bp,order_-two))));
192 hvec_[0] += weight_ * pvalp0;
193 hvec_[1] += weight_ * pvalp1 * (val-threshold_);
194 hvec_[2] += weight_ * pvalp2 * (val-threshold_) * (val-threshold_);
195 hvec_[3] += weight_ * pvalp1 * gv;
196 hvec_[4] += weight_ * pvalp2 * (val-threshold_) * gv;
197 g_->axpy(weight_ * pvalp1, *dualVector_);
198 dualVec1_->axpy(weight_ * pvalp2 * (val-threshold_), *dualVector_);
199 dualVec2_->axpy(weight_ * pvalp2 * gv, *dualVector_);
200 // Hessian only
201 computeHessVec(*dualVector_,obj,v,x,tol);
202 hv_->axpy(weight_ * pvalp1, *dualVector_);
203 }
204 }
205
207 std::vector<Real> &hvstat,
208 const Vector<Real> &v,
209 const std::vector<Real> &vstat,
210 const Vector<Real> &x,
211 const std::vector<Real> &xstat,
212 SampleGenerator<Real> &sampler) {
213 const Real zero(0), one(1), two(2);
214 std::vector<Real> gvals(5);
215 sampler.sumAll(&hvec_[0],&gvals[0],5);
216
217 if ( gvals[0] > zero ) {
218 Real norm0 = ((order_==one) ? one
219 : ((order_==two) ? std::sqrt(gvals[0])
220 : std::pow(gvals[0],(order_-one)/order_)));
221 Real norm1 = ((order_==one) ? gvals[0]
222 : std::pow(gvals[0],(two*order_-one)/order_));
223 hvstat[0] = (order_-one)*((gvals[2]/norm0 - gvals[1]*gvals[1]/norm1)*vstat[0]
224 +xstat[0]*(gvals[4]/norm0 - gvals[3]*gvals[1]/norm1))
225 +(gvals[3]/norm0);
226
227 sampler.sumAll(*hv_,hv);
228 hv.scale(xstat[0]/norm0);
229
230 sampler.sumAll(*g_,*hv_);
231 Real coeff = -(order_-one)*xstat[0]*(xstat[0]*gvals[3]+vstat[0]*gvals[1])/norm1+vstat[0]/norm0;
232 hv.axpy(coeff,*hv_);
233
234 sampler.sumAll(*dualVec1_,*hv_);
235 hv.axpy((order_-one)*vstat[0]*xstat[0]/norm0,*hv_);
236
237 sampler.sumAll(*dualVec2_,*hv_);
238 hv.axpy((order_-one)*xstat[0]*xstat[0]/norm0,*hv_);
239 }
240 }
241};
242
243}
244
245#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the implementation of the buffered probability of exceedance.
Definition ROL_BPOE.hpp:66
ROL::Ptr< Vector< Real > > dualVec1_
Definition ROL_BPOE.hpp:72
std::vector< Real > hvec_
Definition ROL_BPOE.hpp:71
BPOE(ROL::ParameterList &parlist)
Definition ROL_BPOE.hpp:96
Real order_
Definition ROL_BPOE.hpp:69
void updateHessVec(Objective< Real > &obj, const Vector< Real > &v, const std::vector< Real > &vstat, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for Hessian-time-a-vector computation.
Definition ROL_BPOE.hpp:174
bool firstResetBPOE_
Definition ROL_BPOE.hpp:74
BPOE(const Real threshold, const Real order=1)
Definition ROL_BPOE.hpp:91
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
Definition ROL_BPOE.hpp:115
void getHessVec(Vector< Real > &hv, std::vector< Real > &hvstat, const Vector< Real > &v, const std::vector< Real > &vstat, const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.
Definition ROL_BPOE.hpp:206
Real threshold_
Definition ROL_BPOE.hpp:68
ROL::Ptr< Vector< Real > > dualVec2_
Definition ROL_BPOE.hpp:72
void getGradient(Vector< Real > &g, std::vector< Real > &gstat, const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
Definition ROL_BPOE.hpp:153
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
Definition ROL_BPOE.hpp:127
void initialize(const Vector< Real > &x)
Initialize temporary variables.
Definition ROL_BPOE.hpp:103
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
Definition ROL_BPOE.hpp:136
Provides the interface to evaluate objective functions.
Provides the interface to implement any functional that maps a random variable to a (extended) real n...
Real computeValue(Objective< Real > &obj, const Vector< Real > &x, Real &tol)
virtual void initialize(const Vector< Real > &x)
Initialize temporary variables.
void computeHessVec(Vector< Real > &hv, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
void computeGradient(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > dualVector_
Real computeGradVec(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
void sumAll(Real *input, Real *output, int dim) const
Defines the linear algebra or vector space interface.
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .