Pure Julia implementation of SobolSeq #3

Ken-B · 2014-08-20T22:33:25Z

As discussed before my holiday, a pure Julia translation of your code, @stevengj. The motivation came for usage in Optim.jl, which aims to become a pure Julia environment.

The test passed when I last tested it, but today when creating this PR it does not. I tried finding the issue for some time without success.
I noticed that the skip does not produce the same result as the C implementation, but a simply looping over next does.
For many dimensions, some have the same output (usually the first 5 or so) but higher dimensions differ. I'm at a loss where this comes from, but I thought I'd already put it here while I continue investigating, maybe it's something simple you will spot immediately.

stevengj · 2014-08-20T23:20:47Z

src/Sobol.jl

+
+type SobolSeq
+ sdim::Integer #dimension of sequence being generated
+ m::Array{Int,2} #array of size (sdim, WORD_SIZE)


I think these should be Uint32, not Int; I'm not sure the algorithm will work for 64-bit integers.

Unfortunately, you then have to be very careful about type-stability, because operations on smaller integer types currently yield Int results in Julia. This may change in 0.4 with JuliaLang/julia#8028

Actually, we should check the original papers on the 64-bit issue. My worry is that the "seed" data needs to be different for 64-bit integers. If that's not the case, Int is fine.

Ken-B · 2014-08-23T10:09:19Z

Thank you very much for the detailed constructive feedback. I've forced Uint32 everywhere (also ndims for no particular reason), devectorized those assignments and changed the skip.

The test still does not pass. Starting for ndims = 4 the results change.
With the C code:

julia> using Sobol
julia> s=SobolSeq(4)
#4-dimensional Sobol sequence on [0,1]^4
julia> show(next(s))
#[0.5,0.5,0.5,0.5]
julia> show(next(s))
#[0.75,0.25,0.75,0.25]

With the pure Julia translation:

julia> using Sobol
julia> s=SobolSeq(4)
#4-dimensional Sobol sequence on [0,1]^4
julia> show(next(s))
#[0.5,0.5,0.5,0.5]
julia> show(next(s))
#[0.75,0.25,0.75,0.75]

The last element is different, but the first 3 elements continue to be equal between the two versions. That makes me think the issue is with generating the initial m matrix.
I'll continue to investigate.

Ken-B · 2014-08-23T12:18:55Z

Even worse, for high dimensions the number generator goes out of bound. This starts around dimenions 90:

brkpnt = 0
for d = 2:1000
    s=SobolSeq(d)
    for unused = 1:100;next(s);end
    a=next(s)
    if any([a[i]>1 for i=1:length(a)])
        brkpnt = d
        break
    end
end
brkpnt
#90

stevengj · 2014-08-23T12:27:33Z

src/Sobol.jl

+ ac = a
+ m[i,j] = m[i,j-d]
+ for k = 0:d-1
+ m[i,j] $= ((ac & 1) * m[i, j-d+k]) << (d-k)


This * is problematic. Remember that Julia currently promotes the result of every operation to an Int (which is usually 64-bit), whereas we need 32-bit multiplication here. (Nor is it correct to recast the result after the multiplication.)

Ken-B · 2014-08-23T12:41:36Z

I didn't fully understand your original comment and this point goes beyond my knowledge.
Do you mean that the multiplication will not overflow as expected on an Uint32?

Is there a workaround you can recommend? Anything I can do or look up from my side?

stevengj · 2014-08-25T19:07:33Z

Yes, I mean that multiplication will not overflow as expected if you multiply two uint32s.

Actually, now that I look back at the original papers, it seems at first glance like it shouldn't be restricted to the 32-bit case. So, the easiest thing to do is probably to change all of your Uint32 types into Uint.

Ken-B · 2014-08-26T16:00:25Z

Thanks for your comments.
I changed all Uin32 to Uint. The good thing is that seems to produce the same results on julia 32bit and julia 64bit (windows binaries), but still the wrong results and the test still fail.

How it seems now is that it works up to dimension 3, but starts diverging from your C code results at dimension 4 and goes out of bounds at least at dimension 90.

I'm at a loss how to continue. Suggestions?

Ken-B · 2015-09-17T14:45:15Z

@stevengj
I thought with v0.4 coming out I'd have another go at this. With a fresh pair of eyes I found the bug, and it's really embarrassing. Oh, boy. It was on the last line in soboldata.jl (c5032e9).

Now the tests pass (travis fails because it's on v0.3).

It also works for UInt on 64-bit, but I had to add an Int(d) for array indexing (here), while that was not necessary when d was UInt32. Very strange.

There are still some depreciation warnings for array concatenation, shall I change the , to ; in soboldata.jl?

stevengj · 2015-09-17T15:19:00Z

src/Sobol.jl

@@ -1,67 +1,93 @@
 module Sobol



You'll need a using Compat for Julia 0.3 compatibility

stevengj · 2015-09-17T15:41:45Z

Thanks for keeping at it, and I'm glad you found the bug! It would be very useful to have a performance benchmark against the C implementation as well.

stevengj · 2015-09-17T15:45:32Z

src/Sobol.jl

+
+ #set initial values of m from table
+ # TODO transpose sobol_minit for faster column access?
+ m[i, 1:Int(d)] = sobol_minit[i-1, 1:Int(d)]


Just make d = floor(Int, log2(a)) and then you don't need the subsequent Int(d) calls.

Ken-B · 2015-09-17T16:35:11Z

Thanks for all the comments.

I'll try to make it work for v0.3 with Compat, but it might be that operators on UInt32 won't be stable or won't overflow as expected. We'll see.

Good idea also for a speed comparison, I'll just use the tests.

Ken-B · 2015-09-17T16:44:36Z

The latest changes follow your feedback and now it works on v0.3!

Ken-B · 2015-09-17T16:48:43Z

Don't merge yet! With the latest changes the tests fail on v0.4. Let me check.

Ken-B · 2015-09-17T17:06:36Z

Actually, the tests pass if I run Pkg.test("Sobol"), but they fail if I run include("runtests.jl"). I'm investigating.

stevengj · 2015-09-17T17:09:42Z

src/Sobol.jl

- (Ptr{Void}, Ptr{Float64}), s.p, x)
- x
+# Special case N == 0
+next(s::SobolSeq{0}) = Float64[]


This one seems redundant with next(s::SobolSeq) = next!(s, Array(Float64,ndims(s)))

Ken-B · 2015-09-17T17:37:42Z

Oops, nevermind about the tests failing, there was some previous stuff in my workspace. Restarting julia fixed it. Sorry about that.

For the benchmark: running the tests (second run) on my macbook pro with julia v0.4, both pure julia
and your current master (in C) take about 175.1s. I guess they compile to the same code.

I've just added a precompile line. This PR seems finished.

stevengj · 2015-09-17T17:56:19Z

Yes, it makes sense that they would be the same performance, but it's great to see it confirmed.

stevengj · 2015-09-17T17:57:54Z

src/Sobol.jl

+ return nothing
+ end
+
+ s.n += 1


probably += one(s.n) to avoid promoting, adding, then truncating.

Well spotted!

stevengj · 2015-09-17T18:01:54Z

Note that you should also:

Update the README to indicated that it was ported to pure Julia (by you) from the C implementation in NLopt.
squash the commits into a single commit (with git rebase -i followed by git push -f)

@stevengj

adaptation after @stevengj feedback: force Uint32; devectorize some assignments; change skip number Uint32 -> Uint tryout make it work for julia v0.4 and UInt32 support 64bit, make it work for general UInt Revert "support 64bit, make it work for general UInt" This reverts commit 94099ad. make sdims a type parameter, allow sdims==0, compatibility with v0.3 and minor changes from comments add recompilation remove special case next for N==0 update README and small incrementer fix

Ken-B · 2015-09-17T22:02:19Z

Wow, that was my first rebase and force push, pretty cool :)

Did it all arrive well?

stevengj · 2015-09-18T00:01:47Z

README.md

-and takes advantage of the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl)
-Julia module to ensure that the NLopt library is installed when you
-`Pkg.add("Sobol")` in Julia.
+[NLopt library](http:https://ab-initio.mit.edu/nlopt) for nonlinear optimization, it has been ported by [Ken-B](https://github.com/Ken-B) into pure Julia code with the same performance.


has been -> was. (Past tense, not present perfect.)

stevengj · 2015-09-18T00:02:12Z

Looking good!

Pure Julia implementation of SobolSeq

stevengj · 2015-09-18T00:09:29Z

README.md

@@ -28,13 +28,9 @@ described in:
 **14** (1), pp. 88-100 (1988).
 * S. Joe and F. Y. Kuo, *ACM Trans. Math. Soft* **29** (1), 49-57 (2003).

-Our implementation is in C, and was implemented in 2007 by SGJ as
+Originally implementated in C in 2007 as


implementated -> implemented

stevengj · 2015-09-18T00:14:27Z

Now that I look at it, the runtests script is not the best benchmark because it includes the time to evaluate the integrand. Probably something like this would be better:

function bench(N, n)
    s = SobolSeq(N)
    x = Array(Float64, N)
    skip(s, n)
    return @elapsed for j = 1:n
            next!(s, x)
    end
end

stevengj · 2015-09-18T00:21:14Z

Okay, looks like the performance is suboptimal for larger N. C implementation:

julia> bench(3, 2^18)
0.008534276

julia> bench(30, 2^18)
0.031220722

julia> bench(300, 2^18)
0.28640663

Pure-Julia implementation:

julia> bench(3, 2^18)
0.008328307

julia> bench(30, 2^18)
0.040762087

julia> bench(300, 2^18)
0.372826312

stevengj · 2015-09-18T00:31:14Z

Improved the performance somewhat. After f8fbe05, I now get:

julia> bench(3, 2^18)
0.008544011

julia> bench(30, 2^18)
0.03974193

julia> bench(300, 2^18)
0.279880779

which is comparable to the C version. However:

julia> @time bench(300, 2^18)
  0.725053 seconds (263.65 k allocations: 620.109 MB, 11.53% gc time)
0.273653861

involves a disturbing fraction of gc time. Something is doing an allocation somewhere that it shouldn't. I get the same allocations even if I take out the allocation of s and x with:

function bench(s, x, N, n)
           skip(s, n)
           return @elapsed for j = 1:n
                   next!(s, x)
           end
end

stevengj · 2015-09-18T00:53:58Z

Oh, all the allocations are in skip, which calls next repeatedly instead of next!.

Better to benchmark without the skip, in which case I get:

julia> @time bench(s, x, 300, 2^18)
  0.283078 seconds (6 allocations: 192 bytes)
0.283066224

which looks good.

Ken-B · 2015-09-18T10:24:26Z

Thanks for merging this and all the careful tweaks afterwards. Especially thanks for the helpful comments, this has been very instructive for me!
Let me know if I can do anything else.

stevengj reviewed Aug 20, 2014
View reviewed changes

stevengj reviewed Aug 23, 2014
View reviewed changes

stevengj reviewed Sep 17, 2015
View reviewed changes

src/Sobol.jl

@@ -1,67 +1,93 @@

module Sobol

Copy link

Member

stevengj Sep 17, 2015

Choose a reason for hiding this comment

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

You'll need a using Compat for Julia 0.3 compatibility

stevengj reviewed Sep 17, 2015
View reviewed changes

Ken-B force-pushed the SobolJulia branch from 216bdd9 to b17815e Compare September 17, 2015 22:00

stevengj reviewed Sep 18, 2015
View reviewed changes

stevengj added a commit that referenced this pull request Sep 18, 2015

Merge pull request #3 from Ken-B/SobolJulia

28b5fe6

Pure Julia implementation of SobolSeq

stevengj merged commit 28b5fe6 into JuliaMath:master Sep 18, 2015

stevengj reviewed Sep 18, 2015
View reviewed changes

stevengj added a commit that referenced this pull request Sep 18, 2015

make skip faster, as discussed in #3

d4f1402

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Pure Julia implementation of SobolSeq #3

Pure Julia implementation of SobolSeq #3

Ken-B commented Aug 20, 2014

stevengj Aug 20, 2014

stevengj Aug 20, 2014

Ken-B commented Aug 23, 2014

Ken-B commented Aug 23, 2014

stevengj Aug 23, 2014

Ken-B commented Aug 23, 2014

stevengj commented Aug 25, 2014

Ken-B commented Aug 26, 2014

Ken-B commented Sep 17, 2015

stevengj Sep 17, 2015

stevengj commented Sep 17, 2015

stevengj Sep 17, 2015

Ken-B commented Sep 17, 2015

Ken-B commented Sep 17, 2015

Ken-B commented Sep 17, 2015

Ken-B commented Sep 17, 2015

stevengj Sep 17, 2015

Ken-B commented Sep 17, 2015

stevengj commented Sep 17, 2015

stevengj Sep 17, 2015

Ken-B Sep 17, 2015

stevengj commented Sep 17, 2015

Ken-B commented Sep 17, 2015

stevengj Sep 18, 2015

stevengj commented Sep 18, 2015

stevengj Sep 18, 2015

stevengj commented Sep 18, 2015

stevengj commented Sep 18, 2015

stevengj commented Sep 18, 2015

stevengj commented Sep 18, 2015

Ken-B commented Sep 18, 2015

Pure Julia implementation of SobolSeq #3

Pure Julia implementation of SobolSeq #3

Conversation

Ken-B commented Aug 20, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Ken-B commented Aug 23, 2014

Ken-B commented Aug 23, 2014

Choose a reason for hiding this comment

Ken-B commented Aug 23, 2014

stevengj commented Aug 25, 2014

Ken-B commented Aug 26, 2014

Ken-B commented Sep 17, 2015

Choose a reason for hiding this comment

stevengj commented Sep 17, 2015

Choose a reason for hiding this comment

Ken-B commented Sep 17, 2015

Ken-B commented Sep 17, 2015

Ken-B commented Sep 17, 2015

Ken-B commented Sep 17, 2015

Choose a reason for hiding this comment

Ken-B commented Sep 17, 2015

stevengj commented Sep 17, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj commented Sep 17, 2015

Ken-B commented Sep 17, 2015

Choose a reason for hiding this comment

stevengj commented Sep 18, 2015

Choose a reason for hiding this comment

stevengj commented Sep 18, 2015

stevengj commented Sep 18, 2015

stevengj commented Sep 18, 2015

stevengj commented Sep 18, 2015

Ken-B commented Sep 18, 2015