ROL
ROL_MeanDeviationFromTarget.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_MEANDEVIATIONFROMTARGET_HPP
45#define ROL_MEANDEVIATIONFROMTARGET_HPP
46
48#include "ROL_ParameterList.hpp"
50#include "ROL_PlusFunction.hpp"
51#include "ROL_AbsoluteValue.hpp"
52
73namespace ROL {
74
75template<class Real>
77 typedef typename std::vector<Real>::size_type uint;
78private:
79 Ptr<PositiveFunction<Real> > positiveFunction_;
80 std::vector<Real> target_;
81 std::vector<Real> order_;
82 std::vector<Real> coeff_;
84
85 std::vector<Real> pval_;
86 std::vector<Real> pgv_;
87
88 std::vector<Ptr<Vector<Real> > > pg0_;
89 std::vector<Ptr<Vector<Real> > > pg_;
90 std::vector<Ptr<Vector<Real> > > phv_;
91
93
94 using RandVarFunctional<Real>::val_;
95 using RandVarFunctional<Real>::gv_;
96 using RandVarFunctional<Real>::g_;
97 using RandVarFunctional<Real>::hv_;
99
100 using RandVarFunctional<Real>::point_;
101 using RandVarFunctional<Real>::weight_;
102
107
108 void initializeMDT(void) {
109 // Initialize additional storage
110 pg_.clear(); pg0_.clear(); phv_.clear(); pval_.clear(); pgv_.clear();
111 pg_.resize(NumMoments_);
112 pg0_.resize(NumMoments_);
113 phv_.resize(NumMoments_);
114 pval_.resize(NumMoments_);
115 pgv_.resize(NumMoments_);
116 }
117
118 void checkInputs(void) {
119 int oSize = order_.size(), cSize = coeff_.size(), tSize = target_.size();
120 ROL_TEST_FOR_EXCEPTION((oSize!=cSize),std::invalid_argument,
121 ">>> ERROR (ROL::MeanDeviationFromTarget): Order and coefficient arrays have different sizes!");
122 ROL_TEST_FOR_EXCEPTION((oSize!=tSize),std::invalid_argument,
123 ">>> ERROR (ROL::MeanDeviationFromTarget): Order and target arrays have different sizes!");
124 Real zero(0), two(2);
125 for (int i = 0; i < oSize; i++) {
126 ROL_TEST_FOR_EXCEPTION((order_[i] < two), std::invalid_argument,
127 ">>> ERROR (ROL::MeanDeviationFromTarget): Element of order array out of range!");
128 ROL_TEST_FOR_EXCEPTION((coeff_[i] < zero), std::invalid_argument,
129 ">>> ERROR (ROL::MeanDeviationFromTarget): Element of coefficient array out of range!");
130 }
131 ROL_TEST_FOR_EXCEPTION(positiveFunction_ == nullPtr, std::invalid_argument,
132 ">>> ERROR (ROL::MeanDeviationFromTarget): PositiveFunction pointer is null!");
134 }
135
136public:
147 MeanDeviationFromTarget( const Real target, const Real order, const Real coeff,
148 const Ptr<PositiveFunction<Real> > &pf )
150 order_.clear(); order_.push_back(order);
151 coeff_.clear(); coeff_.push_back(coeff);
152 target_.clear(); target_.push_back(target);
153 NumMoments_ = order_.size();
154 checkInputs();
155 }
156
167 MeanDeviationFromTarget( const std::vector<Real> &target,
168 const std::vector<Real> &order,
169 const std::vector<Real> &coeff,
170 const Ptr<PositiveFunction<Real> > &pf )
172 target_.clear(); order_.clear(); coeff_.clear();
173 for ( uint i = 0; i < target.size(); i++ ) {
174 target_.push_back(target[i]);
175 }
176 for ( uint i = 0; i < order.size(); i++ ) {
177 order_.push_back(order[i]);
178 }
179 for ( uint i = 0; i < coeff.size(); i++ ) {
180 coeff_.push_back(coeff[i]);
181 }
182 NumMoments_ = order_.size();
183 checkInputs();
184 }
185
198 MeanDeviationFromTarget( ROL::ParameterList &parlist )
199 : RandVarFunctional<Real>(), firstResetMDT_(true) {
200 ROL::ParameterList &list
201 = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation From Target");
202
203 // Get data from parameter list
204 target_ = ROL::getArrayFromStringParameter<double>(list,"Targets");
205
206 order_ = ROL::getArrayFromStringParameter<double>(list,"Orders");
207
208 coeff_ = ROL::getArrayFromStringParameter<double>(list,"Coefficients");
209
210 // Build (approximate) positive function
211 std::string type = list.get<std::string>("Deviation Type");
212 if ( type == "Upper" ) {
213 positiveFunction_ = makePtr<PlusFunction<Real>>(list);
214 }
215 else if ( type == "Absolute" ) {
216 positiveFunction_ = makePtr<AbsoluteValue<Real>>(list);
217 }
218 else {
219 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
220 ">>> (ROL::MeanDeviation): Deviation type is not recoginized!");
221 }
222 // Check inputs
223 NumMoments_ = order_.size();
224 checkInputs();
225 }
226
227 void initialize(const Vector<Real> &x) {
229 if (firstResetMDT_) {
230 for ( uint p = 0; p < NumMoments_; p++ ) {
231 pg0_[p] = x.dual().clone();
232 pg_[p] = x.dual().clone();
233 phv_[p] = x.dual().clone();
234 }
235 firstResetMDT_ = false;
236 }
237 Real zero(0);
238 for ( uint p = 0; p < NumMoments_; p++ ) {
239 pg0_[p]->zero(); pg_[p]->zero(); phv_[p]->zero();
240 pval_[p] = zero; pgv_[p] = zero;
241 }
242 }
243
245 const Vector<Real> &x,
246 const std::vector<Real> &xstat,
247 Real &tol) {
248 Real diff(0), pf0(0);
249 Real val = computeValue(obj,x,tol);
250 val_ += weight_ * val;
251 for ( uint p = 0; p < NumMoments_; p++ ) {
252 diff = val-target_[p];
253 pf0 = positiveFunction_->evaluate(diff,0);
254 pval_[p] += weight_ * std::pow(pf0,order_[p]);
255 }
256 }
257
258 Real getValue(const Vector<Real> &x,
259 const std::vector<Real> &xstat,
260 SampleGenerator<Real> &sampler) {
261 const Real one(1);
262 Real dev(0);
263 sampler.sumAll(&val_,&dev,1);
264 std::vector<Real> pval_sum(NumMoments_);
265 sampler.sumAll(&pval_[0],&pval_sum[0],NumMoments_);
266 for ( uint p = 0; p < NumMoments_; p++ ) {
267 dev += coeff_[p] * std::pow(pval_sum[p],one/order_[p]);
268 }
269 return dev;
270 }
271
273 const Vector<Real> &x,
274 const std::vector<Real> &xstat,
275 Real &tol) {
276 Real diff(0), pf0(0), pf1(0), c(0), one(1);
277 Real val = computeValue(obj,x,tol);
278 computeGradient(*dualVector_,obj,x,tol);
279 for ( uint p = 0; p < NumMoments_; p++ ) {
280 diff = val-target_[p];
281 pf0 = positiveFunction_->evaluate(diff,0);
282 pf1 = positiveFunction_->evaluate(diff,1);
283 c = std::pow(pf0,order_[p]-one) * pf1;
284 (pg_[p])->axpy(weight_ * c,*dualVector_);
285 pval_[p] += weight_ * std::pow(pf0,order_[p]);
286 }
287 g_->axpy(weight_,*dualVector_);
288 }
289
291 std::vector<Real> &gstat,
292 const Vector<Real> &x,
293 const std::vector<Real> &xstat,
294 SampleGenerator<Real> &sampler) {
295 const Real zero(0), one(1);
296 sampler.sumAll(*g_,g);
297 std::vector<Real> pval_sum(NumMoments_);
298 sampler.sumAll(&pval_[0],&pval_sum[0],NumMoments_);
299 for ( uint p = 0; p < NumMoments_; p++ ) {
300 if ( pval_sum[p] > zero ) {
301 dualVector_->zero();
302 sampler.sumAll(*(pg_[p]),*dualVector_);
303 g.axpy(coeff_[p]/std::pow(pval_sum[p],one-one/order_[p]),*dualVector_);
304 }
305 }
306 }
307
309 const Vector<Real> &v,
310 const std::vector<Real> &vstat,
311 const Vector<Real> &x,
312 const std::vector<Real> &xstat,
313 Real &tol) {
314 const Real one(1), two(2);
315 Real diff(0), pf0(0), pf1(0), pf2(0), p0(0), p1(0), p2(0), c(0);
316 Real val = computeValue(obj,x,tol);
317 Real gv = computeGradVec(*g_,obj,v,x,tol);
318 computeHessVec(*dualVector_,obj,v,x,tol);
319 for ( uint p = 0; p < NumMoments_; p++ ) {
320 diff = val - target_[p];
321 pf0 = positiveFunction_->evaluate(diff,0);
322 pf1 = positiveFunction_->evaluate(diff,1);
323 pf2 = positiveFunction_->evaluate(diff,2);
324 p0 = std::pow(pf0,order_[p]);
325 p1 = std::pow(pf0,order_[p]-one);
326 p2 = std::pow(pf0,order_[p]-two);
327 c = -(order_[p]-one)*p1*pf1;
328 pg0_[p]->axpy(weight_*c,*g_);
329 c = gv*((order_[p]-one)*p2*pf1*pf1 + p1*pf2);
330 pg_[p]->axpy(weight_*c,*g_);
331 c = p1*pf1;
332 phv_[p]->axpy(weight_*c,*dualVector_);
333 pval_[p] += weight_*p0;
334 pgv_[p] += weight_*p1*pf1*gv;
335 }
336 hv_->axpy(weight_,*dualVector_);
337 }
338
340 std::vector<Real> &hvstat,
341 const Vector<Real> &v,
342 const std::vector<Real> &vstat,
343 const Vector<Real> &x,
344 const std::vector<Real> &xstat,
345 SampleGenerator<Real> &sampler) {
346 const Real zero(0), one(1), two(2);
347 sampler.sumAll(*hv_,hv);
348 std::vector<Real> pval_sum(NumMoments_);
349 sampler.sumAll(&(pval_)[0],&pval_sum[0],NumMoments_);
350 std::vector<Real> pgv_sum(NumMoments_);
351 sampler.sumAll(&(pgv_)[0],&pgv_sum[0],NumMoments_);
352 Real c(0);
353 for ( uint p = 0; p < NumMoments_; p++ ) {
354 if ( pval_sum[p] > zero ) {
355 dualVector_->zero();
356 sampler.sumAll(*(pg0_[p]),*dualVector_);
357 c = coeff_[p]*(pgv_sum[p]/std::pow(pval_sum[p],two-one/order_[p]));
358 hv.axpy(c,*dualVector_);
359
360 dualVector_->zero();
361 sampler.sumAll(*(pg_[p]),*dualVector_);
362 c = coeff_[p]/std::pow(pval_sum[p],one-one/order_[p]);
363 hv.axpy(c,*dualVector_);
364
365 dualVector_->zero();
366 sampler.sumAll(*(phv_[p]),*dualVector_);
367 hv.axpy(c,*dualVector_);
368 }
369 }
370 }
371};
372
373}
374
375#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides an interface for the mean plus a sum of arbitrary order deviations from targets.
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
MeanDeviationFromTarget(const Real target, const Real order, const Real coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
MeanDeviationFromTarget(ROL::ParameterList &parlist)
Constructor.
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
std::vector< Ptr< Vector< Real > > > pg_
std::vector< Ptr< Vector< Real > > > pg0_
std::vector< Ptr< Vector< Real > > > phv_
Ptr< PositiveFunction< Real > > positiveFunction_
MeanDeviationFromTarget(const std::vector< Real > &target, const std::vector< Real > &order, const std::vector< Real > &coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
void initialize(const Vector< Real > &x)
Initialize temporary variables.
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.
std::vector< Real >::size_type uint
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.
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.
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 const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
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 .