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  }
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 
99 protected:
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 
115 private:
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
ROL::Distribution::test_moment
void test_moment(const size_t order, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:304
ROL::Distribution::test_centered
void test_centered(const Real x, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:209
ROL::Distribution::upperBound
virtual Real upperBound(void) const
Definition: ROL_Distribution.hpp:89
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:69
ROL::Distribution::lowerBound
virtual Real lowerBound(void) const
Definition: ROL_Distribution.hpp:84
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:74
ROL::Distribution::test_onesided
void test_onesided(const Real x, std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:116
ROL::Distribution::test
virtual void test(std::ostream &outStream=std::cout) const
Definition: ROL_Distribution.hpp:94
ROL::Distribution::evaluateCDF
virtual Real evaluateCDF(const Real input) const
Definition: ROL_Distribution.hpp:64
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:100
ROL::Distribution::moment
virtual Real moment(const size_t m) const
Definition: ROL_Distribution.hpp:79
Generated by   doxygen 1.8.5