Randomized qmc

README.md

@@ -107,3 +107,98 @@ function sample(n,lb,ub,::NewAmazingSamplingAlgorithm)
     end
 end
 ```
+
+## Randomization
+
+API is meant to evolve.
+
+Given a matrix `x` of size `n×d` and `xᵢₛ∈[0,1]ᵈ` one obtain a randomized version `y` using one the following methods
+* `nested_uniform_scramble(x, b; M = M)` where `b` is the base used to scramble and `M` the number of bits in base `b` used to represent digits.
+* `linear_matrix_scramble(x, base; M = M)`.
+* `digital_shift(x, base; M = M)`.
+* `shift(x)`.
+
+All these functions guarantee that the resulting array will have its components uniformly distributed `yᵢₛ∼𝐔([0,1]ᵈ)` (but not independent).
+
+### Example 
+
+Randomization of a Faure sequence with various methods.
+
+```julia
+    m = 4
+    d = 3
+    b = QuasiMonteCarlo.nextprime(d)
+    N = b^m # Number of points
+    M = m
+
+    # Unrandomized low discrepency sequence
+    x_faure = permutedims(QuasiMonteCarlo.sample(N, d, FaureSample()))
+
+    # Randomized version
+    x_nus = nested_uniform_scramble(x_faure, b; M = M)
+    x_lms = linear_matrix_scramble(x_faure, b; M = M)
+    x_digital_shift = digital_shift(x_faure, b; M = M)
+    x_shift = shift(x_faure)
+    x_uniform = rand(N, d) # plain i.i.d. uniform
+```
+
+```julia
+using Plots
+# Setting I like for plotting
+default(fontfamily="Computer Modern", linewidth=1, label=nothing, grid=true, framestyle=:default)
+```
+
+Plot the resulting sequences along dimensions `1` and `3`.
+
+```julia
+begin
+    d1 = 1
+    d2 = 3
+    sequences = [x_uniform, x_faure, x_shift, x_digital_shift, x_lms, x_nus]
+    names = ["Uniform", "Faure (unrandomized)", "Shift", "Digital Shift", "Linear Matrix Scrambling", "Nested Uniform Scrambling"]
+    p = [plot(thickness_scaling=1.5, aspect_ratio=:equal) for i in sequences]
+    for (i, x) in enumerate(sequences)
+        scatter!(p[i], x[:, d1], x[:, d2], ms=0.9, c=1, grid=false)
+        title!(names[i])
+        xlims!(p[i], (0, 1))
+        ylims!(p[i], (0, 1))
+        yticks!(p[i], [0, 1])
+        xticks!(p[i], [0, 1])
+        hline!(p[i], range(0, 1, step=1 / 4), c=:gray, alpha=0.2)
+        vline!(p[i], range(0, 1, step=1 / 4), c=:gray, alpha=0.2)
+        hline!(p[i], range(0, 1, step=1 / 2), c=:gray, alpha=0.8)
+        vline!(p[i], range(0, 1, step=1 / 2), c=:gray, alpha=0.8)
+    end
+    plot(p..., size=(1500, 900))
+end
+```
+
+![Different randomization methods of the same initial set of points](img/various_randomization.svg)
+
+Faure nets and scrambled versions of Faure nets are digital $(t,d,m)$-net ([see this nice reference by A. Owen](https://artowen.su.domains/mc/qmcstuff.pdf)). It basically means that they have strong equipartition properties.
+Here we can (visually) verify that with Nested Uniform Scrambling (it also works with Linear Matrix Scrambling and Digital Shift).
+You must see one point per rectangle of volume $1/b^m$.
+
+```julia
+begin
+    d1 = 1 
+    d2 = 3
+    x = x_nus
+    p = [plot(thickness_scaling=1.5, aspect_ratio=:equal) for i in 0:m]
+    for i in 0:m
+        j = m - i
+        xᵢ = range(0, 1, step=1 / b^(i))
+        xⱼ = range(0, 1, step=1 / b^(j))
+        scatter!(p[i+1], x[:, d1], x[:, d2], ms=2, c=1, grid=false)
+        xlims!(p[i+1], (0, 1.01))
+        ylims!(p[i+1], (0, 1.01))
+        yticks!(p[i+1], [0, 1])
+        xticks!(p[i+1], [0, 1])
+        hline!(p[i+1], xᵢ, c=:gray, alpha=0.2)
+        vline!(p[i+1], xⱼ, c=:gray, alpha=0.2)
+    end
+    plot(p..., size=(1500, 900))
+end
+```
+
+![All the different elementary rectangle contain only one points](img/equidistribution.svg)

src/QuasiMonteCarlo.jl

@@ -7,6 +7,9 @@ abstract type SamplingAlgorithm end
 
 include("Faure.jl")
 include("Kronecker.jl")
+include("RandomizedQuasiMonteCarlo/conversion.jl")
+include("RandomizedQuasiMonteCarlo/shifting.jl")
+include("RandomizedQuasiMonteCarlo/scrambling_base_b.jl")
 
 check_bounds(lb, ub) = all(x -> x[1] <= x[2], zip(lb, ub))
 check_bounds(lb::Number, ub::Number) = lb <= ub
@@ -407,6 +410,13 @@ function generate_design_matrices(n, lb, ub, sampler, num_mats = 2)
     [out[(j * d + 1):((j + 1) * d), :] for j in 0:(num_mats - 1)]
 end
 
+# See https://discourse.julialang.org/t/is-there-a-dedicated-function-computing-m-int-log-b-b-m/89776/10
+function logi(b::Int, n::Int)
+    m = round(Int, log(b, n))
+    b^m == n || throw(ArgumentError("$n is not a power of $b"))
+    return m
+end
+
 export GridSample,
        UniformSample,
        SobolSample,
@@ -417,6 +427,10 @@ export GridSample,
        GoldenSample,
        KroneckerSample,
        SectionSample,
-       FaureSample
+       FaureSample,
+       nested_uniform_scramble,
+       linear_matrix_scramble,
+       shift,
+       digital_shift
 
 end # module

src/RandomizedQuasiMonteCarlo/conversion.jl

@@ -0,0 +1,107 @@
+"""
+    unif2bits(x<:AbstractMatrix, b::Integer; M=32)
+Return the b-adic decomposition of all element y of an array y = ∑ₖ yₖ/bᵏ a number yₖ∈[0,1[ -> [y₁, ⋯, yₘ] 
+"""
+function unif2bits(x::AbstractArray, b::Integer; M = 32)
+    bits = zeros(Int, M, size(x)...)
+    unif2bits!(bits, x, b)
+    return bits
+end
+
+function unif2bits!(bits::AbstractArray, x::AbstractArray, b::Integer)
+    @assert all([size(bits, d + 1) == size(x, d) for d in ndims(x)]) "Bit array of size $(size(bits)) instead of (M, $(size(x)))"
+    for i in CartesianIndices(x)
+        unif2bits!(@view(bits[:, i]), x[i], b)
+    end
+end
+
+function unif2bits(x::AbstractArray; M = 32)
+    bits = BitArray(undef, M, size(x)...)
+    unif2bits!(bits, x, 2)
+    return bits
+end
+"""
+    unif2bits(y<:Real, b::Integer; M=32)
+Return the b-adic decomposition y = ∑ₖ yₖ/bᵏ a number y∈[0,1[ -> [y₁, ⋯, yₘ] 
+"""
+function unif2bits(y::Real, b::Integer; M = 32)
+    bits = zeros(Int, M)
+    unif2bits!(bits, y, b)
+    return bits
+end
+"""
+    unif2bits(y; M=32)
+Return the binary decomposition y = ∑ₖ yₖ/2ᵏ a number y∈[0,1[ -> [y₁, ⋯, yₘ] 
+"""
+function unif2bits(y; M = 32)
+    bits = BitArray(undef, M)
+    unif2bits!(bits, y, 2)
+    return bits
+end
+
+check_zero(a::Rational; kwargs...) = a == 0
+#! isequal(0.18518518518518515... -1/3^2-2/3^3 , 0) Should be true but is not! It fouls the binary expansion! Hence the isapprox with a hand tuned tolerence
+check_zero(a::AbstractFloat; atol = 3e-16) = isapprox(a, 0, atol = atol)
+
+function unif2bits!(bits::AbstractVector{<:Integer}, y, b::Integer; kwargs...)
+    bits .= 0
+    for j in eachindex(bits), bb in (b - 1):-1:1
+        a = y - bb // b^j
+        if check_zero(a; kwargs...)
+            bits[j] = bb
+            break # it breaks from the nested loop (see here)[https://stackoverflow.com/questions/39796234/how-to-break-out-of-nested-for-loops-in-julia]
+        elseif a > 0
+            bits[j] = bb
+            y = a
+        end
+    end
+end
+
+"""
+    bits2unif(bits::AbstractVector{<:Integer}, b::Integer)
+Convert a vector of M "bits" in base b into a number y∈[0,1[.
+"""
+function bits2unif(::Type{T}, bits::AbstractVector{<:Integer},
+                   b::Integer) where {T <: Rational}
+    # Turn sequence of bits into a point in [0,1)
+    # First bits are highest order
+    y = zero(T)
+    for j in lastindex(bits):-1:1
+        y = (y + bits[j]) // b
+    end
+    return y
+end
+
+function bits2unif(::Type{T}, bits::AbstractVector{<:Integer},
+                   b::Integer) where {T <: AbstractFloat}
+    # Turn sequence of bits into a point in [0,1)
+    # First bits are highest order
+    y = zero(T)
+    for j in lastindex(bits):-1:1
+        y = (y + bits[j]) / b
+    end
+    return y
+end
+
+function bits2unif(bits::AbstractVector{<:Integer}, b::Integer)
+    bits2unif(Float64, bits::AbstractVector{<:Integer}, b::Integer)
+end
+
+#? This seems faster than @evalpoly(b, $bit...)
+#?  bi = rand(0:2,32);
+#? @btime @evalpoly(3, $bi...)
+#?   500.515 ns (1 allocation: 272 bytes)
+#? @btime QuasiMonteCarlo.bits2int($bi, 3)
+#?   13.113 ns (0 allocations: 0 bytes)  
+"""
+    bits2int(bit::AbstractMatrix{<:Integer}, b::Integer)
+Convert a vector of M "bits" in base b into an integer.
+"""
+function bits2int(bit::AbstractVector, b::Integer)
+    m = length(bit)
+    y = 0
+    for k in m:-1:1
+        y = y * b + bit[k]
+    end
+    return y
+end