Randomized qmc

Project.toml

@@ -4,8 +4,11 @@ authors = ["ludoro <[email protected]>, Chris Rackauckas <accounts@chrisra
 version = "0.3.0"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LatticeRules = "73f95e8e-ec14-4e6a-8b18-0d2e271c4e55"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"

README.md

@@ -107,3 +107,99 @@ function sample(n,lb,ub,::NewAmazingSamplingAlgorithm)
     end
 end
 ```
+
+## Randomization
+
+Note that this feature is currently experimental and is thus subject to interface changes in
+non-breaking (minor) releases.
+
+Given a matrix `x` of size `d×n` and `xᵢₛ∈[0,1]ᵈ` one obtain a randomized version `y` using one the following methods
+* `owen_scramble(x, b; pad = pad)` where `b` is the base used to scramble and `pad` the number of bits in base `b` used to represent digits.
+* `matousek_scramble(x, base; pad = pad)`.
+* `digital_shift(x, base; pad = pad)`.
+* `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
+    pad = m
+
+    # Unrandomized low discrepency sequence
+    x_faure = QuasiMonteCarlo.sample(N, d, FaureSample())
+
+    # Randomized version
+    x_nus = randomize(x_faure, OwenScramble(base = b, pad = pad))
+    x_lms = randomize(x_faure, MatousekScramble(base = b, pad = pad))
+    x_digital_shift = randomize(x_faure, DigitalShift(base = b, pad = pad))
+    x_shift = randomize(x_faure, Shift())
+    x_uniform = rand(d, N) # 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 randomize 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 book 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/Faure.jl

@@ -20,7 +20,7 @@ end
 end
 
 """
-    FaureSample()
+    FaureSample(R::RandomizationMethod)
 
 A Faure low-discrepancy sequence.
 
@@ -40,8 +40,11 @@ References:
 Faure, H. (1982). Discrepance de suites associees a un systeme de numeration (en dimension s). *Acta Arith.*, 41, 337-351.
 Owen, A. B. (1997). Monte Carlo variance of scrambled net quadrature. *SIAM Journal on Numerical Analysis*, 34(5), 1884-1910.
 """
-struct FaureSample <: SamplingAlgorithm end
-@fastmath function sample(n::Integer, dimension::Integer, ::FaureSample, T = Float64;
+Base.@kwdef struct FaureSample <: DeterministicSamplingAlgorithm
+    R::RandomizationMethod = NoRand()
+end
+
+@fastmath function sample(n::Integer, dimension::Integer, S::FaureSample, T = Float64;
                           skipchecks = false)
     base = nextprime(dimension)
     n_digits = ceil(Int, log(base, n))
@@ -53,7 +56,7 @@ struct FaureSample <: SamplingAlgorithm end
                             "Try $n or $(n+base^power) instead."))
     end
 
-    return _faure_samples(n, n_digits, dimension, T)
+    return randomize(_faure_samples(n, n_digits, dimension, T), S.R)
 end
 
 @fastmath @views function _faure_samples(n_samples::I, n_digits::I, dimension::I,

src/Halton.jl

@@ -3,17 +3,19 @@
 
 Create a Halton sequence.
 """
-struct HaltonSample <: SamplingAlgorithm end
+Base.@kwdef @concrete struct HaltonSample <: DeterministicSamplingAlgorithm
+    R::RandomizationMethod = NoRand()
+end
 
-@views function sample(n::I, d::I, ::HaltonSample, T::Type = Float64) where {I <: Integer}
+@views function sample(n::I, d::I, S::HaltonSample, T::Type = Float64) where {I <: Integer}
     bases = nextprimes(one(n), d)
     n_digits = ceil.(I, log.(bases, n))
     λ = n .÷ bases .^ n_digits
     halton_seq = Matrix{T}(undef, d, n)
     for i in axes(halton_seq, 1)
         halton_seq[i, :] .= _vdc(λ[i], n_digits[i], bases[i], T; n)
     end
-    return halton_seq
+    return randomize(halton_seq, S.R)
 end
 
 # @fastmath @views function vdc(n::I, base::Vector{I}, F=Float64) where I <: Integer

src/Kronecker.jl

@@ -1,5 +1,5 @@
 """
-    KroneckerSample(generator::AbstractVector) <: SamplingAlgorithm
+    KroneckerSample(generator::AbstractVector) <: DeterministicSamplingAlgorithm
 
 A Kronecker sequence is a point set generated using a vector and equation:
 `x[i] = i * generator .% 1`
@@ -19,25 +19,26 @@ References:
 Leobacher, G., & Pillichshammer, F. (2014). *Introduction to quasi-Monte Carlo integration and applications.* Switzerland: Springer International Publishing.
 https://link.springer.com/content/pdf/10.1007/978-3-319-03425-6.pdf
 """
-struct KroneckerSample{V <: Union{AbstractVector, Missing}} <: SamplingAlgorithm
-    generator::V
+Base.@kwdef struct KroneckerSample{V <: Union{AbstractVector, Missing}} <:
+                   DeterministicSamplingAlgorithm
+    generator::V = missing
+    R::RandomizationMethod = NoRand()
 end
 
-KroneckerSample() = KroneckerSample(missing)
-function KroneckerSample(d::Integer, T = Float64)
+function KroneckerSample(d::Integer, R = NoRand(), T = Float64)
     ratio = harmonious(d, 2eps(T))
     generator = [ratio^i for i in 1:d]
-    return KroneckerSample(generator)
+    return KroneckerSample(generator, R)
 end
 
-function sample(n::Integer, d::Integer, ::KroneckerSample{Missing}, T = Float64)
-    return sample(n, d, KroneckerSample(d, T))
+function sample(n::Integer, d::Integer, S::KroneckerSample{Missing}, T = Float64)
+    return sample(n, d, KroneckerSample(d, S.R, T))
 end
 
 function sample(n::Integer, d::Integer, k::KroneckerSample{V},
                 T) where {V <: AbstractVector}
     @assert eltype(V)==T "Sample type must match generator type."
-    return sample(n, d, k)
+    return randomize(sample(n, d, k), k.R)
 end
 
 function sample(n::Integer, d::Integer, k::KroneckerSample{V}) where {V <: AbstractVector}

src/LatinHypercube.jl

@@ -1,9 +1,9 @@
 """
-    LatinHypercubeSample <: SamplingAlgorithm
+    LatinHypercubeSample <: RandomSamplingAlgorithm
 
 A Latin Hypercube is a point set with the property that every one-dimensional interval `(i/n, i+1/n)` contains exactly one point. This is a good way to sample a high-dimensional space, as it is more uniform than a random sample but does not require as many points as a full net.
 """
-Base.@kwdef @concrete struct LatinHypercubeSample <: SamplingAlgorithm
+Base.@kwdef @concrete struct LatinHypercubeSample <: RandomSamplingAlgorithm
     rng::AbstractRNG = Random.GLOBAL_RNG
 end
 

src/Lattices.jl

@@ -7,24 +7,25 @@ In more than 2 dimensions, grids have worse discrepancy than simple Monte Carlo.
 result, they should almost never be used for multivariate integration; their use is as
 a starting point for other algorithms.
 """
-struct GridSample <: SamplingAlgorithm end
+Base.@kwdef @concrete struct GridSample <: DeterministicSamplingAlgorithm
+    R::RandomizationMethod = NoRand()
+end
 
-function sample(n::Integer, d::Integer, ::GridSample, T = Float64)
+function sample(n::Integer, d::Integer, S::GridSample, T = Float64)
     n = convert(T, n)
-    return [(i - convert(T, 0.5)) / n for _ in 1:d, i in 1:n]
+    return randomize([(i - convert(T, 0.5)) / n for _ in 1:d, i in 1:n], S.R)
 end
 
 """
     LatticeRuleSample()
 
 Generate a point set using a lattice rule.
 """
-Base.@kwdef @concrete struct LatticeRuleSample <: SamplingAlgorithm
-    rng::AbstractRNG = Random.GLOBAL_RNG
+Base.@kwdef @concrete struct LatticeRuleSample <: DeterministicSamplingAlgorithm
+    R::RandomizationMethod = NoRand()
 end
 
 function sample(n::Integer, d::Integer, S::LatticeRuleSample, T = Float64)
-    rng = S.rng
     lat = LatticeRules.LatticeRule(d)
-    return reduce(hcat, lat[0:(n - 1)])
+    return randomize(reduce(hcat, lat[0:(n - 1)]), S.R)
 end

src/QuasiMonteCarlo.jl

@@ -1,9 +1,14 @@
 module QuasiMonteCarlo
 
 using Sobol, LatticeRules, Distributions, Primes, LinearAlgebra, Random
+using IntervalArithmetic, Combinatorics
 using ConcreteStructs
+using Accessors
 
 abstract type SamplingAlgorithm end
+abstract type RandomSamplingAlgorithm <: SamplingAlgorithm end
+abstract type DeterministicSamplingAlgorithm <: SamplingAlgorithm end
+abstract type RandomizationMethod end
 
 include("VanDerCorput.jl")
 include("Faure.jl")
@@ -27,11 +32,11 @@ end
 _check_sequence(n::Integer) = @assert n>0 ZERO_SAMPLES_MESSAGE
 
 """
-    RandomSample <: SamplingAlgorithm
+    RandomSample <: RandomSamplingAlgorithm
 
 Randomly distributed random numbers.
 """
-Base.@kwdef @concrete struct RandomSample <: SamplingAlgorithm
+Base.@kwdef @concrete struct RandomSample <: RandomSamplingAlgorithm
     rng::AbstractRNG = Random.GLOBAL_RNG
 end
 
@@ -72,7 +77,28 @@ function sample(n::Integer, d::Integer, D::Distributions.Sampleable)
     return reduce(hcat, x)
 end
 
-function generate_design_matrices(n, lb, ub, sampler, num_mats = 2)
+# 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
+
+"""
+```julia
+generate_design_matrices(n, lb, ub, sample_method, num_mats)
+```
+
+Create `num_mats` matrices each containing a QMC point set, where:
+- `n` is the number of points to sample.
+- `lb` is the lower bound for each variable. The length determines the dimensionality.
+- `ub` is the upper bound.
+- `d` is the dimensionality of the point set.
+- `sample_method` is the quasi-Monte Carlo sampling strategy. Note that any Distributions.jl
+  distribution can be used in addition to any `SamplingAlgorithm`.
+- You can specify a `RandomizationMethod` for `sample_method<:DeterministicSamplingAlgorithm`
+"""
+function generate_design_matrices(n, lb, ub, sampler::RandomSamplingAlgorithm, num_mats = 2)
     if n <= 0
         throw(ZeroSamplesError())
     end
@@ -82,7 +108,49 @@ 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
 
-export GridSample,
+function generate_design_matrices(n, lb, ub, sampler::DeterministicSamplingAlgorithm,
+                                  num_mats = 2)
+    if n <= 0
+        throw(ZeroSamplesError())
+    end
+    @assert length(lb) == length(ub)
+
+    # Generate a vector of num_mats independent "randomized" version of the QMC sequence
+    out = generate_design_matrices(n, length(lb), sampler, sampler.R, num_mats, eltype(lb))
+
+    # Scaling
+    for j in eachindex(out)
+        out[j] = (ub .- lb) .* out[j] .+ lb
+    end
+    return out
+end
+
+"""
+```julia
+NoRand
+```
+
+No Randomization is performed on the sampled sequence.
+"""
+struct NoRand <: RandomizationMethod end
+randomize(x, S::NoRand) = x
+
+"""
+    generate_design_matrices(n, d, sampler, R::NoRand, num_mats, T = Float64)
+`R = NoRand()` produces `num_mats` matrices each containing a different deterministic point set in `[0, 1)ᵈ`.
+Note that this is an ad hoc way to produce different `generate_design_matrices` as it create a deterministic point in dimension `d × num_mats` and split it in `num_mats` point set in dimension `d`.
+"""
+function generate_design_matrices(n, d, sampler, R::NoRand, num_mats, T = Float64)
+    out = sample(n, num_mats * d, sampler, T)
+    return [out[(j * d + 1):((j + 1) * d), :] for j in 0:(num_mats - 1)]
+end
+
+include("RandomizedQuasiMonteCarlo/shifting.jl")
+include("RandomizedQuasiMonteCarlo/net_utilities.jl")
+include("RandomizedQuasiMonteCarlo/scrambling_base_b.jl")
+
+export SamplingAlgorithm,
+       GridSample,
        UniformSample,
        SobolSample,
        LatinHypercubeSample,
@@ -94,6 +162,13 @@ export GridSample,
        KroneckerSample,
        SectionSample,
        FaureSample,
-       SamplingAlgorithm
-
+       randomize,
+       RandomizationMethod,
+       NoRand,
+       Shift,
+       ScrambleMethod,
+       OwenScramble,
+       MatousekScramble,
+       DigitalShift,
+       istmsnet
 end # module