~ubuntu-branches/debian/sid/boost1.49/sid

« back to all changes in this revision

Viewing changes to libs/random/test/test_piecewise_constant.cpp

  • Committer: Package Import Robot
  • Author(s): Steve M. Robbins
  • Date: 2012-02-26 00:31:44 UTC
  • Revision ID: package-import@ubuntu.com-20120226003144-eaytp12cbf6ubpms
Tags: upstream-1.49.0
ImportĀ upstreamĀ versionĀ 1.49.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* test_piecewise_constant.cpp
 
2
 *
 
3
 * Copyright Steven Watanabe 2011
 
4
 * Distributed under the Boost Software License, Version 1.0. (See
 
5
 * accompanying file LICENSE_1_0.txt or copy at
 
6
 * http://www.boost.org/LICENSE_1_0.txt)
 
7
 *
 
8
 * $Id: test_piecewise_constant.cpp 71018 2011-04-05 21:27:52Z steven_watanabe $
 
9
 *
 
10
 */
 
11
 
 
12
#include <boost/random/piecewise_constant_distribution.hpp>
 
13
#include <boost/random/uniform_int.hpp>
 
14
#include <boost/random/mersenne_twister.hpp>
 
15
#include <boost/lexical_cast.hpp>
 
16
#include <boost/exception/diagnostic_information.hpp>
 
17
#include <boost/range/algorithm/lower_bound.hpp>
 
18
#include <boost/range/numeric.hpp>
 
19
#include <vector>
 
20
#include <iostream>
 
21
#include <iomanip>
 
22
 
 
23
#include "statistic_tests.hpp"
 
24
 
 
25
class piecewise_constant
 
26
{
 
27
public:
 
28
    piecewise_constant(const std::vector<double>& intervals, const std::vector<double>& weights)
 
29
      : intervals(intervals),
 
30
        cumulative(1, 0.0)
 
31
    {
 
32
        boost::partial_sum(weights, std::back_inserter(cumulative));
 
33
        for(std::vector<double>::iterator iter = cumulative.begin(), end = cumulative.end();
 
34
            iter != end; ++iter)
 
35
        {
 
36
            *iter /= cumulative.back();
 
37
        }
 
38
    }
 
39
 
 
40
    double cdf(double x) const
 
41
    {
 
42
        std::size_t index = boost::lower_bound(intervals, x) - intervals.begin();
 
43
        if(index == 0) return 0;
 
44
        else if(index == intervals.size()) return 1;
 
45
        else {
 
46
            double lower_weight = cumulative[index - 1];
 
47
            double upper_weight = cumulative[index];
 
48
            double lower = intervals[index - 1];
 
49
            double upper = intervals[index];
 
50
            return lower_weight + (x - lower) / (upper - lower) * (upper_weight - lower_weight);
 
51
        }
 
52
    }
 
53
private:
 
54
    std::vector<double> intervals;
 
55
    std::vector<double> cumulative;
 
56
};
 
57
 
 
58
double cdf(const piecewise_constant& dist, double x)
 
59
{
 
60
    return dist.cdf(x);
 
61
}
 
62
 
 
63
bool do_test(int n, int max) {
 
64
    std::cout << "running piecewise_constant(p0, p1, ..., p" << n-1 << ")" << " " << max << " times: " << std::flush;
 
65
 
 
66
    std::vector<double> weights;
 
67
    {
 
68
        boost::mt19937 egen;
 
69
        for(int i = 0; i < n; ++i) {
 
70
            weights.push_back(egen());
 
71
        }
 
72
    }
 
73
    std::vector<double> intervals;
 
74
    for(int i = 0; i <= n; ++i) {
 
75
        intervals.push_back(i);
 
76
    }
 
77
 
 
78
    piecewise_constant expected(intervals, weights);
 
79
    
 
80
    boost::random::piecewise_constant_distribution<> dist(intervals, weights);
 
81
    boost::mt19937 gen;
 
82
    kolmogorov_experiment test(max);
 
83
    boost::variate_generator<boost::mt19937&, boost::random::piecewise_constant_distribution<> > vgen(gen, dist);
 
84
 
 
85
    double prob = test.probability(test.run(vgen, expected));
 
86
 
 
87
    bool result = prob < 0.99;
 
88
    const char* err = result? "" : "*";
 
89
    std::cout << std::setprecision(17) << prob << err << std::endl;
 
90
 
 
91
    std::cout << std::setprecision(6);
 
92
 
 
93
    return result;
 
94
}
 
95
 
 
96
bool do_tests(int repeat, int max_n, int trials) {
 
97
    boost::mt19937 gen;
 
98
    boost::uniform_int<> idist(1, max_n);
 
99
    int errors = 0;
 
100
    for(int i = 0; i < repeat; ++i) {
 
101
        if(!do_test(idist(gen), trials)) {
 
102
            ++errors;
 
103
        }
 
104
    }
 
105
    if(errors != 0) {
 
106
        std::cout << "*** " << errors << " errors detected ***" << std::endl;
 
107
    }
 
108
    return errors == 0;
 
109
}
 
110
 
 
111
int usage() {
 
112
    std::cerr << "Usage: test_piecewise_constant -r <repeat> -n <max n> -t <trials>" << std::endl;
 
113
    return 2;
 
114
}
 
115
 
 
116
template<class T>
 
117
bool handle_option(int& argc, char**& argv, char opt, T& value) {
 
118
    if(argv[0][1] == opt && argc > 1) {
 
119
        --argc;
 
120
        ++argv;
 
121
        value = boost::lexical_cast<T>(argv[0]);
 
122
        return true;
 
123
    } else {
 
124
        return false;
 
125
    }
 
126
}
 
127
 
 
128
int main(int argc, char** argv) {
 
129
    int repeat = 10;
 
130
    int max_n = 10;
 
131
    int trials = 1000000;
 
132
 
 
133
    if(argc > 0) {
 
134
        --argc;
 
135
        ++argv;
 
136
    }
 
137
    while(argc > 0) {
 
138
        if(argv[0][0] != '-') return usage();
 
139
        else if(!handle_option(argc, argv, 'r', repeat)
 
140
             && !handle_option(argc, argv, 'n', max_n)
 
141
             && !handle_option(argc, argv, 't', trials)) {
 
142
            return usage();
 
143
        }
 
144
        --argc;
 
145
        ++argv;
 
146
    }
 
147
 
 
148
    try {
 
149
        if(do_tests(repeat, max_n, trials)) {
 
150
            return 0;
 
151
        } else {
 
152
            return EXIT_FAILURE;
 
153
        }
 
154
    } catch(...) {
 
155
        std::cerr << boost::current_exception_diagnostic_information() << std::endl;
 
156
        return EXIT_FAILURE;
 
157
    }
 
158
}