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 
55 namespace ROL {
56 
57 template<class Real>
58 class MeanDeviation : public RiskMeasure<Real> {
59  typedef typename std::vector<Real>::size_type uint;
60 private:
61 
62  Teuchos::RCP<PositiveFunction<Real> > positiveFunction_;
63 
64  Teuchos::RCP<Vector<Real> > dualVector1_;
65  Teuchos::RCP<Vector<Real> > dualVector2_;
66 
67  std::vector<Real> order_;
68  std::vector<Real> coeff_;
70 
71  std::vector<Real> weights_;
72  std::vector<Real> value_storage_;
73  std::vector<Real> gradvec_storage_;
74  std::vector<Teuchos::RCP<Vector<Real> > > gradient_storage_;
75  std::vector<Teuchos::RCP<Vector<Real> > > hessvec_storage_;
76 
77  std::vector<Real> dev0_;
78  std::vector<Real> dev1_;
79  std::vector<Real> dev2_;
80  std::vector<Real> dev3_;
81  std::vector<Real> des0_;
82  std::vector<Real> des1_;
83  std::vector<Real> des2_;
84  std::vector<Real> des3_;
85  std::vector<Real> devp_;
86  std::vector<Real> gvp1_;
87  std::vector<Real> gvp2_;
88  std::vector<Real> gvp3_;
89  std::vector<Real> gvs1_;
90  std::vector<Real> gvs2_;
91  std::vector<Real> gvs3_;
92 
94 
95  void initialize(void) {
96  dev0_.clear(); dev1_.clear(); dev2_.clear(); dev3_.clear();
97  des0_.clear(); des1_.clear(); des2_.clear(); des3_.clear();
98  devp_.clear();
99  gvp1_.clear(); gvp2_.clear(); gvp3_.clear();
100  gvs1_.clear(); gvs2_.clear(); gvs3_.clear();
101 
102  dev0_.resize(NumMoments_,0.0); dev1_.resize(NumMoments_,0.0);
103  dev2_.resize(NumMoments_,0.0); dev3_.resize(NumMoments_,0.0);
104  des0_.resize(NumMoments_,0.0); des1_.resize(NumMoments_,0.0);
105  des2_.resize(NumMoments_,0.0); des3_.resize(NumMoments_,0.0);
106  devp_.resize(NumMoments_,0.0);
107  gvp1_.resize(NumMoments_,0.0); gvp2_.resize(NumMoments_,0.0);
108  gvp3_.resize(NumMoments_,0.0);
109  gvs1_.resize(NumMoments_,0.0); gvs2_.resize(NumMoments_,0.0);
110  gvs3_.resize(NumMoments_,0.0);
111  }
112 
113  void clear(void) {
114  dev0_.assign(NumMoments_,0.0); dev1_.assign(NumMoments_,0.0);
115  dev2_.assign(NumMoments_,0.0); dev3_.assign(NumMoments_,0.0);
116  des0_.assign(NumMoments_,0.0); des1_.assign(NumMoments_,0.0);
117  des2_.assign(NumMoments_,0.0); des3_.assign(NumMoments_,0.0);
118  devp_.assign(NumMoments_,0.0);
119  gvp1_.assign(NumMoments_,0.0); gvp2_.assign(NumMoments_,0.0);
120  gvp3_.assign(NumMoments_,0.0);
121  gvs1_.assign(NumMoments_,0.0); gvs2_.assign(NumMoments_,0.0);
122  gvs3_.assign(NumMoments_,0.0);
123 
124  value_storage_.clear();
125  gradient_storage_.clear();
126  gradvec_storage_.clear();
127  hessvec_storage_.clear();
128  weights_.clear();
129  }
130 
131 public:
132 
133  MeanDeviation( Real order, Real coeff,
134  Teuchos::RCP<PositiveFunction<Real> > &pf )
135  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
136  order_.clear(); coeff_.clear();
137  order_.push_back((order < 2.0) ? 2.0 : order);
138  coeff_.push_back((coeff < 0.0) ? 1.0 : coeff);
139  NumMoments_ = order_.size();
140  initialize();
141  }
142 
143  MeanDeviation( std::vector<Real> &order, std::vector<Real> &coeff,
144  Teuchos::RCP<PositiveFunction<Real> > &pf )
145  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
146  order_.clear(); coeff_.clear();
147  NumMoments_ = order.size();
148  if ( NumMoments_ != coeff.size() ) {
149  coeff.resize(NumMoments_,1.0);
150  }
151  for ( uint i = 0; i < NumMoments_; i++ ) {
152  order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
153  coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
154  }
155  initialize();
156  }
157 
158  MeanDeviation( Teuchos::ParameterList &parlist )
159  : RiskMeasure<Real>(), firstReset_(true) {
160  Teuchos::ParameterList &list
161  = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation");
162  // Get data from parameter list
163  Teuchos::Array<Real> order
164  = Teuchos::getArrayFromStringParameter<double>(list,"Orders");
165  Teuchos::Array<Real> coeff
166  = Teuchos::getArrayFromStringParameter<double>(list,"Coefficients");
167  // Check inputs
168  NumMoments_ = order.size();
169  order_.clear(); coeff_.clear();
170  if ( NumMoments_ != static_cast<uint>(coeff.size()) ) {
171  coeff.resize(NumMoments_,1.0);
172  }
173  for ( uint i = 0; i < NumMoments_; i++ ) {
174  order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
175  coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
176  }
177  // Build (approximate) positive function
178  if ( list.get("Deviation Type","Upper") == "Upper" ) {
179  positiveFunction_ = Teuchos::rcp(new PlusFunction<Real>(list));
180  }
181  else {
182  positiveFunction_ = Teuchos::rcp(new AbsoluteValue<Real>(list));
183  }
184  initialize();
185  }
186 
187  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
189  if ( firstReset_ ) {
190  dualVector1_ = (x0->dual()).clone();
191  dualVector2_ = (x0->dual()).clone();
192  firstReset_ = false;
193  }
194  dualVector1_->zero(); dualVector2_->zero();
195  clear();
196  }
197 
198  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
199  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
200  reset(x0,x);
201  v0 = Teuchos::rcp_const_cast<Vector<Real> >(Teuchos::dyn_cast<const RiskVector<Real> >(
202  Teuchos::dyn_cast<const Vector<Real> >(v)).getVector());
203  }
204 
205  void update(const Real val, const Real weight) {
206  RiskMeasure<Real>::val_ += weight * val;
207  value_storage_.push_back(val);
208  weights_.push_back(weight);
209  }
210 
211  void update(const Real val, const Vector<Real> &g, const Real weight) {
212  RiskMeasure<Real>::val_ += weight * val;
213  RiskMeasure<Real>::g_->axpy(weight,g);
214  value_storage_.push_back(val);
215  gradient_storage_.push_back(g.clone());
216  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
217  it--;
218  (*it)->set(g);
219  weights_.push_back(weight);
220  }
221 
222  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
223  const Real weight) {
224  RiskMeasure<Real>::val_ += weight * val;
225  RiskMeasure<Real>::gv_ += weight * gv;
226  RiskMeasure<Real>::g_->axpy(weight,g);
227  RiskMeasure<Real>::hv_->axpy(weight,hv);
228  value_storage_.push_back(val);
229  gradient_storage_.push_back(g.clone());
230  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
231  it--;
232  (*it)->set(g);
233  gradvec_storage_.push_back(gv);
234  hessvec_storage_.push_back(hv.clone());
235  it = hessvec_storage_.end();
236  it--;
237  (*it)->set(hv);
238  weights_.push_back(weight);
239  }
240 
242  // Compute expected value
243  Real val = RiskMeasure<Real>::val_, ev = 0.0;
244  sampler.sumAll(&val,&ev,1);
245  // Compute deviation
246  Real diff = 0.0, pf0 = 0.0, dev = 0.0;
247  for ( uint i = 0; i < weights_.size(); i++ ) {
248  diff = value_storage_[i]-ev;
249  pf0 = positiveFunction_->evaluate(diff,0);
250  for ( uint p = 0; p < NumMoments_; p++ ) {
251  dev0_[p] += std::pow(pf0,order_[p]) * weights_[i];
252  }
253  }
254  sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
255  for ( uint p = 0; p < NumMoments_; p++ ) {
256  dev += coeff_[p]*std::pow(des0_[p],1.0/order_[p]);
257  }
258  // Return mean plus deviation
259  return ev + dev;
260  }
261 
263  // Compute expected value
264  Real val = RiskMeasure<Real>::val_, ev = 0.0;
265  sampler.sumAll(&val,&ev,1);
266  // Compute deviation
267  Real diff = 0.0, pf0 = 0.0, pf1 = 0.0, c = 0.0;
268  for ( uint i = 0; i < weights_.size(); i++ ) {
269  diff = value_storage_[i]-ev;
270  pf0 = positiveFunction_->evaluate(diff,0);
271  pf1 = positiveFunction_->evaluate(diff,1);
272  for ( uint p = 0; p < NumMoments_; p++ ) {
273  dev0_[p] += weights_[i] * std::pow(pf0,order_[p]);
274  dev1_[p] += weights_[i] * std::pow(pf0,order_[p]-1.0) * pf1;
275  }
276  }
277  sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
278  sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
279  for ( uint p = 0; p < NumMoments_; p++ ) {
280  dev0_[p] = std::pow(des0_[p],1.0-1.0/order_[p]);
281  }
282  // Compute derivative
283  for ( uint i = 0; i < weights_.size(); i++ ) {
284  c = 0.0;
285  diff = value_storage_[i]-ev;
286  pf0 = positiveFunction_->evaluate(diff,0);
287  pf1 = positiveFunction_->evaluate(diff,1);
288  for ( uint p = 0; p < NumMoments_; p++ ) {
289  if ( dev0_[p] > 0.0 ) {
290  c += coeff_[p]/dev0_[p] * (std::pow(pf0,order_[p]-1.0)*pf1 - des1_[p]);
291  }
292  }
293  dualVector1_->axpy(weights_[i]*c,*(gradient_storage_[i]));
294  }
295  dualVector1_->axpy(1.0,*(RiskMeasure<Real>::g_));
296  sampler.sumAll(*dualVector1_,*dualVector2_);
297  // Set RiskVector
298  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setVector(*dualVector2_);
299  }
300 
302  // Compute expected value
303  Real val = RiskMeasure<Real>::val_, ev = 0.0;
304  sampler.sumAll(&val,&ev,1);
305  Real gv = RiskMeasure<Real>::gv_, egv = 0.0;
306  sampler.sumAll(&gv,&egv,1);
307  // Compute deviation
308  Real diff = 0.0, pf0 = 0.0, pf1 = 0.0, pf2 = 0.0;
309  Real cg = 0.0, ch = 0.0, diff1 = 0.0, diff2 = 0.0, diff3 = 0.0;
310  for ( uint i = 0; i < weights_.size(); i++ ) {
311  diff = value_storage_[i]-ev;
312  pf0 = positiveFunction_->evaluate(diff,0);
313  pf1 = positiveFunction_->evaluate(diff,1);
314  pf2 = positiveFunction_->evaluate(diff,2);
315  for ( uint p = 0; p < NumMoments_; p++ ) {
316  dev0_[p] += weights_[i] * std::pow(pf0,order_[p]);
317  dev1_[p] += weights_[i] * std::pow(pf0,order_[p]-1.0) * pf1;
318  dev2_[p] += weights_[i] * std::pow(pf0,order_[p]-2.0) * pf1 * pf1;
319  dev3_[p] += weights_[i] * std::pow(pf0,order_[p]-1.0) * pf2;
320  }
321  }
322  sampler.sumAll(&dev0_[0],&des0_[0],NumMoments_);
323  sampler.sumAll(&dev1_[0],&des1_[0],NumMoments_);
324  sampler.sumAll(&dev2_[0],&des2_[0],NumMoments_);
325  sampler.sumAll(&dev3_[0],&des3_[0],NumMoments_);
326  for ( uint p = 0; p < NumMoments_; p++ ) {
327  devp_[p] = std::pow(des0_[p],2.0-1.0/order_[p]);
328  dev0_[p] = std::pow(des0_[p],1.0-1.0/order_[p]);
329  }
330  for ( uint i = 0; i < value_storage_.size(); i++ ) {
331  diff = value_storage_[i]-ev;
332  pf0 = positiveFunction_->evaluate(diff,0);
333  pf1 = positiveFunction_->evaluate(diff,1);
334  pf2 = positiveFunction_->evaluate(diff,2);
335  for ( uint p = 0; p < NumMoments_; p++ ) {
336  gvp1_[p] += weights_[i] * (std::pow(pf0,order_[p]-1.0)*pf1-des1_[p]) *
337  (gradvec_storage_[i] - egv);
338  gvp2_[p] += weights_[i] * (std::pow(pf0,order_[p]-2.0)*pf1*pf1-des2_[p]) *
339  (gradvec_storage_[i] - egv);
340  gvp3_[p] += weights_[i] * (std::pow(pf0,order_[p]-1.0)*pf2-des3_[p]) *
341  (gradvec_storage_[i] - egv);
342  }
343  }
344  sampler.sumAll(&gvp1_[0],&gvs1_[0],NumMoments_);
345  sampler.sumAll(&gvp2_[0],&gvs2_[0],NumMoments_);
346  sampler.sumAll(&gvp3_[0],&gvs3_[0],NumMoments_);
347  // Compute derivative
348  for ( uint i = 0; i < weights_.size(); i++ ) {
349  cg = 1.0;
350  ch = 0.0;
351  diff = value_storage_[i]-ev;
352  pf0 = positiveFunction_->evaluate(diff,0);
353  pf1 = positiveFunction_->evaluate(diff,1);
354  pf2 = positiveFunction_->evaluate(diff,2);
355  for ( uint p = 0; p < NumMoments_; p++ ) {
356  if ( dev0_[p] > 0.0 ) {
357  diff1 = std::pow(pf0,order_[p]-1.0)*pf1-des1_[p];
358  diff2 = std::pow(pf0,order_[p]-2.0)*pf1*pf1*(gradvec_storage_[i]-egv)-gvs2_[p];
359  diff3 = std::pow(pf0,order_[p]-1.0)*pf2*(gradvec_storage_[i]-egv)-gvs3_[p];
360  cg += coeff_[p]*diff1/dev0_[p];
361  ch += coeff_[p]*(((order_[p]-1.0)*diff2+diff3)/dev0_[p] -
362  (order_[p]-1.0)*gvs1_[p]*diff1/devp_[p]);
363  }
364  }
365  dualVector1_->axpy(weights_[i]*ch,*(gradient_storage_[i]));
366  dualVector1_->axpy(weights_[i]*cg,*(hessvec_storage_[i]));
367  }
368  sampler.sumAll(*dualVector1_,*dualVector2_);
369  // Fill RiskVector
370  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setVector(*dualVector2_);
371  }
372 };
373 
374 }
375 
376 #endif
Real getValue(SampleGenerator< Real > &sampler)
MeanDeviation(Real order, Real coeff, Teuchos::RCP< PositiveFunction< Real > > &pf)
std::vector< Real > dev2_
std::vector< Real > gvs3_
std::vector< Real > dev0_
std::vector< Real >::size_type uint
std::vector< Real > dev3_
Teuchos::RCP< Vector< Real > > dualVector2_
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
std::vector< Real > des0_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void sumAll(Real *input, Real *output, int dim) const
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x, Teuchos::RCP< Vector< Real > > &v0, const Vector< Real > &v)
std::vector< Real > gvp3_
Teuchos::RCP< const Vector< Real > > getVector() const
std::vector< Real > order_
MeanDeviation(Teuchos::ParameterList &parlist)
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)
std::vector< Real > gvs2_
void update(const Real val, const Vector< Real > &g, const Real gv, const Vector< Real > &hv, const Real weight)
void update(const Real val, const Vector< Real > &g, const Real weight)
virtual void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
MeanDeviation(std::vector< Real > &order, std::vector< Real > &coeff, Teuchos::RCP< PositiveFunction< Real > > &pf)
Teuchos::RCP< PositiveFunction< Real > > positiveFunction_
std::vector< Real > des1_
std::vector< Real > gvp1_
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
std::vector< Real > weights_
std::vector< Real > coeff_
std::vector< Teuchos::RCP< Vector< Real > > > hessvec_storage_
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
std::vector< Real > gvs1_