1() = evalfile("./test.sl"); 2require("rand.sl"); 3 4private variable CLOSED_UPPER = 1; 5private variable CLOSED_LOWER = 2; 6 7define test_generic (func, ret_type, parms, exp_mean, exp_variance, 8 min_r, max_r, interval_flags) 9{ 10 variable a, b, r; 11 12 if (typeof ((@func) (__push_list(parms))) != ret_type) 13 failed ("${func} did not produce a ${ret_type}"$); 14 variable num = 10000; 15 16 a = (@func) (__push_list(parms), num); 17 if ((length (a) != num) 18 || (_typeof (a) != ret_type)) 19 failed ("${func}(${num}) did not produce a ${num} ${ret_type}"$); 20 21 if (any(isnan (a))) 22 { 23 failed ("${func} produced NaN values"$); 24 } 25 26 if ((min_r != NULL) && (max_r != NULL)) 27 { 28 variable range_error = 0; 29 range_error += ((interval_flags & CLOSED_LOWER) && any (a < min_r)); 30 range_error += (((interval_flags & CLOSED_LOWER) == 0) && any (a <= min_r)); 31 range_error += ((interval_flags & CLOSED_UPPER) && any (a > max_r)); 32 range_error += (((interval_flags & CLOSED_UPPER) == 0) && any (a >= max_r)); 33 34 if (range_error) 35 failed ("${func} produced values outside the ${min_r}-${max_r} range"$); 36 } 37 38 if (exp_variance != NULL) 39 { 40 variable mean = sum (a)/num; 41 variable w = 1.0/sqrt(num); 42 variable exp_stddev = sqrt(exp_variance); 43 variable mean_lo = exp_mean - 3*exp_stddev*w; 44 variable mean_hi = exp_mean + 3*exp_stddev*w; 45 46 ifnot (mean_lo <= mean <= mean_hi) 47 failed ("${func}'s mean ${mean} outside the expected range: ${mean_lo} - ${mean_hi}"$); 48 49 variable stddev = sqrt(sum((a-mean)^2)/num); 50 ifnot (feqs (stddev, exp_stddev, 0.1, 1e-4)) 51 failed ("${func}'s stddev ${stddev} differs from expected value ${exp_stddev} (var=${exp_variance})"$); 52 } 53 54 b = Int_Type[3,2,1]; 55 a = (@func) (__push_list(parms), b); 56 ifnot (all (array_shape (a) == array_shape(b))) 57 failed ("${func}(a) failed to produce an array with the dimensions of a"); 58} 59 60define test_rand_uniform () 61{ 62 test_generic (&rand_uniform, Double_Type, {}, 0.5, 1.0/12, 0, 1, CLOSED_LOWER); 63} 64 65define test_rand_uniform_pos () 66{ 67 test_generic (&rand_uniform_pos, Double_Type, {}, 0.5, 1.0/12, 0, 1, 0); 68} 69 70define test_rand_gauss () 71{ 72 foreach ([0, 0.5, 1.0, 10.0]) 73 { 74 variable sigma = (); 75 test_generic (&rand_gauss, Double_Type, {sigma}, 0, sigma^2, NULL, NULL, CLOSED_LOWER); 76 } 77} 78 79define test_rand_poisson () 80{ 81 foreach ([0.1, 1, 10, 100]) 82 { 83 variable lam = (); 84 test_generic (&rand_poisson, UInt_Type, {lam}, lam, lam, 85 0, _Inf, CLOSED_LOWER); 86 } 87} 88 89define test_rand_gamma () 90{ 91 foreach ([0.1, 0.5, 10, 20]) 92 { 93 variable theta = (); 94 foreach ([1, 2, 4, 16]) 95 { 96 variable k = (); 97 test_generic (&rand_gamma, Double_Type, {k, theta}, 98 k*theta, k*theta^2, 99 0, _Inf, CLOSED_LOWER); 100 } 101 } 102} 103 104define test_rand_binomial() 105{ 106 variable p, n; 107 foreach p ([0, 0.25, 0.75, 1.0]) 108 { 109 foreach n ([1, 2, 10, 20, 50]) % make sure n*MIN(p,1-p) > 10 for one case 110 { 111 test_generic (&rand_binomial, UInt_Type, {p, n}, 112 n*p, n*p*(1-p), 113 0, n, CLOSED_LOWER|CLOSED_UPPER); 114 } 115 } 116} 117 118define test_rand_beta() 119{ 120 variable a, b; 121 foreach a ([0.1, 0.5, 1, 2, 8]) 122 { 123 foreach b ([0.1, 0.5, 1, 4, 16]) 124 { 125 test_generic (&rand_beta, Double_Type, {a, b}, 126 a/(a+b), (a*b)/(a+b)^2/(a+b+1), 127 0, 1, CLOSED_LOWER|CLOSED_UPPER); 128 } 129 } 130} 131 132define test_rand_cauchy () 133{ 134 variable gamma; 135 foreach gamma ([0.01, 0.1, 1, 10, 100]) 136 { 137 test_generic (&rand_cauchy, Double_Type, {gamma}, 138 0.0, NULL, 139 -_Inf, _Inf, 0); 140 } 141} 142 143define test_rand_geometric () 144{ 145 variable p; 146 foreach p ([0.01, 0.2, 0.6, 1.0]) 147 { 148 test_generic (&rand_geometric, UInt_Type, {p}, 149 1/p, (1-p)/p^2, 1, _Inf, CLOSED_LOWER); 150 } 151} 152 153define test_rand_flat () 154{ 155 variable x0 = 2-0.5, x1 = 2+0.5; 156 test_generic (&rand_flat, Double_Type, {x0, x1}, 157 0.5*(x0+x1), (x1-x0)^2/12.0, 158 x0, x1, CLOSED_LOWER); 159} 160 161define test_rand_fdist () 162{ 163 variable nu1, nu2; 164 foreach nu1 ([1,2,5,20]) 165 { 166 foreach nu2 ([1,3,12,20]) 167 { 168 variable mean = NULL, variance = NULL; 169 if (nu2 > 4) 170 { 171 mean = nu2/(nu2-2.0); 172 variance = 2.0*mean^2*(nu1+nu2-2.0)/nu1/(nu2-4.0); 173 } 174 test_generic (&rand_fdist, Double_Type, {nu1, nu2}, 175 mean, variance, 0, _Inf, CLOSED_LOWER); 176 } 177 } 178} 179 180define test_rand_chisq () 181{ 182 variable nu; 183 foreach nu ([1,2,5,20]) 184 { 185 test_generic (&rand_chisq, Double_Type, {nu}, 186 nu, 2.0*nu, 0, _Inf, CLOSED_LOWER); 187 } 188} 189 190define test_rand_tdist () 191{ 192 variable nu; 193 foreach nu ([0.5,1,2,5,20]) 194 { 195 variable mean = 0.0, variance = NULL; 196 if (nu > 2) 197 variance = nu/(nu-2.0); 198 199 test_generic (&rand_tdist, Double_Type, {nu}, 200 mean, variance, -_Inf, _Inf, 0); 201 } 202} 203 204define test_rand_int () 205{ 206 variable x0 = 1, x1 = 10; 207 test_generic (&rand_int, Int_Type, {x0, x1}, 208 0.5*(x0+x1), ((x1-x0+1)^2-1)/12.0, x0, x1, 209 CLOSED_LOWER|CLOSED_UPPER); 210 variable n = 10000, s = sqrt(n); 211 x0 = -2, x1 = 2; 212 variable r = rand_int (x0, x1, (x1-x0+1)*n); 213 variable failed = 0; 214 _for (x0, x1, 1) 215 { 216 variable x = (); 217 variable i = where (x == r); 218 if (n-3*s < length(i) < n+3*s) 219 continue; 220 failed++; 221 } 222 223 x0 = INT_MIN; x1 = INT_MAX; 224 variable nbins = 10, binsize = typecast (x1 - x0, UInt_Type); 225 binsize = typecast (binsize/nbins, Int_Type); 226 227 r = rand_int (x0, x1, nbins * n); 228 loop (nbins) 229 { 230 x1 = x0 + binsize; 231 if (x1 < x0) 232 break; 233 234 i = where (x0 <= r < x1); 235 ifnot (n - 3*s < length(i) < n + 3*s) 236 failed++; 237 238 x0 = x1; 239 } 240 241 if (failed) 242 { 243 () = fputs ("\ 244***WARNING: rand_int produced random numbers outside the 3 sigma range.\n\ 245 Run the test again, and if it fails with this warning then\n\ 246 file a possible bug report. This is a statistical test and can\n\ 247 produce false positives a small fraction of the time.\n", 248 stderr); 249 } 250} 251 252define test_rand_exp () 253{ 254 variable beta; 255 foreach beta ([0.1, 0.5, 1, 2, 10, 20]) 256 { 257 test_generic (&rand_exp, Double_Type, {beta}, 258 beta, beta^2, 0, _Inf, CLOSED_LOWER); 259 } 260} 261 262private define test_rand () 263{ 264 variable r = rand_new (); 265 srand (r, 12345); 266 srand (12345); 267 268 variable expected_values = 269 [746892674U, 935820662U, 3317285904U, 3160065947U, 888929593U, 316432806U, 270 3891177748U, 66504584U, 827237220U, 2731412032U, 105892519U, 1105593792U, 271 4257164826U, 2953826281U, 477842505U, 3161051103U, 654741546U, 2422625584U, 272 3232900523U, 2360188805U, 70104872U, 2440288176U, 1468162482U, 3428658486U, 273 1893960569U, 3842583023U, 4119673423U, 2214061288U, 1769001282U, 3996933411U, 274 3342430755U, 201679192U, 431385446U, 2648942401U, 2718561501U, 1689889009U, 275 403183793U, 662574206U, 2167963286U, 2166423399U, 112978312U, 3881586706U, 276 2111007051U, 2589350153U, 3519959692U, 838415425U, 1148338613U, 1576844827U, 277 1939688263U, 1896225294U, 2843247453U, 179614524U, 594767376U, 4056452573U, 278 2301108737U, 844660562U, 1079103954U, 3239244907U, 3213734172U, 3924276547U 279 ]; 280 variable a = rand (expected_values); 281 variable b = rand (r, a); 282 283 if (any (a != expected_values)) 284 failed ("The generator failed to produce the expected values"); 285 286 if (any (a != b)) 287 failed ("Seeding of the new generator failed"); 288} 289 290define test_rand_permutation () 291{ 292 loop (20) 293 { 294 variable p = 1+rand_permutation (10); 295 if ((prod(p) != 3628800) || (sum(p) != (10*11)/2)) 296 { 297 failed ("rand_permutation failed"); 298 } 299 } 300} 301 302private define check_sampled_array (rt, a, b, n) 303{ 304 variable dims = array_shape (a); 305 dims[0] = n; 306 ifnot (_eqs (dims, array_shape(b))) 307 failed ("rand_sample failed for %S array using rt=%S, n=%S: found %S", a, rt, n,b); 308} 309 310private define sample_and_check_array (rt, a, n) 311{ 312 variable b; 313 if (rt == NULL) 314 b = rand_sample (a, n); 315 else 316 b = rand_sample (rt, a, n); 317 318 check_sampled_array (rt, a, b, n); 319} 320 321define test_rand_sample () 322{ 323 variable a0 = [1:20]; 324 variable a1 = _reshape ([1:20*30], [20, 30]); 325 variable a2 = _reshape ([1:20*30*4], [20, 30, 4]); 326 327 variable dims, a, b, n, b0, b1, b2; 328 variable rt = rand_new (); 329 foreach n ([5, 1, 20, 0]) 330 { 331 foreach a ({a0,a1,a2}) 332 { 333 sample_and_check_array (NULL, a, n); 334 sample_and_check_array (rt, a, n); 335 } 336 (b0, b1, b2) = rand_sample (a0, a1, a2, n); 337 check_sampled_array (NULL, a0, b0, n); 338 check_sampled_array (NULL, a1, b1, n); 339 check_sampled_array (NULL, a2, b2, n); 340 341 (b0, b1, b2) = rand_sample (rt, a0, a1, a2, n); 342 check_sampled_array (rt, a0, b0, n); 343 check_sampled_array (rt, a1, b1, n); 344 check_sampled_array (rt, a2, b2, n); 345 } 346} 347 348define slsh_main () 349{ 350 testing_module ("rand"); 351 352 test_rand (); 353 test_rand_permutation (); 354 test_rand_sample (); 355 356 test_rand_uniform (); 357 test_rand_uniform_pos (); 358 test_rand_gauss (); 359 test_rand_poisson (); 360 test_rand_gamma (); 361 test_rand_binomial (); 362 test_rand_beta (); 363 test_rand_cauchy (); 364 test_rand_geometric (); 365 test_rand_flat (); 366 test_rand_fdist (); 367 test_rand_chisq (); 368 test_rand_tdist (); 369 test_rand_int (); 370 test_rand_exp (); 371 372 end_test (); 373} 374