From 9ea980cbdfea3db04e1e67018392d984fcb738ec Mon Sep 17 00:00:00 2001 From: mlsuh Date: Mon, 29 Jan 2024 15:24:58 +0100 Subject: [PATCH 1/3] Add randomized halton Implement algorithm by A. Owen to sample a randomized halton sequence. Created RandomizedHaltonSample <: RandomSamplingAlgorithm. --- src/QuasiMonteCarlo.jl | 6 ++++-- src/RandomizedHalton.jl | 36 ++++++++++++++++++++++++++++++++++++ test/runtests.jl | 38 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 76 insertions(+), 4 deletions(-) create mode 100644 src/RandomizedHalton.jl diff --git a/src/QuasiMonteCarlo.jl b/src/QuasiMonteCarlo.jl index f8e78a5..74697ed 100644 --- a/src/QuasiMonteCarlo.jl +++ b/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, diff --git a/src/RandomizedHalton.jl b/src/RandomizedHalton.jl new file mode 100644 index 0000000..83efb9f --- /dev/null +++ b/src/RandomizedHalton.jl @@ -0,0 +1,36 @@ +""" + RandomizedHalton(rng::AbstractRNG = Random.GLOBAL_RNG) <: RandomSamplingAlgorithm + + A Halton-sequence randomized using the algorithm presented by A. Owen in "A randomized Halton algorithm in R" (2017). +""" +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) + rng = S.rng + bases = nextprimes(one(n), d) + halton_seq = Matrix{T}(undef, d, n) + + for i in 1:d + halton_seq[i, :] = randradinv(collect(1:n), 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 + + while (1 - b2r < 1) + dig = res .% b + perm = shuffle(rng, collect(1:b)) .- 1 + pdig = perm[convert.(Int, dig .+ 1)] + ans = ans .+ pdig .* b2r + b2r = b2r / b + res = (res .- dig) ./ b + end + return ans +end diff --git a/test/runtests.jl b/test/runtests.jl index ab2d014..db94523 100644 --- a/test/runtests.jl +++ b/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(μ) @@ -194,6 +194,40 @@ end @test istmsnet(s, λ = 1, t = 0, m = 1, s = d, base = n) end +@testset "RandomizedHaltonSample" begin + d = 4 + lb = zeros(d) + ub = ones(d) + bases = nextprimes(1, d) + n = prod(bases)^2 + s = QuasiMonteCarlo.sample(n, lb, ub, RandomizedHaltonSample()) + @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 + μ = 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 + @testset "Van der Corput Sequence" begin lb = 0 ub = 1 From 1d4be78b9e10b09828be3596d5d989df634adbc0 Mon Sep 17 00:00:00 2001 From: mlsuh Date: Tue, 30 Jan 2024 13:36:33 +0100 Subject: [PATCH 2/3] Refactor RandomizedHaltonSample and test --- src/RandomizedHalton.jl | 13 +++++-- test/runtests.jl | 84 +++++++++++++---------------------------- 2 files changed, 35 insertions(+), 62 deletions(-) diff --git a/src/RandomizedHalton.jl b/src/RandomizedHalton.jl index 83efb9f..b7b359a 100644 --- a/src/RandomizedHalton.jl +++ b/src/RandomizedHalton.jl @@ -1,7 +1,10 @@ """ RandomizedHalton(rng::AbstractRNG = Random.GLOBAL_RNG) <: RandomSamplingAlgorithm - A Halton-sequence randomized using the algorithm presented by A. Owen in "A randomized Halton algorithm in R" (2017). + 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 @@ -9,12 +12,12 @@ end function sample(n::Integer, d::Integer, S::RandomizedHaltonSample, T = Float64) _check_sequence(n) - rng = S.rng 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(collect(1:n), S.rng, bases[i]) + halton_seq[i, :] = randradinv(ind, S.rng, bases[i]) end return halton_seq end @@ -24,9 +27,11 @@ function randradinv(ind::Vector{Int}, rng::AbstractRNG, b::Int = 2) ans = ind .* 0 res = ind + base_ind = 1:b + while (1 - b2r < 1) dig = res .% b - perm = shuffle(rng, collect(1:b)) .- 1 + perm = shuffle(rng, base_ind) .- 1 pdig = perm[convert.(Int, dig .+ 1)] ans = ans .+ pdig .* b2r b2r = b2r / b diff --git a/test/runtests.jl b/test/runtests.jl index db94523..f18abd8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -194,40 +194,6 @@ end @test istmsnet(s, λ = 1, t = 0, m = 1, s = d, base = n) end -@testset "RandomizedHaltonSample" begin - d = 4 - lb = zeros(d) - ub = ones(d) - bases = nextprimes(1, d) - n = prod(bases)^2 - s = QuasiMonteCarlo.sample(n, lb, ub, RandomizedHaltonSample()) - @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 - μ = 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 - @testset "Van der Corput Sequence" begin lb = 0 ub = 1 @@ -319,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 From 75cbc8dee8d17bb9a8567a0abf530dc7216d42a5 Mon Sep 17 00:00:00 2001 From: mlsuh Date: Tue, 30 Jan 2024 13:44:14 +0100 Subject: [PATCH 3/3] Add RandomizedHaltonSample to docs --- docs/src/index.md | 1 + docs/src/samplers.md | 1 + 2 files changed, 2 insertions(+) diff --git a/docs/src/index.md b/docs/src/index.md index 6bfeca6..a2a6c8e 100644 --- a/docs/src/index.md +++ b/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 diff --git a/docs/src/samplers.md b/docs/src/samplers.md index 06d2368..883648a 100644 --- a/docs/src/samplers.md +++ b/docs/src/samplers.md @@ -44,4 +44,5 @@ KroneckerSample ```@docs LatinHypercubeSample +RandomizedHaltonSample ```