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_BETA_HPP 45 #define ROL_BETA_HPP 46 47 #include "ROL_Distribution.hpp" 48 #include "ROL_ParameterList.hpp" 49 50 #include <math.h> 51 52 namespace ROL { 53 54 template<class Real> 55 class Beta : public Distribution<Real> { 56 private: 57 Real shape1_; 58 Real shape2_; 59 Real coeff_; 60 61 std::vector<Real> pts_; 62 std::vector<Real> wts_; 63 initializeQuadrature(void)64 void initializeQuadrature(void) { 65 pts_.clear(); pts_.resize(20,0.); wts_.clear(); wts_.resize(20,0.); 66 wts_[0] = 0.1527533871307258; pts_[0] = -0.0765265211334973; 67 wts_[1] = 0.1527533871307258; pts_[1] = 0.0765265211334973; 68 wts_[2] = 0.1491729864726037; pts_[2] = -0.2277858511416451; 69 wts_[3] = 0.1491729864726037; pts_[3] = 0.2277858511416451; 70 wts_[4] = 0.1420961093183820; pts_[4] = -0.3737060887154195; 71 wts_[5] = 0.1420961093183820; pts_[5] = 0.3737060887154195; 72 wts_[6] = 0.1316886384491766; pts_[6] = -0.5108670019508271; 73 wts_[7] = 0.1316886384491766; pts_[7] = 0.5108670019508271; 74 wts_[8] = 0.1181945319615184; pts_[8] = -0.6360536807265150; 75 wts_[9] = 0.1181945319615184; pts_[9] = 0.6360536807265150; 76 wts_[10] = 0.1019301198172404; pts_[10] = -0.7463319064601508; 77 wts_[11] = 0.1019301198172404; pts_[11] = 0.7463319064601508; 78 wts_[12] = 0.0832767415767048; pts_[12] = -0.8391169718222188; 79 wts_[13] = 0.0832767415767048; pts_[13] = 0.8391169718222188; 80 wts_[14] = 0.0626720483341091; pts_[14] = -0.9122344282513259; 81 wts_[15] = 0.0626720483341091; pts_[15] = 0.9122344282513259; 82 wts_[16] = 0.0406014298003869; pts_[16] = -0.9639719272779138; 83 wts_[17] = 0.0406014298003869; pts_[17] = 0.9639719272779138; 84 wts_[18] = 0.0176140071391521; pts_[18] = -0.9931285991850949; 85 wts_[19] = 0.0176140071391521; pts_[19] = 0.9931285991850949; 86 for (size_t i = 0; i < 20; i++) { 87 wts_[i] *= 0.5; 88 pts_[i] += 1.; pts_[i] *= 0.5; 89 } 90 } 91 ibeta(const Real x) const92 Real ibeta(const Real x) const { 93 Real pt = 0., wt = 0., sum = 0.; 94 for (size_t i = 0; i < pts_.size(); i++) { 95 wt = x*wts_[i]; 96 pt = x*pts_[i]; 97 sum += wt*std::pow(pt,shape1_-1)*std::pow(1.-pt,shape2_-1); 98 } 99 return sum; 100 } 101 102 public: Beta(const Real shape1=2.,const Real shape2=2.)103 Beta(const Real shape1 = 2., const Real shape2 = 2.) 104 : shape1_((shape1 > 0.) ? shape1 : 2.), shape2_((shape2 > 0.) ? shape2 : 2.) { 105 coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_)); 106 initializeQuadrature(); 107 } 108 Beta(ROL::ParameterList & parlist)109 Beta(ROL::ParameterList &parlist) { 110 shape1_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 1",2.); 111 shape2_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 2",2.); 112 shape1_ = (shape1_ > 0.) ? shape1_ : 2.; 113 shape2_ = (shape2_ > 0.) ? shape2_ : 2.; 114 coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_)); 115 initializeQuadrature(); 116 } 117 evaluatePDF(const Real input) const118 Real evaluatePDF(const Real input) const { 119 return ((input > 0.) ? ((input < 1.) ? 120 coeff_*std::pow(input,shape1_-1.)*std::pow(1.-input,shape2_-1) : 0.) : 0.); 121 } 122 evaluateCDF(const Real input) const123 Real evaluateCDF(const Real input) const { 124 return ((input > 0.) ? ((input < 1.) ? coeff_*ibeta(input) : 1.) : 0.); 125 } 126 integrateCDF(const Real input) const127 Real integrateCDF(const Real input) const { 128 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument, 129 ">>> ERROR (ROL::Beta): Beta integrateCDF not implemented!"); 130 } 131 invertCDF(const Real input) const132 Real invertCDF(const Real input) const { 133 if ( input <= ROL_EPSILON<Real>() ) { 134 return 0.; 135 } 136 if ( input >= 1.-ROL_EPSILON<Real>() ) { 137 return 1.; 138 } 139 Real a = ROL_EPSILON<Real>(), b = 1.-ROL_EPSILON<Real>(), c = 0.; 140 Real fa = evaluateCDF(a) - input; 141 Real fc = 0.; 142 Real sa = ((fa < 0.) ? -1. : ((fa > 0.) ? 1. : 0.)); 143 Real sc = 0.; 144 for (size_t i = 0; i < 100; i++) { 145 c = (a+b)*0.5; 146 fc = evaluateCDF(c) - input; 147 sc = ((fc < 0.) ? -1. : ((fc > 0.) ? 1. : 0.)); 148 if ( fc == 0. || (b-a)*0.5 < ROL_EPSILON<Real>() ) { 149 break; 150 } 151 if ( sc == sa ) { a = c; fa = fc; sa = sc; } 152 else { b = c; } 153 } 154 return c; 155 } 156 moment(const size_t m) const157 Real moment(const size_t m) const { 158 if ( m == 1 ) { 159 return shape1_/(shape1_ + shape2_); 160 } 161 if ( m == 2 ) { 162 return shape1_*(shape2_/(shape1_+shape2_+1.) + shape1_)/std::pow(shape1_+shape2_,2); 163 } 164 Real val = 1.; 165 for (size_t i = 0; i < m; i++) { 166 val *= (shape1_ + (Real)i)/(shape1_ + shape2_ + (Real)i); 167 } 168 return val; 169 } 170 lowerBound(void) const171 Real lowerBound(void) const { 172 return 0.; 173 } 174 upperBound(void) const175 Real upperBound(void) const { 176 return 1.; 177 } 178 test(std::ostream & outStream=std::cout) const179 void test(std::ostream &outStream = std::cout ) const { 180 size_t size = 5; 181 std::vector<Real> X(size,0.); 182 std::vector<int> T(size,0); 183 X[0] = -4.*(Real)rand()/(Real)RAND_MAX; 184 T[0] = 0; 185 X[1] = 0.; 186 T[1] = 1; 187 X[2] = 0.5*(Real)rand()/(Real)RAND_MAX; 188 T[2] = 0; 189 X[3] = 1.; 190 T[3] = 1; 191 X[4] = 1.+4.*(Real)rand()/(Real)RAND_MAX; 192 T[4] = 0; 193 Distribution<Real>::test(X,T,outStream); 194 } 195 }; 196 197 } 198 199 #endif 200