A Vector-of-vectors view of Sobol quasirandom Uniform and StdNormal draws

sobolrand.jl

using Sobol, ArraysOfArrays

function sobolrand(n,k, extraskip::Int=0)
    ω = SobolSeq(k)
    Sobol.skip(ω, n) # recommended in the Sobol.jl README
    for _ in 1:extraskip
        Sobol.next!(ω)
    end
    x = Matrix{Float64}(undef, k, n)
    for xcol in eachcol(x)
        Sobol.next!(ω, xcol)
    end
    nestedview(x)
end

function boxmuller!(x::AbstractVector)
    n = length(x)
    @assert iseven(n)
    @views for j in 1:n÷2
        twoj = 2j
        u = x[twoj-1]
        v = 2 * x[twoj]
        r = sqrt(-2 * log(u))
        x[twoj - 1] = r * cospi(v)
        x[twoj] = r * sinpi(v)
    end
    x
end

function sobolrandn(n,k, extraskip::Int=0)
    k2 = 2 * (k - k ÷ 2)
    x = sobolrand(n, k2, extraskip)
    boxmuller!.(x)
    nestedview(view(x.data, 1:k, :))
end

# julia> sobolrand(10,3)
# 10-element ArrayOfSimilarArrays{Float64, 1, 1, 2, Matrix{Float64}}:
#  [0.1875, 0.3125, 0.9375]
#  [0.6875, 0.8125, 0.4375]
#  [0.9375, 0.0625, 0.6875]
#  [0.4375, 0.5625, 0.1875]
#  [0.3125, 0.1875, 0.3125]
#  [0.8125, 0.6875, 0.8125]
#  [0.5625, 0.4375, 0.0625]
#  [0.0625, 0.9375, 0.5625]
#  [0.09375, 0.46875, 0.46875]
#  [0.59375, 0.96875, 0.96875]

# julia> sobolrandn(10,3)
# 10-element ArrayOfSimilarArrays{Float64, 1, 1, 2, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}}:
#  [-0.700211643609752, 1.6904604465342195, -0.33192491181244455]
#  [0.33127808631904765, -0.7997760489084498, 1.1879514292027367]
#  [0.33192491181244455, 0.13748780016220807, 0.33127808631904765]
#  [-1.1879514292027367, -0.4920655934162752, -0.700211643609752]
#  [0.5836771236304021, 1.4091212279154346, -1.4091212279154346]
#  [-0.24660933052559192, -0.5953675903626333, 0.5953675903626333]
#  [-0.991064091477675, 0.4105121878710228, 0.9011506174345589]
#  [2.1755700423514006, -0.9011506174345589, -0.4105121878710228]
#  [-2.13402452473067, 0.42448387026042717, -0.6839093044534189]
#  [1.0014551820921285, -0.1992018210945718, 0.13999641949147473]

# julia> using UnicodePlots

# julia> scatterplot(eachrow(sobolrand(1000,2).data)...)
#      ┌────────────────────────────────────────┐ 
#    1 │⠠⠂⣂⠅⠌⠠⠂⣁⠅⠨⠡⠨⢐⠐⡀⠡⠘⢐⠀⡀⠠⠂⣁⠡⠌⠠⠂⣂⠡⠌⠡⠊⢐⠐⡀⠡⠌⠐⠐⡀│ 
#      │⠂⠢⠀⠔⢁⡁⠢⠀⠔⢁⡈⠤⠁⠔⢈⡈⠢⠁⠔⢈⡁⠢⠀⠄⢃⡂⠢⠀⠄⢃⡘⠀⠅⠔⢈⡨⢠⠁⠔⢈│ 
#      │⠠⠂⢅⠑⡄⠠⠁⢅⠑⡄⠠⠊⡈⠐⠄⠐⠊⣈⠐⠄⠠⠁⢅⠃⡠⠠⠂⢅⠃⡄⠐⠘⣈⠀⠄⢐⠐⣈⠐⠄│ 
#      │⡅⢢⠀⠂⠂⡅⢡⠀⠂⠂⠘⠠⡀⡌⢨⠘⠀⠂⡌⢨⡅⢡⠀⠒⠀⡁⢢⠀⠒⠀⠈⠒⠀⡌⢨⠀⠒⡀⡌⢨│ 
#      │⢐⠁⡨⠈⢄⢐⠂⡈⠈⢄⡐⠁⠤⠐⡂⠐⠁⠥⠐⡂⢐⠐⡨⠈⢄⢐⠈⡨⠀⢄⠐⠁⠥⠂⡂⡐⠁⠤⠂⡐│ 
#      │⠌⢀⢂⠊⠠⠔⠰⢀⠊⠠⠄⡑⡀⡂⠡⠅⠑⡀⡂⠡⠄⢒⢀⠊⠠⠄⢑⢀⠊⠠⠁⠑⡀⡊⠠⠄⡑⡀⡊⠠│ 
#      │⢈⠈⠤⢀⢂⠨⠈⠤⢈⢂⡐⡀⠢⠁⠔⡐⠁⠢⠁⠆⢐⠁⠄⢈⢂⢈⠁⠤⢈⢂⡐⠁⠢⠈⠆⡐⡀⠢⠈⠆│ 
#      │⠀⡉⠄⠆⠰⠀⡉⡀⠆⠰⠆⡰⢀⢉⠀⠄⠱⢀⢉⠀⠈⡐⡀⠆⠰⠈⡀⠅⠆⠰⠆⠱⢀⢁⠁⠆⢰⢀⠁⠁│ 
#      │⠌⡀⠒⠈⡅⡈⡀⠒⠈⡅⢨⠀⠔⢀⠢⠨⠁⠄⢀⠢⡈⡀⠒⢁⠅⠌⡀⠒⢁⠌⠨⠈⠔⢀⠢⠨⡀⠔⢀⠢│ 
#      │⠂⠌⠄⡅⠐⠂⡈⠄⡅⠐⢂⢠⠡⠁⠐⢊⢘⠠⢁⠐⠀⡈⠄⢅⠐⠂⠌⠄⢅⠐⠂⡩⠠⢁⠐⢂⡨⠠⢁⠐│ 
#      │⠌⠄⠑⣀⠊⠌⡀⠑⣀⠃⠠⣀⠒⠠⠡⠐⣀⠒⠠⠡⠌⡀⠑⡀⡃⠌⠄⠑⡀⡃⢨⠀⠂⠠⠡⢠⠀⡒⠠⠡│ 
#      │⢃⢘⠠⠠⠀⢂⠸⠠⠀⠀⢀⠄⠄⡃⡘⢀⠄⠂⡃⢘⢃⠜⠠⠠⠀⢃⡘⠠⠠⠀⢀⠌⠂⢃⡘⢀⠄⠄⢃⡘│ 
#      │⠐⣀⠊⡠⠑⠐⢄⠊⡠⠑⠄⢄⠉⡠⠂⠂⢄⠉⡠⠂⠰⠀⡊⡠⠑⢐⢀⠂⡠⠑⠂⢄⠉⡀⠆⠄⢄⠉⡀⠆│ 
#      │⢁⠔⠐⠠⡈⠡⠔⠐⠠⡈⢀⠄⡂⠢⡈⢁⠂⡂⠢⡈⠡⠰⠐⠀⡈⠡⠨⠐⠠⡈⢁⠂⡂⠆⢈⢁⠄⡂⠆⡈│ 
#    0 │⠰⢀⠁⡐⢐⠰⠀⠍⡐⢐⡂⢄⠈⠄⠅⡂⢂⠈⠄⠅⠐⠤⠉⡐⢐⠈⠤⠉⡐⢐⡂⢂⠈⠤⠁⠂⢄⠈⠤⠁│ 
#      └────────────────────────────────────────┘ 
#      ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀1⠀ 

# julia> scatterplot(eachrow(sobolrandn(1000,2).data)...)
#       ┌────────────────────────────────────────┐ 
#     4 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#       │⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠈⠀⠄⠀⠀⡀⣇⠀⠀⢀⠀⠁⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⢂⠀⠀⡀⠐⠂⠄⢀⡔⠀⡀⣇⠀⠡⡀⠄⠐⠐⠀⠀⠀⡐⠀⠀⠀⠐⠀⠀⠀⠀⠀│ 
#       │⠀⠀⠄⠀⠀⡀⠀⢂⠀⠀⡢⠀⡄⠩⡰⠂⡢⣐⢥⢄⡧⠬⣈⣔⠐⢆⠍⢬⠀⢆⠀⠀⡂⠀⡀⠀⢀⠀⠀⠀│ 
#       │⠀⠀⢀⠀⠀⣀⠀⠄⠆⠄⠱⠤⠇⡖⠥⣎⣾⡵⣣⣷⣯⣟⢮⣷⣱⢮⣲⠨⢄⠣⠠⠌⠠⠀⢀⠀⠀⠀⡀⠀│ 
#       │⠀⠀⢀⠀⠀⠠⠐⡐⠪⢌⠱⡘⣌⣾⣧⣹⣷⣻⡏⣷⣿⢹⣏⡾⣯⣜⣶⣱⢣⠃⡡⠕⢐⠂⠄⠀⠀⠠⠀⠀│ 
#       │⠤⠤⠤⠬⠤⡥⠦⢬⡵⠶⡵⢮⢼⣯⢾⡿⣿⣿⣼⣿⣿⣿⣿⣿⢾⡷⣿⡯⣵⢯⠦⢮⡥⠴⢬⠤⠥⠤⠤⠤│ 
#       │⠀⠀⠂⠈⠀⠐⠠⠡⢌⢊⡱⢡⢮⠽⢝⢽⡿⣟⣷⡿⣿⣺⣽⢿⣫⡻⠽⡵⡌⢎⡑⡅⠌⠄⠂⠀⠁⠐⠀⠀│ 
#       │⠀⠀⠁⠀⠈⠈⠀⠁⠐⠂⢆⠂⢆⠿⢓⢏⠹⡻⣝⡻⣏⣫⢞⠯⡹⣚⡻⢰⠒⡰⠐⢐⠈⠀⠉⠀⠀⠈⠀⠀│ 
#       │⠀⠀⠀⠂⠀⠈⠀⠨⠀⠀⠕⠀⠃⡐⡹⠈⠕⡩⡒⠊⡗⢓⠍⠪⠄⢕⢂⠘⠀⠪⠀⠈⠠⠀⠈⠀⠀⠂⠀⠀│ 
#       │⠀⠀⠀⠀⠠⠀⠀⠀⠀⠌⠀⠀⠁⠄⠄⠂⠈⠄⠄⠑⡏⠠⠘⠁⠈⠠⠠⠈⠀⠀⠅⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⡀⠈⠀⠈⠀⡇⠁⠀⠂⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⡇⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#    -4 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
#       └────────────────────────────────────────┘ 
#       ⠀-3⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀3⠀