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 
76 namespace ROL {
77 
78 template<class Real>
79 class MeanVariance : public RiskMeasure<Real> {
80  typedef typename std::vector<Real>::size_type uint;
81 private:
82 
83  Teuchos::RCP<PositiveFunction<Real> > positiveFunction_;
84 
85  Teuchos::RCP<Vector<Real> > dualVector1_;
86  Teuchos::RCP<Vector<Real> > dualVector2_;
87  Teuchos::RCP<Vector<Real> > dualVector3_;
88  Teuchos::RCP<Vector<Real> > dualVector4_;
89 
90  std::vector<Real> order_;
91  std::vector<Real> coeff_;
93 
94  std::vector<Real> weights_;
95  std::vector<Real> value_storage_;
96  std::vector<Teuchos::RCP<Vector<Real> > > gradient_storage_;
97  std::vector<Teuchos::RCP<Vector<Real> > > hessvec_storage_;
98  std::vector<Real> gradvec_storage_;
99 
101 
102  void checkInputs(void) const {
103  int oSize = order_.size(), cSize = coeff_.size();
104  TEUCHOS_TEST_FOR_EXCEPTION((oSize!=cSize),std::invalid_argument,
105  ">>> ERROR (ROL::MeanVariance): Order and coefficient arrays have different sizes!");
106  Real zero(0), two(2);
107  for (int i = 0; i < oSize; i++) {
108  TEUCHOS_TEST_FOR_EXCEPTION((order_[i] < two), std::invalid_argument,
109  ">>> ERROR (ROL::MeanVariance): Element of order array out of range!");
110  TEUCHOS_TEST_FOR_EXCEPTION((coeff_[i] < zero), std::invalid_argument,
111  ">>> ERROR (ROL::MeanVariance): Element of coefficient array out of range!");
112  }
113  TEUCHOS_TEST_FOR_EXCEPTION(positiveFunction_ == Teuchos::null, std::invalid_argument,
114  ">>> ERROR (ROL::MeanVariance): PositiveFunction pointer is null!");
115  }
116 
117 public:
127  MeanVariance( const Real order, const Real coeff,
128  const Teuchos::RCP<PositiveFunction<Real> > &pf )
129  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
130  order_.clear(); order_.push_back(order);
131  coeff_.clear(); coeff_.push_back(coeff);
132  checkInputs();
133  NumMoments_ = order_.size();
134  }
135 
145  MeanVariance( const std::vector<Real> &order,
146  const std::vector<Real> &coeff,
147  const Teuchos::RCP<PositiveFunction<Real> > &pf )
148  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
149  order_.clear(); coeff_.clear();
150  for ( uint i = 0; i < order.size(); i++ ) {
151  order_.push_back(order[i]);
152  }
153  for ( uint i = 0; i < coeff.size(); i++ ) {
154  coeff_.push_back(coeff[i]);
155  }
156  checkInputs();
157  NumMoments_ = order_.size();
158  }
159 
171  MeanVariance( Teuchos::ParameterList &parlist )
172  : RiskMeasure<Real>(), firstReset_(true) {
173  Teuchos::ParameterList &list
174  = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Variance");
175  // Get data from parameter list
176  Teuchos::Array<Real> order
177  = Teuchos::getArrayFromStringParameter<double>(list,"Orders");
178  order_ = order.toVector();
179  Teuchos::Array<Real> coeff
180  = Teuchos::getArrayFromStringParameter<double>(list,"Coefficients");
181  coeff_ = coeff.toVector();
182  // Build (approximate) positive function
183  std::string type = list.get<std::string>("Deviation Type");
184  if ( type == "Upper" ) {
185  positiveFunction_ = Teuchos::rcp(new PlusFunction<Real>(list));
186  }
187  else if ( type == "Absolute" ) {
188  positiveFunction_ = Teuchos::rcp(new AbsoluteValue<Real>(list));
189  }
190  else {
191  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
192  ">>> (ROL::MeanDeviation): Deviation type is not recoginized!");
193  }
194  // Check inputs
195  checkInputs();
196  NumMoments_ = order.size();
197  }
198 
199  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
201  if ( firstReset_ ) {
202  dualVector1_ = (x0->dual()).clone();
203  dualVector2_ = (x0->dual()).clone();
204  dualVector3_ = (x0->dual()).clone();
205  dualVector4_ = (x0->dual()).clone();
206  firstReset_ = false;
207  }
208  dualVector1_->zero(); dualVector2_->zero();
209  dualVector3_->zero(); dualVector4_->zero();
210  value_storage_.clear();
211  gradient_storage_.clear();
212  gradvec_storage_.clear();
213  hessvec_storage_.clear();
214  weights_.clear();
215  }
216 
217  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
218  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
219  reset(x0,x);
220  v0 = Teuchos::rcp_const_cast<Vector<Real> >(Teuchos::dyn_cast<const RiskVector<Real> >(
221  Teuchos::dyn_cast<const Vector<Real> >(v)).getVector());
222  }
223 
224  void update(const Real val, const Real weight) {
225  RiskMeasure<Real>::val_ += weight * val;
226  value_storage_.push_back(val);
227  weights_.push_back(weight);
228  }
229 
231  // Compute expected value
232  Real val = RiskMeasure<Real>::val_, ev(0), zero(0);
233  sampler.sumAll(&val,&ev,1);
234  // Compute deviation
235  val = zero;
236  Real diff(0), pf0(0), var(0);
237  for ( uint i = 0; i < weights_.size(); i++ ) {
238  diff = value_storage_[i]-ev;
239  pf0 = positiveFunction_->evaluate(diff,0);
240  for ( uint p = 0; p < NumMoments_; p++ ) {
241  val += coeff_[p] * std::pow(pf0,order_[p]) * weights_[i];
242  }
243  }
244  sampler.sumAll(&val,&var,1);
245  // Return mean plus deviation
246  return ev + var;
247  }
248 
249  void update(const Real val, const Vector<Real> &g, const Real weight) {
250  RiskMeasure<Real>::val_ += weight * val;
251  RiskMeasure<Real>::g_->axpy(weight,g);
252  value_storage_.push_back(val);
253  gradient_storage_.push_back(g.clone());
254  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
255  it--;
256  (*it)->set(g);
257  weights_.push_back(weight);
258  }
259 
261  // Compute expected value
262  Real val = RiskMeasure<Real>::val_, ev(0), zero(0), one(1);
263  sampler.sumAll(&val,&ev,1);
265  // Compute deviation
266  Real diff(0), pf0(0), pf1(0), c(0), ec(0), ecs(0);
267  for ( uint i = 0; i < weights_.size(); i++ ) {
268  c = zero;
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  c += coeff_[p]*order_[p]*std::pow(pf0,order_[p]-one)*pf1;
274  }
275  ec += weights_[i]*c;
276  dualVector1_->axpy(weights_[i]*c,*(gradient_storage_[i]));
277  }
278  sampler.sumAll(&ec,&ecs,1);
279  dualVector3_->scale(one-ecs);
280  sampler.sumAll(*dualVector1_,*dualVector2_);
281  dualVector3_->plus(*dualVector2_);
282  // Set RiskVector
283  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setVector(*dualVector3_);
284  }
285 
286  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
287  const Real weight) {
288  RiskMeasure<Real>::val_ += weight * val;
289  RiskMeasure<Real>::gv_ += weight * gv;
290  RiskMeasure<Real>::g_->axpy(weight,g);
291  RiskMeasure<Real>::hv_->axpy(weight,hv);
292  value_storage_.push_back(val);
293  gradient_storage_.push_back(g.clone());
294  typename std::vector<Teuchos::RCP<Vector<Real> > >::iterator it = gradient_storage_.end();
295  it--;
296  (*it)->set(g);
297  gradvec_storage_.push_back(gv);
298  hessvec_storage_.push_back(hv.clone());
299  it = hessvec_storage_.end();
300  it--;
301  (*it)->set(hv);
302  weights_.push_back(weight);
303  }
304 
306  hv.zero();
307  // Compute expected value
308  std::vector<Real> myval(2), val(2);
309  myval[0] = RiskMeasure<Real>::val_;
310  myval[1] = RiskMeasure<Real>::gv_;
311  sampler.sumAll(&myval[0],&val[0],2);
312  Real ev = myval[0], egv = myval[1];
315  // Compute deviation
316  Real diff(0), pf0(0), pf1(0), pf2(0), zero(0), one(1), two(2);
317  Real cg(0), ecg(0), ecgs(0), ch(0), ech(0), echs(0);
318  for ( uint i = 0; i < weights_.size(); i++ ) {
319  cg = zero;
320  ch = zero;
321  diff = value_storage_[i]-ev;
322  pf0 = positiveFunction_->evaluate(diff,0);
323  pf1 = positiveFunction_->evaluate(diff,1);
324  pf2 = positiveFunction_->evaluate(diff,2);
325  for ( uint p = 0; p < NumMoments_; p++ ) {
326  cg += coeff_[p]*order_[p]*(gradvec_storage_[i]-egv)*
327  ((order_[p]-one)*std::pow(pf0,order_[p]-two)*pf1*pf1+
328  std::pow(pf0,order_[p]-one)*pf2);
329  ch += coeff_[p]*order_[p]*std::pow(pf0,order_[p]-one)*pf1;
330  }
331  ecg += weights_[i]*cg;
332  ech += weights_[i]*ch;
333  dualVector1_->axpy(weights_[i]*cg,*(gradient_storage_[i]));
334  dualVector1_->axpy(weights_[i]*ch,*(hessvec_storage_[i]));
335  }
336  sampler.sumAll(&ech,&echs,1);
337  dualVector4_->scale(one-echs);
338  sampler.sumAll(&ecg,&ecgs,1);
339  dualVector4_->axpy(-ecgs,*dualVector3_);
340  sampler.sumAll(*dualVector1_,*dualVector2_);
341  dualVector4_->plus(*dualVector2_);
342  // Set RiskVector
343  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setVector(*dualVector4_);
344  }
345 };
346 
347 }
348 
349 #endif
Real getValue(SampleGenerator< Real > &sampler)
Return risk measure value.
std::vector< Teuchos::RCP< Vector< Real > > > gradient_storage_
std::vector< Real > gradvec_storage_
Teuchos::RCP< PositiveFunction< Real > > positiveFunction_
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
void update(const Real val, const Real weight)
Update internal risk measure storage for value computation.
Teuchos::RCP< Vector< Real > > dualVector4_
std::vector< Real > order_
void sumAll(Real *input, Real *output, int dim) const
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
std::vector< Real > weights_
Teuchos::RCP< Vector< Real > > dualVector2_
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)
Update internal risk measure storage for Hessian-time-a-vector computation.
Teuchos::RCP< const Vector< Real > > getVector(void) const
void update(const Real val, const Vector< Real > &g, const Real weight)
Update internal risk measure storage for gradient computation.
MeanVariance(Teuchos::ParameterList &parlist)
Constructor.
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
MeanVariance(const Real order, const Real coeff, const Teuchos::RCP< PositiveFunction< Real > > &pf)
Constructor.
void checkInputs(void) const
MeanVariance(const std::vector< Real > &order, const std::vector< Real > &coeff, const Teuchos::RCP< PositiveFunction< Real > > &pf)
Constructor.
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)
Return risk measure Hessian-times-a-vector.
virtual void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
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.
Provides an interface for the mean plus a sum of arbitrary order variances.
Provides the interface to implement risk measures.
std::vector< Real > value_storage_
Teuchos::RCP< Vector< Real > > dualVector1_