ROL
ROL_MeanVariance.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_MEANVARIANCE_HPP
45 #define ROL_MEANVARIANCE_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 MeanVariance : 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  Teuchos::RCP<Vector<Real> > dualVector3_;
67  Teuchos::RCP<Vector<Real> > dualVector4_;
68 
69  std::vector<Real> order_;
70  std::vector<Real> coeff_;
72 
73  std::vector<Real> weights_;
74  std::vector<Real> value_storage_;
75  std::vector<Teuchos::RCP<Vector<Real> > > gradient_storage_;
76  std::vector<Teuchos::RCP<Vector<Real> > > hessvec_storage_;
77  std::vector<Real> gradvec_storage_;
78 
80 
81 public:
82 
83  MeanVariance( Real order, Real coeff,
84  Teuchos::RCP<PositiveFunction<Real> > &pf )
85  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
86  order_.clear(); coeff_.clear();
87  order_.push_back((order < 2.0) ? 2.0 : order);
88  coeff_.push_back((coeff < 0.0) ? 1.0 : coeff);
89  NumMoments_ = order_.size();
90  }
91  MeanVariance( std::vector<Real> &order, std::vector<Real> &coeff,
92  Teuchos::RCP<PositiveFunction<Real> > &pf )
93  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
94  order_.clear(); coeff_.clear();
95  NumMoments_ = order.size();
96  if ( NumMoments_ != coeff.size() ) {
97  coeff.resize(NumMoments_,1.0);
98  }
99  for ( uint i = 0; i < NumMoments_; i++ ) {
100  order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
101  coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
102  }
103  }
104  MeanVariance( Teuchos::ParameterList &parlist )
105  : RiskMeasure<Real>(), firstReset_(true) {
106  Teuchos::ParameterList &list
107  = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Variance");
108  // Get data from parameter list
109  Teuchos::Array<Real> order
110  = Teuchos::getArrayFromStringParameter<double>(list,"Orders");
111  Teuchos::Array<Real> coeff
112  = Teuchos::getArrayFromStringParameter<double>(list,"Coefficients");
113  // Check inputs
114  order_.clear(); coeff_.clear();
115  NumMoments_ = order.size();
116  if ( NumMoments_ != static_cast<uint>(coeff.size()) ) {
117  coeff.resize(NumMoments_,1.0);
118  }
119  for ( uint i = 0; i < NumMoments_; i++ ) {
120  order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
121  coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
122  }
123  // Build (approximate) positive function
124  if ( list.get("Deviation Type","Upper") == "Upper" ) {
125  positiveFunction_ = Teuchos::rcp(new PlusFunction<Real>(list));
126  }
127  else {
128  positiveFunction_ = Teuchos::rcp(new AbsoluteValue<Real>(list));
129  }
130  }
131 
132  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
134  if ( firstReset_ ) {
135  dualVector1_ = (x0->dual()).clone();
136  dualVector2_ = (x0->dual()).clone();
137  dualVector3_ = (x0->dual()).clone();
138  dualVector4_ = (x0->dual()).clone();
139  firstReset_ = false;
140  }
141  dualVector1_->zero(); dualVector2_->zero();
142  dualVector3_->zero(); dualVector4_->zero();
143  value_storage_.clear();
144  gradient_storage_.clear();
145  gradvec_storage_.clear();
146  hessvec_storage_.clear();
147  weights_.clear();
148  }
149 
150  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
151  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
152  reset(x0,x);
153  v0 = Teuchos::rcp_const_cast<Vector<Real> >(Teuchos::dyn_cast<const RiskVector<Real> >(
154  Teuchos::dyn_cast<const Vector<Real> >(v)).getVector());
155  }
156 
157  void update(const Real val, const Real weight) {
158  RiskMeasure<Real>::val_ += weight * val;
159  value_storage_.push_back(val);
160  weights_.push_back(weight);
161  }
162 
163  void update(const Real val, const Vector<Real> &g, const Real weight) {
164  RiskMeasure<Real>::val_ += weight * val;
165  RiskMeasure<Real>::g_->axpy(weight,g);
166  value_storage_.push_back(val);
167  gradient_storage_.push_back(g.clone());
168  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
169  it--;
170  (*it)->set(g);
171  weights_.push_back(weight);
172  }
173 
174  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
175  const Real weight) {
176  RiskMeasure<Real>::val_ += weight * val;
177  RiskMeasure<Real>::gv_ += weight * gv;
178  RiskMeasure<Real>::g_->axpy(weight,g);
179  RiskMeasure<Real>::hv_->axpy(weight,hv);
180  value_storage_.push_back(val);
181  gradient_storage_.push_back(g.clone());
182  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
183  it--;
184  (*it)->set(g);
185  gradvec_storage_.push_back(gv);
186  hessvec_storage_.push_back(hv.clone());
187  it = hessvec_storage_.end();
188  it--;
189  (*it)->set(hv);
190  weights_.push_back(weight);
191  }
192 
194  // Compute expected value
195  Real val = RiskMeasure<Real>::val_;
196  Real ev = 0.0;
197  sampler.sumAll(&val,&ev,1);
198  // Compute deviation
199  val = 0.0;
200  Real diff = 0.0, pf0 = 0.0, var = 0.0;
201  for ( uint i = 0; i < weights_.size(); i++ ) {
202  diff = value_storage_[i]-ev;
203  pf0 = positiveFunction_->evaluate(diff,0);
204  for ( uint p = 0; p < NumMoments_; p++ ) {
205  val += coeff_[p] * std::pow(pf0,order_[p]) * weights_[i];
206  }
207  }
208  sampler.sumAll(&val,&var,1);
209  // Return mean plus deviation
210  return ev + var;
211  }
212 
214  // Compute expected value
215  Real val = RiskMeasure<Real>::val_, ev = 0.0;
216  sampler.sumAll(&val,&ev,1);
217  sampler.sumAll(*(RiskMeasure<Real>::g_),*dualVector3_);
218  // Compute deviation
219  Real diff = 0.0, pf0 = 0.0, pf1 = 0.0, c = 0.0, ec = 0.0, ecs = 0.0;
220  for ( uint i = 0; i < weights_.size(); i++ ) {
221  c = 0.0;
222  diff = value_storage_[i]-ev;
223  pf0 = positiveFunction_->evaluate(diff,0);
224  pf1 = positiveFunction_->evaluate(diff,1);
225  for ( uint p = 0; p < NumMoments_; p++ ) {
226  c += coeff_[p]*order_[p]*std::pow(pf0,order_[p]-1.0)*pf1;
227  }
228  ec += weights_[i]*c;
229  dualVector1_->axpy(weights_[i]*c,*(gradient_storage_[i]));
230  }
231  sampler.sumAll(&ec,&ecs,1);
232  dualVector3_->scale(1.0-ecs);
233  sampler.sumAll(*dualVector1_,*dualVector2_);
234  dualVector3_->plus(*dualVector2_);
235  // Set RiskVector
236  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setVector(*dualVector3_);
237  }
238 
240  hv.zero();
241  // Compute expected value
242  Real val = RiskMeasure<Real>::val_, ev = 0.0;
243  sampler.sumAll(&val,&ev,1);
244  Real gv = RiskMeasure<Real>::gv_, egv = 0.0;
245  sampler.sumAll(&gv,&egv,1);
246  sampler.sumAll(*(RiskMeasure<Real>::g_),*dualVector3_);
247  sampler.sumAll(*(RiskMeasure<Real>::hv_),*dualVector4_);
248  // Compute deviation
249  Real diff = 0.0, pf0 = 0.0, pf1 = 0.0, pf2 = 0.0;
250  Real cg = 0.0, ecg = 0.0, ecgs = 0.0, ch = 0.0, ech = 0.0, echs = 0.0;
251  for ( uint i = 0; i < weights_.size(); i++ ) {
252  cg = 0.0;
253  ch = 0.0;
254  diff = value_storage_[i]-ev;
255  pf0 = positiveFunction_->evaluate(diff,0);
256  pf1 = positiveFunction_->evaluate(diff,1);
257  pf2 = positiveFunction_->evaluate(diff,2);
258  for ( uint p = 0; p < NumMoments_; p++ ) {
259  cg += coeff_[p]*order_[p]*(gradvec_storage_[i]-egv)*
260  ((order_[p]-1.0)*std::pow(pf0,order_[p]-2.0)*pf1*pf1+
261  std::pow(pf0,order_[p]-1.0)*pf2);
262  ch += coeff_[p]*order_[p]*std::pow(pf0,order_[p]-1.0)*pf1;
263  }
264  ecg += weights_[i]*cg;
265  ech += weights_[i]*ch;
266  dualVector1_->axpy(weights_[i]*cg,*(gradient_storage_[i]));
267  dualVector1_->axpy(weights_[i]*ch,*(hessvec_storage_[i]));
268  }
269  sampler.sumAll(&ech,&echs,1);
270  dualVector4_->scale(1.0-echs);
271  sampler.sumAll(&ecg,&ecgs,1);
272  dualVector4_->axpy(-ecgs,*dualVector3_);
273  sampler.sumAll(*dualVector1_,*dualVector2_);
274  dualVector4_->plus(*dualVector2_);
275  // Set RiskVector
276  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setVector(*dualVector4_);
277  }
278 };
279 
280 }
281 
282 #endif
Real getValue(SampleGenerator< Real > &sampler)
std::vector< Teuchos::RCP< Vector< Real > > > gradient_storage_
MeanVariance(Real order, Real coeff, Teuchos::RCP< PositiveFunction< Real > > &pf)
std::vector< Real > gradvec_storage_
Teuchos::RCP< PositiveFunction< Real > > positiveFunction_
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
void update(const Real val, const Real weight)
Teuchos::RCP< Vector< Real > > dualVector4_
std::vector< Real > order_
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
std::vector< Real > weights_
Teuchos::RCP< Vector< Real > > dualVector2_
MeanVariance(std::vector< Real > &order, std::vector< Real > &coeff, Teuchos::RCP< PositiveFunction< Real > > &pf)
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
std::vector< Real > coeff_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void update(const Real val, const Vector< Real > &g, const Real gv, const Vector< Real > &hv, const Real weight)
void sumAll(Real *input, Real *output, int dim) const
void update(const Real val, const Vector< Real > &g, const Real weight)
Teuchos::RCP< const Vector< Real > > getVector() const
MeanVariance(Teuchos::ParameterList &parlist)
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
Teuchos::RCP< Vector< Real > > dualVector3_
std::vector< Real >::size_type uint
std::vector< Teuchos::RCP< Vector< Real > > > hessvec_storage_
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
virtual void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x, Teuchos::RCP< Vector< Real > > &v0, const Vector< Real > &v)
std::vector< Real > value_storage_
Teuchos::RCP< Vector< Real > > dualVector1_