ROL
ROL_SpectralRisk.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_SPECTRALRISK_HPP
45 #define ROL_SPECTRALRISK_HPP
46 
49 
86 namespace ROL {
87 
88 template<class Real>
89 class SpectralRisk : public RiskMeasure<Real> {
90 private:
91  Teuchos::RCP<MixedQuantileQuadrangle<Real> > mqq_;
92  Teuchos::RCP<PlusFunction<Real> > plusFunction_;
93 
94  std::vector<Real> wts_;
95  std::vector<Real> pts_;
96 
97  void checkInputs(Teuchos::RCP<Distribution<Real> > &dist = Teuchos::null) const {
98  TEUCHOS_TEST_FOR_EXCEPTION(plusFunction_ == Teuchos::null, std::invalid_argument,
99  ">>> ERROR (ROL::SpectralRisk): PlusFunction pointer is null!");
100  if ( dist != Teuchos::null) {
101  Real lb = dist->lowerBound();
102  Real ub = dist->upperBound();
103  TEUCHOS_TEST_FOR_EXCEPTION(lb < static_cast<Real>(0), std::invalid_argument,
104  ">>> ERROR (ROL::SpectralRisk): Distribution lower bound less than zero!");
105  TEUCHOS_TEST_FOR_EXCEPTION(ub > static_cast<Real>(1), std::invalid_argument,
106  ">>> ERROR (ROL::SpectralRisk): Distribution upper bound greater than one!");
107  }
108  }
109 
110 protected:
111  void buildMixedQuantile(const std::vector<Real> &pts, const std::vector<Real> &wts,
112  const Teuchos::RCP<PlusFunction<Real> > &pf) {
113  pts_.clear(); pts_.assign(pts.begin(),pts.end());
114  wts_.clear(); wts_.assign(wts.begin(),wts.end());
115  plusFunction_ = pf;
116  mqq_ = Teuchos::rcp(new MixedQuantileQuadrangle<Real>(pts,wts,pf));
117  }
118 
119  void buildQuadFromDist(std::vector<Real> &pts, std::vector<Real> &wts,
120  const int nQuad, const Teuchos::RCP<Distribution<Real> > &dist) const {
121  const Real lo = dist->lowerBound(), hi = dist->upperBound();
122  const Real half(0.5), one(1), N(nQuad);
123  wts.clear(); wts.resize(nQuad);
124  pts.clear(); pts.resize(nQuad);
125  if ( hi >= one ) {
126  wts[0] = half/(N-half);
127  pts[0] = lo;
128  for (int i = 1; i < nQuad; ++i) {
129  wts[i] = one/(N-half);
130  pts[i] = dist->invertCDF(static_cast<Real>(i)/N);
131  }
132  }
133  else {
134  wts[0] = half/(N-one);
135  pts[0] = lo;
136  for (int i = 1; i < nQuad-1; ++i) {
137  wts[i] = one/(N-one);
138  pts[i] = dist->invertCDF(static_cast<Real>(i)/N);
139  }
140  wts[nQuad-1] = half/(N-one);
141  pts[nQuad-1] = hi;
142  }
143  }
144 
145  void printQuad(const std::vector<Real> &pts,
146  const std::vector<Real> &wts,
147  const bool print = false) const {
148  if ( print ) {
149  const int nQuad = wts.size();
150  std::cout << std::endl;
151  std::cout << std::scientific << std::setprecision(15);
152  std::cout << std::setw(25) << std::left << "Points"
153  << std::setw(25) << std::left << "Weights"
154  << std::endl;
155  for (int i = 0; i < nQuad; ++i) {
156  std::cout << std::setw(25) << std::left << pts[i]
157  << std::setw(25) << std::left << wts[i]
158  << std::endl;
159  }
160  std::cout << std::endl;
161  }
162  }
163 
164 public:
165  SpectralRisk(void) : RiskMeasure<Real>() {}
166 
167  SpectralRisk( const Teuchos::RCP<Distribution<Real> > &dist,
168  const int nQuad,
169  const Teuchos::RCP<PlusFunction<Real> > &pf)
170  : RiskMeasure<Real>() {
171  // Build generalized trapezoidal rule
172  std::vector<Real> wts(nQuad), pts(nQuad);
173  buildQuadFromDist(pts,wts,nQuad,dist);
174  // Build mixed quantile quadrangle risk measure
175  buildMixedQuantile(pts,wts,pf);
176  // Check inputs
177  checkInputs(dist);
178  }
179 
180  SpectralRisk(Teuchos::ParameterList &parlist)
181  : RiskMeasure<Real>() {
182  // Parse parameter list
183  Teuchos::ParameterList &list
184  = parlist.sublist("SOL").sublist("Risk Measure").sublist("Spectral Risk");
185  int nQuad = list.get("Number of Quadrature Points",5);
186  bool print = list.get("Print Quadrature to Screen",false);
187  // Build distribution
188  Teuchos::RCP<Distribution<Real> > dist = DistributionFactory<Real>(list);
189  // Build plus function approximation
190  Teuchos::RCP<PlusFunction<Real> > pf = Teuchos::rcp(new PlusFunction<Real>(list));
191  // Build generalized trapezoidal rule
192  std::vector<Real> wts(nQuad), pts(nQuad);
193  buildQuadFromDist(pts,wts,nQuad,dist);
194  printQuad(pts,wts,print);
195  // Build mixed quantile quadrangle risk measure
196  buildMixedQuantile(pts,wts,pf);
197  // Check inputs
198  checkInputs(dist);
199  }
200 
201  SpectralRisk( const std::vector<Real> &pts, const std::vector<Real> &wts,
202  const Teuchos::RCP<PlusFunction<Real> > &pf)
203  : RiskMeasure<Real>() {
204  buildMixedQuantile(pts,wts,pf);
205  // Check inputs
206  checkInputs();
207  }
208 
209  Real computeStatistic(const Vector<Real> &x) const {
210  std::vector<Real> xstat;
211  Teuchos::dyn_cast<const RiskVector<Real> >(x).getStatistic(xstat);
212  Real stat(0);
213  int nQuad = static_cast<int>(wts_.size());
214  for (int i = 0; i < nQuad; ++i) {
215  stat += wts_[i] * xstat[i];
216  }
217  return stat;
218  }
219 
220  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
221  mqq_->reset(x0,x);
222  }
223 
224  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
225  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
226  mqq_->reset(x0,x,v0,v);
227  }
228 
229  void update(const Real val, const Real weight) {
230  mqq_->update(val,weight);
231  }
232 
233  void update(const Real val, const Vector<Real> &g, const Real weight) {
234  mqq_->update(val,g,weight);
235  }
236 
237  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
238  const Real weight) {
239  mqq_->update(val,g,gv,hv,weight);
240  }
241 
243  return mqq_->getValue(sampler);
244  }
245 
247  mqq_->getGradient(g,sampler);
248  }
249 
251  mqq_->getHessVec(hv,sampler);
252  }
253 };
254 
255 }
256 
257 #endif
void printQuad(const std::vector< Real > &pts, const std::vector< Real > &wts, const bool print=false) const
Provides an interface for a convex combination of conditional value-at-risks.
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
Reset internal risk measure storage. Called for value and gradient computation.
std::vector< Real > pts_
SpectralRisk(Teuchos::ParameterList &parlist)
Real computeStatistic(const Vector< Real > &x) const
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void buildQuadFromDist(std::vector< Real > &pts, std::vector< Real > &wts, const int nQuad, const Teuchos::RCP< Distribution< Real > > &dist) const
SpectralRisk(const Teuchos::RCP< Distribution< Real > > &dist, const int nQuad, const Teuchos::RCP< PlusFunction< Real > > &pf)
Provides an interface for spectral 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, Teuchos::RCP< Vector< Real > > &v0, const Vector< Real > &v)
Reset internal risk measure storage. Called for Hessian-times-a-vector computation.
Teuchos::RCP< PlusFunction< Real > > plusFunction_
void checkInputs(Teuchos::RCP< Distribution< Real > > &dist=Teuchos::null) const
std::vector< Real > wts_
void buildMixedQuantile(const std::vector< Real > &pts, const std::vector< Real > &wts, const Teuchos::RCP< PlusFunction< Real > > &pf)
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< MixedQuantileQuadrangle< Real > > mqq_
void update(const Real val, const Vector< Real > &g, const Real weight)
Update internal risk measure storage for gradient computation.
SpectralRisk(const std::vector< Real > &pts, const std::vector< Real > &wts, const Teuchos::RCP< PlusFunction< Real > > &pf)
Provides the interface to implement risk measures.
Real getValue(SampleGenerator< Real > &sampler)
Return risk measure value.
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.