ROL
ROL_Gaussian.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_GAUSSIAN_HPP
45 #define ROL_GAUSSIAN_HPP
46 
47 #include "ROL_Distribution.hpp"
48 #include "Teuchos_ParameterList.hpp"
49 
50 namespace ROL {
51 
52 template<class Real>
53 class Gaussian : public Distribution<Real> {
54 private:
55  Real mean_;
56  Real variance_;
57 
58  std::vector<Real> a_;
59  std::vector<Real> b_;
60  std::vector<Real> c_;
61  std::vector<Real> d_;
62 
63  Real erfi(const Real p) const {
64  Real val = 0., z = 0.;
65  if ( p < -0.7 ) {
66  z = std::sqrt(-std::log((1.+p)*0.5));
67  val = -(((c_[3]*z+c_[2])*z+c_[1])*z + c_[0])/((d_[1]*z+d_[0])*z + 1.);
68  }
69  else {
70  if ( p < 0.7 ) {
71  z = p*p;
72  val = p*(((a_[3]*z+a_[2])*z+a_[1])*z + a_[0])/((((b_[3]*z+b_[2])*z+b_[1])*z+b_[0])*z+1.);
73  }
74  else {
75  z = std::sqrt(-std::log((1.-p)*0.5));
76  val = (((c_[3]*z+c_[2])*z+c_[1])*z+c_[0])/((d_[1]*z+d_[0])*z+1.);
77  }
78  }
79  val -= (erf(val)-p)/(2.0/std::sqrt(M_PI) * std::exp(-val*val));
80  val -= (erf(val)-p)/(2.0/std::sqrt(M_PI) * std::exp(-val*val));
81  return val;
82  }
83 
84 public:
85 
86  Gaussian(const Real mean = 0., const Real variance = 1.)
87  : mean_(mean), variance_((variance>0.) ? variance : 1.) {
88  a_.clear(); a_.resize(4,0.); b_.clear(); b_.resize(4,0.);
89  c_.clear(); c_.resize(4,0.); d_.clear(); d_.resize(2,0.);
90  a_[0] = 0.886226899; a_[1] = -1.645349621; a_[2] = 0.914624893; a_[3] = -0.140543331;
91  b_[0] = -2.118377725; b_[1] = 1.442710462; b_[2] = -0.329097515; b_[3] = 0.012229801;
92  c_[0] = -1.970840454; c_[1] = -1.624906493; c_[2] = 3.429567803; c_[3] = 1.641345311;
93  d_[0] = 3.543889200; d_[1] = 1.637067800;
94  }
95 
96  Gaussian(Teuchos::ParameterList &parlist) {
97  mean_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gaussian").get("Mean",0.);
98  variance_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gaussian").get("Variance",1.);
99  variance_ = (variance_ > 0.) ? variance_ : 1.;
100  a_.clear(); a_.resize(4,0.); b_.clear(); b_.resize(4,0.);
101  c_.clear(); c_.resize(4,0.); d_.clear(); d_.resize(2,0.);
102  a_[0] = 0.886226899; a_[1] = -1.645349621; a_[2] = 0.914624893; a_[3] = -0.140543331;
103  b_[0] = -2.118377725; b_[1] = 1.442710462; b_[2] = -0.329097515; b_[3] = 0.012229801;
104  c_[0] = -1.970840454; c_[1] = -1.624906493; c_[2] = 3.429567803; c_[3] = 1.641345311;
105  d_[0] = 3.543889200; d_[1] = 1.637067800;
106  }
107 
108  Real evaluatePDF(const Real input) const {
109  return std::exp(-std::pow(input-mean_,2)/(2.*variance_))/(std::sqrt(2.*M_PI*variance_));
110  }
111 
112  Real evaluateCDF(const Real input) const {
113  return 0.5*(1.+erf((input-mean_)/std::sqrt(2.*variance_)));
114  }
115 
116  Real integrateCDF(const Real input) const {
117  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
118  ">>> ERROR (ROL::Gaussian): Gaussian integrateCDF not implemented!");
119  return ((input < mean_) ? 0.0 : input);
120  }
121 
122  Real invertCDF(const Real input) const {
123  return std::sqrt(2.*variance_)*erfi(2.*input-1.) + mean_;
124  }
125 
126  Real moment(const size_t m) const {
127  Real val = 0.;
128  switch(m) {
129  case 1: val = mean_; break;
130  case 2: val = std::pow(mean_,2) + variance_; break;
131  case 3: val = std::pow(mean_,3)
132  + 3.*mean_*variance_; break;
133  case 4: val = std::pow(mean_,4)
134  + 6.*std::pow(mean_,2)*variance_
135  + 3.*std::pow(variance_,2); break;
136  case 5: val = std::pow(mean_,5)
137  + 10.*std::pow(mean_,3)*variance_
138  + 15.*mean_*std::pow(variance_,2); break;
139  case 6: val = std::pow(mean_,6)
140  + 15.*std::pow(mean_,4)*variance_
141  + 45.*std::pow(mean_*variance_,2)
142  + 15.*std::pow(variance_,3); break;
143  case 7: val = std::pow(mean_,7)
144  + 21.*std::pow(mean_,5)*variance_
145  + 105.*std::pow(mean_,3)*std::pow(variance_,2)
146  + 105.*mean_*std::pow(variance_,3); break;
147  case 8: val = std::pow(mean_,8)
148  + 28.*std::pow(mean_,6)*variance_
149  + 210.*std::pow(mean_,4)*std::pow(variance_,2)
150  + 420.*std::pow(mean_,2)*std::pow(variance_,3)
151  + 105.*std::pow(variance_,4); break;
152  default:
153  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
154  ">>> ERROR (ROL::Distribution): Gaussian moment not implemented for m > 8!");
155  }
156  return val;
157  }
158 
159  Real lowerBound(void) const {
160  return ROL_NINF<Real>();
161  }
162 
163  Real upperBound(void) const {
164  return ROL_INF<Real>();
165  }
166 
167  void test(std::ostream &outStream = std::cout ) const {
168  size_t size = 1;
169  std::vector<Real> X(size,4.*(Real)rand()/(Real)RAND_MAX - 2.);
170  std::vector<int> T(size,0);
171  Distribution<Real>::test(X,T,outStream);
172  }
173 };
174 
175 }
176 
177 #endif
Real evaluateCDF(const Real input) const
Real evaluatePDF(const Real input) const
virtual void test(std::ostream &outStream=std::cout) const
Real upperBound(void) const
std::vector< Real > b_
void test(std::ostream &outStream=std::cout) const
Gaussian(Teuchos::ParameterList &parlist)
Real moment(const size_t m) const
Gaussian(const Real mean=0., const Real variance=1.)
Real lowerBound(void) const
std::vector< Real > a_
Real invertCDF(const Real input) const
std::vector< Real > d_
std::vector< Real > c_
Real integrateCDF(const Real input) const
Real erfi(const Real p) const