ROL
ROL_ExpectationQuad.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_EXPECTATIONQUAD_HPP
45 #define ROL_EXPECTATIONQUAD_HPP
46 
47 #include "ROL_RiskVector.hpp"
48 #include "ROL_RiskMeasure.hpp"
49 #include "ROL_Types.hpp"
50 
86 namespace ROL {
87 
88 template<class Real>
89 class ExpectationQuad : public RiskMeasure<Real> {
90 private:
91 
92  Teuchos::RCP<Vector<Real> > dualVector_;
93 
94  Real xstat_;
95  Real vstat_;
96 
98 
99 public:
100  ExpectationQuad(void) : RiskMeasure<Real>(), xstat_(0), vstat_(0), firstReset_(true) {}
101 
109  virtual Real regret(Real x, int deriv = 0) = 0;
110 
113  virtual void checkRegret(void) {
114  Real zero(0), half(0.5), two(2), one(1), oem3(1.e-3), fem4(5.e-4), p1(0.1);
115  // Check v(0) = 0
116  Real x = zero;
117  Real vx = regret(x,0);
118  std::cout << std::right << std::setw(20) << "CHECK REGRET: v(0) = 0? \n";
119  std::cout << std::right << std::setw(20) << "v(0)" << "\n";
120  std::cout << std::scientific << std::setprecision(11) << std::right
121  << std::setw(20) << std::abs(vx)
122  << "\n";
123  std::cout << "\n";
124  // Check v(x) > x
125  Real scale = two;
126  std::cout << std::right << std::setw(20) << "CHECK REGRET: x < v(x) for |x| > 0? \n";
127  std::cout << std::right << std::setw(20) << "x"
128  << std::right << std::setw(20) << "v(x)"
129  << "\n";
130  for (int i = 0; i < 10; i++) {
131  x = scale*(Real)rand()/(Real)RAND_MAX - scale*half;
132  vx = regret(x,0);
133  std::cout << std::scientific << std::setprecision(11) << std::right
134  << std::setw(20) << x
135  << std::setw(20) << vx
136  << "\n";
137  scale *= two;
138  }
139  std::cout << "\n";
140  // Check v(x) is convex
141  Real y = zero;
142  Real vy = zero;
143  Real z = zero;
144  Real vz = zero;
145  Real l = zero;
146  scale = two;
147  std::cout << std::right << std::setw(20) << "CHECK REGRET: v(x) is convex? \n";
148  std::cout << std::right << std::setw(20) << "v(l*x+(1-l)*y)"
149  << std::setw(20) << "l*v(x)+(1-l)*v(y)"
150  << "\n";
151  for (int i = 0; i < 10; i++) {
152  x = scale*(Real)rand()/(Real)RAND_MAX - scale*half;
153  vx = regret(x,0);
154  y = scale*(Real)rand()/(Real)RAND_MAX - scale*half;
155  vy = regret(y,0);
156  l = (Real)rand()/(Real)RAND_MAX;
157  z = l*x + (one-l)*y;
158  vz = regret(z,0);
159  std::cout << std::scientific << std::setprecision(11) << std::right
160  << std::setw(20) << vz
161  << std::setw(20) << l*vx + (one-l)*vy
162  << "\n";
163  scale *= two;
164  }
165  std::cout << "\n";
166  // Check v'(x)
167  x = oem3*(Real)rand()/(Real)RAND_MAX - fem4;
168  vx = regret(x,0);
169  Real dv = regret(x,1);
170  Real t = one;
171  Real diff = zero;
172  Real err = zero;
173  std::cout << std::right << std::setw(20) << "CHECK REGRET: v'(x) is correct? \n";
174  std::cout << std::right << std::setw(20) << "t"
175  << std::setw(20) << "v'(x)"
176  << std::setw(20) << "(v(x+t)-v(x))/t"
177  << std::setw(20) << "Error"
178  << "\n";
179  for (int i = 0; i < 13; i++) {
180  y = x + t;
181  vy = regret(y,0);
182  diff = (vy-vx)/t;
183  err = std::abs(diff-dv);
184  std::cout << std::scientific << std::setprecision(11) << std::right
185  << std::setw(20) << t
186  << std::setw(20) << dv
187  << std::setw(20) << diff
188  << std::setw(20) << err
189  << "\n";
190  t *= p1;
191  }
192  std::cout << "\n";
193  // Check v''(x)
194  x = oem3*(Real)rand()/(Real)RAND_MAX - fem4;
195  vx = regret(x,1);
196  dv = regret(x,2);
197  t = one;
198  diff = zero;
199  err = zero;
200  std::cout << std::right << std::setw(20) << "CHECK REGRET: v''(x) is correct? \n";
201  std::cout << std::right << std::setw(20) << "t"
202  << std::setw(20) << "v''(x)"
203  << std::setw(20) << "(v'(x+t)-v'(x))/t"
204  << std::setw(20) << "Error"
205  << "\n";
206  for (int i = 0; i < 13; i++) {
207  y = x + t;
208  vy = regret(y,1);
209  diff = (vy-vx)/t;
210  err = std::abs(diff-dv);
211  std::cout << std::scientific << std::setprecision(11) << std::right
212  << std::setw(20) << t
213  << std::setw(20) << dv
214  << std::setw(20) << diff
215  << std::setw(20) << err
216  << "\n";
217  t *= p1;
218  }
219  std::cout << "\n";
220  }
221 
222  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
224  xstat_ = Teuchos::dyn_cast<const RiskVector<Real> >(
225  Teuchos::dyn_cast<const Vector<Real> >(x)).getStatistic(0);
226  if (firstReset_) {
227  dualVector_ = (x0->dual()).clone();
228  firstReset_ = false;
229  }
230  dualVector_->zero();
231  }
232 
233  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
234  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
235  reset(x0,x);
236  v0 = Teuchos::rcp_const_cast<Vector<Real> >(Teuchos::dyn_cast<const RiskVector<Real> >(
237  Teuchos::dyn_cast<const Vector<Real> >(v)).getVector());
238  vstat_ = Teuchos::dyn_cast<const RiskVector<Real> >(
239  Teuchos::dyn_cast<const Vector<Real> >(v)).getStatistic(0);
240  }
241 
242  void update(const Real val, const Real weight) {
243  Real r = regret(val-xstat_,0);
244  RiskMeasure<Real>::val_ += weight * r;
245  }
246 
247  void update(const Real val, const Vector<Real> &g, const Real weight) {
248  Real r = regret(val-xstat_,1);
249  RiskMeasure<Real>::val_ -= weight * r;
250  RiskMeasure<Real>::g_->axpy(weight*r,g);
251  }
252 
253  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
254  const Real weight) {
255  Real r1 = regret(val-xstat_,1);
256  Real r2 = regret(val-xstat_,2);
257  RiskMeasure<Real>::val_ += weight * r2 * (vstat_ - gv);
258  RiskMeasure<Real>::hv_->axpy(weight*r2*(gv-vstat_),g);
259  RiskMeasure<Real>::hv_->axpy(weight*r1,hv);
260  }
261 
263  Real val = RiskMeasure<Real>::val_, gval(0);
264  sampler.sumAll(&val,&gval,1);
265  gval += xstat_;
266  return gval;
267  }
268 
270  RiskVector<Real> &gs = Teuchos::dyn_cast<RiskVector<Real> >(Teuchos::dyn_cast<Vector<Real> >(g));
271  Real stat = RiskMeasure<Real>::val_, gstat(0), one(1);
272  sampler.sumAll(&stat,&gstat,1);
273  gstat += one;
274  gs.setStatistic(gstat);
275 
277  gs.setVector(*dualVector_);
278  }
279 
281  RiskVector<Real> &hs = Teuchos::dyn_cast<RiskVector<Real> >(Teuchos::dyn_cast<Vector<Real> >(hv));
282  Real stat = RiskMeasure<Real>::val_, gstat(0);
283  sampler.sumAll(&stat,&gstat,1);
284  hs.setStatistic(gstat);
285 
287  hs.setVector(*dualVector_);
288  }
289 };
290 
291 }
292 
293 #endif
void update(const Real val, const Vector< Real > &g, const Real weight)
Update internal risk measure storage for gradient computation.
Provides a general interface for risk measures generated through the expectation risk quadrangle...
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
virtual void checkRegret(void)
Run default derivative tests for the scalar regret function.
Real getValue(SampleGenerator< Real > &sampler)
Return risk measure value.
Contains definitions of custom data types in ROL.
void sumAll(Real *input, Real *output, int dim) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< const Vector< Real > > getVector(void) const
void setVector(const Vector< Real > &vec)
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.
void setStatistic(const Real stat)
virtual Real regret(Real x, int deriv=0)=0
Evaluate the scalar regret function at x.
Teuchos::RCP< const StdVector< Real > > getStatistic(void) const
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.
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.
Provides the interface to implement risk measures.
void update(const Real val, const Real weight)
Update internal risk measure storage for value computation.
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
Teuchos::RCP< Vector< Real > > dualVector_