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
18namespace ROL {
19
20template<class Real>
22public:
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
65protected:
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
81private:
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
Contains definitions of custom data types in ROL.
virtual Real evaluatePDF(const Real input) const
void test_onesided(const Real x, std::ostream &outStream=std::cout) const
virtual Real evaluateCDF(const Real input) const
virtual Real invertCDF(const Real input) const
virtual Real integrateCDF(const Real input) const
virtual void test(std::ostream &outStream=std::cout) const
virtual Real moment(const size_t m) const
virtual Real upperBound(void) const
virtual ~Distribution(void)
void test(const std::vector< Real > &X, const std::vector< int > &T, std::ostream &outStream=std::cout) const
void test_moment(const size_t order, std::ostream &outStream=std::cout) const
void test_centered(const Real x, std::ostream &outStream=std::cout) const
virtual Real lowerBound(void) const