ROL
ROL_Distribution.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_DISTRIBUTION_HPP
11 #define ROL_DISTRIBUTION_HPP
12 
13 #include "ROL_Types.hpp"
14 
15 #include <cmath>
16 #include <iostream>
17 
18 namespace ROL {
19 
20 template<class Real>
21 class Distribution {
22 public:
23  virtual ~Distribution(void) {}
24 
25  virtual Real evaluatePDF(const Real input) const {
26  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
27  ">>> ERROR (ROL::Distribution): evaluatePDF not implemented!");
28  }
29 
30  virtual Real evaluateCDF(const Real input) const {
31  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
32  ">>> ERROR (ROL::Distribution): evaluateCDF not implemented!");
33  }
34 
35  virtual Real integrateCDF(const Real input) const {
36  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
37  ">>> ERROR (ROL::Distribution): integrateCDF not implemented!");
38  }
39 
40  virtual Real invertCDF(const Real input) const {
41  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
42  ">>> ERROR (ROL::Distribution): invertCDF not implemented!");
43  }
44 
45  virtual Real moment(const size_t m) const {
46  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
47  ">>> ERROR (ROL::Distribution): moment not implemented!");
48  }
49 
50  virtual Real lowerBound(void) const {
51  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
52  ">>> ERROR (ROL::Distribution): lowerBound not implemented!");
53  }
54 
55  virtual Real upperBound(void) const {
56  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
57  ">>> ERROR (ROL::Distribution): upperBound not implemented!");
58  }
59 
60  virtual void test(std::ostream &outStream = std::cout) const {
61  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
62  ">>> ERROR (ROL::Distribution): test not implemented!");
63  }
64 
65 protected:
66  void test(const std::vector<Real> &X, const std::vector<int> &T,
67  std::ostream &outStream = std::cout ) const {
68  size_t size = X.size();
69  for ( size_t k = 0; k < size; k++ ) {
70  if ( T[k] == 0 ) {
71  test_onesided(X[k],outStream);
72  }
73  else {
74  test_centered(X[k],outStream);
75  }
76  }
77  size_t order = 2;
78  test_moment(order,outStream);
79  }
80 
81 private:
82  void test_onesided(const Real x, std::ostream &outStream = std::cout) const {
83  Real X = x, vx = 0., vy = 0., dv = 0., t = 1., diff = 0., err = 0.;
84  try {
85  vx = evaluateCDF(X);
86  vy = 0.0;
87  dv = evaluatePDF(X);
88  t = 1.0;
89  diff = 0.0;
90  err = 0.0;
91  outStream << std::scientific << std::setprecision(11);
92  outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = cdf(x) with x = "
93  << X << " is correct?" << std::endl;
94  outStream << std::right << std::setw(20) << "t"
95  << std::setw(20) << "f'(x)"
96  << std::setw(20) << "(f(x+t)-f(x))/t"
97  << std::setw(20) << "Error"
98  << std::endl;
99  for (int i = 0; i < 13; i++) {
100  vy = evaluateCDF(X+t);
101  diff = (vy-vx)/t;
102  err = std::abs(diff-dv);
103  outStream << std::scientific << std::setprecision(11) << std::right
104  << std::setw(20) << t
105  << std::setw(20) << dv
106  << std::setw(20) << diff
107  << std::setw(20) << err
108  << std::endl;
109  t *= 0.1;
110  }
111  outStream << std::endl;
112  }
113  catch(std::exception &e) {
114  outStream << "Either evaluateCDF or evaluatePDF is not implemented!"
115  << std::endl << std::endl;
116  }
117  // CHECK INTCDF
118  try {
119  vx = integrateCDF(X);
120  vy = 0.0;
121  dv = evaluateCDF(X);
122  t = 1.0;
123  diff = 0.0;
124  err = 0.0;
125  outStream << std::scientific << std::setprecision(11);
126  outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = intcdf(x) with x = "
127  << X << " is correct?" << std::endl;
128  outStream << std::right << std::setw(20) << "t"
129  << std::setw(20) << "f'(x)"
130  << std::setw(20) << "(f(x+t)-f(x))/t"
131  << std::setw(20) << "Error"
132  << std::endl;
133  for (int i = 0; i < 13; i++) {
134  vy = integrateCDF(X+t);
135  diff = (vy-vx)/t;
136  err = std::abs(diff-dv);
137  outStream << std::scientific << std::setprecision(11) << std::right
138  << std::setw(20) << t
139  << std::setw(20) << dv
140  << std::setw(20) << diff
141  << std::setw(20) << err
142  << std::endl;
143  t *= 0.1;
144  }
145  outStream << std::endl;
146  }
147  catch(std::exception &e) {
148  outStream << "Either evaluateCDF or integrateCDF is not implemented!"
149  << std::endl << std::endl;
150  }
151  // CHECK INVCDF
152  try {
153  vx = evaluateCDF(X);
154  vy = invertCDF(vx);
155  err = std::abs(x-vy);
156  outStream << std::scientific << std::setprecision(11);
157  outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = invcdf(x) with x = "
158  << X << " is correct?" << std::endl;
159  outStream << std::right << std::setw(20) << "cdf(x)"
160  << std::setw(20) << "invcdf(cdf(x))"
161  << std::setw(20) << "Error"
162  << std::endl;
163  outStream << std::scientific << std::setprecision(11) << std::right
164  << std::setw(20) << vx
165  << std::setw(20) << vy
166  << std::setw(20) << err
167  << std::endl << std::endl;
168  }
169  catch(std::exception &e) {
170  outStream << "Either evaluateCDF or invertCDF is not implemented!"
171  << std::endl << std::endl;
172  }
173  }
174 
175  void test_centered(const Real x, std::ostream &outStream = std::cout) const {
176  Real X = x, vx = 0., vy = 0., dv = 0., t = 1., diff = 0., err = 0.;
177  try {
178  vx = 0.0;
179  vy = 0.0;
180  dv = evaluatePDF(X);
181  t = 1.0;
182  diff = 0.0;
183  err = 0.0;
184  outStream << std::scientific << std::setprecision(11);
185  outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = cdf(x) with x = "
186  << X << " is correct?" << std::endl;
187  outStream << std::right << std::setw(20) << "t"
188  << std::setw(20) << "f'(x)"
189  << std::setw(20) << "(f(x+t)-f(x-t))/2t"
190  << std::setw(20) << "Error"
191  << std::endl;
192  for (int i = 0; i < 13; i++) {
193  vx = evaluateCDF(X+t);
194  vy = evaluateCDF(X-t);
195  diff = 0.5*(vx-vy)/t;
196  err = std::abs(diff-dv);
197  outStream << std::scientific << std::setprecision(11) << std::right
198  << std::setw(20) << t
199  << std::setw(20) << dv
200  << std::setw(20) << diff
201  << std::setw(20) << err
202  << std::endl;
203  t *= 0.1;
204  }
205  outStream << "\n";
206  }
207  catch(std::exception &e) {
208  outStream << "Either evaluateCDF or evaluatePDF is not implemented!"
209  << std::endl << std::endl;
210  }
211  // CHECK INTCDF
212  try {
213  vx = 0.0;
214  vy = 0.0;
215  dv = evaluateCDF(X);
216  t = 1.0;
217  diff = 0.0;
218  err = 0.0;
219  outStream << std::scientific << std::setprecision(11);
220  outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = intcdf(x) with x = "
221  << X << " is correct?" << std::endl;
222  outStream << std::right << std::setw(20) << "t"
223  << std::setw(20) << "f'(x)"
224  << std::setw(20) << "(f(x+t)-f(x-t))/2t"
225  << std::setw(20) << "Error"
226  << std::endl;
227  for (int i = 0; i < 13; i++) {
228  vx = integrateCDF(X+t);
229  vy = integrateCDF(X-t);
230  diff = 0.5*(vx-vy)/t;
231  err = std::abs(diff-dv);
232  outStream << std::scientific << std::setprecision(11) << std::right
233  << std::setw(20) << t
234  << std::setw(20) << dv
235  << std::setw(20) << diff
236  << std::setw(20) << err
237  << std::endl;
238  t *= 0.1;
239  }
240  outStream << std::endl;
241  }
242  catch(std::exception &e) {
243  outStream << "Either evaluateCDF or integrateCDF is not implemented!"
244  << std::endl << std::endl;
245  }
246  // CHECK INVCDF
247  try {
248  vx = evaluateCDF(X);
249  vy = invertCDF(vx);
250  err = std::abs(X-vy);
251  outStream << std::scientific << std::setprecision(11);
252  outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = invcdf(x) with x = "
253  << X << " is correct?" << std::endl;
254  outStream << std::right << std::setw(20) << "cdf(x)"
255  << std::setw(20) << "invcdf(cdf(x))"
256  << std::setw(20) << "Error"
257  << std::endl;
258  outStream << std::scientific << std::setprecision(11) << std::right
259  << std::setw(20) << vx
260  << std::setw(20) << vy
261  << std::setw(20) << err
262  << std::endl << std::endl;
263  }
264  catch(std::exception &e) {
265  outStream << "Either evaluateCDF or invertCDF is not implemented!"
266  << std::endl << std::endl;
267  }
268  }
269 
270  void test_moment(const size_t order, std::ostream &outStream = std::cout) const {
271  try {
272  const size_t numPts = 10000;
273  Real pt = 0., wt = 1./(Real)numPts;
274  std::vector<Real> mVec(order,0.);
275  for (size_t i = 0; i < numPts; i++) {
276  pt = invertCDF((Real)rand()/(Real)RAND_MAX);
277  mVec[0] += wt*pt;
278  for (size_t q = 1; q < order; q++) {
279  mVec[q] += wt*std::pow(pt,q+1);
280  }
281  }
282  outStream << std::scientific << std::setprecision(0);
283  outStream << std::right << std::setw(20) << "CHECK DENSITY: Check first " << order
284  << " moments against Monte Carlo using " << numPts << " samples"
285  << std::endl;
286  outStream << std::setw(20) << "Error should be O(" << 1./std::sqrt(numPts) << ")" << std::endl;
287  outStream << std::scientific << std::setprecision(11);
288  for (size_t q = 0; q < order; q++) {
289  outStream << std::setw(20) << "Error in " << q+1 << " moment: "
290  << std::abs(mVec[q]-moment(q+1)) << std::endl;
291  }
292  outStream << std::endl;
293  }
294  catch(std::exception &e) {
295  outStream << "moment is not implemented!"
296  << std::endl << std::endl;
297  }
298  }
299 };
300 
301 }
302 
303 #endif
ROL::Distribution::test_moment
void test_moment(const size_t order, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:270
ROL::Distribution::test_centered
void test_centered(const Real x, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:175
ROL::Distribution::upperBound
virtual Real upperBound(void) const
Definition: ROL_Distribution.hpp:55
ROL::Distribution::evaluatePDF
virtual Real evaluatePDF(const Real input) const
Definition: ROL_Distribution.hpp:25
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::Distribution::integrateCDF
virtual Real integrateCDF(const Real input) const
Definition: ROL_Distribution.hpp:35
ROL::Distribution::lowerBound
virtual Real lowerBound(void) const
Definition: ROL_Distribution.hpp:50
ROL::Distribution::~Distribution
virtual ~Distribution(void)
Definition: ROL_Distribution.hpp:23
ROL::Distribution::invertCDF
virtual Real invertCDF(const Real input) const
Definition: ROL_Distribution.hpp:40
ROL::Distribution::test_onesided
void test_onesided(const Real x, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:82
ROL::Distribution::test
virtual void test(std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:60
ROL::Distribution::evaluateCDF
virtual Real evaluateCDF(const Real input) const
Definition: ROL_Distribution.hpp:30
ROL::Distribution::test
void test(const std::vector< Real > &X, const std::vector< int > &T, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:66
ROL::Distribution::moment
virtual Real moment(const size_t m) const
Definition: ROL_Distribution.hpp:45
Generated by   doxygen 1.8.5