ROL
ROL_CoherentExpUtility.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_COHERENTEXPUTILITY_HPP
45 #define ROL_COHERENTEXPUTILITY_HPP
46 
47 #include "ROL_RiskMeasure.hpp"
48 
64 namespace ROL {
65 
66 template<class Real>
67 class CoherentExpUtility : public RiskMeasure<Real> {
68 private:
70 
71  Teuchos::RCP<Vector<Real> > scaledGradient1_;
72  Teuchos::RCP<Vector<Real> > scaledGradient2_;
73  Real dval1_;
74  Real dval2_;
75  Real dval3_;
76 
77  Teuchos::RCP<Vector<Real> > dualVector1_;
78  Teuchos::RCP<Vector<Real> > dualVector2_;
79 
80  Real xstat_;
81  Real vstat_;
82 
83 public:
84  CoherentExpUtility(void) : RiskMeasure<Real>(), firstReset_(true),
85  dval1_(0), dval2_(0), dval3_(0), xstat_(0), vstat_(0) {}
86 
87  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
88  Real zero(0);
90  xstat_ = Teuchos::dyn_cast<const RiskVector<Real> >(x).getStatistic(0);
91  if ( firstReset_ ) {
92  scaledGradient1_ = (x0->dual()).clone();
93  scaledGradient2_ = (x0->dual()).clone();
94  dualVector1_ = (x0->dual()).clone();
95  dualVector2_ = (x0->dual()).clone();
96  firstReset_ = false;
97  }
98  scaledGradient1_->zero(); scaledGradient2_->zero();
99  dualVector1_->zero(); dualVector2_->zero();
100  dval1_ = zero; dval2_ = zero; dval3_ = zero;
101  }
102 
103  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
104  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
105  reset(x0,x);
106  v0 = Teuchos::rcp_const_cast<Vector<Real> >(
107  Teuchos::dyn_cast<const RiskVector<Real> >(v).getVector());
108  vstat_ = Teuchos::dyn_cast<const RiskVector<Real> >(v).getStatistic(0);
109  }
110 
111  void update(const Real val, const Real weight) {
112  RiskMeasure<Real>::val_ += weight * std::exp(val/xstat_);
113  }
114 
116  Real val = RiskMeasure<Real>::val_, ev(0);
117  sampler.sumAll(&val,&ev,1);
118  return xstat_*std::log(ev);
119  }
120 
121  void update(const Real val, const Vector<Real> &g, const Real weight) {
122  Real ev = std::exp(val/xstat_);
123  RiskMeasure<Real>::val_ += weight * ev;
124  RiskMeasure<Real>::gv_ += weight * ev * val;
125  RiskMeasure<Real>::g_->axpy(weight*ev,g);
126  }
127 
129  Real one(1);
130  // Perform sum over batches
131  std::vector<Real> myval(2,0), val(2,0);
132  myval[0] = RiskMeasure<Real>::val_;
133  myval[1] = RiskMeasure<Real>::gv_;
134  sampler.sumAll(&myval[0],&val[0],2);
136  // Compute partial derivatives
137  Real gstat = std::log(myval[0]) - myval[1]/(myval[0]*xstat_);
138  dualVector1_->scale(one/myval[0]);
139  // Set partial derivatives in g vector
140  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setVector(*dualVector1_);
141  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setStatistic(gstat);
142  }
143 
144  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
145  const Real weight) {
146  Real ev = std::exp(val/xstat_);
147  RiskMeasure<Real>::val_ += weight * ev;
148  RiskMeasure<Real>::gv_ += weight * ev * gv;
149  dval1_ += weight * ev * val;
150  dval2_ += weight * ev * val * val;
151  dval3_ += weight * ev * val * gv;
152  RiskMeasure<Real>::g_->axpy(weight*ev,g);
153  RiskMeasure<Real>::hv_->axpy(weight*ev,hv);
154  scaledGradient1_->axpy(weight*ev*gv,g);
155  scaledGradient2_->axpy(weight*ev*val,g);
156  }
157 
159  Real one(1);
160  std::vector<Real> myval(5,0), val(5,0);
161  myval[0] = RiskMeasure<Real>::val_;
162  myval[1] = RiskMeasure<Real>::gv_;
163  myval[2] = dval1_;
164  myval[3] = dval2_;
165  myval[4] = dval3_;
166  sampler.sumAll(&myval[0],&val[0],5);
167 
168  Real xs2 = xstat_*xstat_;
169  Real xs3 = xs2*xstat_;
170  Real v02 = val[0]*val[0];
171  Real h11 = (val[3]*val[0] - val[2]*val[2])/(v02*xs3) * vstat_;
172  Real h12 = (val[1]*val[2] - val[4]*val[0])/(v02*xs2);
173 
176  dualVector1_->axpy(one/xstat_,*dualVector2_);
177  dualVector1_->scale(one/val[0]);
178  dualVector2_->zero();
180  dualVector1_->axpy(vstat_*val[2]/(xs2*v02)-val[1]/(v02*xstat_),*dualVector2_);
181  dualVector2_->zero();
183  dualVector1_->axpy(-vstat_/val[0],*dualVector2_);
184 
185  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setVector(*dualVector1_);
186  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setStatistic(h11+h12);
187  }
188 };
189 
190 }
191 
192 #endif
Teuchos::RCP< Vector< Real > > dualVector2_
void sumAll(Real *input, Real *output, int dim) const
Teuchos::RCP< Vector< Real > > scaledGradient1_
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void update(const Real val, const Real weight)
Update internal risk measure storage for value computation.
void update(const Real val, const Vector< Real > &g, const Real weight)
Update internal risk measure storage for gradient computation.
Teuchos::RCP< Vector< Real > > dualVector1_
Provides the interface for the coherent entropic risk measure.
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
Real getValue(SampleGenerator< Real > &sampler)
Return risk measure value.
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
virtual 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 Vector< Real > &g, const Real gv, const Vector< Real > &hv, const Real weight)
Update internal risk measure storage for Hessian-time-a-vector computation.
Provides the interface to implement risk measures.
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.
Teuchos::RCP< Vector< Real > > scaledGradient2_