ROL
ROL_MeanDeviation.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_MEANDEVIATION_HPP
45#define ROL_MEANDEVIATION_HPP
46
48#include "ROL_ParameterList.hpp"
50#include "ROL_PlusFunction.hpp"
51#include "ROL_AbsoluteValue.hpp"
52
77namespace ROL {
78
79template<class Real>
80class MeanDeviation : public RandVarFunctional<Real> {
81 typedef typename std::vector<Real>::size_type uint;
82private:
83 Ptr<PositiveFunction<Real> > positiveFunction_;
84 std::vector<Real> order_;
85 std::vector<Real> coeff_;
87
88 std::vector<Real> dev0_;
89 std::vector<Real> dev1_;
90 std::vector<Real> dev2_;
91 std::vector<Real> dev3_;
92 std::vector<Real> des0_;
93 std::vector<Real> des1_;
94 std::vector<Real> des2_;
95 std::vector<Real> des3_;
96 std::vector<Real> devp_;
97 std::vector<Real> gvp1_;
98 std::vector<Real> gvp2_;
99 std::vector<Real> gvp3_;
100 std::vector<Real> gvs1_;
101 std::vector<Real> gvs2_;
102 std::vector<Real> gvs3_;
103
104 Ptr<ScalarController<Real>> values_;
105 Ptr<ScalarController<Real>> gradvecs_;
106 Ptr<VectorController<Real>> gradients_;
107 Ptr<VectorController<Real>> hessvecs_;
108
109 using RandVarFunctional<Real>::val_;
110 using RandVarFunctional<Real>::gv_;
111 using RandVarFunctional<Real>::g_;
112 using RandVarFunctional<Real>::hv_;
114
115 using RandVarFunctional<Real>::point_;
116 using RandVarFunctional<Real>::weight_;
117
122
123 void initializeStorage(void) {
124 NumMoments_ = order_.size();
125
126 dev0_.clear(); dev1_.clear(); dev2_.clear(); dev3_.clear();
127 des0_.clear(); des1_.clear(); des2_.clear(); des3_.clear();
128 devp_.clear();
129 gvp1_.clear(); gvp2_.clear(); gvp3_.clear();
130 gvs1_.clear(); gvs2_.clear(); gvs3_.clear();
131
132 dev0_.resize(NumMoments_); dev1_.resize(NumMoments_);
133 dev2_.resize(NumMoments_); dev3_.resize(NumMoments_);
134 des0_.resize(NumMoments_); des1_.resize(NumMoments_);
135 des2_.resize(NumMoments_); des3_.resize(NumMoments_);
136 devp_.resize(NumMoments_);
137 gvp1_.resize(NumMoments_); gvp2_.resize(NumMoments_);
138 gvp3_.resize(NumMoments_);
139 gvs1_.resize(NumMoments_); gvs2_.resize(NumMoments_);
140 gvs3_.resize(NumMoments_);
141
142 values_ = makePtr<ScalarController<Real>>();
143 gradvecs_ = makePtr<ScalarController<Real>>();
144 gradients_ = makePtr<VectorController<Real>>();
145 hessvecs_ = makePtr<VectorController<Real>>();
146
149 }
150
151 void clear(void) {
152 const Real zero(0);
153 dev0_.assign(NumMoments_,zero); dev1_.assign(NumMoments_,zero);
154 dev2_.assign(NumMoments_,zero); dev3_.assign(NumMoments_,zero);
155 des0_.assign(NumMoments_,zero); des1_.assign(NumMoments_,zero);
156 des2_.assign(NumMoments_,zero); des3_.assign(NumMoments_,zero);
157 devp_.assign(NumMoments_,zero);
158 gvp1_.assign(NumMoments_,zero); gvp2_.assign(NumMoments_,zero);
159 gvp3_.assign(NumMoments_,zero);
160 gvs1_.assign(NumMoments_,zero); gvs2_.assign(NumMoments_,zero);
161 gvs3_.assign(NumMoments_,zero);
162
163 gradvecs_->reset();
164 hessvecs_->reset();
165 }
166
167 void checkInputs(void) {
168 int oSize = order_.size(), cSize = coeff_.size();
169 ROL_TEST_FOR_EXCEPTION((oSize!=cSize),std::invalid_argument,
170 ">>> ERROR (ROL::MeanDeviation): Order and coefficient arrays have different sizes!");
171 Real zero(0), two(2);
172 for (int i = 0; i < oSize; i++) {
173 ROL_TEST_FOR_EXCEPTION((order_[i] < two), std::invalid_argument,
174 ">>> ERROR (ROL::MeanDeviation): Element of order array out of range!");
175 ROL_TEST_FOR_EXCEPTION((coeff_[i] < zero), std::invalid_argument,
176 ">>> ERROR (ROL::MeanDeviation): Element of coefficient array out of range!");
177 }
178 ROL_TEST_FOR_EXCEPTION(positiveFunction_ == nullPtr, std::invalid_argument,
179 ">>> ERROR (ROL::MeanDeviation): PositiveFunction pointer is null!");
181 }
182
183public:
193 MeanDeviation( const Real order, const Real coeff,
194 const Ptr<PositiveFunction<Real> > &pf )
195 : RandVarFunctional<Real>(), positiveFunction_(pf) {
196 order_.clear(); order_.push_back(order);
197 coeff_.clear(); coeff_.push_back(coeff);
198 checkInputs();
199 }
200
210 MeanDeviation( const std::vector<Real> &order,
211 const std::vector<Real> &coeff,
212 const Ptr<PositiveFunction<Real> > &pf )
213 : RandVarFunctional<Real>(), positiveFunction_(pf) {
214 order_.clear(); coeff_.clear();
215 for ( uint i = 0; i < order.size(); i++ ) {
216 order_.push_back(order[i]);
217 }
218 for ( uint i = 0; i < coeff.size(); i++ ) {
219 coeff_.push_back(coeff[i]);
220 }
221 checkInputs();
222 }
223
235 MeanDeviation( ROL::ParameterList &parlist )
236 : RandVarFunctional<Real>() {
237 ROL::ParameterList &list
238 = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation");
239
240 // Get data from parameter list
241 order_ = ROL::getArrayFromStringParameter<double>(list,"Orders");
242
243 coeff_ = ROL::getArrayFromStringParameter<double>(list,"Coefficients");
244
245 // Build (approximate) positive function
246 std::string type = list.get<std::string>("Deviation Type");
247 if ( type == "Upper" ) {
248 positiveFunction_ = makePtr<PlusFunction<Real>>(list);
249 }
250 else if ( type == "Absolute" ) {
251 positiveFunction_ = makePtr<AbsoluteValue<Real>>(list);
252 }
253 else {
254 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
255 ">>> (ROL::MeanDeviation): Deviation type is not recoginized!");
256 }
257 // Check inputs
258 checkInputs();
259 }
260
261 void setStorage(const Ptr<ScalarController<Real>> &value_storage,
262 const Ptr<VectorController<Real>> &gradient_storage) {
263 values_ = value_storage;
264 gradients_ = gradient_storage;
266 }
267
268 void setHessVecStorage(const Ptr<ScalarController<Real>> &gradvec_storage,
269 const Ptr<VectorController<Real>> &hessvec_storage) {
270 gradvecs_ = gradvec_storage;
271 hessvecs_ = hessvec_storage;
273 }
274
275 void initialize(const Vector<Real> &x) {
277 clear();
278 }
279
281 const Vector<Real> &x,
282 const std::vector<Real> &xstat,
283 Real &tol) {
284 Real val = computeValue(obj,x,tol);
285 val_ += weight_ * val;
286 }
287
289 const Vector<Real> &x,
290 const std::vector<Real> &xstat,
291 Real &tol) {
292 Real val = computeValue(obj,x,tol);
293 val_ += weight_ * val;
294 computeGradient(*dualVector_,obj,x,tol);
295 g_->axpy(weight_,*dualVector_);
296 }
297
299 const Vector<Real> &v,
300 const std::vector<Real> &vstat,
301 const Vector<Real> &x,
302 const std::vector<Real> &xstat,
303 Real &tol) {
304 Real val = computeValue(obj,x,tol);
305 val_ += weight_ * val;
306 Real gv = computeGradVec(*dualVector_,obj,v,x,tol);
307 gv_ += weight_ * gv;
308 //g_->axpy(weight_,*dualVector_);
309 computeHessVec(*dualVector_,obj,v,x,tol);
310 //hv_->axpy(weight_,hv);
311 }
312
313 Real getValue(const Vector<Real> &x,
314 const std::vector<Real> &xstat,
315 SampleGenerator<Real> &sampler) {
316 // Compute expected value
317 Real ev(0);
318 sampler.sumAll(&val_,&ev,1);
319 // Compute deviation
320 Real diff(0), pf0(0), dev(0), one(1), weight(0);
321 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
322 values_->get(diff,sampler.getMyPoint(i));
323 weight = sampler.getMyWeight(i);
324 diff -= ev;
325 pf0 = positiveFunction_->evaluate(diff,0);
326 for ( uint p = 0; p < NumMoments_; p++ ) {
327 dev0_[p] += std::pow(pf0,order_[p]) * weight;
328 }
329 }
330 sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
331 for ( uint p = 0; p < NumMoments_; p++ ) {
332 dev += coeff_[p]*std::pow(des0_[p],one/order_[p]);
333 }
334 // Return mean plus deviation
335 return ev + dev;
336 }
337
339 std::vector<Real> &gstat,
340 const Vector<Real> &x,
341 const std::vector<Real> &xstat,
342 SampleGenerator<Real> &sampler) {
343 // Compute expected value
344 Real ev(0);
345 sampler.sumAll(&val_,&ev,1);
346 // Compute deviation
347 Real diff(0), pf0(0), pf1(0), c(0), one(1), zero(0), weight(0);
348 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
349 values_->get(diff,sampler.getMyPoint(i));
350 weight = sampler.getMyWeight(i);
351 diff -= ev;
352 pf0 = positiveFunction_->evaluate(diff,0);
353 pf1 = positiveFunction_->evaluate(diff,1);
354 for ( uint p = 0; p < NumMoments_; p++ ) {
355 dev0_[p] += weight * std::pow(pf0,order_[p]);
356 dev1_[p] += weight * std::pow(pf0,order_[p]-one) * pf1;
357 }
358 }
359 sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
360 sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
361 for ( uint p = 0; p < NumMoments_; p++ ) {
362 dev0_[p] = std::pow(des0_[p],one-one/order_[p]);
363 }
364 // Compute derivative
365 dualVector_->zero();
366 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
367 values_->get(diff,sampler.getMyPoint(i));
368 weight = sampler.getMyWeight(i);
369 diff -= ev;
370 pf0 = positiveFunction_->evaluate(diff,0);
371 pf1 = positiveFunction_->evaluate(diff,1);
372 c = zero;
373 for ( uint p = 0; p < NumMoments_; p++ ) {
374 if ( dev0_[p] > zero ) {
375 c += coeff_[p]/dev0_[p] * (std::pow(pf0,order_[p]-one)*pf1 - des1_[p]);
376 }
377 }
378 gradients_->get(*dualVector_,sampler.getMyPoint(i));
379 g_->axpy(weight*c,*dualVector_);
380 }
381 sampler.sumAll(*g_,g);
382 }
383
385 std::vector<Real> &hvstat,
386 const Vector<Real> &v,
387 const std::vector<Real> &vstat,
388 const Vector<Real> &x,
389 const std::vector<Real> &xstat,
390 SampleGenerator<Real> &sampler) {
391 // Compute expected value
392 std::vector<Real> myval(2), val(2);
393 myval[0] = val_;
394 myval[1] = gv_;
395 sampler.sumAll(&myval[0],&val[0],2);
396 Real ev = val[0], egv = val[1];
397 // Compute deviation
398 Real diff(0), pf0(0), pf1(0), pf2(0), zero(0), one(1), two(2);
399 Real cg(0), ch(0), diff1(0), diff2(0), diff3(0), weight(0), gv(0);
400 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
401 values_->get(diff,sampler.getMyPoint(i));
402 weight = sampler.getMyWeight(i);
403 diff -= ev;
404 pf0 = positiveFunction_->evaluate(diff,0);
405 pf1 = positiveFunction_->evaluate(diff,1);
406 pf2 = positiveFunction_->evaluate(diff,2);
407 for ( uint p = 0; p < NumMoments_; p++ ) {
408 dev0_[p] += weight * std::pow(pf0,order_[p]);
409 dev1_[p] += weight * std::pow(pf0,order_[p]-one) * pf1;
410 dev2_[p] += weight * std::pow(pf0,order_[p]-two) * pf1 * pf1;
411 dev3_[p] += weight * std::pow(pf0,order_[p]-one) * pf2;
412 }
413 }
414 sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
415 sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
416 sampler.sumAll(&dev2_[0],&des2_[0],NumMoments_);
417 sampler.sumAll(&dev3_[0],&des3_[0],NumMoments_);
418 for ( uint p = 0; p < NumMoments_; p++ ) {
419 devp_[p] = std::pow(des0_[p],two-one/order_[p]);
420 dev0_[p] = std::pow(des0_[p],one-one/order_[p]);
421 }
422 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
423 values_->get(diff,sampler.getMyPoint(i));
424 weight = sampler.getMyWeight(i);
425 diff -= ev;
426 gradvecs_->get(gv,sampler.getMyPoint(i));
427 pf0 = positiveFunction_->evaluate(diff,0);
428 pf1 = positiveFunction_->evaluate(diff,1);
429 pf2 = positiveFunction_->evaluate(diff,2);
430 for ( uint p = 0; p < NumMoments_; p++ ) {
431 gvp1_[p] += weight * (std::pow(pf0,order_[p]-one)*pf1-des1_[p]) *
432 (gv - egv);
433 gvp2_[p] += weight * (std::pow(pf0,order_[p]-two)*pf1*pf1-des2_[p]) *
434 (gv - egv);
435 gvp3_[p] += weight * (std::pow(pf0,order_[p]-one)*pf2-des3_[p]) *
436 (gv - egv);
437 }
438 }
439 sampler.sumAll(&gvp1_[0],&gvs1_[0],NumMoments_);
440 sampler.sumAll(&gvp2_[0],&gvs2_[0],NumMoments_);
441 sampler.sumAll(&gvp3_[0],&gvs3_[0],NumMoments_);
442 // Compute derivative
443 dualVector_->zero();
444 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
445 values_->get(diff,sampler.getMyPoint(i));
446 weight = sampler.getMyWeight(i);
447 diff -= ev;
448 gradvecs_->get(gv,sampler.getMyPoint(i));
449 pf0 = positiveFunction_->evaluate(diff,0);
450 pf1 = positiveFunction_->evaluate(diff,1);
451 pf2 = positiveFunction_->evaluate(diff,2);
452 cg = one;
453 ch = zero;
454 for ( uint p = 0; p < NumMoments_; p++ ) {
455 if ( dev0_[p] > zero ) {
456 diff1 = std::pow(pf0,order_[p]-one)*pf1-des1_[p];
457 diff2 = std::pow(pf0,order_[p]-two)*pf1*pf1*(gv-egv)-gvs2_[p];
458 diff3 = std::pow(pf0,order_[p]-one)*pf2*(gv-egv)-gvs3_[p];
459 cg += coeff_[p]*diff1/dev0_[p];
460 ch += coeff_[p]*(((order_[p]-one)*diff2+diff3)/dev0_[p] -
461 (order_[p]-one)*gvs1_[p]*diff1/devp_[p]);
462 }
463 }
464 gradients_->get(*g_,sampler.getMyPoint(i));
465 dualVector_->axpy(weight*ch,*g_);
466 hessvecs_->get(*hv_,sampler.getMyPoint(i));
467 dualVector_->axpy(weight*cg,*hv_);
468 }
469 sampler.sumAll(*dualVector_,hv);
470 }
471};
472
473}
474
475#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.
MeanDeviation(const Real order, const Real coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
std::vector< Real > gvs2_
std::vector< Real > des0_
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 initialize(const Vector< Real > &x)
Initialize temporary variables.
MeanDeviation(ROL::ParameterList &parlist)
Constructor.
Ptr< ScalarController< Real > > gradvecs_
std::vector< Real > gvp2_
Ptr< PositiveFunction< Real > > positiveFunction_
std::vector< Real > coeff_
std::vector< Real > dev0_
void setHessVecStorage(const Ptr< ScalarController< Real > > &gradvec_storage, const Ptr< VectorController< Real > > &hessvec_storage)
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
Ptr< ScalarController< Real > > values_
std::vector< Real > des1_
std::vector< Real > des3_
std::vector< Real > gvs1_
void setStorage(const Ptr< ScalarController< Real > > &value_storage, const Ptr< VectorController< Real > > &gradient_storage)
std::vector< Real > gvp3_
std::vector< Real >::size_type uint
Ptr< VectorController< Real > > gradients_
std::vector< Real > dev2_
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 > gvp1_
std::vector< Real > dev1_
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
std::vector< Real > order_
MeanDeviation(const std::vector< Real > &order, const std::vector< Real > &coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
std::vector< Real > gvs3_
std::vector< Real > devp_
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
Ptr< VectorController< Real > > hessvecs_
std::vector< Real > dev3_
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.
std::vector< Real > des2_
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)
Ptr< Vector< Real > > g_
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)
virtual void setStorage(const Ptr< ScalarController< Real > > &value_storage, const Ptr< VectorController< Real > > &gradient_storage)
Ptr< Vector< Real > > hv_
void computeGradient(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > dualVector_
virtual void setHessVecStorage(const Ptr< ScalarController< Real > > &gradvec_storage, const Ptr< VectorController< Real > > &hessvec_storage)
Real computeGradVec(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
virtual int numMySamples(void) const
virtual std::vector< Real > getMyPoint(const int i) const
void sumAll(Real *input, Real *output, int dim) const
virtual Real getMyWeight(const int i) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84