Add optional rng arguments to random methods

src/Gabs.jl

@@ -6,7 +6,7 @@ using LinearAlgebra: I, det, mul!, diag, qr, eigvals, Diagonal, cholesky, Symmet
 import QuantumInterface: StateVector, AbstractOperator, apply!, tensor, ⊗, directsum, ⊕, entropy_vn, fidelity, logarithmic_negativity, ptrace, embed
 
 import Random
-using Random: randn!
+using Random: randn!, AbstractRNG
 
 import SymplecticMatrices: williamson, Williamson, polar, Polar, blochmessiah, BlochMessiah, randsymplectic, symplecticform, issymplectic
 using SymplecticMatrices: williamson, Williamson, polar, Polar, blochmessiah, BlochMessiah, BlockForm, PairForm

src/homodyne.jl

@@ -74,15 +74,16 @@ true
 function homodyne(
     state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
     indices::R, 
-    angles::G
+    angles::G;
+    rng::AbstractRNG = Random.default_rng()
 ) where {Tm,Tc,R,G}
     basis = state.basis
     nmodes = basis.nmodes
     indlength = length(indices)
     indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR))
     indlength == length(angles) || throw(ArgumentError(Gabs.GENERALDYNE_ERROR))
     # perform conditional mapping of Gaussian quantum state
-    result′, a, A = _homodyne_filter(state, indices, angles)
+    result′, a, A = _homodyne_filter(rng, state, indices, angles)
     mean′ = zeros(eltype(Tm), 2*nmodes)
     covar′ = Matrix{eltype(Tc)}((state.ħ/2)*I, 2*nmodes, 2*nmodes)
     # fill in measured modes with vacuum states 
@@ -111,15 +112,16 @@ end
 function homodyne(
     state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
     indices::R, 
-    angles::G
+    angles::G;
+    rng::AbstractRNG = Random.default_rng()
 ) where {Tm,Tc,R,G}
     basis = state.basis
     nmodes = basis.nmodes
     indlength = length(indices)
     indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR))
     indlength == length(angles) || throw(ArgumentError(Gabs.GENERALDYNE_ERROR))
     # perform conditional mapping of Gaussian quantum state
-    result′, a, A = _homodyne_filter(state, indices, angles)
+    result′, a, A = _homodyne_filter(rng, state, indices, angles)
     mean′ = zeros(eltype(Tm), 2*nmodes)
     covar′ = Matrix{eltype(Tc)}((state.ħ/2)*I, 2*nmodes, 2*nmodes)
     nmodes′ = nmodes - length(indices)
@@ -147,9 +149,11 @@ function homodyne(
     state′ = GaussianState(basis, mean′′, covar′′, ħ = state.ħ)
     return Homodyne(result′, state′)
 end
+homodyne(rng::AbstractRNG, state::GaussianState{<:QuadPairBasis,Tm,Tc}, indices::R, angles::G) where {Tm,Tc,R,G} = homodyne(state, indices, angles; rng = rng)
+homodyne(rng::AbstractRNG, state::GaussianState{<:QuadBlockBasis,Tm,Tc}, indices::R, angles::G) where {Tm,Tc,R,G} = homodyne(state, indices, angles; rng = rng)
 
 """
-    rand(::Type{Homodyne}, state::GaussianState, indices::Vector, angles::Vector; shots = 1) -> Array
+    rand([rng::AbstractRNG], ::Type{Homodyne}, state::GaussianState, indices::Vector, angles::Vector; shots = 1) -> Array
 
 Compute the projection of the subsystem of a Gaussian state `state` indicated by `indices`
 on rotated quadrature states with homodyne phases given by `angles` and return an array of measured modes.
@@ -172,6 +176,17 @@ julia> rand(Homodyne, st, [1, 3], [pi/2, 0], shots = 5)
 ```
 """
 function Base.rand(
+    ::Type{Homodyne}, 
+    state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
+    indices::R, 
+    angles::G; 
+    shots::Int = 1,
+    rng::AbstractRNG = Random.default_rng()
+) where {Tm,Tc,R,G}
+    return Base.rand(rng, Homodyne, state, indices, angles; shots = shots)
+end
+function Base.rand(
+    rng::AbstractRNG,
     ::Type{Homodyne}, 
     state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
     indices::R, 
@@ -216,12 +231,23 @@ function Base.rand(
     buf = zeros(2*indlength)
     results = zeros(2*indlength, shots)
     @inbounds for i in Base.OneTo(shots)
-        mul!(@view(results[:,i]), L, randn!(buf))
+        mul!(@view(results[:,i]), L, randn!(rng, buf))
         @view(results[:,i]) .+= b
     end
     return results
 end
 function Base.rand(
+    ::Type{Homodyne}, 
+    state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
+    indices::R, 
+    angles::G;
+    shots::Int = 1,
+    rng::AbstractRNG = Random.default_rng()
+) where {Tm,Tc,R,G}
+    return Base.rand(rng, Homodyne, state, indices, angles; shots = shots)
+end
+function Base.rand(
+    rng::AbstractRNG,
     ::Type{Homodyne}, 
     state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
     indices::R, 
@@ -278,13 +304,14 @@ function Base.rand(
     buf = zeros(2*indlength)
     results = zeros(2*indlength, shots)
     @inbounds for i in Base.OneTo(shots)
-        mul!(@view(results[:,i]), L, randn!(buf))
+        mul!(@view(results[:,i]), L, randn!(rng, buf))
         @view(results[:,i]) .+= b
     end
     return results
 end
 
 function _homodyne_filter(
+    rng::AbstractRNG,
     state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
     indices::R, 
     angles::G
@@ -307,7 +334,7 @@ function _homodyne_filter(
     # Cholesky decomposition of the covariance matrix
     symB = Symmetric(B)
     L = cholesky(symB).L
-    resultmean = L * randn(2*indlength) + b
+    resultmean = L * randn(rng, 2*indlength) + b
     meandiff = resultmean - b
     # conditional mapping (see Serafini's Quantum Continuous Variables textbook for reference)
     buf = C * inv(symB)
@@ -318,6 +345,7 @@ function _homodyne_filter(
     return result′, a, A
 end
 function _homodyne_filter(
+    rng::AbstractRNG,
     state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
     indices::R, 
     angles::G
@@ -340,7 +368,7 @@ function _homodyne_filter(
     # Cholesky decomposition of the covariance matrix
     symB = Symmetric(B)
     L = cholesky(symB).L
-    resultmean = L * randn(2*indlength) + b
+    resultmean = L * randn(rng, 2*indlength) + b
     meandiff = resultmean - b
     # conditional mapping (see Serafini's Quantum Continuous Variables textbook for reference)
     buf = C * inv(symB)

src/randoms.jl

@@ -1,22 +1,25 @@
 """
-    randstate([Tm=Vector{Float64}, Tc=Matrix{Float64},] basis::SymplecticBasis; pure=false)
+    randstate([Tm=Vector{Float64}, Tc=Matrix{Float64},] basis::SymplecticBasis; pure=false, rng=Random.default_rng())
 
 Calculate a random Gaussian state in symplectic representation defined by `basis`.
 """
-function randstate(::Type{Tm}, ::Type{Tc}, basis::SymplecticBasis{N}; pure = false, ħ = 2) where {Tm,Tc,N<:Int}
-    mean, covar = _randstate(basis; pure = pure, ħ = ħ)
+function randstate(::Type{Tm}, ::Type{Tc}, basis::SymplecticBasis{N}; pure = false, ħ = 2, rng::AbstractRNG = Random.default_rng()) where {Tm,Tc,N<:Int}
+    mean, covar = _randstate(rng, basis; pure = pure, ħ = ħ)
     return GaussianState(basis, Tm(mean), Tc(covar), ħ = ħ)
 end
-randstate(::Type{T}, basis::SymplecticBasis{N}; pure = false, ħ = 2) where {T,N<:Int} = randstate(T,T,basis; pure = pure, ħ = ħ)
-function randstate(basis::SymplecticBasis{N}; pure = false, ħ = 2) where {N<:Int}
-    mean, covar = _randstate(basis; pure = pure, ħ = ħ)
+randstate(::Type{T}, basis::SymplecticBasis{N}; pure = false, ħ = 2, rng::AbstractRNG = Random.default_rng()) where {T,N<:Int} = randstate(T,T,basis; pure = pure, ħ = ħ, rng = rng)
+function randstate(basis::SymplecticBasis{N}; pure = false, ħ = 2, rng::AbstractRNG = Random.default_rng()) where {N<:Int}
+    mean, covar = _randstate(rng, basis; pure = pure, ħ = ħ)
     return GaussianState(basis, mean, covar, ħ = ħ)
 end
-function _randstate(basis::QuadPairBasis{N}; pure = false, ħ = 2) where {N<:Int}
+randstate(rng::AbstractRNG, ::Type{Tm}, ::Type{Tc}, basis::SymplecticBasis{N}; pure = false, ħ = 2) where {Tm,Tc,N<:Int} = randstate(Tm, Tc, basis; pure = pure, ħ = ħ, rng = rng)
+randstate(rng::AbstractRNG, ::Type{T}, basis::SymplecticBasis{N}; pure = false, ħ = 2) where {T,N<:Int} = randstate(T, basis; pure = pure, ħ = ħ, rng = rng)
+randstate(rng::AbstractRNG, basis::SymplecticBasis{N}; pure = false, ħ = 2) where {N<:Int} = randstate(basis; pure = pure, ħ = ħ, rng = rng)
+function _randstate(rng::AbstractRNG, basis::QuadPairBasis{N}; pure = false, ħ = 2) where {N<:Int}
     nmodes = basis.nmodes
-    mean = randn(2*nmodes)
+    mean = randn(rng, 2*nmodes)
     covar = zeros(2*nmodes, 2*nmodes)
-    symp = randsymplectic(basis)
+    symp = randsymplectic(rng, basis)
     # generate pure Gaussian state
     if pure
         mul!(covar, symp, symp')
@@ -26,16 +29,16 @@ function _randstate(basis::QuadPairBasis{N}; pure = false, ħ = 2) where {N<:Int
     # create buffer for matrix multiplication
     buf = zeros(2*nmodes, 2*nmodes)
     # William decomposition for mixed Gaussian states
-    sympeigs = (ħ/2)*(abs.(rand(nmodes)) .+ 1.0)
+    sympeigs = (ħ/2)*(abs.(rand(rng, nmodes)) .+ 1.0)
     will = Diagonal(repeat(sympeigs, inner = 2))
     mul!(covar, symp, mul!(buf, will, symp'))
     return mean, covar
 end
-function _randstate(basis::QuadBlockBasis{N}; pure = false, ħ = 2) where {N<:Int}
+function _randstate(rng::AbstractRNG, basis::QuadBlockBasis{N}; pure = false, ħ = 2) where {N<:Int}
     nmodes = basis.nmodes
-    mean = randn(2*nmodes)
+    mean = randn(rng, 2*nmodes)
     covar = zeros(2*nmodes, 2*nmodes)
-    symp = randsymplectic(basis)
+    symp = randsymplectic(rng, basis)
     # generate pure Gaussian state
     if pure
         mul!(covar, symp, symp')
@@ -45,53 +48,59 @@ function _randstate(basis::QuadBlockBasis{N}; pure = false, ħ = 2) where {N<:In
     # create buffer for matrix multiplication
     buf = zeros(2*nmodes, 2*nmodes)
     # William decomposition for mixed Gaussian states
-    sympeigs = (ħ/2)*(abs.(rand(nmodes)) .+ 1.0)
+    sympeigs = (ħ/2)*(abs.(rand(rng, nmodes)) .+ 1.0)
     will = Diagonal(repeat(sympeigs, outer = 2))
     mul!(covar, symp, mul!(buf, will, symp'))
     return mean, covar
 end
 
 """
-    randunitary([Td=Vector{Float64}, Ts=Matrix{Float64},] basis::SymplecticBasis; passive=false)
+    randunitary([Td=Vector{Float64}, Ts=Matrix{Float64},] basis::SymplecticBasis; passive=false, rng=Random.default_rng())
 
 Calculate a random Gaussian unitary operator in symplectic representation defined by `basis`.
 """
-function randunitary(::Type{Td}, ::Type{Ts}, basis::SymplecticBasis{N}; passive = false, ħ = 2) where {Td,Ts,N<:Int}
-    disp, symp = _randunitary(basis, passive = passive)
+function randunitary(::Type{Td}, ::Type{Ts}, basis::SymplecticBasis{N}; passive = false, ħ = 2, rng::AbstractRNG = Random.default_rng()) where {Td,Ts,N<:Int}
+    disp, symp = _randunitary(rng, basis, passive = passive)
     return GaussianUnitary(basis, Td(disp), Ts(symp), ħ = ħ)
 end
-randunitary(::Type{T}, basis::SymplecticBasis{N}; passive = false, ħ = 2) where {T,N<:Int} = randunitary(T,T,basis; passive = passive, ħ = ħ)
-function randunitary(basis::SymplecticBasis{N}; passive = false, ħ = 2) where {N<:Int}
-    disp, symp = _randunitary(basis, passive = passive)
+randunitary(::Type{T}, basis::SymplecticBasis{N}; passive = false, ħ = 2, rng::AbstractRNG = Random.default_rng()) where {T,N<:Int} = randunitary(T,T,basis; passive = passive, ħ = ħ, rng = rng)
+function randunitary(basis::SymplecticBasis{N}; passive = false, ħ = 2, rng::AbstractRNG = Random.default_rng()) where {N<:Int}
+    disp, symp = _randunitary(rng, basis, passive = passive)
     return GaussianUnitary(basis, disp, symp, ħ = ħ)
 end
-function _randunitary(basis::SymplecticBasis{N}; passive = false) where {N<:Int}
+randunitary(rng::AbstractRNG, ::Type{Td}, ::Type{Ts}, basis::SymplecticBasis{N}; passive = false, ħ = 2) where {Td,Ts,N<:Int} = randunitary(Td, Ts, basis; passive = passive, ħ = ħ, rng = rng)
+randunitary(rng::AbstractRNG, ::Type{T}, basis::SymplecticBasis{N}; passive = false, ħ = 2) where {T,N<:Int} = randunitary(T, basis; passive = passive, ħ = ħ, rng = rng)
+randunitary(rng::AbstractRNG, basis::SymplecticBasis{N}; passive = false, ħ = 2) where {N<:Int} = randunitary(basis; passive = passive, ħ = ħ, rng = rng)
+function _randunitary(rng::AbstractRNG, basis::SymplecticBasis{N}; passive = false) where {N<:Int}
     nmodes = basis.nmodes
-    disp = rand(2*nmodes)
-    symp = randsymplectic(basis, passive = passive)
+    disp = rand(rng, 2*nmodes)
+    symp = randsymplectic(rng, basis, passive = passive)
     return disp, symp
 end
 
 """
-    randchannel([Td=Vector{Float64}, Tt=Matrix{Float64},] basis::SymplecticBasis)
+    randchannel([Td=Vector{Float64}, Tt=Matrix{Float64},] basis::SymplecticBasis; rng=Random.default_rng())
 
 Calculate a random Gaussian channel in symplectic representation defined by `basis`.
 """
-function randchannel(::Type{Td}, ::Type{Tt}, basis::SymplecticBasis{N}; ħ = 2) where {Td,Tt,N<:Int}
-    disp, transform, noise = _randchannel(basis)
+function randchannel(::Type{Td}, ::Type{Tt}, basis::SymplecticBasis{N}; ħ = 2, rng::AbstractRNG = Random.default_rng()) where {Td,Tt,N<:Int}
+    disp, transform, noise = _randchannel(rng, basis)
     return GaussianChannel(basis, Td(disp), Tt(transform), Tt(noise), ħ = ħ)
 end
-randchannel(::Type{T}, basis::SymplecticBasis{N}; ħ = 2) where {T,N<:Int} = randchannel(T,T,basis; ħ = ħ)
-function randchannel(basis::SymplecticBasis{N}; ħ = 2) where {N<:Int}
-    disp, transform, noise = _randchannel(basis)
+randchannel(::Type{T}, basis::SymplecticBasis{N}; ħ = 2, rng::AbstractRNG = Random.default_rng()) where {T,N<:Int} = randchannel(T,T,basis; ħ = ħ, rng = rng)
+function randchannel(basis::SymplecticBasis{N}; ħ = 2, rng::AbstractRNG = Random.default_rng()) where {N<:Int}
+    disp, transform, noise = _randchannel(rng, basis)
     return GaussianChannel(basis, disp, transform, noise, ħ = ħ)
 end
-function _randchannel(basis::SymplecticBasis{N}) where {N<:Int}
+randchannel(rng::AbstractRNG, ::Type{Td}, ::Type{Tt}, basis::SymplecticBasis{N}; ħ = 2) where {Td,Tt,N<:Int} = randchannel(Td, Tt, basis; ħ = ħ, rng = rng)
+randchannel(rng::AbstractRNG, ::Type{T}, basis::SymplecticBasis{N}; ħ = 2) where {T,N<:Int} = randchannel(T, basis; ħ = ħ, rng = rng)
+randchannel(rng::AbstractRNG, basis::SymplecticBasis{N}; ħ = 2) where {N<:Int} = randchannel(basis; ħ = ħ, rng = rng)
+function _randchannel(rng::AbstractRNG, basis::SymplecticBasis{N}) where {N<:Int}
     nmodes = basis.nmodes
-    disp = randn(2*nmodes)
+    disp = randn(rng, 2*nmodes)
     # generate symplectic matrix describing the evolution of the system with N modes
     # and environment with 2N modes
-    symp = randsymplectic(typeof(basis)(3*nmodes))
+    symp = randsymplectic(rng, typeof(basis)(3*nmodes))
     transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:6*nmodes])
     # generate noise matrix from evolution of environment
     noise = zeros(2*nmodes, 2*nmodes)
@@ -100,26 +109,31 @@ function _randchannel(basis::SymplecticBasis{N}) where {N<:Int}
 end
 
 """
-    randsymplectic([T=Matrix{Float64},] basis::SymplecticBasis, passive=false)
+    randsymplectic([T=Matrix{Float64},] basis::SymplecticBasis, passive=false, rng=Random.default_rng())
 
 Calculate a random symplectic matrix in symplectic representation defined by `basis`.
 """
-function randsymplectic(::Type{T}, basis::SymplecticBasis{N}; passive = false) where {T, N<:Int} 
-    symp = randsymplectic(basis, passive = passive)
+function randsymplectic(::Type{T}, basis::SymplecticBasis{N}; passive = false, rng::AbstractRNG = Random.default_rng()) where {T, N<:Int}
+    symp = randsymplectic(rng, basis; passive = passive)
     return T(symp)
 end
-randsymplectic(basis::SymplecticBasis{N}; passive = false) where {N<:Int} = _randsymplectic(basis, passive = passive)
-function _randsymplectic(basis::QuadPairBasis{N}; passive = false) where {N<:Int}
+function randsymplectic(rng::AbstractRNG, ::Type{T}, basis::SymplecticBasis{N}; passive = false) where {T, N<:Int}
+    symp = randsymplectic(rng, basis; passive = passive)
+    return T(symp)
+end
+randsymplectic(basis::SymplecticBasis{N}; passive = false, rng::AbstractRNG = Random.default_rng()) where {N<:Int} = _randsymplectic(rng, basis, passive = passive)
+randsymplectic(rng::AbstractRNG, basis::SymplecticBasis{N}; passive = false) where {N<:Int} = _randsymplectic(rng, basis, passive = passive)
+function _randsymplectic(rng::AbstractRNG, basis::QuadPairBasis{N}; passive = false) where {N<:Int}
     nmodes = basis.nmodes
     # generate random orthogonal symplectic matrix
-    O = _rand_orthogonal_symplectic(basis)
+    O = _rand_orthogonal_symplectic(rng, basis)
     if passive
         return O
     end
     # instead generate random active symplectic matrix
-    O′ = _rand_orthogonal_symplectic(basis)
+    O′ = _rand_orthogonal_symplectic(rng, basis)
     # direct sum of symplectic matrices for single-mode squeeze transformations
-    rs = rand(nmodes)
+    rs = rand(rng, nmodes)
     squeezes = zeros(2*nmodes, 2*nmodes)
     @inbounds for i in Base.OneTo(nmodes)
         val = rs[i]
@@ -128,17 +142,17 @@ function _randsymplectic(basis::QuadPairBasis{N}; passive = false) where {N<:Int
     end
     return O * squeezes * O′
 end
-function _randsymplectic(basis::QuadBlockBasis{N}; passive = false) where {N<:Int}
+function _randsymplectic(rng::AbstractRNG, basis::QuadBlockBasis{N}; passive = false) where {N<:Int}
     nmodes = basis.nmodes
     # generate random orthogonal symplectic matrix
-    O = _rand_orthogonal_symplectic(basis)
+    O = _rand_orthogonal_symplectic(rng, basis)
     if passive
         return O
     end
     # instead generate random active symplectic matrix
-    O′ = _rand_orthogonal_symplectic(basis)
+    O′ = _rand_orthogonal_symplectic(rng, basis)
     # direct sum of symplectic matrices for single-mode squeeze transformations
-    rs = rand(nmodes)
+    rs = rand(rng, nmodes)
     squeezes = zeros(2*nmodes, 2*nmodes)
     @inbounds for i in Base.OneTo(nmodes)
         val = rs[i]
@@ -151,9 +165,9 @@ end
 
 # Generates random orthogonal symplectic matrix by blocking real
 # and imaginary parts of a random unitary matrix
-function _rand_orthogonal_symplectic(basis::QuadPairBasis{N}) where {N<:Int}
+function _rand_orthogonal_symplectic(rng::AbstractRNG, basis::QuadPairBasis{N}) where {N<:Int}
     nmodes = basis.nmodes
-    U = _rand_unitary(basis)
+    U = _rand_unitary(rng, basis)
     O = zeros(2*nmodes, 2*nmodes)
     @inbounds for i in Base.OneTo(nmodes), j in Base.OneTo(nmodes)
         val = U[i,j]
@@ -164,9 +178,9 @@ function _rand_orthogonal_symplectic(basis::QuadPairBasis{N}) where {N<:Int}
     end
     return O
 end
-function _rand_orthogonal_symplectic(basis::QuadBlockBasis{N}) where {N<:Int}
+function _rand_orthogonal_symplectic(rng::AbstractRNG, basis::QuadBlockBasis{N}) where {N<:Int}
     nmodes = basis.nmodes
-    U = _rand_unitary(basis)
+    U = _rand_unitary(rng, basis)
     O = zeros(2*nmodes, 2*nmodes)
     @inbounds for i in Base.OneTo(nmodes), j in Base.OneTo(nmodes)
         val = U[i,j]
@@ -177,12 +191,14 @@ function _rand_orthogonal_symplectic(basis::QuadBlockBasis{N}) where {N<:Int}
     end
     return O
 end
-# Generates unitary matrix randomly distributed over Haar measure;
-# see https://arxiv.org/abs/math-ph/0609050 for algorithm.
-# This approach is faster and creates less allocations than rand(Haar(2), nmodes) from RandomMatrices.jl
-function _rand_unitary(basis::SymplecticBasis{N}) where {N<:Int}
+"""
+Generates unitary matrix randomly distributed over Haar measure;
+see https://arxiv.org/abs/math-ph/0609050 for algorithm.
+This approach is faster and creates less allocations than rand(Haar(2), nmodes) from RandomMatrices.jl
+"""
+function _rand_unitary(rng::AbstractRNG, basis::SymplecticBasis{N}) where {N<:Int}
     nmodes = basis.nmodes
-    M = rand(ComplexF64, nmodes, nmodes) ./ sqrt(2.0)
+    M = rand(rng, ComplexF64, nmodes, nmodes) ./ sqrt(2.0)
     q, r = qr(M)
     d = Diagonal([r[i, i] / abs(r[i, i]) for i in Base.OneTo(nmodes)])
     return q * d

test/Project.toml

@@ -5,6 +5,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"