Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
76 changes: 42 additions & 34 deletions src/Sobol.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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())
Expand Down Expand Up @@ -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
Expand All @@ -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))")
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down
9 changes: 5 additions & 4 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)]
Expand All @@ -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

Expand Down
Loading