Add rng argument, fixes #44 (#58)

* Add rng parameter, fixes #44

* Remove now-unused method

src/optimize.jl

@@ -1,30 +1,31 @@
 ## Optimization =======================================
 
-function perturb(p, Cp)
+function perturb(rng::AbstractRNG, p, Cp)
     d = MvNormal(mean(p), 1.1Cp + 1e-12I)
-    Particles(length(p[1]), d)
+    Particles(rng, length(p[1]), d)
 end
 
 
 
 """
-    res = optimize(f,p,τ=1,iters=10000)
+    res = optimize([rng::AbstractRNG,] f,p,τ=1,iters=10000)
 
 Find the minimum of Function `f`, starting with initial distribution described by `p::Vector{Particles}`. `τ` is the initial temperature.
 """
-function optimize(f,p,τ=1; τi=1.005, iters=10000, tol=1e-8)
+function optimize(rng::AbstractRNG,f,p,τ=1; τi=1.005, iters=10000, tol=1e-8)
     p = deepcopy(p)
     N = length(p[1])
     we = zeros(N)
     for i = 1:iters
         y  = -(f(p).particles)
         we .= exp.(τ.*y)
-        j  = sample(1:N, ProbabilityWeights(we), N, ordered=true)
+        j  = sample(rng,1:N, ProbabilityWeights(we), N, ordered=true)
         foreach(x->(x.particles .= x.particles[j]), p); # @test length(unique(p[1].particles)) == length(unique(j))
         Cp = cov(p)
         tr(Cp) < tol && (@info "Converged at iteration $i"; return p)
-        p = perturb(p, Cp)
+        p = perturb(rng, p, Cp)
         τ *= τi
     end
     p
 end
+optimize(f,p,τ=1; kwargs...) = optimize(Random.GLOBAL_RNG,f,p,τ; kwargs...)

src/particles.jl

@@ -38,24 +38,27 @@ See also [`outer_product`](@ref), [`±`](@ref)
 ⊗(μ,σ) = outer_product(Normal.(μ,σ))
 
 """
-    p = outer_product(dists::Vector{<:Distribution}, N=100_000)
+    p = outer_product([rng::AbstractRNG,] dists::Vector{<:Distribution}, N=100_000)
 
 Creates a multivariate systematic sample where each dimension is sampled according to the corresponding univariate distribution in `dists`. Returns `p::Vector{Particles}` where each Particles has a length approximately equal to `N`.
 The particles form the outer product between `d` systematically sampled vectors with length given by the d:th root of N, where `d` is the length of `dists`, All particles will be independent and have marginal distributions given by `dists`.
 
 See also `MonteCarloMeasurements.⊗`
 """
-function outer_product(dists::AbstractVector{<:Distribution}, N=100_000)
+function outer_product(rng::AbstractRNG, dists::AbstractVector{<:Distribution}, N=100_000)
     d = length(dists)
     N = floor(Int,N^(1/d))
     dims = map(dists) do dist
-        v = systematic_sample(N,dist; permute=true)
+        v = systematic_sample(rng,N,dist; permute=true)
     end
     cart_prod = vec(collect(Iterators.product(dims...)))
     p = map(1:d) do i
         Particles(getindex.(cart_prod,i))
     end
 end
+function outer_product(dists::AbstractVector{<:Distribution}, N=100_000)
+    return outer_product(Random.GLOBAL_RNG, dists, N)
+end
 
 # StaticParticles(N::Integer = DEFAULT_NUM_PARTICLES; permute=true) = StaticParticles{Float64,N}(SVector{N,Float64}(systematic_sample(N, permute=permute)))
 

src/resampling.jl

@@ -27,7 +27,7 @@
 #     w[i] += 1
 #     s
 # end
-# 
+#
 # """
 #     loglik = resample!(p::WeightedParticles)
 # Resample the particles based on the `p.logweights`. After a call to this function, weights will be reset to sum to one. Returns log-likelihood.
@@ -67,10 +67,11 @@
 # end
 
 """
-    bootstrap(p::Particles)
+    bootstrap([rng::AbstractRNG,] p::Particles)
 
 Return Particles resampled with replacement.
 """
-function bootstrap(p::T) where T <: AbstractParticles
-    T(p.particles[rand(1:length(p))])
+function bootstrap(rng::AbstractRNG, p::T) where T <: AbstractParticles
+    T(p.particles[rand(rng, 1:length(p))])
 end
+bootstrap(p::T) where T <: AbstractParticles = bootstrap(Random.GLOBAL_RNG, p)

src/sampling.jl

@@ -1,13 +1,17 @@
 """
-    systematic_sample(N, d=Normal(0,1); permute=true)
+    systematic_sample([rng::AbstractRNG,] N, d=Normal(0,1); permute=true)
+
 returns a `Vector` of length `N` sampled systematically from the distribution `d`. If `permute=false`, this vector will be sorted.
 """
-function systematic_sample(N, d=Normal(0,1); permute=true)
+function systematic_sample(rng::AbstractRNG, N, d=Normal(0,1); permute=true)
     e   = 0.5/N # rand()/N
     y   = e:1/N:1
     o = map(y) do y
         quantile(d,y)
     end
-    permute && permute!(o, randperm(N))
+    permute && permute!(o, randperm(rng, N))
     o
 end
+function systematic_sample(N, d=Normal(0,1); kwargs...)
+    return systematic_sample(Random.GLOBAL_RNG, N, d; kwargs...)
+end

src/types.jl

@@ -9,8 +9,8 @@ This type represents uncertainty using a cloud of particles.
 # Constructors:
 - `Particles()`
 - `Particles(N::Integer)`
-- `Particles(d::Distribution)`
-- `Particles(N::Integer, d::Distribution; permute=true, systematic=true)`
+- `Particles([rng::AbstractRNG,] d::Distribution)`
+- `Particles([rng::AbstractRNG,] N::Integer, d::Distribution; permute=true, systematic=true)`
 - `Particles(v::Vector{T} where T)`
 - `Particles(m::Matrix{T} where T)`: Creates multivariate particles (Vector{Particles})
 """
@@ -56,33 +56,38 @@ for PT in (:Particles, :StaticParticles)
 
         $PT{T,N}(p::$PT{T,N}) where {T,N} = p
 
-        function $PT(N::Integer=DEFAULT_NUM_PARTICLES, d::Distribution{<:Any,VS}=Normal(0,1); permute=true, systematic=VS==Continuous) where VS
+        function $PT(rng::AbstractRNG, N::Integer=DEFAULT_NUM_PARTICLES, d::Distribution{<:Any,VS}=Normal(0,1); permute=true, systematic=VS==Continuous) where VS
             if systematic
-                v = systematic_sample(N,d; permute=permute)
+                v = systematic_sample(rng,N,d; permute=permute)
             else
-                v = rand(d, N)
+                v = rand(rng, d, N)
             end
             $PT{eltype(v),N}(v)
         end
+        function $PT(N::Integer=DEFAULT_NUM_PARTICLES, d::Distribution{<:Any,VS}=Normal(0,1); kwargs...) where VS
+            return $PT(Random.GLOBAL_RNG, N, d; kwargs...)
+        end
 
-
-        function $PT(N::Integer, d::MultivariateDistribution)
-            v = rand(d,N)' |> copy # For cache locality
+        function $PT(rng::AbstractRNG, N::Integer, d::MultivariateDistribution)
+            v = rand(rng,d,N)' |> copy # For cache locality
             $PT(v)
         end
+        $PT(N::Integer, d::MultivariateDistribution) = $PT(Random.GLOBAL_RNG, N, d)
 
         nakedtypeof(p::$PT{T,N}) where {T,N} = $PT
         nakedtypeof(::Type{$PT{T,N}}) where {T,N} = $PT
     end
 end
 
-function Particles(d::Distribution;kwargs...)
-    Particles(DEFAULT_NUM_PARTICLES, d; kwargs...)
+function Particles(rng::AbstractRNG, d::Distribution;kwargs...)
+    Particles(rng, DEFAULT_NUM_PARTICLES, d; kwargs...)
 end
+Particles(d::Distribution; kwargs...) = Particles(Random.GLOBAL_RNG, d; kwargs...)
 
-function StaticParticles(d::Distribution;kwargs...)
-    StaticParticles(DEFAULT_STATIC_NUM_PARTICLES, d; kwargs...)
+function StaticParticles(rng::AbstractRNG, d::Distribution;kwargs...)
+    StaticParticles(rng, DEFAULT_STATIC_NUM_PARTICLES, d; kwargs...)
 end
+StaticParticles(d::Distribution;kwargs...) = StaticParticles(Random.GLOBAL_RNG, d; kwargs...)
 
 
 const MvParticles = Vector{<:AbstractParticles} # This can not be AbstractVector since it causes some methods below to be less specific than desired

test/runtests.jl

@@ -158,6 +158,24 @@ Random.seed!(0)
                 @test pn ≈ 2
                 @test issorted(pn.particles)
 
+                rng = MersenneTwister(657)
+                pn1 = Particles(rng, 100, Normal(2,1), systematic=true, permute=true)
+                rng = MersenneTwister(657)
+                pn2 = Particles(rng, 100, Normal(2,1), systematic=true, permute=true)
+                @test pn1 == pn2
+
+                rng = MersenneTwister(27)
+                pn1 = Particles(rng, Normal(2,1), systematic=true, permute=true)
+                rng = MersenneTwister(27)
+                pn2 = Particles(rng, Normal(2,1), systematic=true, permute=true)
+                @test pn1 == pn2
+
+                rng = MersenneTwister(932)
+                pn1 = Particles(rng, 100, systematic=true, permute=true)
+                rng = MersenneTwister(932)
+                pn2 = Particles(rng, 100, systematic=true, permute=true)
+                @test pn1 == pn2
+
                 @info "Tests for $PT done"
 
                 p = PT{Float64,10}(2)
@@ -403,6 +421,11 @@ Random.seed!(0)
         @test wasserstein(p,p,1) == 0
         @test wasserstein(p,q,1) >= 0
         @test bootstrap(p) ≈ p
+        rng = MersenneTwister(453)
+        p1 = bootstrap(rng,p)
+        rng = MersenneTwister(453)
+        p2 = bootstrap(rng,p)
+        @test p1 == p2
     end
 
     @time @testset "mutation" begin
@@ -435,6 +458,12 @@ Random.seed!(0)
         @test length(p) == 2
         @test length(p[1]) <= 100_000
         @test cov(p) ≈ I atol=1e-1
+
+        rng = MersenneTwister(38)
+        p1 = outer_product(rng, Normal.(μ,σ))
+        rng = MersenneTwister(38)
+        p2 = outer_product(rng, Normal.(μ,σ))
+        @test p1 == p2
     end
 
     @time @testset "plotting" begin
@@ -467,6 +496,13 @@ Random.seed!(0)
             popt = optimize(rosenbrock2d, deepcopy(p))
             popt ≈ [1,1]
         end
+
+        p = -1ones(2) .+ 2 .*Particles.(200) # Optimum is in [1,1]
+        rng = MersenneTwister(876)
+        popt1 = optimize(rng, rosenbrock2d, deepcopy(p))
+        rng = MersenneTwister(876)
+        popt2 = optimize(rng, rosenbrock2d, deepcopy(p))
+        @test popt1 == popt2
     end