ROL
ROL_Distribution.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_DISTRIBUTION_HPP
45#define ROL_DISTRIBUTION_HPP
46
47#include "ROL_Types.hpp"
48
49#include <cmath>
50#include <iostream>
51
52namespace ROL {
53
54template<class Real>
56public:
57 virtual ~Distribution(void) {}
58
59 virtual Real evaluatePDF(const Real input) const {
60 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
61 ">>> ERROR (ROL::Distribution): evaluatePDF not implemented!");
62 }
63
64 virtual Real evaluateCDF(const Real input) const {
65 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
66 ">>> ERROR (ROL::Distribution): evaluateCDF not implemented!");
67 }
68
69 virtual Real integrateCDF(const Real input) const {
70 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
71 ">>> ERROR (ROL::Distribution): integrateCDF not implemented!");
72 }
73
74 virtual Real invertCDF(const Real input) const {
75 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
76 ">>> ERROR (ROL::Distribution): invertCDF not implemented!");
77 }
78
79 virtual Real moment(const size_t m) const {
80 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
81 ">>> ERROR (ROL::Distribution): moment not implemented!");
82 }
83
84 virtual Real lowerBound(void) const {
85 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
86 ">>> ERROR (ROL::Distribution): lowerBound not implemented!");
87 }
88
89 virtual Real upperBound(void) const {
90 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
91 ">>> ERROR (ROL::Distribution): upperBound not implemented!");
92 }
93
94 virtual void test(std::ostream &outStream = std::cout) const {
95 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
96 ">>> ERROR (ROL::Distribution): test not implemented!");
97 }
98
99protected:
100 void test(const std::vector<Real> &X, const std::vector<int> &T,
101 std::ostream &outStream = std::cout ) const {
102 size_t size = X.size();
103 for ( size_t k = 0; k < size; k++ ) {
104 if ( T[k] == 0 ) {
105 test_onesided(X[k],outStream);
106 }
107 else {
108 test_centered(X[k],outStream);
109 }
110 }
111 size_t order = 2;
112 test_moment(order,outStream);
113 }
114
115private:
116 void test_onesided(const Real x, std::ostream &outStream = std::cout) const {
117 Real X = x, vx = 0., vy = 0., dv = 0., t = 1., diff = 0., err = 0.;
118 try {
119 vx = evaluateCDF(X);
120 vy = 0.0;
121 dv = evaluatePDF(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) = cdf(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 = evaluateCDF(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 evaluatePDF is not implemented!"
149 << std::endl << std::endl;
150 }
151 // CHECK INTCDF
152 try {
153 vx = integrateCDF(X);
154 vy = 0.0;
155 dv = evaluateCDF(X);
156 t = 1.0;
157 diff = 0.0;
158 err = 0.0;
159 outStream << std::scientific << std::setprecision(11);
160 outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = intcdf(x) with x = "
161 << X << " is correct?" << std::endl;
162 outStream << std::right << std::setw(20) << "t"
163 << std::setw(20) << "f'(x)"
164 << std::setw(20) << "(f(x+t)-f(x))/t"
165 << std::setw(20) << "Error"
166 << std::endl;
167 for (int i = 0; i < 13; i++) {
168 vy = integrateCDF(X+t);
169 diff = (vy-vx)/t;
170 err = std::abs(diff-dv);
171 outStream << std::scientific << std::setprecision(11) << std::right
172 << std::setw(20) << t
173 << std::setw(20) << dv
174 << std::setw(20) << diff
175 << std::setw(20) << err
176 << std::endl;
177 t *= 0.1;
178 }
179 outStream << std::endl;
180 }
181 catch(std::exception &e) {
182 outStream << "Either evaluateCDF or integrateCDF is not implemented!"
183 << std::endl << std::endl;
184 }
185 // CHECK INVCDF
186 try {
187 vx = evaluateCDF(X);
188 vy = invertCDF(vx);
189 err = std::abs(x-vy);
190 outStream << std::scientific << std::setprecision(11);
191 outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = invcdf(x) with x = "
192 << X << " is correct?" << std::endl;
193 outStream << std::right << std::setw(20) << "cdf(x)"
194 << std::setw(20) << "invcdf(cdf(x))"
195 << std::setw(20) << "Error"
196 << std::endl;
197 outStream << std::scientific << std::setprecision(11) << std::right
198 << std::setw(20) << vx
199 << std::setw(20) << vy
200 << std::setw(20) << err
201 << std::endl << std::endl;
202 }
203 catch(std::exception &e) {
204 outStream << "Either evaluateCDF or invertCDF is not implemented!"
205 << std::endl << std::endl;
206 }
207 }
208
209 void test_centered(const Real x, std::ostream &outStream = std::cout) const {
210 Real X = x, vx = 0., vy = 0., dv = 0., t = 1., diff = 0., err = 0.;
211 try {
212 vx = 0.0;
213 vy = 0.0;
214 dv = evaluatePDF(X);
215 t = 1.0;
216 diff = 0.0;
217 err = 0.0;
218 outStream << std::scientific << std::setprecision(11);
219 outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = cdf(x) with x = "
220 << X << " is correct?" << std::endl;
221 outStream << std::right << std::setw(20) << "t"
222 << std::setw(20) << "f'(x)"
223 << std::setw(20) << "(f(x+t)-f(x-t))/2t"
224 << std::setw(20) << "Error"
225 << std::endl;
226 for (int i = 0; i < 13; i++) {
227 vx = evaluateCDF(X+t);
228 vy = evaluateCDF(X-t);
229 diff = 0.5*(vx-vy)/t;
230 err = std::abs(diff-dv);
231 outStream << std::scientific << std::setprecision(11) << std::right
232 << std::setw(20) << t
233 << std::setw(20) << dv
234 << std::setw(20) << diff
235 << std::setw(20) << err
236 << std::endl;
237 t *= 0.1;
238 }
239 outStream << "\n";
240 }
241 catch(std::exception &e) {
242 outStream << "Either evaluateCDF or evaluatePDF is not implemented!"
243 << std::endl << std::endl;
244 }
245 // CHECK INTCDF
246 try {
247 vx = 0.0;
248 vy = 0.0;
249 dv = evaluateCDF(X);
250 t = 1.0;
251 diff = 0.0;
252 err = 0.0;
253 outStream << std::scientific << std::setprecision(11);
254 outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = intcdf(x) with x = "
255 << X << " is correct?" << std::endl;
256 outStream << std::right << std::setw(20) << "t"
257 << std::setw(20) << "f'(x)"
258 << std::setw(20) << "(f(x+t)-f(x-t))/2t"
259 << std::setw(20) << "Error"
260 << std::endl;
261 for (int i = 0; i < 13; i++) {
262 vx = integrateCDF(X+t);
263 vy = integrateCDF(X-t);
264 diff = 0.5*(vx-vy)/t;
265 err = std::abs(diff-dv);
266 outStream << std::scientific << std::setprecision(11) << std::right
267 << std::setw(20) << t
268 << std::setw(20) << dv
269 << std::setw(20) << diff
270 << std::setw(20) << err
271 << std::endl;
272 t *= 0.1;
273 }
274 outStream << std::endl;
275 }
276 catch(std::exception &e) {
277 outStream << "Either evaluateCDF or integrateCDF is not implemented!"
278 << std::endl << std::endl;
279 }
280 // CHECK INVCDF
281 try {
282 vx = evaluateCDF(X);
283 vy = invertCDF(vx);
284 err = std::abs(X-vy);
285 outStream << std::scientific << std::setprecision(11);
286 outStream << std::right << std::setw(20) << "CHECK DENSITY: f(x) = invcdf(x) with x = "
287 << X << " is correct?" << std::endl;
288 outStream << std::right << std::setw(20) << "cdf(x)"
289 << std::setw(20) << "invcdf(cdf(x))"
290 << std::setw(20) << "Error"
291 << std::endl;
292 outStream << std::scientific << std::setprecision(11) << std::right
293 << std::setw(20) << vx
294 << std::setw(20) << vy
295 << std::setw(20) << err
296 << std::endl << std::endl;
297 }
298 catch(std::exception &e) {
299 outStream << "Either evaluateCDF or invertCDF is not implemented!"
300 << std::endl << std::endl;
301 }
302 }
303
304 void test_moment(const size_t order, std::ostream &outStream = std::cout) const {
305 try {
306 const size_t numPts = 10000;
307 Real pt = 0., wt = 1./(Real)numPts;
308 std::vector<Real> mVec(order,0.);
309 for (size_t i = 0; i < numPts; i++) {
310 pt = invertCDF((Real)rand()/(Real)RAND_MAX);
311 mVec[0] += wt*pt;
312 for (size_t q = 1; q < order; q++) {
313 mVec[q] += wt*std::pow(pt,q+1);
314 }
315 }
316 outStream << std::scientific << std::setprecision(0);
317 outStream << std::right << std::setw(20) << "CHECK DENSITY: Check first " << order
318 << " moments against Monte Carlo using " << numPts << " samples"
319 << std::endl;
320 outStream << std::setw(20) << "Error should be O(" << 1./std::sqrt(numPts) << ")" << std::endl;
321 outStream << std::scientific << std::setprecision(11);
322 for (size_t q = 0; q < order; q++) {
323 outStream << std::setw(20) << "Error in " << q+1 << " moment: "
324 << std::abs(mVec[q]-moment(q+1)) << std::endl;
325 }
326 outStream << std::endl;
327 }
328 catch(std::exception &e) {
329 outStream << "moment is not implemented!"
330 << std::endl << std::endl;
331 }
332 }
333};
334
335}
336
337#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