-
-
Notifications
You must be signed in to change notification settings - Fork 5.4k
/
misc.jl
374 lines (312 loc) · 9.18 KB
/
misc.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
# This file is a part of Julia. License is MIT: https://julialang.org/license
## rand!(::BitArray) && bitrand
function rand!(rng::AbstractRNG, B::BitArray, ::SamplerType{Bool})
isempty(B) && return B
Bc = B.chunks
rand!(rng, Bc)
Bc[end] &= Base._msk_end(B)
return B
end
"""
bitrand([rng=GLOBAL_RNG], [dims...])
Generate a `BitArray` of random boolean values.
# Examples
```jldoctest
julia> rng = MersenneTwister(1234);
julia> bitrand(rng, 10)
10-element BitArray{1}:
false
true
true
true
true
false
true
false
false
true
```
"""
bitrand(r::AbstractRNG, dims::Dims) = rand!(r, BitArray(undef, dims))
bitrand(r::AbstractRNG, dims::Integer...) = rand!(r, BitArray(undef, convert(Dims, dims)))
bitrand(dims::Dims) = rand!(BitArray(undef, dims))
bitrand(dims::Integer...) = rand!(BitArray(undef, convert(Dims, dims)))
## randstring (often useful for temporary filenames/dirnames)
"""
randstring([rng=GLOBAL_RNG], [chars], [len=8])
Create a random string of length `len`, consisting of characters from
`chars`, which defaults to the set of upper- and lower-case letters
and the digits 0-9. The optional `rng` argument specifies a random
number generator, see [Random Numbers](@ref).
# Examples
```jldoctest
julia> Random.seed!(0); randstring()
"0IPrGg0J"
julia> randstring(MersenneTwister(0), 'a':'z', 6)
"aszvqk"
julia> randstring("ACGT")
"TATCGGTC"
```
!!! note
`chars` can be any collection of characters, of type `Char` or
`UInt8` (more efficient), provided [`rand`](@ref) can randomly
pick characters from it.
"""
function randstring end
let b = UInt8['0':'9';'A':'Z';'a':'z']
global randstring
randstring(r::AbstractRNG, chars=b, n::Integer=8) = String(rand(r, chars, n))
randstring(r::AbstractRNG, n::Integer) = randstring(r, b, n)
randstring(chars=b, n::Integer=8) = randstring(GLOBAL_RNG, chars, n)
randstring(n::Integer) = randstring(GLOBAL_RNG, b, n)
end
## randsubseq & randsubseq!
# Fill S (resized as needed) with a random subsequence of A, where
# each element of A is included in S with independent probability p.
# (Note that this is different from the problem of finding a random
# size-m subset of A where m is fixed!)
function randsubseq!(r::AbstractRNG, S::AbstractArray, A::AbstractArray, p::Real)
@assert !has_offset_axes(S, A)
0 <= p <= 1 || throw(ArgumentError("probability $p not in [0,1]"))
n = length(A)
p == 1 && return copyto!(resize!(S, n), A)
empty!(S)
p == 0 && return S
nexpected = p * length(A)
sizehint!(S, round(Int,nexpected + 5*sqrt(nexpected)))
if p > 0.15 # empirical threshold for trivial O(n) algorithm to be better
for i = 1:n
rand(r) <= p && push!(S, A[i])
end
else
# Skip through A, in order, from each element i to the next element i+s
# included in S. The probability that the next included element is
# s==k (k > 0) is (1-p)^(k-1) * p, and hence the probability (CDF) that
# s is in {1,...,k} is 1-(1-p)^k = F(k). Thus, we can draw the skip s
# from this probability distribution via the discrete inverse-transform
# method: s = ceil(F^{-1}(u)) where u = rand(), which is simply
# s = ceil(log(rand()) / log1p(-p)).
# -log(rand()) is an exponential variate, so can use randexp().
L = -1 / log1p(-p) # L > 0
i = 0
while true
s = randexp(r) * L
s >= n - i && return S # compare before ceil to avoid overflow
push!(S, A[i += ceil(Int,s)])
end
# [This algorithm is similar in spirit to, but much simpler than,
# the one by Vitter for a related problem in "Faster methods for
# random sampling," Comm. ACM Magazine 7, 703-718 (1984).]
end
return S
end
"""
randsubseq!([rng=GLOBAL_RNG,] S, A, p)
Like [`randsubseq`](@ref), but the results are stored in `S`
(which is resized as needed).
# Examples
```jldoctest
julia> rng = MersenneTwister(1234);
julia> S = Int64[];
julia> randsubseq!(rng, S, collect(1:8), 0.3);
julia> S
2-element Array{Int64,1}:
7
8
```
"""
randsubseq!(S::AbstractArray, A::AbstractArray, p::Real) = randsubseq!(GLOBAL_RNG, S, A, p)
randsubseq(r::AbstractRNG, A::AbstractArray{T}, p::Real) where {T} =
randsubseq!(r, T[], A, p)
"""
randsubseq([rng=GLOBAL_RNG,] A, p) -> Vector
Return a vector consisting of a random subsequence of the given array `A`, where each
element of `A` is included (in order) with independent probability `p`. (Complexity is
linear in `p*length(A)`, so this function is efficient even if `p` is small and `A` is
large.) Technically, this process is known as "Bernoulli sampling" of `A`.
# Examples
```jldoctest
julia> rng = MersenneTwister(1234);
julia> randsubseq(rng, collect(1:8), 0.3)
2-element Array{Int64,1}:
7
8
```
"""
randsubseq(A::AbstractArray, p::Real) = randsubseq(GLOBAL_RNG, A, p)
## rand Less Than Masked 52 bits (helper function)
"Return a sampler generating a random `Int` (masked with `mask`) in ``[0, n)``, when `n <= 2^52`."
ltm52(n::Int, mask::Int=nextpow(2, n)-1) = LessThan(n-1, Masked(mask, UInt52Raw(Int)))
## shuffle & shuffle!
"""
shuffle!([rng=GLOBAL_RNG,] v::AbstractArray)
In-place version of [`shuffle`](@ref): randomly permute `v` in-place,
optionally supplying the random-number generator `rng`.
# Examples
```jldoctest
julia> rng = MersenneTwister(1234);
julia> shuffle!(rng, Vector(1:16))
16-element Array{Int64,1}:
2
15
5
14
1
9
10
6
11
3
16
7
4
12
8
13
```
"""
function shuffle!(r::AbstractRNG, a::AbstractArray)
@assert !has_offset_axes(a)
n = length(a)
n <= 1 && return a # nextpow below won't work with n == 0
@assert n <= Int64(2)^52
mask = nextpow(2, n) - 1
for i = n:-1:2
(mask >> 1) == i && (mask >>= 1)
j = 1 + rand(r, ltm52(i, mask))
a[i], a[j] = a[j], a[i]
end
return a
end
shuffle!(a::AbstractArray) = shuffle!(GLOBAL_RNG, a)
"""
shuffle([rng=GLOBAL_RNG,] v::AbstractArray)
Return a randomly permuted copy of `v`. The optional `rng` argument specifies a random
number generator (see [Random Numbers](@ref)).
To permute `v` in-place, see [`shuffle!`](@ref). To obtain randomly permuted
indices, see [`randperm`](@ref).
# Examples
```jldoctest
julia> rng = MersenneTwister(1234);
julia> shuffle(rng, Vector(1:10))
10-element Array{Int64,1}:
6
1
10
2
3
9
5
7
4
8
```
"""
shuffle(r::AbstractRNG, a::AbstractArray) = shuffle!(r, copymutable(a))
shuffle(a::AbstractArray) = shuffle(GLOBAL_RNG, a)
## randperm & randperm!
"""
randperm([rng=GLOBAL_RNG,] n::Integer)
Construct a random permutation of length `n`. The optional `rng`
argument specifies a random number generator (see [Random
Numbers](@ref)). The element type of the result is the same as the type
of `n`.
To randomly permute an arbitrary vector, see [`shuffle`](@ref) or
[`shuffle!`](@ref).
# Examples
```jldoctest
julia> randperm(MersenneTwister(1234), 4)
4-element Array{Int64,1}:
2
1
4
3
```
"""
randperm(r::AbstractRNG, n::T) where {T <: Integer} = randperm!(r, Vector{T}(undef, n))
randperm(n::Integer) = randperm(GLOBAL_RNG, n)
"""
randperm!([rng=GLOBAL_RNG,] A::Array{<:Integer})
Construct in `A` a random permutation of length `length(A)`. The
optional `rng` argument specifies a random number generator (see
[Random Numbers](@ref)). To randomly permute an arbitrary vector, see
[`shuffle`](@ref) or [`shuffle!`](@ref).
# Examples
```jldoctest
julia> randperm!(MersenneTwister(1234), Vector{Int}(undef, 4))
4-element Array{Int64,1}:
2
1
4
3
```
"""
function randperm!(r::AbstractRNG, a::Array{<:Integer})
n = length(a)
@assert n <= Int64(2)^52
n == 0 && return a
a[1] = 1
mask = 3
@inbounds for i = 2:n
j = 1 + rand(r, ltm52(i, mask))
if i != j # a[i] is undef (and could be #undef)
a[i] = a[j]
end
a[j] = i
i == 1+mask && (mask = 2mask + 1)
end
return a
end
randperm!(a::Array{<:Integer}) = randperm!(GLOBAL_RNG, a)
## randcycle & randcycle!
"""
randcycle([rng=GLOBAL_RNG,] n::Integer)
Construct a random cyclic permutation of length `n`. The optional `rng`
argument specifies a random number generator, see [Random Numbers](@ref).
The element type of the result is the same as the type of `n`.
# Examples
```jldoctest
julia> randcycle(MersenneTwister(1234), 6)
6-element Array{Int64,1}:
3
5
4
6
1
2
```
"""
randcycle(r::AbstractRNG, n::T) where {T <: Integer} = randcycle!(r, Vector{T}(undef, n))
randcycle(n::Integer) = randcycle(GLOBAL_RNG, n)
"""
randcycle!([rng=GLOBAL_RNG,] A::Array{<:Integer})
Construct in `A` a random cyclic permutation of length `length(A)`.
The optional `rng` argument specifies a random number generator, see
[Random Numbers](@ref).
# Examples
```jldoctest
julia> randcycle!(MersenneTwister(1234), Vector{Int}(undef, 6))
6-element Array{Int64,1}:
3
5
4
6
1
2
```
"""
function randcycle!(r::AbstractRNG, a::Array{<:Integer})
n = length(a)
n == 0 && return a
@assert n <= Int64(2)^52
a[1] = 1
mask = 3
@inbounds for i = 2:n
j = 1 + rand(r, ltm52(i-1, mask))
a[i] = a[j]
a[j] = i
i == 1+mask && (mask = 2mask + 1)
end
return a
end
randcycle!(a::Array{<:Integer}) = randcycle!(GLOBAL_RNG, a)