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