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
49#include "ROL_PlusFunction.hpp"
50#include "ROL_AbsoluteValue.hpp"
51
52#include "ROL_ParameterList.hpp"
53
75namespace ROL {
76
77template<class Real>
78class MeanVariance : public RandVarFunctional<Real> {
79 typedef typename std::vector<Real>::size_type uint;
80private:
81 Ptr<PositiveFunction<Real> > positiveFunction_;
82 std::vector<Real> order_;
83 std::vector<Real> coeff_;
85
86 Ptr<ScalarController<Real>> values_;
87 Ptr<ScalarController<Real>> gradvecs_;
88 Ptr<VectorController<Real>> gradients_;
89 Ptr<VectorController<Real>> hessvecs_;
90
91 using RandVarFunctional<Real>::val_;
92 using RandVarFunctional<Real>::gv_;
93 using RandVarFunctional<Real>::g_;
94 using RandVarFunctional<Real>::hv_;
96
97 using RandVarFunctional<Real>::point_;
99
104
105 void initializeMV(void) {
106 values_ = makePtr<ScalarController<Real>>();
107 gradvecs_ = makePtr<ScalarController<Real>>();
108 gradients_ = makePtr<VectorController<Real>>();
109 hessvecs_ = makePtr<VectorController<Real>>();
110
113 }
114
115 void checkInputs(void) {
116 int oSize = order_.size(), cSize = coeff_.size();
117 ROL_TEST_FOR_EXCEPTION((oSize!=cSize),std::invalid_argument,
118 ">>> ERROR (ROL::MeanVariance): Order and coefficient arrays have different sizes!");
119 Real zero(0), two(2);
120 for (int i = 0; i < oSize; i++) {
121 ROL_TEST_FOR_EXCEPTION((order_[i] < two), std::invalid_argument,
122 ">>> ERROR (ROL::MeanVariance): Element of order array out of range!");
123 ROL_TEST_FOR_EXCEPTION((coeff_[i] < zero), std::invalid_argument,
124 ">>> ERROR (ROL::MeanVariance): Element of coefficient array out of range!");
125 }
126 ROL_TEST_FOR_EXCEPTION(positiveFunction_ == nullPtr, std::invalid_argument,
127 ">>> ERROR (ROL::MeanVariance): PositiveFunction pointer is null!");
128 initializeMV();
129 }
130
131public:
141 MeanVariance( const Real order, const Real coeff,
142 const Ptr<PositiveFunction<Real> > &pf )
143 : RandVarFunctional<Real>(), positiveFunction_(pf) {
144 order_.clear(); order_.push_back(order);
145 coeff_.clear(); coeff_.push_back(coeff);
146 checkInputs();
147 NumMoments_ = order_.size();
148 }
149
159 MeanVariance( const std::vector<Real> &order,
160 const std::vector<Real> &coeff,
161 const Ptr<PositiveFunction<Real> > &pf )
162 : RandVarFunctional<Real>(), positiveFunction_(pf) {
163 order_.clear(); coeff_.clear();
164 for ( uint i = 0; i < order.size(); i++ ) {
165 order_.push_back(order[i]);
166 }
167 for ( uint i = 0; i < coeff.size(); i++ ) {
168 coeff_.push_back(coeff[i]);
169 }
170 checkInputs();
171 NumMoments_ = order_.size();
172 }
173
185 MeanVariance( ROL::ParameterList &parlist )
186 : RandVarFunctional<Real>() {
187 ROL::ParameterList &list
188 = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Variance");
189 // Get data from parameter list
190 order_ = ROL::getArrayFromStringParameter<double>(list,"Orders");
191 coeff_ = ROL::getArrayFromStringParameter<double>(list,"Coefficients");
192 // Build (approximate) positive function
193 std::string type = list.get<std::string>("Deviation Type");
194 if ( type == "Upper" ) {
195 positiveFunction_ = makePtr<PlusFunction<Real>>(list);
196 }
197 else if ( type == "Absolute" ) {
198 positiveFunction_ = makePtr<AbsoluteValue<Real>>(list);
199 }
200 else {
201 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
202 ">>> (ROL::MeanVariance): Variance type is not recoginized!");
203 }
204 // Check inputs
205 checkInputs();
206 NumMoments_ = order_.size();
207 }
208
209 void setStorage(const Ptr<ScalarController<Real>> &value_storage,
210 const Ptr<VectorController<Real>> &gradient_storage) {
211 values_ = value_storage;
212 gradients_ = gradient_storage;
214 }
215
216 void setHessVecStorage(const Ptr<ScalarController<Real>> &gradvec_storage,
217 const Ptr<VectorController<Real>> &hessvec_storage) {
218 gradvecs_ = gradvec_storage;
219 hessvecs_ = hessvec_storage;
221 }
222
224 const Vector<Real> &x,
225 const std::vector<Real> &xstat,
226 Real &tol) {
227 Real val = computeValue(obj,x,tol);
228 val_ += weight_ * val;
229 }
230
231 Real getValue(const Vector<Real> &x,
232 const std::vector<Real> &xstat,
233 SampleGenerator<Real> &sampler) {
234 // Compute expected value
235 Real ev(0);
236 sampler.sumAll(&val_,&ev,1);
237 // Compute deviation
238 Real val(0), diff(0), pf0(0), var(0), weight(0);
239 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
240 values_->get(diff,sampler.getMyPoint(i));
241 weight = sampler.getMyWeight(i);
242 diff -= ev;
243 pf0 = positiveFunction_->evaluate(diff,0);
244 for ( uint p = 0; p < NumMoments_; p++ ) {
245 val += coeff_[p] * std::pow(pf0,order_[p]) * weight;
246 }
247 }
248 sampler.sumAll(&val,&var,1);
249 // Return mean plus deviation
250 return ev + var;
251 }
252
254 const Vector<Real> &x,
255 const std::vector<Real> &xstat,
256 Real &tol) {
257 Real val = computeValue(obj,x,tol);
258 val_ += weight_ * val;
259 computeGradient(*dualVector_,obj,x,tol);
260 g_->axpy(weight_,*dualVector_);
261 }
262
264 std::vector<Real> &gstat,
265 const Vector<Real> &x,
266 const std::vector<Real> &xstat,
267 SampleGenerator<Real> &sampler) {
268 // Compute expected value
269 Real ev(0), zero(0), one(1);
270 sampler.sumAll(&val_,&ev,1);
271 sampler.sumAll(*g_,g);
272 // Compute deviation
273 g_->zero(); dualVector_->zero();
274 Real diff(0), pf0(0), pf1(0), c(0), ec(0), ecs(0), weight(0);
275 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
276 values_->get(diff,sampler.getMyPoint(i));
277 weight = sampler.getMyWeight(i);
278 diff -= ev;
279 pf0 = positiveFunction_->evaluate(diff,0);
280 pf1 = positiveFunction_->evaluate(diff,1);
281 c = zero;
282 for ( uint p = 0; p < NumMoments_; p++ ) {
283 c += coeff_[p]*order_[p]*std::pow(pf0,order_[p]-one)*pf1;
284 }
285 ec += weight*c;
286 gradients_->get(*dualVector_,sampler.getMyPoint(i));
287 g_->axpy(weight*c,*dualVector_);
288 }
289 dualVector_->zero();
290 sampler.sumAll(&ec,&ecs,1);
291 g.scale(one-ecs);
292 sampler.sumAll(*g_,*dualVector_);
293 g.plus(*dualVector_);
294 }
295
297 const Vector<Real> &v,
298 const std::vector<Real> &vstat,
299 const Vector<Real> &x,
300 const std::vector<Real> &xstat,
301 Real &tol) {
302 Real val = computeValue(obj,x,tol);
303 val_ += weight_ * val;
304 Real gv = computeGradVec(*dualVector_,obj,v,x,tol);
305 gv_ += weight_ * gv;
306 g_->axpy(weight_,*dualVector_);
307 computeHessVec(*dualVector_,obj,v,x,tol);
308 hv_->axpy(weight_,*dualVector_);
309 }
310
312 std::vector<Real> &hvstat,
313 const Vector<Real> &v,
314 const std::vector<Real> &vstat,
315 const Vector<Real> &x,
316 const std::vector<Real> &xstat,
317 SampleGenerator<Real> &sampler) {
318 // Compute expected value
319 std::vector<Real> myval(2), val(2);
320 myval[0] = val_;
321 myval[1] = gv_;
322 sampler.sumAll(&myval[0],&val[0],2);
323 Real ev = myval[0], egv = myval[1];
324 dualVector_->zero();
325 sampler.sumAll(*g_,*dualVector_);
326 sampler.sumAll(*hv_,hv);
327 // Compute deviation
328 g_->zero(); hv_->zero();
329 Real diff(0), pf0(0), pf1(0), pf2(0), zero(0), one(1), two(2);
330 Real cg(0), ecg(0), ecgs(0), ch(0), ech(0), echs(0), weight(0), gv(0);
331 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) {
332 values_->get(diff,sampler.getMyPoint(i));
333 gradvecs_->get(gv,sampler.getMyPoint(i));
334 weight = sampler.getMyWeight(i);
335 diff -= ev;
336 pf0 = positiveFunction_->evaluate(diff,0);
337 pf1 = positiveFunction_->evaluate(diff,1);
338 pf2 = positiveFunction_->evaluate(diff,2);
339 cg = zero;
340 ch = zero;
341 for ( uint p = 0; p < NumMoments_; p++ ) {
342 cg += coeff_[p]*order_[p]*(gv-egv)*
343 ((order_[p]-one)*std::pow(pf0,order_[p]-two)*pf1*pf1+
344 std::pow(pf0,order_[p]-one)*pf2);
345 ch += coeff_[p]*order_[p]*std::pow(pf0,order_[p]-one)*pf1;
346 }
347 ecg += weight*cg;
348 ech += weight*ch;
349 gradients_->get(*hv_,sampler.getMyPoint(i));
350 g_->axpy(weight*cg,*hv_);
351 hessvecs_->get(*hv_,sampler.getMyPoint(i));
352 g_->axpy(weight*ch,*hv_);
353 }
354 sampler.sumAll(&ech,&echs,1);
355 hv.scale(one-echs);
356 sampler.sumAll(&ecg,&ecgs,1);
357 hv.axpy(-ecgs,*dualVector_);
358 dualVector_->zero();
359 sampler.sumAll(*g_,*dualVector_);
360 hv.plus(*dualVector_);
361 }
362};
363
364}
365
366#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides an interface for the mean plus a sum of arbitrary order variances.
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
MeanVariance(const Real order, const Real coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
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.
MeanVariance(ROL::ParameterList &parlist)
Constructor.
Ptr< ScalarController< Real > > gradvecs_
Ptr< ScalarController< Real > > values_
std::vector< Real >::size_type uint
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.
Ptr< VectorController< Real > > gradients_
Ptr< VectorController< Real > > hessvecs_
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
MeanVariance(const std::vector< Real > &order, const std::vector< Real > &coeff, const Ptr< PositiveFunction< Real > > &pf)
Constructor.
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< PositiveFunction< Real > > positiveFunction_
void setHessVecStorage(const Ptr< ScalarController< Real > > &gradvec_storage, const Ptr< VectorController< Real > > &hessvec_storage)
void setStorage(const Ptr< ScalarController< Real > > &value_storage, const Ptr< VectorController< Real > > &gradient_storage)
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
std::vector< Real > order_
std::vector< Real > coeff_
Provides the interface to evaluate objective functions.
Provides the interface to implement any functional that maps a random variable to a (extended) real n...
Real computeValue(Objective< Real > &obj, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > g_
void computeHessVec(Vector< Real > &hv, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
virtual void setStorage(const Ptr< ScalarController< Real > > &value_storage, const Ptr< VectorController< Real > > &gradient_storage)
Ptr< Vector< Real > > hv_
void computeGradient(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > dualVector_
virtual void setHessVecStorage(const Ptr< ScalarController< Real > > &gradvec_storage, const Ptr< VectorController< Real > > &hessvec_storage)
Real computeGradVec(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
virtual int numMySamples(void) const
virtual std::vector< Real > getMyPoint(const int i) const
void sumAll(Real *input, Real *output, int dim) const
virtual Real getMyWeight(const int i) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void scale(const Real alpha)=0
Compute where .
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153