Merge pull request #108 from mlsuh/randomized-halton

Add randomized halton sequence

docs/src/index.md

@@ -29,6 +29,7 @@ s = QuasiMonteCarlo.sample(n, lb, ub, SobolSample())
 s = QuasiMonteCarlo.sample(n, lb, ub, LatinHypercubeSample())
 s = QuasiMonteCarlo.sample(n, lb, ub, LatticeRuleSample())
 s = QuasiMonteCarlo.sample(n, lb, ub, HaltonSample())
+s = QuasiMonteCarlo.sample(n, lb, ub, RandomizedHaltonSample())
 ```
 
 The output `s` is a matrix, so one can use things like `@uview` from

docs/src/samplers.md

@@ -44,4 +44,5 @@ KroneckerSample
 
 ```@docs
 LatinHypercubeSample
+RandomizedHaltonSample
 ```

src/QuasiMonteCarlo.jl

@@ -55,8 +55,8 @@ Return a QMC point set where:
 In the first method the type of the point set is specified by `T` while in the second method the output type is inferred from the bound types.
 """
 function sample(n::Integer, lb::T, ub::T,
-    S::D) where {T <: Union{Base.AbstractVecOrTuple, Number},
-    D <: SamplingAlgorithm}
+        S::D) where {T <: Union{Base.AbstractVecOrTuple, Number},
+        D <: SamplingAlgorithm}
     _check_sequence(lb, ub, n)
     lb = float.(lb)
     ub = float.(ub)
@@ -78,6 +78,7 @@ include("Kronecker.jl")
 include("Halton.jl")
 include("Sobol.jl")
 include("LatinHypercube.jl")
+include("RandomizedHalton.jl")
 include("Lattices.jl")
 include("Section.jl")
 
@@ -99,6 +100,7 @@ export SamplingAlgorithm,
     GridSample,
     SobolSample,
     LatinHypercubeSample,
+    RandomizedHaltonSample,
     LatticeRuleSample,
     RandomSample,
     HaltonSample,

src/RandomizedHalton.jl

@@ -0,0 +1,41 @@
+"""
+    RandomizedHalton(rng::AbstractRNG = Random.GLOBAL_RNG) <: RandomSamplingAlgorithm
+
+    Create a randomized Halton sequence.
+
+    References:
+    Owen, A. (2017). *A randomized Halton algorithm in R*. https://doi.org/10.48550/arXiv.1706.02808
+"""
+Base.@kwdef @concrete struct RandomizedHaltonSample <: RandomSamplingAlgorithm
+    rng::AbstractRNG = Random.GLOBAL_RNG
+end
+
+function sample(n::Integer, d::Integer, S::RandomizedHaltonSample, T = Float64)
+    _check_sequence(n)
+    bases = nextprimes(one(n), d)
+    halton_seq = Matrix{T}(undef, d, n)
+
+    ind = collect(1:n)
+    for i in 1:d
+        halton_seq[i, :] = randradinv(ind, S.rng, bases[i])
+    end
+    return halton_seq
+end
+
+function randradinv(ind::Vector{Int}, rng::AbstractRNG, b::Int = 2)
+    b2r = 1 / b
+    ans = ind .* 0
+    res = ind
+
+    base_ind = 1:b
+
+    while (1 - b2r < 1)
+        dig = res .% b
+        perm = shuffle(rng, base_ind) .- 1
+        pdig = perm[convert.(Int, dig .+ 1)]
+        ans = ans .+ pdig .* b2r
+        b2r = b2r / b
+        res = (res .- dig) ./ b
+    end
+    return ans
+end

test/runtests.jl

@@ -43,7 +43,7 @@ The test is exact if the element of `net` are of type `Rational`. Indeed, in tha
 The conversion `Float` to `Rational` is usually possible with usual nets, e.g., Sobol, Faure (may require `Rational{BigInt}.(net)`).
 """
 function istmsnet(net::AbstractMatrix{T}; λ::I, t::I, m::I, s::I,
-    base::I) where {I <: Integer, T <: Real}
+        base::I) where {I <: Integer, T <: Real}
     pass = true
 
     @assert size(net, 2)==λ * (base^m) "Number of points must be as size(net, 2) = $(size(net, 2)) == λ * (base^m) = $(λ * (base^m))"
@@ -145,7 +145,7 @@ end
     s = sortslices(s; dims = 2)
     differences = diff(s; dims = 2)
     # test for grid
-    @test all(x -> all(y -> y in 0.0:1/n:1.0, x), eachrow(s))
+    @test all(x -> all(y -> y in 0.0:(1 / n):1.0, x), eachrow(s))
     μ = mean(s; dims = 2)
     variance = var(s; corrected = false, dims = 2)
     for i in eachindex(μ)
@@ -285,37 +285,39 @@ end
     end
 end
 
-@testset "Halton Sequence" begin
+@testset "(Randomized) Halton Sequence" begin
     d = 4
     lb = zeros(d)
     ub = ones(d)
     bases = nextprimes(1, d)
     n = prod(bases)^2
-    s = QuasiMonteCarlo.sample(n, lb, ub, HaltonSample())
-    @test isa(s, Matrix)
-    @test size(s) == (d, n)
-    sorted = reduce(vcat, sort.(eachslice(s; dims = 1))')
-    each_dim = eachrow(sorted)
-
-    # each 1d sequence should have base b stratification property
-    # (property inherited from van der Corput)
-    @test all(zip(each_dim, bases)) do (seq, base)
-        theoretical_count = n ÷ base
-        part = Iterators.partition(seq, theoretical_count)
-        all(enumerate(part)) do (i, subseq)
-            all(subseq) do x
-                i - 1 ≤ base * x ≤ i
+    for sample_type in [HaltonSample(), RandomizedHaltonSample()]
+        s = QuasiMonteCarlo.sample(n, lb, ub, sample_type)
+        @test isa(s, Matrix)
+        @test size(s) == (d, n)
+        sorted = reduce(vcat, sort.(eachslice(s; dims = 1))')
+        each_dim = eachrow(sorted)
+
+        # each 1d sequence should have base b stratification property
+        # (property inherited from van der Corput)
+        @test all(zip(each_dim, bases)) do (seq, base)
+            theoretical_count = n ÷ base
+            part = Iterators.partition(seq, theoretical_count)
+            all(enumerate(part)) do (i, subseq)
+                all(subseq) do x
+                    i - 1 ≤ base * x ≤ i
+                end
             end
         end
-    end
-    μ = mean(s; dims = 2)
-    variance = var(s; dims = 2)
-    for i in eachindex(μ)
-        @test μ[i]≈0.5 atol=1 / sqrt(n)
-        @test variance[i]≈1 / 12 rtol=1 / sqrt(n)
-    end
-    for (i, j) in combinations(1:d, 2)
-        @test pvalue(SignedRankTest(s[i, :], s[j, :])) > 0.0001
+        μ = mean(s; dims = 2)
+        variance = var(s; dims = 2)
+        for i in eachindex(μ)
+            @test μ[i]≈0.5 atol=1 / sqrt(n)
+            @test variance[i]≈1 / 12 rtol=1 / sqrt(n)
+        end
+        for (i, j) in combinations(1:d, 2)
+            @test pvalue(SignedRankTest(s[i, :], s[j, :])) > 0.0001
+        end
     end
 end