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_QUANTILERADIUSQUADRANGLE_HPP
45 #define ROL_QUANTILERADIUSQUADRANGLE_HPP
46 
47 #include "ROL_RandVarFunctional.hpp"
48 #include "ROL_PlusFunction.hpp"
49 
50 #include "ROL_ParameterList.hpp"
51 
52 namespace ROL {
53 
54 template<class Real>
55 class QuantileRadius : public RandVarFunctional<Real> {
56 private:
57   Ptr<PlusFunction<Real> > plusFunction_;
58   Real prob_;
59   Real coeff_;
60   std::vector<Real> vec_;
61 
62   using RandVarFunctional<Real>::val_;
63   using RandVarFunctional<Real>::gv_;
64   using RandVarFunctional<Real>::g_;
65   using RandVarFunctional<Real>::hv_;
66   using RandVarFunctional<Real>::dualVector_;
67 
68   using RandVarFunctional<Real>::point_;
69   using RandVarFunctional<Real>::weight_;
70 
71   using RandVarFunctional<Real>::computeValue;
72   using RandVarFunctional<Real>::computeGradient;
73   using RandVarFunctional<Real>::computeGradVec;
74   using RandVarFunctional<Real>::computeHessVec;
75 
initializeQR(void)76   void initializeQR(void) {
77     Real zero(0);
78     // Initialize temporary storage
79     vec_.clear();  vec_.resize(2,zero);
80   }
81 
checkInputs(void)82   void checkInputs(void) {
83     Real zero(0), one(1);
84     // Check inputs
85     ROL_TEST_FOR_EXCEPTION((prob_>one || prob_<zero), std::invalid_argument,
86       ">>> ERROR (ROL::QuantileRadius): Confidence level out of range!");
87     ROL_TEST_FOR_EXCEPTION((coeff_<zero), std::invalid_argument,
88       ">>> ERROR (ROL::QuantileRadius): Coefficient is negative!");
89      initializeQR();
90   }
91 
92 public:
93 
QuantileRadius(ROL::ParameterList & parlist)94   QuantileRadius( ROL::ParameterList &parlist )
95     : RandVarFunctional<Real>() {
96     ROL::ParameterList &list
97       = parlist.sublist("SOL").sublist("Risk Measure").sublist("Quantile Radius");
98     // Grab probability and coefficient arrays
99     prob_  = list.get<Real>("Confidence Level");
100     coeff_ = list.get<Real>("Coefficient");
101     // Build (approximate) plus function
102     plusFunction_ = makePtr<PlusFunction<Real>>(list);
103     checkInputs();
104   }
105 
QuantileRadius(const Real prob,const Real coeff,const Ptr<PlusFunction<Real>> & pf)106   QuantileRadius(const Real prob, const Real coeff,
107                  const Ptr<PlusFunction<Real> > &pf)
108     : RandVarFunctional<Real>(), plusFunction_(pf), prob_(prob), coeff_(coeff) {
109     checkInputs();
110   }
111 
initialize(const Vector<Real> & x)112   void initialize(const Vector<Real> &x) {
113     RandVarFunctional<Real>::initialize(x);
114     vec_.assign(2,static_cast<Real>(0));
115   }
116 
computeStatistic(const Ptr<std::vector<Real>> & xstat) const117   Real computeStatistic(const Ptr<std::vector<Real>> &xstat) const {
118     Real stat(0), half(0.5);
119     if (xstat != nullPtr) {
120       stat = half*((*xstat)[0] + (*xstat)[1]);
121     }
122     return stat;
123   }
124 
updateValue(Objective<Real> & obj,const Vector<Real> & x,const std::vector<Real> & xstat,Real & tol)125   void updateValue(Objective<Real>         &obj,
126                    const Vector<Real>      &x,
127                    const std::vector<Real> &xstat,
128                    Real                    &tol) {
129     const Real half(0.5), one(1);
130     Real val = computeValue(obj,x,tol);
131     Real pf1 = plusFunction_->evaluate(val-xstat[0],0);
132     Real pf2 = plusFunction_->evaluate(-val-xstat[1],0);
133     RandVarFunctional<Real>::val_ += weight_*(val + half*coeff_/(one-prob_)*(pf1 + pf2));
134   }
135 
getValue(const Vector<Real> & x,const std::vector<Real> & xstat,SampleGenerator<Real> & sampler)136   Real getValue(const Vector<Real>      &x,
137                 const std::vector<Real> &xstat,
138                 SampleGenerator<Real>   &sampler) {
139     const Real half(0.5);
140     Real cvar(0);
141     sampler.sumAll(&val_,&cvar,1);
142     cvar += half*coeff_*(xstat[0] + xstat[1]);
143     return cvar;
144   }
145 
updateGradient(Objective<Real> & obj,const Vector<Real> & x,const std::vector<Real> & xstat,Real & tol)146   void updateGradient(Objective<Real>         &obj,
147                       const Vector<Real>      &x,
148                       const std::vector<Real> &xstat,
149                       Real                    &tol) {
150     const Real half(0.5), one(1);
151     Real val = computeValue(obj,x,tol);
152     Real pf1 = plusFunction_->evaluate(val-xstat[0],1);
153     Real pf2 = plusFunction_->evaluate(-val-xstat[1],1);
154     Real c   = half*weight_*coeff_/(one-prob_);
155     vec_[0] -= c*pf1;
156     vec_[1] -= c*pf2;
157     computeGradient(*dualVector_,obj,x,tol);
158     g_->axpy(weight_ + c * (pf1 - pf2),*dualVector_);
159   }
160 
getGradient(Vector<Real> & g,std::vector<Real> & gstat,const Vector<Real> & x,const std::vector<Real> & xstat,SampleGenerator<Real> & sampler)161   void getGradient(Vector<Real>            &g,
162                    std::vector<Real>       &gstat,
163                    const Vector<Real>      &x,
164                    const std::vector<Real> &xstat,
165                    SampleGenerator<Real>   &sampler) {
166     const Real half(0.5);
167     sampler.sumAll(&vec_[0],&gstat[0],2);
168     sampler.sumAll(*g_,g);
169     gstat[0] += half*coeff_;
170     gstat[1] += half*coeff_;
171   }
172 
updateHessVec(Objective<Real> & obj,const Vector<Real> & v,const std::vector<Real> & vstat,const Vector<Real> & x,const std::vector<Real> & xstat,Real & tol)173   void updateHessVec(Objective<Real>         &obj,
174                      const Vector<Real>      &v,
175                      const std::vector<Real> &vstat,
176                      const Vector<Real>      &x,
177                      const std::vector<Real> &xstat,
178                      Real                    &tol) {
179     const Real half(0.5), one(1);
180     Real val = computeValue(obj,x,tol);
181     Real pf11 = plusFunction_->evaluate(val-xstat[0],1);
182     Real pf12 = plusFunction_->evaluate(val-xstat[0],2);
183     Real pf21 = plusFunction_->evaluate(-val-xstat[1],1);
184     Real pf22 = plusFunction_->evaluate(-val-xstat[1],2);
185     Real c    = half*weight_*coeff_/(one-prob_);
186     Real gv   = computeGradVec(*dualVector_,obj,v,x,tol);
187     vec_[0]  -= c*pf12*(gv-vstat[0]);
188     vec_[1]  += c*pf22*(gv+vstat[1]);
189     hv_->axpy(c*(pf12*(gv-vstat[0]) + pf22*(gv+vstat[1])),*dualVector_);
190     computeHessVec(*dualVector_,obj,v,x,tol);
191     hv_->axpy(weight_ + c * (pf11 - pf21),*dualVector_);
192   }
193 
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)194   void getHessVec(Vector<Real>            &hv,
195                   std::vector<Real>       &hvstat,
196                   const Vector<Real>      &v,
197                   const std::vector<Real> &vstat,
198                   const Vector<Real>      &x,
199                   const std::vector<Real> &xstat,
200                   SampleGenerator<Real>   &sampler) {
201     sampler.sumAll(&vec_[0],&hvstat[0],2);
202     sampler.sumAll(*hv_,hv);
203   }
204 };
205 
206 }
207 
208 #endif
209