Allow general RNG in random_* functions (#106)

* allow general RNG in random_* functions

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* fix missing method

* add another missing method

* add defaults

* remove default

* remove default

* add another default...

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Project.toml

@@ -8,6 +8,7 @@ DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -32,6 +33,7 @@ ForwardDiff = "0.10.36"
 LinearAlgebra = "1"
 Meshes = "0.52.1"
 Printf = "1"
+Random = "1"
 ReadVTK = "0.2"
 RecipesBase = "1.3.4"
 Reexport = "1.2"

src/KernelInterpolation.jl

@@ -15,6 +15,7 @@ using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 using ForwardDiff: ForwardDiff
 using LinearAlgebra: Symmetric, I, norm, tr, muladd, dot, diagind
 using Printf: @sprintf
+using Random: Random
 using ReadVTK: VTKFile, get_points, get_point_data, get_data
 using RecipesBase: RecipesBase, @recipe, @series
 using SciMLBase: ODEFunction, ODEProblem, ODESolution, DiscreteCallback, u_modified!

src/nodes.jl

@@ -220,65 +220,80 @@ end
 
 # Some convenience function to create some specific `NodeSet`s
 """
-    random_hypercube(n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])
+    random_hypercube([rng], n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])
 
 Create a [`NodeSet`](@ref) with `n` random nodes each of dimension `dim` inside a hypercube defined by
 the bounds `x_min` and `x_max`. If the bounds are given as single values, they are applied for
 each dimension. If they are `Tuple`s of size `dim` the hypercube has the according bounds.
 If `dim` is not given explicitly, it is inferred by the lengths of `x_min` and `x_max` if possible.
+Optionally, pass a random number generator `rng`.
 """
-function random_hypercube(n::Int, x_min::Real = 0.0, x_max::Real = 1.0; dim = 1)
-    nodes = x_min .+ (x_max - x_min) .* rand(n, dim)
+function random_hypercube(n, x_min = 0.0, x_max = 1.0; kwargs...)
+    random_hypercube(Random.default_rng(), n, x_min, x_max; kwargs...)
+end
+
+function random_hypercube(rng::Random.AbstractRNG, n::Int, x_min::Real = 0.0,
+                          x_max::Real = 1.0; dim = 1)
+    nodes = x_min .+ (x_max - x_min) .* rand(rng, n, dim)
     return NodeSet(nodes)
 end
 
-function random_hypercube(n::Int, x_min::NTuple{Dim}, x_max::NTuple{Dim};
+function random_hypercube(rng::Random.AbstractRNG, n::Int, x_min::NTuple{Dim},
+                          x_max::NTuple{Dim};
                           dim = Dim) where {Dim}
     @assert dim == Dim
-    nodes = rand(n, dim)
+    nodes = rand(rng, n, dim)
     for i in 1:dim
         nodes[:, i] = x_min[i] .+ (x_max[i] - x_min[i]) .* view(nodes, :, i)
     end
     return NodeSet(nodes)
 end
 
 """
-    random_hypercube_boundary(n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])
+    random_hypercube_boundary([rng], n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])
 
 Create a [`NodeSet`](@ref) with `n` random nodes each of dimension `dim` on the boundary of a hypercube
 defined by the bounds `x_min` and `x_max`. If the bounds are given as single values, they are
 applied for each dimension. If they are `Tuple`s of size `dim` the hypercube has the according bounds.
 If `dim` is not given explicitly, it is inferred by the lengths of `x_min` and `x_max` if possible.
+Optionally, pass a random number generator `rng`.
 """
-function random_hypercube_boundary(n::Int, x_min::Real = 0.0, x_max::Real = 1.0; dim = 1)
-    random_hypercube_boundary(n, ntuple(_ -> x_min, dim), ntuple(_ -> x_max, dim))
+function random_hypercube_boundary(n, x_min = 0.0, x_max = 1.0; kwargs...)
+    random_hypercube_boundary(Random.default_rng(), n, x_min, x_max; kwargs...)
+end
+
+function random_hypercube_boundary(rng::Random.AbstractRNG, n::Int, x_min::Real = 0.0,
+                                   x_max::Real = 1.0; dim = 1)
+    random_hypercube_boundary(rng, n, ntuple(_ -> x_min, dim), ntuple(_ -> x_max, dim))
 end
 
-function project_on_hypercube_boundary!(nodeset::NodeSet{Dim}, x_min::NTuple{Dim},
+function project_on_hypercube_boundary!(rng::Random.AbstractRNG, nodeset::NodeSet{Dim},
+                                        x_min::NTuple{Dim},
                                         x_max::NTuple{Dim}) where {Dim}
     for i in eachindex(nodeset)
         # j = argmin([abs.(nodeset[i] .- x_min); abs.(nodeset[i] .- x_max)])
         # Project to random axis
-        j = rand(1:Dim)
-        if rand([1, 2]) == 1
+        j = rand(rng, 1:Dim)
+        if rand(rng, [1, 2]) == 1
             nodeset[i][j] = x_min[j]
         else
             nodeset[i][j] = x_max[j]
         end
     end
 end
 
-function random_hypercube_boundary(n::Int, x_min::NTuple{Dim}, x_max::NTuple{Dim};
+function random_hypercube_boundary(rng::Random.AbstractRNG, n::Int, x_min::NTuple{Dim},
+                                   x_max::NTuple{Dim};
                                    dim = Dim) where {Dim}
     @assert dim == Dim
     if dim == 1 && n >= 2
         @warn "For one dimension the boundary of the hypercube consists only of 2 points"
         return NodeSet([x_min[1], x_max[1]])
     end
     # First, create random nodes *inside* hypercube
-    nodeset = random_hypercube(n, x_min, x_max)
+    nodeset = random_hypercube(rng, n, x_min, x_max)
     # Then, project all the nodes on the boundary
-    project_on_hypercube_boundary!(nodeset, x_min, x_max)
+    project_on_hypercube_boundary!(rng, nodeset, x_min, x_max)
     return nodeset
 end
 
@@ -406,44 +421,64 @@ function homogeneous_hypercube_boundary(n::NTuple{Dim},
 end
 
 """
-    random_hypersphere(n, r = 1.0, center = zeros(dim); [dim])
+    random_hypersphere([rng], n, r = 1.0, center = zeros(dim); [dim])
 
 Create a [`NodeSet`](@ref) with `n` random nodes each of dimension `dim` inside a hypersphere with
 radius `r` around the center `center`.
 If `dim` is not given explicitly, it is inferred by the length of `center` if possible.
+Optionally, pass a random number generator `rng`.
 """
-function random_hypersphere(n::Int, r = 1.0; dim = 2)
-    random_hypersphere(n, r, zeros(dim))
+function random_hypersphere(n, r = 1.0; kwargs...)
+    random_hypersphere(Random.default_rng(), n, r; kwargs...)
+end
+
+function random_hypersphere(n, r, center; kwargs...)
+    random_hypersphere(Random.default_rng(), n, r, center; kwargs...)
 end
 
-function random_hypersphere(n::Int, r::Real, center::AbstractVector; dim = length(center))
+function random_hypersphere(rng::Random.AbstractRNG, n, r = 1.0; dim = 2)
+    random_hypersphere(rng, n, r, zeros(dim))
+end
+
+function random_hypersphere(rng::Random.AbstractRNG, n, r,
+                            center; dim = length(center))
     @assert length(center) == dim
-    nodes = randn(n, dim)
+    nodes = randn(rng, n, dim)
     for i in 1:n
-        nodes[i, :] .= center .+ r .* nodes[i, :] ./ norm(nodes[i, :]) * rand()^(1 / dim)
+        nodes[i, :] .= center .+ r .* nodes[i, :] ./ norm(nodes[i, :]) * rand(rng)^(1 / dim)
     end
     return NodeSet(nodes)
 end
 
 """
-    random_hypersphere_boundary(n, r = 1.0, center = zeros(dim); [dim])
+    random_hypersphere_boundary([rng], n, r = 1.0, center = zeros(dim); [dim])
 
 Create a [`NodeSet`](@ref) with `n` random nodes each of dimension `dim` at the boundary of a
 hypersphere with radius `r` around the center `center`.
 If `dim` is not given explicitly, it is inferred by the length of `center` if possible.
+Optionally, pass a random number generator `rng`.
 """
-function random_hypersphere_boundary(n::Int, r = 1.0; dim = 2)
-    random_hypersphere_boundary(n, r, zeros(dim))
+function random_hypersphere_boundary(n, r = 1.0; kwargs...)
+    random_hypersphere_boundary(Random.default_rng(), n, r; kwargs...)
+end
+
+function random_hypersphere_boundary(n, r, center; kwargs...)
+    random_hypersphere_boundary(Random.default_rng(), n, r, center; kwargs...)
+end
+
+function random_hypersphere_boundary(rng::Random.AbstractRNG, n, r = 1.0; dim = 2)
+    random_hypersphere_boundary(rng, n, r, zeros(dim))
 end
 
-function random_hypersphere_boundary(n::Int, r::Real, center::AbstractVector;
+function random_hypersphere_boundary(rng::Random.AbstractRNG, n, r,
+                                     center;
                                      dim = length(center))
     @assert length(center) == dim
     if dim == 1 && n >= 2
         @warn "For one dimension the boundary of the hypersphere consists only of 2 points"
         return NodeSet([-r, r])
     end
-    nodes = randn(n, dim)
+    nodes = randn(rng, n, dim)
     for i in 1:n
         nodes[i, :] .= center .+ r .* nodes[i, :] ./ norm(nodes[i, :])
     end