diff --git a/src/Gabs.jl b/src/Gabs.jl index d3ddf80..5dd5f1c 100644 --- a/src/Gabs.jl +++ b/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 diff --git a/src/homodyne.jl b/src/homodyne.jl index dbe4e1f..6bd4f74 100644 --- a/src/homodyne.jl +++ b/src/homodyne.jl @@ -74,7 +74,8 @@ 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 @@ -82,7 +83,7 @@ function homodyne( 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,7 +112,8 @@ 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 @@ -119,7 +121,7 @@ function homodyne( 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) diff --git a/src/randoms.jl b/src/randoms.jl index 2f58551..2a08a3b 100644 --- a/src/randoms.jl +++ b/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 diff --git a/test/Project.toml b/test/Project.toml index 8c69541..3aa6dee 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/test_measurements.jl b/test/test_measurements.jl index 52649c6..3585f95 100644 --- a/test/test_measurements.jl +++ b/test/test_measurements.jl @@ -1,5 +1,6 @@ @testitem "Measurements" begin using Gabs + using Random using StaticArrays using LinearAlgebra: det, I, cholesky, Symmetric @@ -150,6 +151,24 @@ @test (hstatic.state).mean isa SVector && (hstatic.state).covar isa SMatrix @test isequal(hstatic.state.mean[1:2], zeros(2)) + seed = 4242 + base_state = squeezedstate(QuadPairBasis(2), 0.3, π/4) + base_state_block = changebasis(QuadBlockBasis, base_state) + + h1 = homodyne(MersenneTwister(seed), base_state, [1], [0.0]) + h2 = homodyne(MersenneTwister(seed), base_state, [1], [0.0]) + @test h1.result == h2.result + @test h1.state == h2.state + + hb1 = homodyne(MersenneTwister(seed), base_state_block, [1], [0.0]) + hb2 = homodyne(MersenneTwister(seed), base_state_block, [1], [0.0]) + @test hb1.result == hb2.result + @test hb1.state == hb2.state + + samples_kw = rand(Homodyne, base_state, [1], [π/2]; shots = 3, rng = MersenneTwister(seed)) + samples_pos = rand(MersenneTwister(seed), Homodyne, base_state, [1], [π/2]; shots = 3) + @test samples_kw == samples_pos + @test_throws ArgumentError rand(Homodyne, rs_qpair, collect(1:5), [0.0]) @test_throws ArgumentError rand(Homodyne, rs_qblock, collect(1:5), [π/2]) end diff --git a/test/test_randoms.jl b/test/test_randoms.jl index 349a37e..b8d8af0 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -1,5 +1,6 @@ @testitem "Random objects" begin using Gabs + using Random using StaticArrays using LinearAlgebra: eigvals, adjoint @@ -7,13 +8,14 @@ nmodes = rand(1:5) qpairbasis = QuadPairBasis(nmodes) qblockbasis = QuadBlockBasis(nmodes) - U_qpair = Gabs._rand_unitary(qpairbasis) - U_qblock = Gabs._rand_unitary(qblockbasis) + rng = Random.default_rng() + U_qpair = Gabs._rand_unitary(rng, qpairbasis) + U_qblock = Gabs._rand_unitary(rng, qblockbasis) @test isapprox(adjoint(U_qpair), inv(U_qpair), atol = 1e-5) @test isapprox(adjoint(U_qblock), inv(U_qblock), atol = 1e-5) - O_qpair = Gabs._rand_orthogonal_symplectic(qpairbasis) - O_qblock = Gabs._rand_orthogonal_symplectic(qblockbasis) + O_qpair = Gabs._rand_orthogonal_symplectic(rng, qpairbasis) + O_qblock = Gabs._rand_orthogonal_symplectic(rng, qblockbasis) @test isapprox(O_qpair', inv(O_qpair), atol = 1e-5) @test isapprox(O_qblock', inv(O_qblock), atol = 1e-5) @test issymplectic(qpairbasis, O_qpair, atol = 1e-5) @@ -135,4 +137,29 @@ @test rc_static.ħ == 2 @test isgaussian(rc_static, atol = 1e-5) end + + @testset "rng control" begin + seed = 2026 + basis = QuadPairBasis(3) + + rs1 = randstate(MersenneTwister(seed), basis) + rs2 = randstate(MersenneTwister(seed), basis) + @test rs1 == rs2 + + rs_kw = randstate(basis; rng = MersenneTwister(seed)) + rs_pos = randstate(MersenneTwister(seed), basis) + @test rs_kw == rs_pos + + ru1 = randunitary(MersenneTwister(seed), basis) + ru2 = randunitary(MersenneTwister(seed), basis) + @test ru1 == ru2 + + rc1 = randchannel(MersenneTwister(seed), basis) + rc2 = randchannel(MersenneTwister(seed), basis) + @test rc1 == rc2 + + S1 = randsymplectic(MersenneTwister(seed), basis) + S2 = randsymplectic(MersenneTwister(seed), basis) + @test S1 == S2 + end end \ No newline at end of file