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