From 8615b6538b61558650157d91d32b7a6394e4c036 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 1 Feb 2020 02:24:11 -0800 Subject: [PATCH] Add rng argument, fixes #44 (#58) * Add rng parameter, fixes #44 * Remove now-unused method --- src/optimize.jl | 13 +++++++------ src/particles.jl | 9 ++++++--- src/resampling.jl | 9 +++++---- src/sampling.jl | 10 +++++++--- src/types.jl | 29 +++++++++++++++++------------ test/runtests.jl | 36 ++++++++++++++++++++++++++++++++++++ 6 files changed, 78 insertions(+), 28 deletions(-) diff --git a/src/optimize.jl b/src/optimize.jl index 07e7dfd5..39b34069 100644 --- a/src/optimize.jl +++ b/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...) diff --git a/src/particles.jl b/src/particles.jl index cc6b6d8b..3ab5ecba 100644 --- a/src/particles.jl +++ b/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))) diff --git a/src/resampling.jl b/src/resampling.jl index a1538184..3e28b485 100644 --- a/src/resampling.jl +++ b/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) diff --git a/src/sampling.jl b/src/sampling.jl index 50ef4a69..093348c1 100644 --- a/src/sampling.jl +++ b/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 diff --git a/src/types.jl b/src/types.jl index b3c1a634..987e8c68 100644 --- a/src/types.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 4598f9fb..6cd01fe7 100644 --- a/test/runtests.jl +++ b/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