| 1 | #! parrot |
|---|
| 2 | # Copyright (C) 2010, Parrot Foundation. |
|---|
| 3 | # $Id: $ |
|---|
| 4 | |
|---|
| 5 | =head1 NAME |
|---|
| 6 | |
|---|
| 7 | t/dynoplibs/random-range.t - Tests random range behavior of math_ops lib. |
|---|
| 8 | |
|---|
| 9 | =head1 SYNOPSIS |
|---|
| 10 | |
|---|
| 11 | % prove t/dynoblibs/random-range.t |
|---|
| 12 | |
|---|
| 13 | =head1 DESCRIPTION |
|---|
| 14 | |
|---|
| 15 | Tests random range using chi-square |
|---|
| 16 | |
|---|
| 17 | =cut |
|---|
| 18 | |
|---|
| 19 | .loadlib 'math_ops' |
|---|
| 20 | .sub add_chi2_entry |
|---|
| 21 | .param pmc table |
|---|
| 22 | .param int nu |
|---|
| 23 | .param num p_10 |
|---|
| 24 | .param num p_05 |
|---|
| 25 | .param num p_025 |
|---|
| 26 | .param num p_01 |
|---|
| 27 | .param num p_001 |
|---|
| 28 | |
|---|
| 29 | $P0 = new ['FixedFloatArray'] |
|---|
| 30 | $P0 = 5 |
|---|
| 31 | |
|---|
| 32 | $P0[0] = p_10 |
|---|
| 33 | $P0[1] = p_05 |
|---|
| 34 | $P0[2] = p_025 |
|---|
| 35 | $P0[3] = p_01 |
|---|
| 36 | $P0[4] = p_001 |
|---|
| 37 | |
|---|
| 38 | table[nu] = $P0 |
|---|
| 39 | .end |
|---|
| 40 | |
|---|
| 41 | .sub compute_chi_square |
|---|
| 42 | .param pmc histogram |
|---|
| 43 | .param int num_samples |
|---|
| 44 | |
|---|
| 45 | .local int possible_values |
|---|
| 46 | possible_values = elements histogram |
|---|
| 47 | |
|---|
| 48 | .local num expected |
|---|
| 49 | expected = num_samples |
|---|
| 50 | expected /= possible_values |
|---|
| 51 | |
|---|
| 52 | .local num K |
|---|
| 53 | K = 0.0 |
|---|
| 54 | |
|---|
| 55 | $I0 = possible_values |
|---|
| 56 | if $I0 == 0 goto end |
|---|
| 57 | chi2_loop: |
|---|
| 58 | dec $I0 |
|---|
| 59 | |
|---|
| 60 | # Compute: (O - E)^2 / E |
|---|
| 61 | $N0 = histogram[$I0] |
|---|
| 62 | $N0 -= expected |
|---|
| 63 | $N0 *= $N0 |
|---|
| 64 | $N0 /= expected |
|---|
| 65 | |
|---|
| 66 | K += $N0 |
|---|
| 67 | |
|---|
| 68 | unless $I0 == 0 goto chi2_loop |
|---|
| 69 | |
|---|
| 70 | #~ print "Chi-squared K is: " |
|---|
| 71 | #~ say K |
|---|
| 72 | end: |
|---|
| 73 | .return (K) |
|---|
| 74 | .end |
|---|
| 75 | |
|---|
| 76 | .sub make_chi2_table |
|---|
| 77 | # This info comes from http://itl.nist.gov/div898/handbook/eda/section3/eda3674.htm |
|---|
| 78 | # Obtained 26 Feb 2010 -- Austin |
|---|
| 79 | |
|---|
| 80 | # Probability of exceeding the critical value |
|---|
| 81 | # nu 0.10 0.05 0.025 0.01 0.001 |
|---|
| 82 | |
|---|
| 83 | .local pmc table |
|---|
| 84 | table = new ['ResizablePMCArray'] |
|---|
| 85 | |
|---|
| 86 | add_chi2_entry(table, 1, 2.706, 3.841, 5.024, 6.635, 10.828) |
|---|
| 87 | add_chi2_entry(table, 2, 4.605, 5.991, 7.378, 9.210, 13.816) |
|---|
| 88 | add_chi2_entry(table, 3, 6.251, 7.815, 9.348, 11.345, 16.266) |
|---|
| 89 | add_chi2_entry(table, 4, 7.779, 9.488, 11.143, 13.277, 18.467) |
|---|
| 90 | add_chi2_entry(table, 5, 9.236, 11.070, 12.833, 15.086, 20.515) |
|---|
| 91 | add_chi2_entry(table, 6, 10.645, 12.592, 14.449, 16.812, 22.458) |
|---|
| 92 | add_chi2_entry(table, 7, 12.017, 14.067, 16.013, 18.475, 24.322) |
|---|
| 93 | add_chi2_entry(table, 8, 13.362, 15.507, 17.535, 20.090, 26.125) |
|---|
| 94 | add_chi2_entry(table, 9, 14.684, 16.919, 19.023, 21.666, 27.877) |
|---|
| 95 | add_chi2_entry(table, 10, 15.987, 18.307, 20.483, 23.209, 29.588) |
|---|
| 96 | add_chi2_entry(table, 11, 17.275, 19.675, 21.920, 24.725, 31.264) |
|---|
| 97 | add_chi2_entry(table, 12, 18.549, 21.026, 23.337, 26.217, 32.910) |
|---|
| 98 | add_chi2_entry(table, 13, 19.812, 22.362, 24.736, 27.688, 34.528) |
|---|
| 99 | add_chi2_entry(table, 14, 21.064, 23.685, 26.119, 29.141, 36.123) |
|---|
| 100 | add_chi2_entry(table, 15, 22.307, 24.996, 27.488, 30.578, 37.697) |
|---|
| 101 | add_chi2_entry(table, 16, 23.542, 26.296, 28.845, 32.000, 39.252) |
|---|
| 102 | add_chi2_entry(table, 17, 24.769, 27.587, 30.191, 33.409, 40.790) |
|---|
| 103 | add_chi2_entry(table, 18, 25.989, 28.869, 31.526, 34.805, 42.312) |
|---|
| 104 | add_chi2_entry(table, 19, 27.204, 30.144, 32.852, 36.191, 43.820) |
|---|
| 105 | add_chi2_entry(table, 20, 28.412, 31.410, 34.170, 37.566, 45.315) |
|---|
| 106 | add_chi2_entry(table, 21, 29.615, 32.671, 35.479, 38.932, 46.797) |
|---|
| 107 | add_chi2_entry(table, 22, 30.813, 33.924, 36.781, 40.289, 48.268) |
|---|
| 108 | add_chi2_entry(table, 23, 32.007, 35.172, 38.076, 41.638, 49.728) |
|---|
| 109 | add_chi2_entry(table, 24, 33.196, 36.415, 39.364, 42.980, 51.179) |
|---|
| 110 | add_chi2_entry(table, 25, 34.382, 37.652, 40.646, 44.314, 52.620) |
|---|
| 111 | add_chi2_entry(table, 26, 35.563, 38.885, 41.923, 45.642, 54.052) |
|---|
| 112 | add_chi2_entry(table, 27, 36.741, 40.113, 43.195, 46.963, 55.476) |
|---|
| 113 | add_chi2_entry(table, 28, 37.916, 41.337, 44.461, 48.278, 56.892) |
|---|
| 114 | add_chi2_entry(table, 29, 39.087, 42.557, 45.722, 49.588, 58.301) |
|---|
| 115 | add_chi2_entry(table, 30, 40.256, 43.773, 46.979, 50.892, 59.703) |
|---|
| 116 | add_chi2_entry(table, 31, 41.422, 44.985, 48.232, 52.191, 61.098) |
|---|
| 117 | add_chi2_entry(table, 32, 42.585, 46.194, 49.480, 53.486, 62.487) |
|---|
| 118 | add_chi2_entry(table, 33, 43.745, 47.400, 50.725, 54.776, 63.870) |
|---|
| 119 | add_chi2_entry(table, 34, 44.903, 48.602, 51.966, 56.061, 65.247) |
|---|
| 120 | add_chi2_entry(table, 35, 46.059, 49.802, 53.203, 57.342, 66.619) |
|---|
| 121 | add_chi2_entry(table, 36, 47.212, 50.998, 54.437, 58.619, 67.985) |
|---|
| 122 | add_chi2_entry(table, 37, 48.363, 52.192, 55.668, 59.893, 69.347) |
|---|
| 123 | add_chi2_entry(table, 38, 49.513, 53.384, 56.896, 61.162, 70.703) |
|---|
| 124 | add_chi2_entry(table, 39, 50.660, 54.572, 58.120, 62.428, 72.055) |
|---|
| 125 | add_chi2_entry(table, 40, 51.805, 55.758, 59.342, 63.691, 73.402) |
|---|
| 126 | add_chi2_entry(table, 41, 52.949, 56.942, 60.561, 64.950, 74.745) |
|---|
| 127 | add_chi2_entry(table, 42, 54.090, 58.124, 61.777, 66.206, 76.084) |
|---|
| 128 | add_chi2_entry(table, 43, 55.230, 59.304, 62.990, 67.459, 77.419) |
|---|
| 129 | add_chi2_entry(table, 44, 56.369, 60.481, 64.201, 68.710, 78.750) |
|---|
| 130 | add_chi2_entry(table, 45, 57.505, 61.656, 65.410, 69.957, 80.077) |
|---|
| 131 | add_chi2_entry(table, 46, 58.641, 62.830, 66.617, 71.201, 81.400) |
|---|
| 132 | add_chi2_entry(table, 47, 59.774, 64.001, 67.821, 72.443, 82.720) |
|---|
| 133 | add_chi2_entry(table, 48, 60.907, 65.171, 69.023, 73.683, 84.037) |
|---|
| 134 | add_chi2_entry(table, 49, 62.038, 66.339, 70.222, 74.919, 85.351) |
|---|
| 135 | add_chi2_entry(table, 50, 63.167, 67.505, 71.420, 76.154, 86.661) |
|---|
| 136 | add_chi2_entry(table, 51, 64.295, 68.669, 72.616, 77.386, 87.968) |
|---|
| 137 | add_chi2_entry(table, 52, 65.422, 69.832, 73.810, 78.616, 89.272) |
|---|
| 138 | add_chi2_entry(table, 53, 66.548, 70.993, 75.002, 79.843, 90.573) |
|---|
| 139 | add_chi2_entry(table, 54, 67.673, 72.153, 76.192, 81.069, 91.872) |
|---|
| 140 | add_chi2_entry(table, 55, 68.796, 73.311, 77.380, 82.292, 93.168) |
|---|
| 141 | add_chi2_entry(table, 56, 69.919, 74.468, 78.567, 83.513, 94.461) |
|---|
| 142 | add_chi2_entry(table, 57, 71.040, 75.624, 79.752, 84.733, 95.751) |
|---|
| 143 | add_chi2_entry(table, 58, 72.160, 76.778, 80.936, 85.950, 97.039) |
|---|
| 144 | add_chi2_entry(table, 59, 73.279, 77.931, 82.117, 87.166, 98.324) |
|---|
| 145 | add_chi2_entry(table, 60, 74.397, 79.082, 83.298, 88.379, 99.607) |
|---|
| 146 | add_chi2_entry(table, 61, 75.514, 80.232, 84.476, 89.591, 100.888) |
|---|
| 147 | add_chi2_entry(table, 62, 76.630, 81.381, 85.654, 90.802, 102.166) |
|---|
| 148 | add_chi2_entry(table, 63, 77.745, 82.529, 86.830, 92.010, 103.442) |
|---|
| 149 | add_chi2_entry(table, 64, 78.860, 83.675, 88.004, 93.217, 104.716) |
|---|
| 150 | add_chi2_entry(table, 65, 79.973, 84.821, 89.177, 94.422, 105.988) |
|---|
| 151 | add_chi2_entry(table, 66, 81.085, 85.965, 90.349, 95.626, 107.258) |
|---|
| 152 | add_chi2_entry(table, 67, 82.197, 87.108, 91.519, 96.828, 108.526) |
|---|
| 153 | add_chi2_entry(table, 68, 83.308, 88.250, 92.689, 98.028, 109.791) |
|---|
| 154 | add_chi2_entry(table, 69, 84.418, 89.391, 93.856, 99.228, 111.055) |
|---|
| 155 | add_chi2_entry(table, 70, 85.527, 90.531, 95.023, 100.425, 112.317) |
|---|
| 156 | add_chi2_entry(table, 71, 86.635, 91.670, 96.189, 101.621, 113.577) |
|---|
| 157 | add_chi2_entry(table, 72, 87.743, 92.808, 97.353, 102.816, 114.835) |
|---|
| 158 | add_chi2_entry(table, 73, 88.850, 93.945, 98.516, 104.010, 116.092) |
|---|
| 159 | add_chi2_entry(table, 74, 89.956, 95.081, 99.678, 105.202, 117.346) |
|---|
| 160 | add_chi2_entry(table, 75, 91.061, 96.217, 100.839, 106.393, 118.599) |
|---|
| 161 | add_chi2_entry(table, 76, 92.166, 97.351, 101.999, 107.583, 119.850) |
|---|
| 162 | add_chi2_entry(table, 77, 93.270, 98.484, 103.158, 108.771, 121.100) |
|---|
| 163 | add_chi2_entry(table, 78, 94.374, 99.617, 104.316, 109.958, 122.348) |
|---|
| 164 | add_chi2_entry(table, 79, 95.476, 100.749, 105.473, 111.144, 123.594) |
|---|
| 165 | add_chi2_entry(table, 80, 96.578, 101.879, 106.629, 112.329, 124.839) |
|---|
| 166 | add_chi2_entry(table, 81, 97.680, 103.010, 107.783, 113.512, 126.083) |
|---|
| 167 | add_chi2_entry(table, 82, 98.780, 104.139, 108.937, 114.695, 127.324) |
|---|
| 168 | add_chi2_entry(table, 83, 99.880, 105.267, 110.090, 115.876, 128.565) |
|---|
| 169 | add_chi2_entry(table, 84, 100.980, 106.395, 111.242, 117.057, 129.804) |
|---|
| 170 | add_chi2_entry(table, 85, 102.079, 107.522, 112.393, 118.236, 131.041) |
|---|
| 171 | add_chi2_entry(table, 86, 103.177, 108.648, 113.544, 119.414, 132.277) |
|---|
| 172 | add_chi2_entry(table, 87, 104.275, 109.773, 114.693, 120.591, 133.512) |
|---|
| 173 | add_chi2_entry(table, 88, 105.372, 110.898, 115.841, 121.767, 134.746) |
|---|
| 174 | add_chi2_entry(table, 89, 106.469, 112.022, 116.989, 122.942, 135.978) |
|---|
| 175 | add_chi2_entry(table, 90, 107.565, 113.145, 118.136, 124.116, 137.208) |
|---|
| 176 | add_chi2_entry(table, 91, 108.661, 114.268, 119.282, 125.289, 138.438) |
|---|
| 177 | add_chi2_entry(table, 92, 109.756, 115.390, 120.427, 126.462, 139.666) |
|---|
| 178 | add_chi2_entry(table, 93, 110.850, 116.511, 121.571, 127.633, 140.893) |
|---|
| 179 | add_chi2_entry(table, 94, 111.944, 117.632, 122.715, 128.803, 142.119) |
|---|
| 180 | add_chi2_entry(table, 95, 113.038, 118.752, 123.858, 129.973, 143.344) |
|---|
| 181 | add_chi2_entry(table, 96, 114.131, 119.871, 125.000, 131.141, 144.567) |
|---|
| 182 | add_chi2_entry(table, 97, 115.223, 120.990, 126.141, 132.309, 145.789) |
|---|
| 183 | add_chi2_entry(table, 98, 116.315, 122.108, 127.282, 133.476, 147.010) |
|---|
| 184 | add_chi2_entry(table, 99, 117.407, 123.225, 128.422, 134.642, 148.230) |
|---|
| 185 | add_chi2_entry(table, 100, 118.498, 124.342, 129.561, 135.807, 149.449) |
|---|
| 186 | |
|---|
| 187 | .return (table) |
|---|
| 188 | .end |
|---|
| 189 | |
|---|
| 190 | .sub make_random_hist |
|---|
| 191 | .param int min |
|---|
| 192 | .param int max |
|---|
| 193 | .param int num_samples |
|---|
| 194 | |
|---|
| 195 | if max > min goto get_pv |
|---|
| 196 | die "Max must be > min" |
|---|
| 197 | |
|---|
| 198 | get_pv: |
|---|
| 199 | .local int possible_values |
|---|
| 200 | possible_values = max - min |
|---|
| 201 | inc possible_values |
|---|
| 202 | |
|---|
| 203 | .local pmc histogram |
|---|
| 204 | histogram = new ['FixedIntegerArray'] |
|---|
| 205 | histogram = possible_values |
|---|
| 206 | |
|---|
| 207 | $I0 = max - min |
|---|
| 208 | init: |
|---|
| 209 | histogram[$I0] = 0 |
|---|
| 210 | dec $I0 |
|---|
| 211 | unless $I0 < 0 goto init |
|---|
| 212 | |
|---|
| 213 | $I0 = num_samples |
|---|
| 214 | if $I0 == 0 goto loop_done |
|---|
| 215 | loop: |
|---|
| 216 | .local int random |
|---|
| 217 | random = rand min, max |
|---|
| 218 | |
|---|
| 219 | random -= min |
|---|
| 220 | $I1 = histogram[random] |
|---|
| 221 | inc $I1 |
|---|
| 222 | histogram[random] = $I1 |
|---|
| 223 | |
|---|
| 224 | dec $I0 |
|---|
| 225 | unless $I0 <= 0 goto loop |
|---|
| 226 | loop_done: |
|---|
| 227 | |
|---|
| 228 | .return (histogram) |
|---|
| 229 | .end |
|---|
| 230 | |
|---|
| 231 | .sub print_histogram |
|---|
| 232 | .param pmc histogram |
|---|
| 233 | .param int num_samples |
|---|
| 234 | |
|---|
| 235 | .local num expected |
|---|
| 236 | expected = num_samples |
|---|
| 237 | $I0 = elements histogram |
|---|
| 238 | expected /= $I0 |
|---|
| 239 | |
|---|
| 240 | .local pmc sprintf_args |
|---|
| 241 | sprintf_args = new ['FixedPMCArray'] |
|---|
| 242 | sprintf_args = 2 |
|---|
| 243 | |
|---|
| 244 | $I0 = elements histogram |
|---|
| 245 | hist_loop: |
|---|
| 246 | dec $I0 |
|---|
| 247 | |
|---|
| 248 | $N0 = histogram[$I0] |
|---|
| 249 | $N0 /= expected |
|---|
| 250 | $P0 = box $N0 |
|---|
| 251 | sprintf_args[0] = $P0 |
|---|
| 252 | |
|---|
| 253 | $N1 = $N0 * 40 |
|---|
| 254 | $I1 = $N1 |
|---|
| 255 | $I1 -= 6 # width of printed $N0, plus space |
|---|
| 256 | if $I1 > 0 goto make_stars |
|---|
| 257 | $I1 = 0 |
|---|
| 258 | make_stars: |
|---|
| 259 | |
|---|
| 260 | $S0 = repeat '*', $I1 |
|---|
| 261 | $P0 = box $S0 |
|---|
| 262 | sprintf_args[1] = $P0 |
|---|
| 263 | |
|---|
| 264 | $S0 = sprintf "%5.3f %s", sprintf_args |
|---|
| 265 | say $S0 |
|---|
| 266 | |
|---|
| 267 | unless $I0 <= 0 goto hist_loop |
|---|
| 268 | |
|---|
| 269 | end: |
|---|
| 270 | .end |
|---|
| 271 | |
|---|
| 272 | .sub test_histogram |
|---|
| 273 | .param pmc histogram |
|---|
| 274 | .param num K |
|---|
| 275 | .param pmc table |
|---|
| 276 | |
|---|
| 277 | .local int degrees_of_freedom |
|---|
| 278 | $I0 = elements histogram |
|---|
| 279 | degrees_of_freedom = $I0 - 1 |
|---|
| 280 | |
|---|
| 281 | .local pmc sprintf_args |
|---|
| 282 | sprintf_args = new ['ResizablePMCArray'] |
|---|
| 283 | $P0 = box degrees_of_freedom |
|---|
| 284 | push sprintf_args, $P0 |
|---|
| 285 | |
|---|
| 286 | $I0 = exists table[degrees_of_freedom] |
|---|
| 287 | if $I0 goto get_data |
|---|
| 288 | |
|---|
| 289 | $S0 = sprintf "Don't have chi2 data for %d degrees of freedom", sprintf_args |
|---|
| 290 | skip(1, $S0) |
|---|
| 291 | goto end |
|---|
| 292 | |
|---|
| 293 | get_data: |
|---|
| 294 | |
|---|
| 295 | .local pmc chi2_data |
|---|
| 296 | chi2_data = table[degrees_of_freedom] |
|---|
| 297 | $N0 = chi2_data[0] |
|---|
| 298 | |
|---|
| 299 | $I0 = islt K, $N0 |
|---|
| 300 | |
|---|
| 301 | $P0 = box $N0 |
|---|
| 302 | unshift sprintf_args, $P0 |
|---|
| 303 | $P0= box K |
|---|
| 304 | unshift sprintf_args, $P0 |
|---|
| 305 | |
|---|
| 306 | $S0 = sprintf "K (%5.3f) should be less than limit (%5.3f) for %d degrees of freedom", sprintf_args |
|---|
| 307 | ok($I0, $S0) |
|---|
| 308 | |
|---|
| 309 | end: |
|---|
| 310 | .return ($I0) |
|---|
| 311 | .end |
|---|
| 312 | |
|---|
| 313 | .sub test_random_range |
|---|
| 314 | .param int range |
|---|
| 315 | .param int times |
|---|
| 316 | .param pmc table |
|---|
| 317 | |
|---|
| 318 | .local int num_samples |
|---|
| 319 | num_samples = range * times |
|---|
| 320 | |
|---|
| 321 | .local int min, max |
|---|
| 322 | min = 10 |
|---|
| 323 | max = min + range |
|---|
| 324 | |
|---|
| 325 | .local pmc histogram |
|---|
| 326 | histogram = make_random_hist(min, max, num_samples) |
|---|
| 327 | |
|---|
| 328 | .local num K |
|---|
| 329 | K = compute_chi_square(histogram, num_samples) |
|---|
| 330 | |
|---|
| 331 | $I0 = test_histogram(histogram, K, table) |
|---|
| 332 | |
|---|
| 333 | #if $I0 goto end |
|---|
| 334 | #print_histogram(histogram, num_samples) |
|---|
| 335 | end: |
|---|
| 336 | .end |
|---|
| 337 | |
|---|
| 338 | .sub main :main |
|---|
| 339 | .include 'test_more.pir' |
|---|
| 340 | |
|---|
| 341 | .const int TIMES = 10000 |
|---|
| 342 | |
|---|
| 343 | .local int num_tries |
|---|
| 344 | num_tries = 101 |
|---|
| 345 | |
|---|
| 346 | plan(num_tries) |
|---|
| 347 | |
|---|
| 348 | $S0 = "Running tests expecting #TIMES# hits per bucket" |
|---|
| 349 | $P0 = split '#', $S0 |
|---|
| 350 | $S0 = TIMES |
|---|
| 351 | $P0[1] = $S0 |
|---|
| 352 | $S0 = join '', $P0 |
|---|
| 353 | diag($S0) |
|---|
| 354 | |
|---|
| 355 | .local pmc table |
|---|
| 356 | table = make_chi2_table() |
|---|
| 357 | |
|---|
| 358 | $I0 = 0 |
|---|
| 359 | loop: |
|---|
| 360 | inc $I0 |
|---|
| 361 | test_random_range($I0, TIMES, table) |
|---|
| 362 | if $I0 < num_tries goto loop |
|---|
| 363 | |
|---|
| 364 | .end |
|---|
| 365 | |
|---|
