Fix constructor of ScaledSobolSeq (#17)

* fix bug in constructors for ScaledSobolSeq


add next!(s::ScaledSobolSeq) to allow for iteration over ScaledSobolSeq

* use type params in inner constructor, more tests, users don't call ScaledSobolSeq directly, support AbstractVector

* copyto! for Julia 0.7 and 1.0

* don't make extra copies if data is Vector{Float64} already

* support SobolSeq(lb, ub) for iterators

* safer to always make a copy

Author: moritzdrechselgrau

README.md

@@ -45,7 +45,7 @@ facts it is not copyrightable.)  SGJ's implementation in NLopt, along
 with this Julia translation, is free/open-source software under the [MIT
 ("expat") license](http://opensource.org/licenses/MIT).
 
-Direction numbers used were derived from the file 
+Direction numbers used were derived from the file
 http://web.maths.unsw.edu.au/~fkuo/sobol/new-joe-kuo-6.21201
 
 Technically, we implement a 32-bit Sobol sequence.  After
@@ -84,13 +84,13 @@ Note, however, that the loop will *never terminate* unless you explicitly
 call `break` (or similar) in the loop body at some point of your choosing.
 
 We also provide a different `SobolSeq` constructor to provide
-a Sobol sequence rescaled to an arbitrary hypercube:
+an `N`-dimensional Sobol sequence rescaled to an arbitrary hypercube:
 ```
-s = SobolSeq(N, lb, ub)
+s = SobolSeq(lb, ub)
 ```
 where `lb` and `ub` are arrays (or other iterables) of length `N`, giving
 the lower and upper bounds of the hypercube, respectively.   For example,
-`SobolSeq(2, [-1,0,0],[1,3,2])` generates points in the box [-1,1]×[0,3]×[0,2].
+`SobolSeq([-1,0,0],[1,3,2])` generates points in the box [-1,1]×[0,3]×[0,2].
 
 If you know in advance the number `n` of points that you plan to
 generate, some authors suggest that better uniformity can be attained

src/Sobol.jl

@@ -122,11 +122,14 @@ struct ScaledSobolSeq{N} <: AbstractSobolSeq{N}
     s::SobolSeq{N}
     lb::Vector{Float64}
     ub::Vector{Float64}
-    ScaledSobolSeq(lb::Vector{Float64}, ub::Vector{Float64}) =
-        new{N}(SobolSeq(N), lb, ub)
+    function ScaledSobolSeq{N}(lb::Vector{Float64}, ub::Vector{Float64}) where {N}
+        length(lb)==length(ub)==N || throw(DimensionMismatch("lb and ub do not have length $N"))
+        new(SobolSeq(N), lb, ub)
+    end
 end
 SobolSeq(N::Integer, lb, ub) =
-    ScaledSobolSeq{Int(N)}(copy!(Vector{Float64}(undef,N), lb), copy!(Vector{Float64}(undef,N), ub))
+    ScaledSobolSeq{Int(N)}(copyto!(Vector{Float64}(undef,N), lb), copyto!(Vector{Float64}(undef,N), ub))
+SobolSeq(lb, ub) = SobolSeq(length(lb), lb, ub)
 
 function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat},
                lb::AbstractVector, ub::AbstractVector)
@@ -140,6 +143,7 @@ end
 next!(s::SobolSeq{N}, lb::AbstractVector, ub::AbstractVector) where {N} = next!(s, Vector{Float64}(undef, N), lb, ub)
 
 next!(s::ScaledSobolSeq, x::AbstractVector{<:AbstractFloat}) = next!(s.s, x, s.lb, s.ub)
+next!(s::ScaledSobolSeq) = next!(s.s, Array{Float64,1}(undef, ndims(s)), s.lb, s.ub)
 next(s::ScaledSobolSeq) = next!(s, Vector{Float64}(undef, ndims(s)))
 
 Base.skip(s::ScaledSobolSeq, n) = skip(s.s, n)

test/runtests.jl

@@ -1,33 +1,50 @@
 using Sobol, Compat
 using Compat.Test
 
-# compare results with results from C++ code sobol.cc published on 
+# compare results with results from C++ code sobol.cc published on
 # http://web.maths.unsw.edu.au/~fkuo/sobol/
 # with new-joe-kuo-6.21201 file used as input.
 #
-# Command line used to generate output is 
+# Command line used to generate output is
 #  ./sobol $N 1024 new-joe-kuo-6.21201 > exp_results_$N
 #
 # For each of the dimensions below
 # (except 21201 where only 64 samples are generated)
 
-dimensions = [5, 10, 20, 50, 100, 500, 21201]
+@testset "published results" begin
+    dimensions = [5, 10, 20, 50, 100, 500, 21201]
 
-for dim in dimensions
-    println("Testing dimension $(dim)")
-    open(joinpath(dirname(@__FILE__), "results", "exp_results_$(dim)")) do exp_results_file
-        s = SobolSeq(dim)
-        x = zeros(dim)
-        for line in eachline(exp_results_file)
-            values = [parse(Float64, item) for item in split(line)]
-            if length(values) > 0
-                @test x == values
-                next!(s, x)
+    for dim in dimensions
+        println("Testing dimension $(dim)")
+        open(joinpath(dirname(@__FILE__), "results", "exp_results_$(dim)")) do exp_results_file
+            s = SobolSeq(dim)
+            x = zeros(dim)
+            for line in eachline(exp_results_file)
+                values = [parse(Float64, item) for item in split(line)]
+                if length(values) > 0
+                    @test x == values
+                    next!(s, x)
+                end
             end
         end
     end
 end
 
-# issue #8
 using Base.Iterators: take
-@test [x[1] for x in collect(take(Sobol.SobolSeq(1),5))] == [0.5,0.75,0.25,0.375,0.875]
+@testset "iterators" begin
+    # issue #8
+    @test [x[1] for x in collect(take(Sobol.SobolSeq(1),5))] == [0.5,0.75,0.25,0.375,0.875]
+end
+
+@testset "scaled" begin
+    # ScaledSobolSeq constructors
+    lb = [-1,0,0]
+    ub = [1,3,2]
+    N = length(lb)
+    s = SobolSeq(lb,ub)
+    @test s isa ScaledSobolSeq{3}
+    @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}
+    @test_throws BoundsError SobolSeq(2,lb,ub)
+end