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