store RNG instead of seed

allows to customize `DEFAULT_RNG[]` on a global basis to use
stable RNGs

Author: hexaeder

Project.toml

@@ -9,26 +9,28 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
+[weakdeps]
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+
+[extensions]
+NetworkLayoutGraphsExt = "Graphs"
+
 [compat]
 GeometryBasics = "0.4"
 Graphs = "1"
 Requires = "1"
+StableRNGs = "1.0.2"
 StaticArrays = "1"
 julia = "1.6"
 
-[extensions]
-NetworkLayoutGraphsExt = "Graphs"
-
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Graphs", "Test", "DelimitedFiles", "GeometryTypes", "SparseArrays"]
-
-[weakdeps]
-Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+test = ["Graphs", "Test", "DelimitedFiles", "GeometryTypes", "SparseArrays", "StableRNGs"]

src/NetworkLayout.jl

@@ -7,6 +7,11 @@ using StaticArrays
 
 export LayoutIterator
 
+"""
+Default RNG for layouts.
+"""
+const DEFAULT_RNG = Ref{DataType}(MersenneTwister)
+
 """
     AbstractLayout{Dim,Ptype}
 

src/sfdp.jl

@@ -33,36 +33,40 @@ the nodes.
     - `(true, false, false)` : only pin certain coordinates
 
 - `seed=1`: Seed for random initial positions.
+- `rng=DEFAULT_RNG[](seed)`
+
+  Create rng based on seed. Defaults to `MersenneTwister`, can be specified
+  by overwriting `DEFAULT_RNG[]`
 """
-@addcall struct SFDP{Dim,Ptype,T<:AbstractFloat} <: IterativeLayout{Dim,Ptype}
+@addcall struct SFDP{Dim,Ptype,T<:AbstractFloat,RNG} <: IterativeLayout{Dim,Ptype}
     tol::T
     C::T
     K::T
     iterations::Int
     initialpos::Dict{Int,Point{Dim,Ptype}}
     pin::Dict{Int,SVector{Dim,Bool}}
-    seed::UInt
+    rng::RNG
 end
 
 # TODO: check SFDP default parameters
 function SFDP(; dim=2, Ptype=Float64,
               tol=1.0, C=0.2, K=1.0,
               iterations=100,
               initialpos=[], pin=[],
-              seed=1)
+              seed=1, rng=DEFAULT_RNG[](seed))
     if !isempty(initialpos)
         dim, Ptype = infer_pointtype(initialpos)
         Ptype = promote_type(Float32, Ptype) # make sure to get at least f32 if given as int
     end
     _initialpos, _pin = _sanitize_initialpos_pin(dim, Ptype, initialpos, pin)
 
-    return SFDP{dim,Ptype,typeof(tol)}(tol, C, K, iterations, _initialpos, _pin, seed)
+    return SFDP{dim,Ptype,typeof(tol),typeof(rng)}(tol, C, K, iterations, _initialpos, _pin, rng)
 end
 
-function Base.iterate(iter::LayoutIterator{SFDP{Dim,Ptype,T}}) where {Dim,Ptype,T}
+function Base.iterate(iter::LayoutIterator{<:SFDP{Dim,Ptype,T}}) where {Dim,Ptype,T}
     algo, adj_matrix = iter.algorithm, iter.adj_matrix
     N = size(adj_matrix, 1)
-    rng = MersenneTwister(algo.seed)
+    rng = copy(algo.rng)
     startpos = [2 .* rand(rng, Point{Dim,Ptype}) .- 1 for _ in 1:N]
 
     for (k, v) in algo.initialpos
@@ -106,8 +110,8 @@ function Base.iterate(iter::LayoutIterator{<:SFDP}, state)
         if any(isnan, force)
             # if two points are at the exact same location
             # use random force in any direction
-            rng = MersenneTwister(algo.seed + i)
-            force = randn(rng, Ftype)
+            rng = copy(algo.rng)
+            force += randn(rng, Ftype)
         end
         mask = (!).(pin[i]) # where pin=true mask will multiply with 0
         locs[i] = locs[i] .+ (step .* (force ./ norm(force))) .* mask

src/spring.jl

@@ -35,33 +35,37 @@ the nodes.
     - `(true, false, false)` : only pin certain coordinates
 
 - `seed=1`: Seed for random initial positions.
+- `rng=DEFAULT_RNG[](seed)`
+
+  Create rng based on seed. Defaults to `MersenneTwister`, can be specified
+  by overwriting `DEFAULT_RNG[]`
 """
-@addcall struct Spring{Dim,Ptype} <: IterativeLayout{Dim,Ptype}
+@addcall struct Spring{Dim,Ptype,RNG} <: IterativeLayout{Dim,Ptype}
     C::Float64
     iterations::Int
     initialtemp::Float64
     initialpos::Dict{Int,Point{Dim,Ptype}}
     pin::Dict{Int,SVector{Dim,Bool}}
-    seed::UInt
+    rng::RNG
 end
 
 function Spring(; dim=2, Ptype=Float64,
                 C=2.0, iterations=100, initialtemp=2.0,
                 initialpos=[], pin=[],
-                seed=1)
+                seed=1, rng=DEFAULT_RNG[](seed))
     if !isempty(initialpos)
         dim, Ptype = infer_pointtype(initialpos)
         Ptype = promote_type(Float32, Ptype) # make sure to get at least f32 if given as int
     end
     _initialpos, _pin = _sanitize_initialpos_pin(dim, Ptype, initialpos, pin)
 
-    return Spring{dim,Ptype}(C, iterations, initialtemp, _initialpos, _pin, seed)
+    return Spring{dim,Ptype,typeof(rng)}(C, iterations, initialtemp, _initialpos, _pin, rng)
 end
 
 function Base.iterate(iter::LayoutIterator{<:Spring{Dim,Ptype}}) where {Dim,Ptype}
     algo, adj_matrix = iter.algorithm, iter.adj_matrix
     N = size(adj_matrix, 1)
-    rng = MersenneTwister(algo.seed)
+    rng = copy(algo.rng)
     startpos = [2 .* rand(rng, Point{Dim,Ptype}) .- 1 for _ in 1:N]
 
     for (k, v) in algo.initialpos
@@ -110,7 +114,7 @@ function Base.iterate(iter::LayoutIterator{<:Spring}, state)
             else
                 # if two points are at the exact same location
                 # use random force in any direction
-                rng = MersenneTwister(algo.seed + i)
+                rng = copy(algo.rng)
                 force_vec += randn(rng, Ftype)
             end
 

src/stress.jl

@@ -53,8 +53,12 @@ The main equation to solve is (8) in Gansner, Koren and North (2005,
     - `(true, false, false)` : only pin certain coordinates
 
 - `seed=1`: Seed for random initial positions.
+- `rng=DEFAULT_RNG[](seed)`
+
+  Create rng based on seed. Defaults to `MersenneTwister`, can be specified
+  by overwriting `DEFAULT_RNG[]`
 """
-@addcall struct Stress{Dim,Ptype,IT<:Union{Symbol,Int},FT<:AbstractFloat,M<:AbstractMatrix} <:
+@addcall struct Stress{Dim,Ptype,IT<:Union{Symbol,Int},FT<:AbstractFloat,M<:AbstractMatrix,RNG} <:
                 IterativeLayout{Dim,Ptype}
     iterations::IT
     abstols::FT
@@ -63,7 +67,7 @@ The main equation to solve is (8) in Gansner, Koren and North (2005,
     weights::M
     initialpos::Dict{Int,Point{Dim,Ptype}}
     pin::Dict{Int,SVector{Dim,Bool}}
-    seed::UInt
+    rng::RNG
 end
 
 function Stress(; dim=2,
@@ -74,23 +78,23 @@ function Stress(; dim=2,
                 abstolx=10e-6,
                 weights=Array{Float64}(undef, 0, 0),
                 initialpos=[], pin=[],
-                seed=1)
+                seed=1, rng=DEFAULT_RNG[](seed))
     if !isempty(initialpos)
         dim, Ptype = infer_pointtype(initialpos)
         Ptype = promote_type(Float32, Ptype) # make sure to get at least f32 if given as int
     end
 
     _initialpos, _pin = _sanitize_initialpos_pin(dim, Ptype, initialpos, pin)
 
-    IT, FT, WT = typeof(iterations), typeof(abstols), typeof(weights)
-    Stress{dim,Ptype,IT,FT,WT}(iterations, abstols, reltols, abstolx, weights, _initialpos, _pin, seed)
+    IT, FT, WT, RNG = typeof(iterations), typeof(abstols), typeof(weights), typeof(rng)
+    Stress{dim,Ptype,IT,FT,WT,RNG}(iterations, abstols, reltols, abstolx, weights, _initialpos, _pin, rng)
 end
 
 function Base.iterate(iter::LayoutIterator{<:Stress{Dim,Ptype,IT,FT}}) where {Dim,Ptype,IT,FT}
     algo, δ = iter.algorithm, iter.adj_matrix
     N = size(δ, 1)
     M = length(algo.initialpos)
-    rng = MersenneTwister(algo.seed)
+    rng = copy(algo.rng)
     startpos = randn(rng, Point{Dim,Ptype}, N)
 
     for (k, v) in algo.initialpos

test/runtests.jl

@@ -4,7 +4,9 @@ using GeometryBasics
 using DelimitedFiles: readdlm
 using SparseArrays: sparse
 using StaticArrays
+using StableRNGs
 using Test
+using Random
 
 function jagmesh()
     jagmesh_path = joinpath(dirname(@__FILE__), "jagmesh1.mtx")
@@ -63,6 +65,14 @@ jagmesh_adj = jagmesh()
             positions = @time SFDP(; dim=3, Ptype=Float32, tol=0.1, K=1)(adj_matrix)
             @test typeof(positions) == Vector{Point3f}
             @test positions == sfdp(adj_matrix; dim=3, Ptype=Float32, tol=0.1, K=1)
+
+            NetworkLayout.DEFAULT_RNG[] = StableRNG
+            l = SFDP()
+            @test l.rng isa StableRNG
+            li = LayoutIterator(l, g)
+            p1, p2 = iterate(li)[1], iterate(li)[1]
+            @test p1 == p2
+            NetworkLayout.DEFAULT_RNG[] = MersenneTwister
         end
     end
 
@@ -112,6 +122,14 @@ jagmesh_adj = jagmesh()
             positions = @time Stress(; iterations=10, dim=3, Ptype=Float32)(adj_matrix)
             @test typeof(positions) == Vector{Point3f}
             @test positions == stress(adj_matrix; iterations=10, dim=3, Ptype=Float32)
+
+            NetworkLayout.DEFAULT_RNG[] = StableRNG
+            l = Stress()
+            @test l.rng isa StableRNG
+            li = LayoutIterator(l, g)
+            p1, p2 = iterate(li)[1], iterate(li)[1]
+            @test p1 == p2
+            NetworkLayout.DEFAULT_RNG[] = MersenneTwister
         end
 
         @testset "test pairwise_distance" begin
@@ -170,6 +188,14 @@ jagmesh_adj = jagmesh()
             @test typeof(positions) == Vector{Point3f}
             @test positions ==
                   spring(adj_matrix; C=2.0, iterations=100, initialtemp=2.0, Ptype=Float32, dim=3)
+
+            NetworkLayout.DEFAULT_RNG[] = StableRNG
+            l = Spring()
+            @test l.rng isa StableRNG
+            li = LayoutIterator(l, g)
+            p1, p2 = iterate(li)[1], iterate(li)[1]
+            @test p1 == p2
+            NetworkLayout.DEFAULT_RNG[] = MersenneTwister
         end
 
         @testset "test single node graph" begin