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"
15 #include "ROL_PositiveFunction.hpp"
16 #include "ROL_PlusFunction.hpp"
17 #include "ROL_AbsoluteValue.hpp"
18 
43 namespace ROL {
44 
45 template<class Real>
46 class MeanDeviation : public RandVarFunctional<Real> {
48 private:
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 
80 
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 
114  RandVarFunctional<Real>::setHessVecStorage(gradvecs_,hessvecs_);
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 
149 public:
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 
241  void initialize(const Vector<Real> &x) {
243  clear();
244  }
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
Provides the interface to evaluate objective functions.
void computeHessVec(Vector< Real > &hv, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
typename PV< Real >::size_type size_type
std::vector< Real > dev2_
std::vector< Real > gvs3_
Ptr< Vector< Real > > g_
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
MeanDeviation(const std::vector< Real > &order, const std::vector< Real > &coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
Real computeValue(Objective< Real > &obj, const Vector< Real > &x, Real &tol)
void setHessVecStorage(const Ptr< ScalarController< Real >> &gradvec_storage, const Ptr< VectorController< Real >> &hessvec_storage)
std::vector< Real > dev0_
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.
Ptr< Vector< Real > > hv_
Ptr< ScalarController< Real > > values_
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.
std::vector< Real >::size_type uint
std::vector< Real > dev3_
MeanDeviation(const Real order, const Real coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
virtual std::vector< Real > getMyPoint(const int i) const
virtual Real getMyWeight(const int i) const
Ptr< Vector< Real > > dualVector_
std::vector< Real > des0_
virtual void setStorage(const Ptr< ScalarController< Real >> &value_storage, const Ptr< VectorController< Real >> &gradient_storage)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual int numMySamples(void) const
MeanDeviation(ROL::ParameterList &parlist)
Constructor.
void sumAll(Real *input, Real *output, int dim) const
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides an interface for the mean plus a sum of arbitrary order deviations.
std::vector< Real > gvp3_
std::vector< Real > order_
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.
void initialize(const Vector< Real > &x)
Initialize temporary variables.
std::vector< Real > gvp2_
Ptr< ScalarController< Real > > gradvecs_
std::vector< Real > des3_
std::vector< Real > dev1_
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
std::vector< Real > des2_
std::vector< Real > devp_
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
void computeGradient(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &x, Real &tol)
void setStorage(const Ptr< ScalarController< Real >> &value_storage, const Ptr< VectorController< Real >> &gradient_storage)
std::vector< Real > gvs2_
Real computeGradVec(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Provides the interface to implement any functional that maps a random variable to a (extended) real n...
Ptr< VectorController< Real > > gradients_
std::vector< Real > des1_
virtual void setHessVecStorage(const Ptr< ScalarController< Real >> &gradvec_storage, const Ptr< VectorController< Real >> &hessvec_storage)
std::vector< Real > gvp1_
std::vector< Real > coeff_
Ptr< VectorController< Real > > hessvecs_
virtual void initialize(const Vector< Real > &x)
Initialize temporary variables.
std::vector< Real > gvs1_
Ptr< PositiveFunction< Real > > positiveFunction_