Merge pull request #119 from SciML/dw/sobol_allocations

Faster generation of Sobol samples

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "QuasiMonteCarlo"
 uuid = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 authors = ["ludoro <ludovicobessi@gmail.com>, Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -50,4 +50,4 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Combinatorics", "Distributions", "HypothesisTests", "IntervalArithmetic", "LinearAlgebra", "Primes", "Random", "Statistics", "StatsBase", "Test"]
+test = ["Aqua", "Combinatorics", "Distributions", "HypothesisTests", "IntervalArithmetic", "LinearAlgebra", "Primes", "Random", "Statistics", "StatsBase", "Test"]

src/Sobol.jl

@@ -8,7 +8,24 @@ Base.@kwdef @concrete struct SobolSample <: DeterministicSamplingAlgorithm
 end
 
 function sample(n::Integer, d::Integer, S::SobolSample, T = Float64)
-    s = Sobol.SobolSeq(zeros(T, d), ones(T, d))
-    skip(s, n)
-    return randomize(reduce(hcat, [Sobol.next!(s) for i in 1:n]), S.R)
+    if n < 0
+        throw(ArgumentError("number of samples must be non-negative"))
+    end
+
+    seq = Matrix{T}(undef, d, n)
+    if n == 0
+        return seq
+    end
+
+    # Use function barrier since `Sobol.SobolSeq(d)` can't be inferred
+    return _sample!(seq, Sobol.SobolSeq(d), S.R)
+end
+
+function _sample!(seq::AbstractMatrix, s::Sobol.SobolSeq, R::RandomizationMethod)
+    n = size(seq, 2)
+    Sobol.skip!(s, n, @view(seq[:, begin]))
+    for x in eachcol(seq)
+        Sobol.next!(s, x)
+    end
+    return randomize(seq, R)
 end