ROL
ROL_MeanDeviation.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_MEANDEVIATION_HPP
11#define ROL_MEANDEVIATION_HPP
12
14#include "ROL_ParameterList.hpp"
16#include "ROL_PlusFunction.hpp"
17#include "ROL_AbsoluteValue.hpp"
18
43namespace ROL {
44
45template<class Real>
46class MeanDeviation : public RandVarFunctional<Real> {
47 typedef typename std::vector<Real>::size_type uint;
48private:
49 Ptr<PositiveFunction<Real> > positiveFunction_;
50 std::vector<Real> order_;
51 std::vector<Real> coeff_;
53
54 std::vector<Real> dev0_;
55 std::vector<Real> dev1_;
56 std::vector<Real> dev2_;
57 std::vector<Real> dev3_;
58 std::vector<Real> des0_;
59 std::vector<Real> des1_;
60 std::vector<Real> des2_;
61 std::vector<Real> des3_;
62 std::vector<Real> devp_;
63 std::vector<Real> gvp1_;
64 std::vector<Real> gvp2_;
65 std::vector<Real> gvp3_;
66 std::vector<Real> gvs1_;
67 std::vector<Real> gvs2_;
68 std::vector<Real> gvs3_;
69
70 Ptr<ScalarController<Real>> values_;
71 Ptr<ScalarController<Real>> gradvecs_;
72 Ptr<VectorController<Real>> gradients_;
73 Ptr<VectorController<Real>> hessvecs_;
74
75 using RandVarFunctional<Real>::val_;
76 using RandVarFunctional<Real>::gv_;
77 using RandVarFunctional<Real>::g_;
78 using RandVarFunctional<Real>::hv_;
80
81 using RandVarFunctional<Real>::point_;
83
88
89 void initializeStorage(void) {
90 NumMoments_ = order_.size();
91
92 dev0_.clear(); dev1_.clear(); dev2_.clear(); dev3_.clear();
93 des0_.clear(); des1_.clear(); des2_.clear(); des3_.clear();
94 devp_.clear();
95 gvp1_.clear(); gvp2_.clear(); gvp3_.clear();
96 gvs1_.clear(); gvs2_.clear(); gvs3_.clear();
97
98 dev0_.resize(NumMoments_); dev1_.resize(NumMoments_);
99 dev2_.resize(NumMoments_); dev3_.resize(NumMoments_);
100 des0_.resize(NumMoments_); des1_.resize(NumMoments_);
101 des2_.resize(NumMoments_); des3_.resize(NumMoments_);
102 devp_.resize(NumMoments_);
103 gvp1_.resize(NumMoments_); gvp2_.resize(NumMoments_);
104 gvp3_.resize(NumMoments_);
105 gvs1_.resize(NumMoments_); gvs2_.resize(NumMoments_);
106 gvs3_.resize(NumMoments_);
107
108 values_ = makePtr<ScalarController<Real>>();
109 gradvecs_ = makePtr<ScalarController<Real>>();
110 gradients_ = makePtr<VectorController<Real>>();
111 hessvecs_ = makePtr<VectorController<Real>>();
112
115 }
116
117 void clear(void) {
118 const Real zero(0);
119 dev0_.assign(NumMoments_,zero); dev1_.assign(NumMoments_,zero);
120 dev2_.assign(NumMoments_,zero); dev3_.assign(NumMoments_,zero);
121 des0_.assign(NumMoments_,zero); des1_.assign(NumMoments_,zero);
122 des2_.assign(NumMoments_,zero); des3_.assign(NumMoments_,zero);
123 devp_.assign(NumMoments_,zero);
124 gvp1_.assign(NumMoments_,zero); gvp2_.assign(NumMoments_,zero);
125 gvp3_.assign(NumMoments_,zero);
126 gvs1_.assign(NumMoments_,zero); gvs2_.assign(NumMoments_,zero);
127 gvs3_.assign(NumMoments_,zero);
128
129 gradvecs_->reset();
130 hessvecs_->reset();
131 }
132
133 void checkInputs(void) {
134 int oSize = order_.size(), cSize = coeff_.size();
135 ROL_TEST_FOR_EXCEPTION((oSize!=cSize),std::invalid_argument,
136 ">>> ERROR (ROL::MeanDeviation): Order and coefficient arrays have different sizes!");
137 Real zero(0), two(2);
138 for (int i = 0; i < oSize; i++) {
139 ROL_TEST_FOR_EXCEPTION((order_[i] < two), std::invalid_argument,
140 ">>> ERROR (ROL::MeanDeviation): Element of order array out of range!");
141 ROL_TEST_FOR_EXCEPTION((coeff_[i] < zero), std::invalid_argument,
142 ">>> ERROR (ROL::MeanDeviation): Element of coefficient array out of range!");
143 }
144 ROL_TEST_FOR_EXCEPTION(positiveFunction_ == nullPtr, std::invalid_argument,
145 ">>> ERROR (ROL::MeanDeviation): PositiveFunction pointer is null!");
147 }
148
149public:
159 MeanDeviation( const Real order, const Real coeff,
160 const Ptr<PositiveFunction<Real> > &pf )
161 : RandVarFunctional<Real>(), positiveFunction_(pf) {
162 order_.clear(); order_.push_back(order);
163 coeff_.clear(); coeff_.push_back(coeff);
164 checkInputs();
165 }
166
176 MeanDeviation( const std::vector<Real> &order,
177 const std::vector<Real> &coeff,
178 const Ptr<PositiveFunction<Real> > &pf )
179 : RandVarFunctional<Real>(), positiveFunction_(pf) {
180 order_.clear(); coeff_.clear();
181 for ( uint i = 0; i < order.size(); i++ ) {
182 order_.push_back(order[i]);
183 }
184 for ( uint i = 0; i < coeff.size(); i++ ) {
185 coeff_.push_back(coeff[i]);
186 }
187 checkInputs();
188 }
189
201 MeanDeviation( ROL::ParameterList &parlist )
202 : RandVarFunctional<Real>() {
203 ROL::ParameterList &list
204 = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation");
205
206 // Get data from parameter list
207 order_ = ROL::getArrayFromStringParameter<double>(list,"Orders");
208
209 coeff_ = ROL::getArrayFromStringParameter<double>(list,"Coefficients");
210
211 // Build (approximate) positive function
212 std::string type = list.get<std::string>("Deviation Type");
213 if ( type == "Upper" ) {
214 positiveFunction_ = makePtr<PlusFunction<Real>>(list);
215 }
216 else if ( type == "Absolute" ) {
217 positiveFunction_ = makePtr<AbsoluteValue<Real>>(list);
218 }
219 else {
220 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
221 ">>> (ROL::MeanDeviation): Deviation type is not recoginized!");
222 }
223 // Check inputs
224 checkInputs();
225 }
226
227 void setStorage(const Ptr<ScalarController<Real>> &value_storage,
228 const Ptr<VectorController<Real>> &gradient_storage) {
229 values_ = value_storage;
230 gradients_ = gradient_storage;
232 }
233
234 void setHessVecStorage(const Ptr<ScalarController<Real>> &gradvec_storage,
235 const Ptr<VectorController<Real>> &hessvec_storage) {
236 gradvecs_ = gradvec_storage;
237 hessvecs_ = hessvec_storage;
239 }
240
245
247 const Vector<Real> &x,
248 const std::vector<Real> &xstat,
249 Real &tol) {
250 Real val = computeValue(obj,x,tol);
251 val_ += weight_ * val;
252 }
253
255 const Vector<Real> &x,
256 const std::vector<Real> &xstat,
257 Real &tol) {
258 Real val = computeValue(obj,x,tol);
259 val_ += weight_ * val;
260 computeGradient(*dualVector_,obj,x,tol);
261 g_->axpy(weight_,*dualVector_);
262 }
263
265 const Vector<Real> &v,
266 const std::vector<Real> &vstat,
267 const Vector<Real> &x,
268 const std::vector<Real> &xstat,
269 Real &tol) {
270 Real val = computeValue(obj,x,tol);
271 val_ += weight_ * val;
272 Real gv = computeGradVec(*dualVector_,obj,v,x,tol);
273 gv_ += weight_ * gv;
274 //g_->axpy(weight_,*dualVector_);
275 computeHessVec(*dualVector_,obj,v,x,tol);
276 //hv_->axpy(weight_,hv);
277 }
278
279 Real getValue(const Vector<Real> &x,
280 const std::vector<Real> &xstat,
281 SampleGenerator<Real> &sampler) {
282 // Compute expected value
283 Real ev(0);
284 sampler.sumAll(&val_,&ev,1);
285 // Compute deviation
286 Real diff(0), pf0(0), dev(0), one(1), weight(0);
287 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
288 values_->get(diff,sampler.getMyPoint(i));
289 weight = sampler.getMyWeight(i);
290 diff -= ev;
291 pf0 = positiveFunction_->evaluate(diff,0);
292 for ( uint p = 0; p < NumMoments_; p++ ) {
293 dev0_[p] += std::pow(pf0,order_[p]) * weight;
294 }
295 }
296 sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
297 for ( uint p = 0; p < NumMoments_; p++ ) {
298 dev += coeff_[p]*std::pow(des0_[p],one/order_[p]);
299 }
300 // Return mean plus deviation
301 return ev + dev;
302 }
303
305 std::vector<Real> &gstat,
306 const Vector<Real> &x,
307 const std::vector<Real> &xstat,
308 SampleGenerator<Real> &sampler) {
309 // Compute expected value
310 Real ev(0);
311 sampler.sumAll(&val_,&ev,1);
312 // Compute deviation
313 Real diff(0), pf0(0), pf1(0), c(0), one(1), zero(0), weight(0);
314 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
315 values_->get(diff,sampler.getMyPoint(i));
316 weight = sampler.getMyWeight(i);
317 diff -= ev;
318 pf0 = positiveFunction_->evaluate(diff,0);
319 pf1 = positiveFunction_->evaluate(diff,1);
320 for ( uint p = 0; p < NumMoments_; p++ ) {
321 dev0_[p] += weight * std::pow(pf0,order_[p]);
322 dev1_[p] += weight * std::pow(pf0,order_[p]-one) * pf1;
323 }
324 }
325 sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
326 sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
327 for ( uint p = 0; p < NumMoments_; p++ ) {
328 dev0_[p] = std::pow(des0_[p],one-one/order_[p]);
329 }
330 // Compute derivative
331 dualVector_->zero();
332 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
333 values_->get(diff,sampler.getMyPoint(i));
334 weight = sampler.getMyWeight(i);
335 diff -= ev;
336 pf0 = positiveFunction_->evaluate(diff,0);
337 pf1 = positiveFunction_->evaluate(diff,1);
338 c = zero;
339 for ( uint p = 0; p < NumMoments_; p++ ) {
340 if ( dev0_[p] > zero ) {
341 c += coeff_[p]/dev0_[p] * (std::pow(pf0,order_[p]-one)*pf1 - des1_[p]);
342 }
343 }
344 gradients_->get(*dualVector_,sampler.getMyPoint(i));
345 g_->axpy(weight*c,*dualVector_);
346 }
347 sampler.sumAll(*g_,g);
348 }
349
351 std::vector<Real> &hvstat,
352 const Vector<Real> &v,
353 const std::vector<Real> &vstat,
354 const Vector<Real> &x,
355 const std::vector<Real> &xstat,
356 SampleGenerator<Real> &sampler) {
357 // Compute expected value
358 std::vector<Real> myval(2), val(2);
359 myval[0] = val_;
360 myval[1] = gv_;
361 sampler.sumAll(&myval[0],&val[0],2);
362 Real ev = val[0], egv = val[1];
363 // Compute deviation
364 Real diff(0), pf0(0), pf1(0), pf2(0), zero(0), one(1), two(2);
365 Real cg(0), ch(0), diff1(0), diff2(0), diff3(0), weight(0), gv(0);
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 pf2 = positiveFunction_->evaluate(diff,2);
373 for ( uint p = 0; p < NumMoments_; p++ ) {
374 dev0_[p] += weight * std::pow(pf0,order_[p]);
375 dev1_[p] += weight * std::pow(pf0,order_[p]-one) * pf1;
376 dev2_[p] += weight * std::pow(pf0,order_[p]-two) * pf1 * pf1;
377 dev3_[p] += weight * std::pow(pf0,order_[p]-one) * pf2;
378 }
379 }
380 sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
381 sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
382 sampler.sumAll(&dev2_[0],&des2_[0],NumMoments_);
383 sampler.sumAll(&dev3_[0],&des3_[0],NumMoments_);
384 for ( uint p = 0; p < NumMoments_; p++ ) {
385 devp_[p] = std::pow(des0_[p],two-one/order_[p]);
386 dev0_[p] = std::pow(des0_[p],one-one/order_[p]);
387 }
388 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
389 values_->get(diff,sampler.getMyPoint(i));
390 weight = sampler.getMyWeight(i);
391 diff -= ev;
392 gradvecs_->get(gv,sampler.getMyPoint(i));
393 pf0 = positiveFunction_->evaluate(diff,0);
394 pf1 = positiveFunction_->evaluate(diff,1);
395 pf2 = positiveFunction_->evaluate(diff,2);
396 for ( uint p = 0; p < NumMoments_; p++ ) {
397 gvp1_[p] += weight * (std::pow(pf0,order_[p]-one)*pf1-des1_[p]) *
398 (gv - egv);
399 gvp2_[p] += weight * (std::pow(pf0,order_[p]-two)*pf1*pf1-des2_[p]) *
400 (gv - egv);
401 gvp3_[p] += weight * (std::pow(pf0,order_[p]-one)*pf2-des3_[p]) *
402 (gv - egv);
403 }
404 }
405 sampler.sumAll(&gvp1_[0],&gvs1_[0],NumMoments_);
406 sampler.sumAll(&gvp2_[0],&gvs2_[0],NumMoments_);
407 sampler.sumAll(&gvp3_[0],&gvs3_[0],NumMoments_);
408 // Compute derivative
409 dualVector_->zero();
410 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
411 values_->get(diff,sampler.getMyPoint(i));
412 weight = sampler.getMyWeight(i);
413 diff -= ev;
414 gradvecs_->get(gv,sampler.getMyPoint(i));
415 pf0 = positiveFunction_->evaluate(diff,0);
416 pf1 = positiveFunction_->evaluate(diff,1);
417 pf2 = positiveFunction_->evaluate(diff,2);
418 cg = one;
419 ch = zero;
420 for ( uint p = 0; p < NumMoments_; p++ ) {
421 if ( dev0_[p] > zero ) {
422 diff1 = std::pow(pf0,order_[p]-one)*pf1-des1_[p];
423 diff2 = std::pow(pf0,order_[p]-two)*pf1*pf1*(gv-egv)-gvs2_[p];
424 diff3 = std::pow(pf0,order_[p]-one)*pf2*(gv-egv)-gvs3_[p];
425 cg += coeff_[p]*diff1/dev0_[p];
426 ch += coeff_[p]*(((order_[p]-one)*diff2+diff3)/dev0_[p] -
427 (order_[p]-one)*gvs1_[p]*diff1/devp_[p]);
428 }
429 }
430 gradients_->get(*g_,sampler.getMyPoint(i));
431 dualVector_->axpy(weight*ch,*g_);
432 hessvecs_->get(*hv_,sampler.getMyPoint(i));
433 dualVector_->axpy(weight*cg,*hv_);
434 }
435 sampler.sumAll(*dualVector_,hv);
436 }
437};
438
439}
440
441#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)
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)
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.