Line data Source code
1 : % Copyright (C) 2012-2021,2022 John E. Davis
2 : %
3 : % This file is part of the S-Lang Library and may be distributed under the
4 : % terms of the GNU General Public License. See the file COPYING for
5 : % more information.
6 : %---------------------------------------------------------------------------
7 4 : import ("rand");
8 :
9 : private define get_generator_args (nargs, num_parms, parmsp, rtp, nump,
10 : usage_str)
11 : {
12 122 : @rtp = NULL;
13 122 : @nump = NULL;
14 122 : if (nargs == num_parms)
15 : {
16 39 : @parmsp = __pop_list (num_parms);
17 : return;
18 : }
19 :
20 83 : if (nargs == num_parms + 1)
21 : {
22 81 : @nump = ();
23 81 : variable parms = __pop_list (num_parms);
24 81 : if (typeof (parms[0]) == Rand_Type)
25 : {
26 1 : @rtp = list_pop (parms);
27 1 : list_append (parms, @nump);
28 1 : @nump = NULL;
29 : }
30 81 : @parmsp = parms;
31 : return;
32 : }
33 :
34 2 : if (nargs == num_parms + 2)
35 : {
36 1 : @nump = ();
37 1 : @parmsp = __pop_list (num_parms);
38 1 : variable rt = ();
39 1 : if (typeof (rt) == Rand_Type)
40 : {
41 1 : @rtp = rt;
42 : return;
43 : }
44 : }
45 1 : else _pop_n (nargs);
46 :
47 1 : usage (usage_str);
48 : }
49 :
50 : private define call_rand_func ()
51 : {
52 184 : variable num = ();
53 184 : variable args = __pop_list (_NARGS-3);
54 : variable rt, func;
55 184 : (func, rt) = ();
56 :
57 184 : if (rt == NULL)
58 : {
59 182 : if (num == NULL)
60 60 : return (@func) (__push_list(args));
61 :
62 122 : return (@func) (__push_list(args), num);
63 : }
64 :
65 2 : if (num == NULL)
66 1 : return (@func)(rt, __push_list(args));
67 :
68 1 : return (@func)(rt, __push_list(args), num);
69 : }
70 :
71 : define rand_flat ()
72 : {
73 : variable parms, rt, num;
74 :
75 3 : get_generator_args (_NARGS, 2, &parms, &rt, &num,
76 : "r = rand_flat ([Rand_Type,] xmin, xmax [,num])");
77 :
78 3 : variable r = call_rand_func (&rand_uniform, rt, num);
79 :
80 3 : return parms[0] + (parms[1] - parms[0])*__tmp(r);
81 : }
82 :
83 : define rand_chisq ()
84 : {
85 : variable parms, rt, num;
86 :
87 27 : get_generator_args (_NARGS, 1, &parms, &rt, &num,
88 : "r = rand_chisq ([Rand_Type,] nu [,num])");
89 27 : return 2.0 * call_rand_func (&rand_gamma, rt, 0.5*parms[0], 1.0, num);
90 : }
91 :
92 : define rand_fdist ()
93 : {
94 : variable parms, rt, num;
95 :
96 48 : get_generator_args (_NARGS, 2, &parms, &rt, &num,
97 : "r = rand_fdist ([Rand_Type,] nu1, nu2 [,num])");
98 96 : variable nu1 = parms[0], nu2 = parms[1];
99 :
100 48 : return (call_rand_func (&rand_gamma, rt, 0.5*nu1, 1.0, num)/nu1)
101 : / (call_rand_func(&rand_gamma, rt, 0.5*nu2, 1.0, num)/nu2);
102 : }
103 :
104 : define rand_tdist ()
105 : {
106 : variable parms, rt, num;
107 :
108 15 : get_generator_args (_NARGS, 1, &parms, &rt, &num,
109 : "r = rand_tdist ([Rand_Type,] nu, [,num])");
110 15 : variable nu = parms[0];
111 15 : return call_rand_func (&rand_gauss, rt, 1.0, num)
112 : / sqrt(call_rand_func(&rand_chisq, rt, nu, num)/nu);
113 : }
114 :
115 : define rand_int ()
116 : {
117 : variable parms, rt, num;
118 :
119 11 : get_generator_args (_NARGS, 2, &parms, &rt, &num,
120 : "r = rand_int ([Rand_Type,] imin, imax [,num])");
121 :
122 : variable
123 10 : imin = typecast (parms[0], Int32_Type),
124 10 : imax = typecast (parms[1], Int32_Type), di;
125 :
126 10 : if (imin > imax)
127 0 : throw InvalidParmError, "rand_int: imax < imin";
128 :
129 10 : di = typecast (imax - imin, UInt32_Type);
130 :
131 10 : variable r = call_rand_func (&rand, rt, num);
132 :
133 10 : if (di + 1 != 0) % UINT32_MAX does not exist
134 9 : r = __tmp(r) mod (di + 1);
135 :
136 10 : return typecast (imin + r, Int32_Type);
137 : }
138 :
139 : define rand_exp ()
140 : {
141 : variable parms, rt, num;
142 :
143 18 : get_generator_args (_NARGS, 1, &parms, &rt, &num,
144 : "r = rand_exp ([Rand_Type,] beta [,num])");
145 :
146 18 : return (-parms[0]) * log (call_rand_func (&rand_uniform_pos, rt, num));
147 : }
148 :
149 : private define make_indices (a, d, i)
150 : {
151 48 : _for (0, length(array_shape(a))-1, 1)
152 : {
153 96 : variable j = ();
154 96 : if (j == d)
155 48 : i;
156 : else
157 48 : [:];
158 : }
159 : }
160 :
161 : define rand_sample ()
162 : {
163 32 : if (_NARGS < 2)
164 : {
165 0 : _pop_n (_NARGS);
166 0 : usage ("(B1 [,B2,...]) = rand_sample ([Rand_Type,] A1 [A2,...], num)");
167 : }
168 :
169 32 : variable num = ();
170 32 : variable arrays = __pop_list (_NARGS-1);
171 32 : variable rt = NULL;
172 :
173 32 : if (typeof (arrays[0]) == Rand_Type)
174 16 : rt = list_pop (arrays);
175 :
176 32 : variable n0 = NULL, dim0;
177 : variable a, indices;
178 :
179 32 : foreach a (arrays)
180 : {
181 48 : dim0 = array_shape (a)[0];
182 48 : if (n0 == NULL)
183 : {
184 32 : n0 = dim0;
185 32 : continue;
186 : }
187 16 : if (n0 != dim0)
188 0 : throw TypeMismatchError, "The arrays passed to rand_sample must have the same leading dimension";
189 : }
190 :
191 32 : if (num > n0)
192 0 : num = n0;
193 :
194 32 : if (rt == NULL)
195 16 : indices = rand_permutation (n0);
196 : else
197 16 : indices = rand_permutation (rt, n0);
198 32 : if (num < n0)
199 24 : indices = indices[[0:num-1]];
200 :
201 32 : foreach a (arrays)
202 48 : a[make_indices(a, 0, indices)];
203 : }
204 :
|
