Skip to content
Open
Show file tree
Hide file tree
Changes from 3 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
64 changes: 32 additions & 32 deletions src/Sobol.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,27 +4,26 @@ 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}(N::Integer) where {X<:AbstractVector{UInt32}}
(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(undef,0),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(undef,1),zero(UInt32)),zero(UInt32))

for i = 2:N
a = sobol_a[i-1]
Expand All @@ -42,9 +41,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(undef,N),zero(UInt32)),zero(UInt32))
end
SobolSeq(N::Integer) = SobolSeq(Int(N))
SobolSeq(args...) = SobolSeq{Vector{UInt32}}(args...)

function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat})
length(x) != ndims(s) && throw(BoundsError())
Expand Down Expand Up @@ -72,7 +71,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 +94,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 +114,21 @@ Base.IteratorEltype(::Type{<:AbstractSobolSeq}) = Base.HasEltype()

# Convenience wrapper for scaled Sobol sequences

struct ScaledSobolSeq{N,T} <: AbstractSobolSeq{N}
s::SobolSeq{N}
struct ScaledSobolSeq{T,X<:AbstractVector{UInt32}} <: AbstractSobolSeq
s::SobolSeq{X}
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)
function ScaledSobolSeq{T,X}(lb::Vector{T}, ub::Vector{T}) where {T,X<:AbstractVector{UInt32}}
length(lb)==length(ub) || throw(DimensionMismatch("lb and ub do not have same length"))
new{T,X}(SobolSeq{X}(length(lb)), lb, ub)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe

Suggested change
new{T,X}(SobolSeq{X}(length(lb)), lb, ub)
x = similar(lb, T) # determine the type of array from lb
new{T,typeof{x}}(x, lb, ub)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But if you want to have a much greater level of flexibility, then this is overly restrictive - the bounds could very well be immutable SVectors, but x would have to be its mutable complement.

end
end
function SobolSeq(N::Integer, lb, ub)
function SobolSeq{X}(N::Integer, lb, ub) where {X<:AbstractVector{UInt32}}
T = typeof(sum(ub) - sum(lb))
ScaledSobolSeq{Int(N),T}(copyto!(Vector{T}(undef,N), lb), copyto!(Vector{T}(undef,N), ub))
ScaledSobolSeq{T,X}(copyto!(Vector{T}(undef,N), lb), copyto!(Vector{T}(undef,N), ub))
end
SobolSeq(lb, ub) = SobolSeq(length(lb), lb, ub)
SobolSeq{X}(lb, ub) where {X<:AbstractVector{UInt32}} = SobolSeq{X}(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 +139,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 +170,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