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