diff --git a/Project.toml b/Project.toml index dc2faec..cef11b8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Sobol" uuid = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" -version = "1.5.0" +version = "2.0.0" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/src/Sobol.jl b/src/Sobol.jl index 2bbb87b..810d530 100644 --- a/src/Sobol.jl +++ b/src/Sobol.jl @@ -4,27 +4,27 @@ export SobolSeq, ScaledSobolSeq, next! include("soboldata.jl") #loads `sobol_a` and `sobol_minit` -abstract type AbstractSobolSeq{N} end +abstract type AbstractSobolSeq end -# N iis the dimension of sequence being generated -mutable struct SobolSeq{N} <: AbstractSobolSeq{N} - m::Array{UInt32,2} #array of size (sdim, 32) - x::Array{UInt32,1} #previous x = x_n, array of length sdim - b::Array{UInt32,1} #position of fixed point in x[i] is after bit b[i] +mutable struct SobolSeq{X<:AbstractVector{UInt32}} <: AbstractSobolSeq + m::Matrix{UInt32} #array of size (sdim, 32) + x::X #previous x = x_n, array of length sdim + b::Vector{UInt32} #position of fixed point in x[i] is after bit b[i] n::UInt32 #number of x's generated so far end -ndims(s::AbstractSobolSeq{N}) where {N} = N::Int +Base.ndims(s::SobolSeq) = length(s.x) -function SobolSeq(N::Int) +function SobolSeq(x::AbstractVector{UInt32}) + N = length(x) (N < 0 || N > (length(sobol_a) + 1)) && error("invalid Sobol dimension") m = ones(UInt32, (N, 32)) #special cases - N == 0 && return(SobolSeq{0})(m,UInt32[],UInt32[],zero(UInt32)) + N == 0 && return SobolSeq(m,UInt32[],x,zero(UInt32)) #special cases 1 - N == 1 && return(SobolSeq{N}(m,UInt32[0],UInt32[0],zero(UInt32))) + N == 1 && return SobolSeq(m,UInt32[0],fill!(x,zero(UInt32)),zero(UInt32)) for i = 2:N a = sobol_a[i-1] @@ -42,9 +42,9 @@ function SobolSeq(N::Int) end end end - SobolSeq{N}(m,zeros(UInt32,N),zeros(UInt32,N),zero(UInt32)) + SobolSeq(m,zeros(UInt32,N),fill!(x,zero(UInt32)),zero(UInt32)) end -SobolSeq(N::Integer) = SobolSeq(Int(N)) +SobolSeq(N::Integer) = SobolSeq(Vector{UInt32}(undef, N)) function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat}) length(x) != ndims(s) && throw(BoundsError()) @@ -72,7 +72,7 @@ function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat}) end return x end -next!(s::SobolSeq) = next!(s, Array{Float64,1}(undef, ndims(s))) +next!(s::SobolSeq) = next!(s, Vector{Float64}(undef, ndims(s))) # if we know in advance how many points (n) we want to compute, then # adopt a suggestion similar to the Joe and Kuo paper, which in turn @@ -95,7 +95,7 @@ function skip!(s::SobolSeq, n::Integer, x; exact=false) for unused=1:nskip; next!(s,x); end return s end -Base.skip(s::SobolSeq, n::Integer; exact=false) = skip!(s, n, Array{Float64,1}(undef, ndims(s)); exact=exact) +Base.skip(s::SobolSeq, n::Integer; exact=false) = skip!(s, n, Vector{Float64}(undef, ndims(s)); exact=exact) function Base.show(io::IO, s::SobolSeq) print(io, "$(ndims(s))-dimensional Sobol sequence on [0,1]^$(ndims(s))") @@ -115,20 +115,28 @@ Base.IteratorEltype(::Type{<:AbstractSobolSeq}) = Base.HasEltype() # Convenience wrapper for scaled Sobol sequences -struct ScaledSobolSeq{N,T} <: AbstractSobolSeq{N} - s::SobolSeq{N} - lb::Vector{T} - ub::Vector{T} - function ScaledSobolSeq{N,T}(lb::Vector{T}, ub::Vector{T}) where {N,T} - length(lb)==length(ub)==N || throw(DimensionMismatch("lb and ub do not have length $N")) - new(SobolSeq(N), lb, ub) +struct ScaledSobolSeq{T,X<:AbstractVector{UInt32},B<:AbstractVector{T}} <: AbstractSobolSeq + s::SobolSeq{X} + lb::B + ub::B + function ScaledSobolSeq(x::X, lb::B, ub::B) where {T,X<:AbstractVector{UInt32},B<:AbstractVector{T}} + length(x)==length(lb)==length(ub) || throw(DimensionMismatch("x, lb, and ub do not all have same length")) + new{T,X,B}(SobolSeq(x), lb, ub) end end -function SobolSeq(N::Integer, lb, ub) - T = typeof(sum(ub) - sum(lb)) - ScaledSobolSeq{Int(N),T}(copyto!(Vector{T}(undef,N), lb), copyto!(Vector{T}(undef,N), ub)) +function SobolSeq(x::AbstractVector{UInt32}, lb::B, ub::B) where {B} + if Base.IteratorEltype(B) === Base.HasEltype() + T = eltype(B) + else + i = iterate(lb) + isnothing(i) && error("Cannot determine eltype of empty iterator") + T = typeof(i[1]) + end + ScaledSobolSeq(x, copyto!(Vector{T}(undef,length(x)), lb), copyto!(Vector{T}(undef, length(x)), ub)) end +SobolSeq(N::Integer, lb, ub) = SobolSeq(Vector{UInt32}(undef, N), lb, ub) SobolSeq(lb, ub) = SobolSeq(length(lb), lb, ub) +Base.ndims(s::ScaledSobolSeq) = ndims(s.s) function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat}, lb::AbstractVector, ub::AbstractVector) @@ -139,22 +147,22 @@ function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat}, end return x end -function next!(s::SobolSeq{N}, lb::AbstractVector, ub::AbstractVector) where {N} +function next!(s::SobolSeq, lb::AbstractVector, ub::AbstractVector) T = typeof(float((zero(eltype(ub)) - zero(eltype(lb))))) - next!(s, Vector{T}(undef, N), lb, ub) + next!(s, Vector{T}(undef, ndims(s)), lb, ub) end next!(s::ScaledSobolSeq, x::AbstractVector{<:AbstractFloat}) = next!(s.s, x, s.lb, s.ub) -next!(s::ScaledSobolSeq{N,T}) where {N,T} = next!(s.s, Vector{float(T)}(undef, N), s.lb, s.ub) -Base.eltype(::Type{ScaledSobolSeq{N,T}}) where {N,T} = Vector{float(T)} +next!(s::ScaledSobolSeq{T}) where {T} = next!(s.s, Vector{float(T)}(undef, ndims(s)), s.lb, s.ub) +Base.eltype(::Type{<:ScaledSobolSeq{T}}) where {T} = Vector{float(T)} Base.skip(s::ScaledSobolSeq, n; exact = false) = (skip(s.s, n; exact = exact); s) -function Base.show(io::IO, s::ScaledSobolSeq{N,T}) where {N,T} +function Base.show(io::IO, s::ScaledSobolSeq{T}) where {T} lb = s.lb; ub = s.ub - print(io, "$N-dimensional scaled $(float(T)) Sobol sequence on [$(lb[1]),$(ub[1])]") + print(io, "$(ndims(s))-dimensional scaled $(float(T)) Sobol sequence on [$(lb[1]),$(ub[1])]") cnt = 1 - for i = 2:N + for i = 2:ndims(s) if lb[i] == lb[i-1] && ub[i] == ub[i-1] cnt += 1 else @@ -170,11 +178,11 @@ function Base.show(io::IO, s::ScaledSobolSeq{N,T}) where {N,T} end end -function Base.show(io::IO, ::MIME"text/html", s::ScaledSobolSeq{N,T}) where {N,T} +function Base.show(io::IO, ::MIME"text/html", s::ScaledSobolSeq{T}) where {T} lb = s.lb; ub = s.ub - print(io, "$N-dimensional scaled $(float(T)) Sobol sequence on [$(lb[1]),$(ub[1])]") + print(io, "$(ndims(s))-dimensional scaled $(float(T)) Sobol sequence on [$(lb[1]),$(ub[1])]") cnt = 1 - for i = 2:N + for i = 2:ndims(s) if lb[i] == lb[i-1] && ub[i] == ub[i-1] cnt += 1 else diff --git a/test/runtests.jl b/test/runtests.jl index 1555d98..a53b7be 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,7 +16,7 @@ using Sobol, Test for dim in dimensions println("Testing dimension $(dim)") open(joinpath(dirname(@__FILE__), "results", "exp_results_$(dim)")) do exp_results_file - s = SobolSeq(dim) + s = @inferred SobolSeq(dim) x = zeros(dim) for line in eachline(exp_results_file) values = [parse(Float64, item) for item in split(line)] @@ -40,13 +40,14 @@ end lb = [-1,0,0] ub = [1,3,2] N = length(lb) - s = SobolSeq(lb,ub) - @test s isa ScaledSobolSeq{3,Int} + s = @inferred SobolSeq(lb,ub) + @test s isa ScaledSobolSeq{Int} && ndims(s) == 3 @test eltype(s) == Vector{Float64} @test eltype(SobolSeq(Float32.(lb),Float32.(ub))) == Vector{Float32} @test first(s) == [0,1.5,1] @test first(SobolSeq((x for x in lb), (x for x in ub))) == [0,1.5,1] - @test SobolSeq(N,lb,ub) isa ScaledSobolSeq{3,Int} + s = @inferred SobolSeq(N,lb,ub) + @test s isa ScaledSobolSeq{Int} && ndims(s) == 3 @test_throws BoundsError SobolSeq(2,lb,ub) end