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 
47 #include "ROL_RiskMeasure.hpp"
48 #include "ROL_PositiveFunction.hpp"
49 #include "ROL_PlusFunction.hpp"
50 #include "ROL_AbsoluteValue.hpp"
51 
52 #include "Teuchos_ParameterList.hpp"
53 #include "Teuchos_Array.hpp"
54 
79 namespace ROL {
80 
81 template<class Real>
82 class MeanDeviation : public RiskMeasure<Real> {
83  typedef typename std::vector<Real>::size_type uint;
84 private:
85 
86  Teuchos::RCP<PositiveFunction<Real> > positiveFunction_;
87 
88  Teuchos::RCP<Vector<Real> > dualVector1_;
89  Teuchos::RCP<Vector<Real> > dualVector2_;
90 
91  std::vector<Real> order_;
92  std::vector<Real> coeff_;
94 
95  std::vector<Real> weights_;
96  std::vector<Real> value_storage_;
97  std::vector<Real> gradvec_storage_;
98  std::vector<Teuchos::RCP<Vector<Real> > > gradient_storage_;
99  std::vector<Teuchos::RCP<Vector<Real> > > hessvec_storage_;
100 
101  std::vector<Real> dev0_;
102  std::vector<Real> dev1_;
103  std::vector<Real> dev2_;
104  std::vector<Real> dev3_;
105  std::vector<Real> des0_;
106  std::vector<Real> des1_;
107  std::vector<Real> des2_;
108  std::vector<Real> des3_;
109  std::vector<Real> devp_;
110  std::vector<Real> gvp1_;
111  std::vector<Real> gvp2_;
112  std::vector<Real> gvp3_;
113  std::vector<Real> gvs1_;
114  std::vector<Real> gvs2_;
115  std::vector<Real> gvs3_;
116 
118 
119  void initialize(void) {
120  dev0_.clear(); dev1_.clear(); dev2_.clear(); dev3_.clear();
121  des0_.clear(); des1_.clear(); des2_.clear(); des3_.clear();
122  devp_.clear();
123  gvp1_.clear(); gvp2_.clear(); gvp3_.clear();
124  gvs1_.clear(); gvs2_.clear(); gvs3_.clear();
125 
126  dev0_.resize(NumMoments_); dev1_.resize(NumMoments_);
127  dev2_.resize(NumMoments_); dev3_.resize(NumMoments_);
128  des0_.resize(NumMoments_); des1_.resize(NumMoments_);
129  des2_.resize(NumMoments_); des3_.resize(NumMoments_);
130  devp_.resize(NumMoments_);
131  gvp1_.resize(NumMoments_); gvp2_.resize(NumMoments_);
132  gvp3_.resize(NumMoments_);
133  gvs1_.resize(NumMoments_); gvs2_.resize(NumMoments_);
134  gvs3_.resize(NumMoments_);
135  }
136 
137  void clear(void) {
138  Real zero(0);
139  dev0_.assign(NumMoments_,zero); dev1_.assign(NumMoments_,zero);
140  dev2_.assign(NumMoments_,zero); dev3_.assign(NumMoments_,zero);
141  des0_.assign(NumMoments_,zero); des1_.assign(NumMoments_,zero);
142  des2_.assign(NumMoments_,zero); des3_.assign(NumMoments_,zero);
143  devp_.assign(NumMoments_,zero);
144  gvp1_.assign(NumMoments_,zero); gvp2_.assign(NumMoments_,zero);
145  gvp3_.assign(NumMoments_,zero);
146  gvs1_.assign(NumMoments_,zero); gvs2_.assign(NumMoments_,zero);
147  gvs3_.assign(NumMoments_,zero);
148 
149  value_storage_.clear();
150  gradient_storage_.clear();
151  gradvec_storage_.clear();
152  hessvec_storage_.clear();
153  weights_.clear();
154  }
155 
156  void checkInputs(void) const {
157  int oSize = order_.size(), cSize = coeff_.size();
158  TEUCHOS_TEST_FOR_EXCEPTION((oSize!=cSize),std::invalid_argument,
159  ">>> ERROR (ROL::MeanDeviation): Order and coefficient arrays have different sizes!");
160  Real zero(0), two(2);
161  for (int i = 0; i < oSize; i++) {
162  TEUCHOS_TEST_FOR_EXCEPTION((order_[i] < two), std::invalid_argument,
163  ">>> ERROR (ROL::MeanDeviation): Element of order array out of range!");
164  TEUCHOS_TEST_FOR_EXCEPTION((coeff_[i] < zero), std::invalid_argument,
165  ">>> ERROR (ROL::MeanDeviation): Element of coefficient array out of range!");
166  }
167  TEUCHOS_TEST_FOR_EXCEPTION(positiveFunction_ == Teuchos::null, std::invalid_argument,
168  ">>> ERROR (ROL::MeanDeviation): PositiveFunction pointer is null!");
169  }
170 
171 public:
181  MeanDeviation( const Real order, const Real coeff,
182  const Teuchos::RCP<PositiveFunction<Real> > &pf )
183  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
184  order_.clear(); order_.push_back(order);
185  coeff_.clear(); coeff_.push_back(coeff);
186  checkInputs();
187  NumMoments_ = order_.size();
188  initialize();
189  }
190 
200  MeanDeviation( const std::vector<Real> &order,
201  const std::vector<Real> &coeff,
202  const Teuchos::RCP<PositiveFunction<Real> > &pf )
203  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
204  order_.clear(); coeff_.clear();
205  for ( uint i = 0; i < order.size(); i++ ) {
206  order_.push_back(order[i]);
207  }
208  for ( uint i = 0; i < coeff.size(); i++ ) {
209  coeff_.push_back(coeff[i]);
210  }
211  checkInputs();
212  NumMoments_ = order_.size();
213  initialize();
214  }
215 
227  MeanDeviation( Teuchos::ParameterList &parlist )
228  : RiskMeasure<Real>(), firstReset_(true) {
229  Teuchos::ParameterList &list
230  = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation");
231  // Get data from parameter list
232  Teuchos::Array<Real> order
233  = Teuchos::getArrayFromStringParameter<double>(list,"Orders");
234  order_ = order.toVector();
235  Teuchos::Array<Real> coeff
236  = Teuchos::getArrayFromStringParameter<double>(list,"Coefficients");
237  coeff_ = coeff.toVector();
238  // Build (approximate) positive function
239  std::string type = list.get<std::string>("Deviation Type");
240  if ( type == "Upper" ) {
241  positiveFunction_ = Teuchos::rcp(new PlusFunction<Real>(list));
242  }
243  else if ( type == "Absolute" ) {
244  positiveFunction_ = Teuchos::rcp(new AbsoluteValue<Real>(list));
245  }
246  else {
247  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
248  ">>> (ROL::MeanDeviation): Deviation type is not recoginized!");
249  }
250  // Check inputs
251  checkInputs();
252  NumMoments_ = order.size();
253  initialize();
254  }
255 
256  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
258  if ( firstReset_ ) {
259  dualVector1_ = (x0->dual()).clone();
260  dualVector2_ = (x0->dual()).clone();
261  firstReset_ = false;
262  }
263  dualVector1_->zero(); dualVector2_->zero();
264  clear();
265  }
266 
267  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
268  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
269  reset(x0,x);
270  v0 = Teuchos::rcp_const_cast<Vector<Real> >(Teuchos::dyn_cast<const RiskVector<Real> >(
271  Teuchos::dyn_cast<const Vector<Real> >(v)).getVector());
272  }
273 
274  void update(const Real val, const Real weight) {
275  RiskMeasure<Real>::val_ += weight * val;
276  value_storage_.push_back(val);
277  weights_.push_back(weight);
278  }
279 
280  void update(const Real val, const Vector<Real> &g, const Real weight) {
281  RiskMeasure<Real>::val_ += weight * val;
282  RiskMeasure<Real>::g_->axpy(weight,g);
283  value_storage_.push_back(val);
284  gradient_storage_.push_back(g.clone());
285  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
286  it--;
287  (*it)->set(g);
288  weights_.push_back(weight);
289  }
290 
291  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
292  const Real weight) {
293  RiskMeasure<Real>::val_ += weight * val;
294  RiskMeasure<Real>::gv_ += weight * gv;
295  RiskMeasure<Real>::g_->axpy(weight,g);
296  RiskMeasure<Real>::hv_->axpy(weight,hv);
297  value_storage_.push_back(val);
298  gradient_storage_.push_back(g.clone());
299  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
300  it--;
301  (*it)->set(g);
302  gradvec_storage_.push_back(gv);
303  hessvec_storage_.push_back(hv.clone());
304  it = hessvec_storage_.end();
305  it--;
306  (*it)->set(hv);
307  weights_.push_back(weight);
308  }
309 
311  // Compute expected value
312  Real val = RiskMeasure<Real>::val_, ev(0);
313  sampler.sumAll(&val,&ev,1);
314  // Compute deviation
315  Real diff(0), pf0(0), dev(0), one(1);
316  for ( uint i = 0; i < weights_.size(); i++ ) {
317  diff = value_storage_[i]-ev;
318  pf0 = positiveFunction_->evaluate(diff,0);
319  for ( uint p = 0; p < NumMoments_; p++ ) {
320  dev0_[p] += std::pow(pf0,order_[p]) * weights_[i];
321  }
322  }
323  sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
324  for ( uint p = 0; p < NumMoments_; p++ ) {
325  dev += coeff_[p]*std::pow(des0_[p],one/order_[p]);
326  }
327  // Return mean plus deviation
328  return ev + dev;
329  }
330 
332  // Compute expected value
333  Real val = RiskMeasure<Real>::val_, ev(0);
334  sampler.sumAll(&val,&ev,1);
335  // Compute deviation
336  Real diff(0), pf0(0), pf1(0), c(0), one(1), zero(0);
337  for ( uint i = 0; i < weights_.size(); i++ ) {
338  diff = value_storage_[i]-ev;
339  pf0 = positiveFunction_->evaluate(diff,0);
340  pf1 = positiveFunction_->evaluate(diff,1);
341  for ( uint p = 0; p < NumMoments_; p++ ) {
342  dev0_[p] += weights_[i] * std::pow(pf0,order_[p]);
343  dev1_[p] += weights_[i] * std::pow(pf0,order_[p]-one) * pf1;
344  }
345  }
346  sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
347  sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
348  for ( uint p = 0; p < NumMoments_; p++ ) {
349  dev0_[p] = std::pow(des0_[p],one-one/order_[p]);
350  }
351  // Compute derivative
352  for ( uint i = 0; i < weights_.size(); i++ ) {
353  c = zero;
354  diff = value_storage_[i]-ev;
355  pf0 = positiveFunction_->evaluate(diff,0);
356  pf1 = positiveFunction_->evaluate(diff,1);
357  for ( uint p = 0; p < NumMoments_; p++ ) {
358  if ( dev0_[p] > zero ) {
359  c += coeff_[p]/dev0_[p] * (std::pow(pf0,order_[p]-one)*pf1 - des1_[p]);
360  }
361  }
362  dualVector1_->axpy(weights_[i]*c,*(gradient_storage_[i]));
363  }
365  sampler.sumAll(*dualVector1_,*dualVector2_);
366  // Set RiskVector
367  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setVector(*dualVector2_);
368  }
369 
371  // Compute expected value
372  std::vector<Real> myval(2), val(2);
373  myval[0] = RiskMeasure<Real>::val_;
374  myval[1] = RiskMeasure<Real>::gv_;
375  sampler.sumAll(&myval[0],&val[0],2);
376  Real ev = val[0], egv = val[1];
377  // Compute deviation
378  Real diff(0), pf0(0), pf1(0), pf2(0), zero(0), one(1), two(2);
379  Real cg(0), ch(0), diff1(0), diff2(0), diff3(0);
380  for ( uint i = 0; i < weights_.size(); i++ ) {
381  diff = value_storage_[i]-ev;
382  pf0 = positiveFunction_->evaluate(diff,0);
383  pf1 = positiveFunction_->evaluate(diff,1);
384  pf2 = positiveFunction_->evaluate(diff,2);
385  for ( uint p = 0; p < NumMoments_; p++ ) {
386  dev0_[p] += weights_[i] * std::pow(pf0,order_[p]);
387  dev1_[p] += weights_[i] * std::pow(pf0,order_[p]-one) * pf1;
388  dev2_[p] += weights_[i] * std::pow(pf0,order_[p]-two) * pf1 * pf1;
389  dev3_[p] += weights_[i] * std::pow(pf0,order_[p]-one) * pf2;
390  }
391  }
392  sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
393  sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
394  sampler.sumAll(&dev2_[0],&des2_[0],NumMoments_);
395  sampler.sumAll(&dev3_[0],&des3_[0],NumMoments_);
396  for ( uint p = 0; p < NumMoments_; p++ ) {
397  devp_[p] = std::pow(des0_[p],two-one/order_[p]);
398  dev0_[p] = std::pow(des0_[p],one-one/order_[p]);
399  }
400  for ( uint i = 0; i < value_storage_.size(); i++ ) {
401  diff = value_storage_[i]-ev;
402  pf0 = positiveFunction_->evaluate(diff,0);
403  pf1 = positiveFunction_->evaluate(diff,1);
404  pf2 = positiveFunction_->evaluate(diff,2);
405  for ( uint p = 0; p < NumMoments_; p++ ) {
406  gvp1_[p] += weights_[i] * (std::pow(pf0,order_[p]-one)*pf1-des1_[p]) *
407  (gradvec_storage_[i] - egv);
408  gvp2_[p] += weights_[i] * (std::pow(pf0,order_[p]-two)*pf1*pf1-des2_[p]) *
409  (gradvec_storage_[i] - egv);
410  gvp3_[p] += weights_[i] * (std::pow(pf0,order_[p]-one)*pf2-des3_[p]) *
411  (gradvec_storage_[i] - egv);
412  }
413  }
414  sampler.sumAll(&gvp1_[0],&gvs1_[0],NumMoments_);
415  sampler.sumAll(&gvp2_[0],&gvs2_[0],NumMoments_);
416  sampler.sumAll(&gvp3_[0],&gvs3_[0],NumMoments_);
417  // Compute derivative
418  for ( uint i = 0; i < weights_.size(); i++ ) {
419  cg = one;
420  ch = zero;
421  diff = value_storage_[i]-ev;
422  pf0 = positiveFunction_->evaluate(diff,0);
423  pf1 = positiveFunction_->evaluate(diff,1);
424  pf2 = positiveFunction_->evaluate(diff,2);
425  for ( uint p = 0; p < NumMoments_; p++ ) {
426  if ( dev0_[p] > zero ) {
427  diff1 = std::pow(pf0,order_[p]-one)*pf1-des1_[p];
428  diff2 = std::pow(pf0,order_[p]-two)*pf1*pf1*(gradvec_storage_[i]-egv)-gvs2_[p];
429  diff3 = std::pow(pf0,order_[p]-one)*pf2*(gradvec_storage_[i]-egv)-gvs3_[p];
430  cg += coeff_[p]*diff1/dev0_[p];
431  ch += coeff_[p]*(((order_[p]-one)*diff2+diff3)/dev0_[p] -
432  (order_[p]-one)*gvs1_[p]*diff1/devp_[p]);
433  }
434  }
435  dualVector1_->axpy(weights_[i]*ch,*(gradient_storage_[i]));
436  dualVector1_->axpy(weights_[i]*cg,*(hessvec_storage_[i]));
437  }
438  sampler.sumAll(*dualVector1_,*dualVector2_);
439  // Fill RiskVector
440  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setVector(*dualVector2_);
441  }
442 };
443 
444 }
445 
446 #endif
Real getValue(SampleGenerator< Real > &sampler)
Return risk measure value.
std::vector< Real > dev2_
std::vector< Real > gvs3_
std::vector< Real > dev0_
std::vector< Real >::size_type uint
std::vector< Real > dev3_
void sumAll(Real *input, Real *output, int dim) const
Teuchos::RCP< Vector< Real > > dualVector2_
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
std::vector< Real > des0_
MeanDeviation(const std::vector< Real > &order, const std::vector< Real > &coeff, const Teuchos::RCP< PositiveFunction< Real > > &pf)
Constructor.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< const Vector< Real > > getVector(void) const
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
Provides an interface for the mean plus a sum of arbitrary order deviations.
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x, Teuchos::RCP< Vector< Real > > &v0, const Vector< Real > &v)
Reset internal risk measure storage. Called for Hessian-times-a-vector computation.
std::vector< Real > gvp3_
std::vector< Real > order_
MeanDeviation(Teuchos::ParameterList &parlist)
Constructor.
MeanDeviation(const Real order, const Real coeff, const Teuchos::RCP< PositiveFunction< Real > > &pf)
Constructor.
std::vector< Real > gvp2_
std::vector< Real > value_storage_
std::vector< Teuchos::RCP< Vector< Real > > > gradient_storage_
std::vector< Real > des3_
std::vector< Real > dev1_
std::vector< Real > gradvec_storage_
Teuchos::RCP< Vector< Real > > dualVector1_
std::vector< Real > des2_
std::vector< Real > devp_
void update(const Real val, const Real weight)
Update internal risk measure storage for value computation.
std::vector< Real > gvs2_
void update(const Real val, const Vector< Real > &g, const Real gv, const Vector< Real > &hv, const Real weight)
Update internal risk measure storage for Hessian-time-a-vector computation.
void update(const Real val, const Vector< Real > &g, const Real weight)
Update internal risk measure storage for gradient computation.
virtual void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
void checkInputs(void) const
Provides the interface to implement risk measures.
Teuchos::RCP< PositiveFunction< Real > > positiveFunction_
std::vector< Real > des1_
std::vector< Real > gvp1_
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
std::vector< Real > weights_
std::vector< Real > coeff_
std::vector< Teuchos::RCP< Vector< Real > > > hessvec_storage_
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.
std::vector< Real > gvs1_