From b63b282f3ee0c7bd08d924b4209c579d5563bc5f Mon Sep 17 00:00:00 2001 From: apkille Date: Tue, 12 Nov 2024 23:45:13 -0500 Subject: [PATCH 1/8] add random features --- src/Gabs.jl | 14 +++-- src/metrics.jl | 6 ++ src/randoms.jl | 135 ++++++++++++++++++++++++++++++++++++------- test/test_randoms.jl | 18 ++++++ 4 files changed, 147 insertions(+), 26 deletions(-) create mode 100644 src/metrics.jl diff --git a/src/Gabs.jl b/src/Gabs.jl index 8d3a5c4..1a5a45c 100644 --- a/src/Gabs.jl +++ b/src/Gabs.jl @@ -3,7 +3,7 @@ module Gabs using BlockArrays: BlockedArray, BlockArray, Block, mortar import LinearAlgebra -using LinearAlgebra: I, det, mul! +using LinearAlgebra: I, det, mul!, diagm, diag, qr import QuantumInterface: StateVector, AbstractOperator, apply!, tensor, ⊗ @@ -13,14 +13,18 @@ export # operations tensor, ⊗, apply!, ptrace, output, prob, # predefined Gaussian states - vacuumstate, thermalstate, coherentstate, squeezedstate, eprstate, randstate, + vacuumstate, thermalstate, coherentstate, squeezedstate, eprstate, # predefined Gaussian channels displace, squeeze, twosqueeze, phaseshift, beamsplitter, - attenuator, amplifier, randchannel, + attenuator, amplifier, + # random objects + randstate, randunitary, randchannel, randsymplectic, # wigner functions wigner, wignerchar, # symplectic form - symplecticform + symplecticform, + # metrics + purity include("errors.jl") @@ -40,4 +44,6 @@ include("measurements.jl") include("wigner.jl") +include("metrics.jl") + end diff --git a/src/metrics.jl b/src/metrics.jl new file mode 100644 index 0000000..e145493 --- /dev/null +++ b/src/metrics.jl @@ -0,0 +1,6 @@ +""" + purity(state::GaussianState) + +Calculate the purity of a Gaussian state, defined by `1/sqrt(det(V))`. +""" +purity(state::GaussianState) = 1/sqrt(det(state.covar)) \ No newline at end of file diff --git a/src/randoms.jl b/src/randoms.jl index d8f8a6a..c6212e8 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -1,23 +1,54 @@ """ - randstate([Tm=Vector{Float64}, Tc=Matrix{Float64},] nmodes<:Int) + randstate([Tm=Vector{Float64}, Tc=Matrix{Float64},] nmodes<:Int; pure=false) Calculate a random Gaussian state. """ -function randstate(::Type{Tm}, ::Type{Tc}, nmodes::N) where {Tm,Tc,N<:Int} - mean = randn(2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - covar = randmat' * randmat +function randstate(::Type{Tm}, ::Type{Tc}, nmodes::N; pure = false) where {Tm,Tc,N<:Int} + mean, covar = _randstate_fields(nmodes, pure) return GaussianState(Tm(mean), Tc(covar), nmodes) end -randstate(::Type{T}, nmodes::N) where {T,N<:Int} = randstate(T,T,nmodes) -function randstate(nmodes::N) where {N<:Int} - mean = randn(2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - covar = randmat' * randmat +randstate(::Type{T}, nmodes::N; pure = false) where {T,N<:Int} = randstate(T,T,nmodes,pure=pure) +function randstate(nmodes::N; pure = false) where {N<:Int} + mean, covar = _randstate_fields(nmodes, pure) return GaussianState(mean, covar, nmodes) end +function _randstate_fields(nmodes::N, pure) where {N<:Int} + mean = randn(2*nmodes) + covar = zeros(2*nmodes, 2*nmodes) + symp = randsymplectic(nmodes) + # generate pure Gaussian state + if pure + mul!(covar, symp, symp') + return mean, covar + end + # create buffer for matrix multiplication + buf = zeros(2*nmodes, 2*nmodes) + # William decomposition for mixed Gaussian states + sympeigs = abs.(rand(nmodes)) .+ 1.0 + will = diagm(repeat(sympeigs, inner = 2)) + mul!(covar, symp, mul!(buf, will, symp')) + return mean, covar +end + +""" + randunitary([Td=Vector{Float64}, Ts=Matrix{Float64},] nmodes<:Int; passive=false) + +Calculate a random Gaussian unitary operator. +""" +function randunitary(::Type{Td}, ::Type{Ts}, nmodes::N; passive = false) where {Td,Ts,N} + disp, symp = _randunitary_fields(nmodes, passive) + return GaussianUnitary(Td(disp), Ts(symp)) +end +randunitary(::Type{T}, nmodes::N; passive = false) where {T,N<:Int} = randunitary(T,T,nmodes; passive = passive) +function randunitary(nmodes::N; passive = false) where {N<:Int} + disp, symp = _randunitary_fields(nmodes, passive) + return GaussianUnitary(disp, symp) +end +function _randunitary_fields(nmodes::N, passive) where {N<:Int} + disp = rand(2*nmodes) + symp = randsymplectic(nmodes, passive = passive) + return disp, symp +end """ randchannel([Td=Vector{Float64}, Tt=Matrix{Float64},] nmodes<:Int) @@ -25,19 +56,79 @@ end Calculate a random Gaussian channel. """ function randchannel(::Type{Td}, ::Type{Tt}, nmodes::N) where {Td,Tt,N<:Int} - disp = randn(2*nmodes) - transform = randn(2*nmodes, 2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - noise = randmat' * randmat + disp, transform, noise = _randchannel(nmodes) return GaussianChannel(Td(disp), Tt(transform), Tt(noise), nmodes) end randchannel(::Type{T}, nmodes::N) where {T,N<:Int} = randchannel(T,T,nmodes) function randchannel(nmodes::N) where {N<:Int} - disp = randn(2*nmodes) - transform = randn(2*nmodes, 2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - noise = randmat' * randmat + disp, transform, noise = _randchannel(nmodes) return GaussianChannel(disp, transform, noise, nmodes) +end +function _randchannel(nmodes::N) where {N<:Int} + disp = randn(2*nmodes) + # generate symplectic matrix describing the evolution of the system + # and environment + m = rand(1:nmodes) # number of modes in environment + symp = randsymplectic(nmodes + m) + transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:2*(nmodes+m)]) + # generate random covariance matrix of environment + envcovar = zeros(2*m, 2*m) + envsymp = randsymplectic(m) + # create buffers for matrix multiplication + buf1, buf2 = zeros(2*m, 2*m), zeros(2*m, 2*nmodes) + # William decomposition for mixed Gaussian states + sympeigs = abs.(rand(m)) .+ 1.0 + will = diagm(repeat(sympeigs, inner = 2)) + mul!(envcovar, envsymp, mul!(buf1, will, envsymp')) + # generate noise matrix from evolution of environment + noise = zeros(2*nmodes, 2*nmodes) + mul!(noise, B, mul!(buf2, envcovar, B')) + return disp, transform, noise +end + +""" + randsymplectic([T=Matrix{Float64},] nmodes<:Int, passive=false) + +Calculate a random symplectic matrix of size `2*nmodes x 2*nmodes`. +""" +function randsymplectic(nmodes::N; passive = false) where {N<:Int} + # generate random orthogonal symplectic matrix + O = _rand_orthogonal_symplectic(nmodes) + if passive + return O + end + # instead generate random active symplectic matrix + O′ = _rand_orthogonal_symplectic(nmodes) + # direct sum of symplectic matrices for single-mode squeeze transformations + rs = rand(nmodes) + squeezes = diagm(vcat(exp.(-rs), exp.(rs))) + return O * squeezes * O′ +end +function randsymplectic(::Type{T}, nmodes::N; passive = false) where {T, N<:Int} + symp = randsymplectic(nmodes, passive) + return T(symp) +end + +# Generates unitary matrix randomly distributed over Haar measure; +# see https://arxiv.org/abs/math-ph/0609050 for algorithm. +# This approach is much faster than using rand(Haar(2), nmodes) from RandomMatrices.jl +function _rand_unitary(nmodes::N) where {N<:Int} + M = rand(ComplexF64, nmodes, nmodes) ./ sqrt(2.0) + q, r = qr(M) + d = diagm([r[i, i] / abs(r[i, i]) for i in Base.OneTo(nmodes)]) + return q * d +end +# Generates random orthogonal symplectic matrix by blocking real +# and imaginary parts of a random unitary matrix +function _rand_orthogonal_symplectic(nmodes::N) where {N<:Int} + U = _rand_unitary(nmodes) + O = zeros(2*nmodes, 2*nmodes) + @inbounds for i in Base.OneTo(nmodes), j in Base.OneTo(nmodes) + val = U[i,j] + O[i,j] = real(val) + O[i+nmodes, j] = imag(val) + O[i, j+nmodes] = -imag(val) + O[i+nmodes, j+nmodes] = real(val) + end + return O end \ No newline at end of file diff --git a/test/test_randoms.jl b/test/test_randoms.jl index c8a65bc..0a664b7 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -1,6 +1,18 @@ @testitem "Measurements" begin using Gabs using StaticArrays + using LinearAlgebra: eigvals, adjoint + + @testset "random utils" begin + nmodes = rand(1:20) + U = Gabs._rand_unitary(nmodes) + @test isapprox(adjoint(U), inv(U)) + + O = Gabs._rand_orthogonal_symplectic(nmodes) + @test isapprox(O', inv(O)) + Omega = symplecticform(nmodes) + @test isapprox(O' * Omega * O, Omega) + end @testset "random types" begin nmodes = rand(1:20) @@ -8,6 +20,12 @@ rc = randchannel(nmodes) @test rc isa GaussianChannel @test rc * rs isa GaussianState + Omega = symplecticform(nmodes) + @test all(i -> (i >= 0) || isapprox(i, 0.0, atol = 1e-5), real(eigvals(rs.covar .+ im*Omega))) + + rspure = randstate(nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rspure.covar .+ im*Omega))) + @test isapprox(purity(rspure), 1.0, atol = 1e-5) rsfloat32 = randstate(AbstractArray{Float32}, nmodes) rcfloat32 = randchannel(AbstractArray{Float32}, nmodes) From ab9f768100566740e2e5f0bacf59f3154ab62020 Mon Sep 17 00:00:00 2001 From: apkille Date: Wed, 13 Nov 2024 11:32:23 -0500 Subject: [PATCH 2/8] add more tests --- src/randoms.jl | 15 ++++---- test/test_randoms.jl | 83 +++++++++++++++++++++++++++++++++++++------- 2 files changed, 77 insertions(+), 21 deletions(-) diff --git a/src/randoms.jl b/src/randoms.jl index c6212e8..a5987f5 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -66,17 +66,16 @@ function randchannel(nmodes::N) where {N<:Int} end function _randchannel(nmodes::N) where {N<:Int} disp = randn(2*nmodes) - # generate symplectic matrix describing the evolution of the system - # and environment - m = rand(1:nmodes) # number of modes in environment + # generate symplectic matrix describing the evolution of the system with N modes + # and environment with 2N modes + m = 2*nmodes # number of modes in environment symp = randsymplectic(nmodes + m) - transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:2*(nmodes+m)]) + transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:2*(nmodes + m)]) # generate random covariance matrix of environment - envcovar = zeros(2*m, 2*m) envsymp = randsymplectic(m) - # create buffers for matrix multiplication - buf1, buf2 = zeros(2*m, 2*m), zeros(2*m, 2*nmodes) - # William decomposition for mixed Gaussian states + # create buffer for matrix multiplication + buf = zeros(2*m, 2*m) + # William decomposition for Gaussian states sympeigs = abs.(rand(m)) .+ 1.0 will = diagm(repeat(sympeigs, inner = 2)) mul!(envcovar, envsymp, mul!(buf1, will, envsymp')) diff --git a/test/test_randoms.jl b/test/test_randoms.jl index 0a664b7..563af1f 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -6,35 +6,92 @@ @testset "random utils" begin nmodes = rand(1:20) U = Gabs._rand_unitary(nmodes) - @test isapprox(adjoint(U), inv(U)) + @test isapprox(adjoint(U), inv(U), atol = 1e-5) O = Gabs._rand_orthogonal_symplectic(nmodes) - @test isapprox(O', inv(O)) + @test isapprox(O', inv(O), atol = 1e-5) Omega = symplecticform(nmodes) - @test isapprox(O' * Omega * O, Omega) + @test isapprox(O' * Omega * O, Omega, atol = 1e-5) + + Spassive = randsymplectic(nmodes, passive = true) + @test isapprox(Spassive', inv(Spassive), atol = 1e-5) + @test isapprox(Spassive' * Omega * Spassive, Omega, atol = 1e-5) + + S = randsymplectic(nmodes) + @test isapprox(S' * Omega * S, Omega, atol = 1e-5) end - @testset "random types" begin + @testset "random states" begin nmodes = rand(1:20) rs = randstate(nmodes) rc = randchannel(nmodes) @test rc isa GaussianChannel @test rc * rs isa GaussianState Omega = symplecticform(nmodes) - @test all(i -> (i >= 0) || isapprox(i, 0.0, atol = 1e-5), real(eigvals(rs.covar .+ im*Omega))) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs.covar .+ im*Omega))) rspure = randstate(nmodes, pure = true) @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rspure.covar .+ im*Omega))) @test isapprox(purity(rspure), 1.0, atol = 1e-5) - rsfloat32 = randstate(AbstractArray{Float32}, nmodes) - rcfloat32 = randchannel(AbstractArray{Float32}, nmodes) - @test rcfloat32 isa GaussianChannel - @test rcfloat32 * rsfloat32 isa GaussianState + rs_float32 = randstate(AbstractArray{Float32}, nmodes) + rc_float32 = randchannel(AbstractArray{Float32}, nmodes) + @test rc_float32 isa GaussianChannel + @test rc_float32 * rs_float32 isa GaussianState + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs_float32.covar .+ im*Omega))) + + rspure_float32 = randstate(AbstractArray{Float32}, nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rspure_float32.covar .+ im*Omega))) + @test isapprox(purity(rspure_float32), 1.0, atol = 1e-5) + + rs_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) + rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) + @test rc_static isa GaussianChannel + @test rc_static * rs_static isa GaussianState + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rs_static.covar) .+ im*Omega))) + + rspure_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rspure_static.covar) .+ im*Omega))) + @test isapprox(purity(rspure_float32), 1.0, atol = 1e-5) + end + + @testset "random unitaries" begin + nmodes = rand(1:20) + ru = randunitary(nmodes) + Omega = symplecticform(nmodes) + @test isapprox(ru.symplectic' * Omega * ru.symplectic, Omega, atol = 1e-5) + + rupassive = randunitary(nmodes, passive = true) + @test isapprox(rupassive.symplectic', inv(rupassive.symplectic), atol = 1e-5) + @test isapprox(rupassive.symplectic' * Omega * rupassive.symplectic, Omega, atol = 1e-5) + + ru_float32 = randunitary(AbstractArray{Float32}, nmodes) + @test isapprox(ru_float32.symplectic' * Omega * ru_float32.symplectic, Omega, atol = 1e-5) + + rupassive_float32 = randunitary(nmodes, passive = true) + @test isapprox(rupassive_float32.symplectic', inv(rupassive_float32.symplectic), atol = 1e-5) + @test isapprox(rupassive_float32.symplectic' * Omega * rupassive_float32.symplectic, Omega, atol = 1e-5) + + ru_static = randunitary(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) + @test isapprox(ru_static.symplectic' * Omega * ru_static.symplectic, Omega, atol = 1e-5) + + rupassive_static = randunitary(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes, passive = true) + @test isapprox(rupassive_static.symplectic', inv(rupassive_static.symplectic), atol = 1e-5) + @test isapprox(rupassive_static.symplectic' * Omega * rupassive_static.symplectic, Omega, atol = 1e-5) + end + + @testset "random channels" begin + #= + nmodes = rand(1:20) + rc = randchannel(nmodes) + Omega = symplecticform(nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc.noise .+ im*Omega .- im*rc.transform*Omega*rc.transform'))) + + rc_float32 = randchannel(AbstractArray{Float32}, nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc_float32.noise .+ im*Omega .- im*rc_float32.transform'*Omega*rc_float32.transform))) - rsstatic = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) - rcstatic = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) - @test rcstatic isa GaussianChannel - @test rcstatic * rsstatic isa GaussianState + rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes, 2*nmodes}, nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rc_static.noise) .+ im*Omega .- im*Array(rc_static.transform)'*Omega*Array(rc_static.transform)))) + =# end end \ No newline at end of file From d95fc50770cce0da22acde6c2c13d5127ada353a Mon Sep 17 00:00:00 2001 From: apkille Date: Wed, 13 Nov 2024 11:46:15 -0500 Subject: [PATCH 3/8] fix --- src/randoms.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/randoms.jl b/src/randoms.jl index a5987f5..0632cd9 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -73,15 +73,12 @@ function _randchannel(nmodes::N) where {N<:Int} transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:2*(nmodes + m)]) # generate random covariance matrix of environment envsymp = randsymplectic(m) - # create buffer for matrix multiplication - buf = zeros(2*m, 2*m) - # William decomposition for Gaussian states - sympeigs = abs.(rand(m)) .+ 1.0 - will = diagm(repeat(sympeigs, inner = 2)) - mul!(envcovar, envsymp, mul!(buf1, will, envsymp')) + envcovar = envsymp * envsymp' + # create buffers for matrix multiplication + buf = zeros(2*m, 2*nmodes) # generate noise matrix from evolution of environment noise = zeros(2*nmodes, 2*nmodes) - mul!(noise, B, mul!(buf2, envcovar, B')) + mul!(noise, B, mul!(buf, envcovar, B')) return disp, transform, noise end From e7231fc120dca8b1b678ef6c975a9b44ae870cc0 Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 17 Nov 2024 20:20:21 -0500 Subject: [PATCH 4/8] change basis of symplecticform --- src/randoms.jl | 45 ++++++++++++++++++++++---------------------- src/wigner.jl | 20 +++++++------------- test/test_randoms.jl | 16 +++++++--------- test/test_wigner.jl | 6 +++--- 4 files changed, 39 insertions(+), 48 deletions(-) diff --git a/src/randoms.jl b/src/randoms.jl index 0632cd9..e13e4eb 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -68,17 +68,11 @@ function _randchannel(nmodes::N) where {N<:Int} disp = randn(2*nmodes) # generate symplectic matrix describing the evolution of the system with N modes # and environment with 2N modes - m = 2*nmodes # number of modes in environment - symp = randsymplectic(nmodes + m) - transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:2*(nmodes + m)]) - # generate random covariance matrix of environment - envsymp = randsymplectic(m) - envcovar = envsymp * envsymp' - # create buffers for matrix multiplication - buf = zeros(2*m, 2*nmodes) + symp = randsymplectic(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) - mul!(noise, B, mul!(buf, envcovar, B')) + mul!(noise, B, B') return disp, transform, noise end @@ -97,7 +91,12 @@ function randsymplectic(nmodes::N; passive = false) where {N<:Int} O′ = _rand_orthogonal_symplectic(nmodes) # direct sum of symplectic matrices for single-mode squeeze transformations rs = rand(nmodes) - squeezes = diagm(vcat(exp.(-rs), exp.(rs))) + squeezes = zeros(2*nmodes, 2*nmodes) + @inbounds for i in Base.OneTo(nmodes) + val = rs[i] + squeezes[2*i-1, 2*i-1] = val + squeezes[2*i, 2*i] = 1/val + end return O * squeezes * O′ end function randsymplectic(::Type{T}, nmodes::N; passive = false) where {T, N<:Int} @@ -105,15 +104,6 @@ function randsymplectic(::Type{T}, nmodes::N; passive = false) where {T, N<:Int} return T(symp) end -# Generates unitary matrix randomly distributed over Haar measure; -# see https://arxiv.org/abs/math-ph/0609050 for algorithm. -# This approach is much faster than using rand(Haar(2), nmodes) from RandomMatrices.jl -function _rand_unitary(nmodes::N) where {N<:Int} - M = rand(ComplexF64, nmodes, nmodes) ./ sqrt(2.0) - q, r = qr(M) - d = diagm([r[i, i] / abs(r[i, i]) for i in Base.OneTo(nmodes)]) - return q * d -end # Generates random orthogonal symplectic matrix by blocking real # and imaginary parts of a random unitary matrix function _rand_orthogonal_symplectic(nmodes::N) where {N<:Int} @@ -121,10 +111,19 @@ function _rand_orthogonal_symplectic(nmodes::N) where {N<:Int} O = zeros(2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes), j in Base.OneTo(nmodes) val = U[i,j] - O[i,j] = real(val) - O[i+nmodes, j] = imag(val) - O[i, j+nmodes] = -imag(val) - O[i+nmodes, j+nmodes] = real(val) + O[2*i-1,2*j-1] = real(val) + O[2*i, 2*j-1] = -imag(val) + O[2*i-1, 2*j] = imag(val) + O[2*i, 2*j] = real(val) 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), size) from RandomMatrices.jl +function _rand_unitary(size::N) where {N<:Int} + M = rand(ComplexF64, size, size) ./ sqrt(2.0) + q, r = qr(M) + d = diagm([r[i, i] / abs(r[i, i]) for i in Base.OneTo(size)]) + return q * d end \ No newline at end of file diff --git a/src/wigner.jl b/src/wigner.jl index b6e2987..5b95fb0 100644 --- a/src/wigner.jl +++ b/src/wigner.jl @@ -1,19 +1,13 @@ """ - symplecticform([T = Matrix{Float64},] modes<:Int) + symplecticform([T = Matrix{Float64},] nmodes<:Int) -Compute the symplectic form matrix of size 2N x 2N, where N is given by `modes`. +Compute the symplectic form matrix of size 2N x 2N, where N is given by `nmodes`. """ -function symplecticform(modes::N) where {N<:Int} - Omega = zeros(2modes, 2modes) - @inbounds for i in 1:modes, j in modes:2modes - if isequal(i, j-modes) - Omega[i,j] = 1.0 - end - end - @inbounds for i in modes:2modes, j in 1:modes - if isequal(i-modes,j) - Omega[i, j] = -1.0 - end +function symplecticform(nmodes::N) where {N<:Int} + Omega = zeros(2*nmodes, 2*nmodes) + @inbounds for i in Base.OneTo(nmodes) + Omega[2*i-1, 2*i] = 1.0 + Omega[2*i, 2*i-1] = -1.0 end return Omega end diff --git a/test/test_randoms.jl b/test/test_randoms.jl index 563af1f..bf618df 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -41,8 +41,8 @@ @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs_float32.covar .+ im*Omega))) rspure_float32 = randstate(AbstractArray{Float32}, nmodes, pure = true) - @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rspure_float32.covar .+ im*Omega))) - @test isapprox(purity(rspure_float32), 1.0, atol = 1e-5) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-3)), real(eigvals(rspure_float32.covar .+ im*Omega))) + @test isapprox(purity(rspure_float32), 1.0, atol = 1e-3) rs_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) @@ -81,17 +81,15 @@ end @testset "random channels" begin - #= - nmodes = rand(1:20) - rc = randchannel(nmodes) - Omega = symplecticform(nmodes) + nmodes = rand(1:20); + rc = randchannel(nmodes); + Omega = symplecticform(nmodes); @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc.noise .+ im*Omega .- im*rc.transform*Omega*rc.transform'))) rc_float32 = randchannel(AbstractArray{Float32}, nmodes) - @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc_float32.noise .+ im*Omega .- im*rc_float32.transform'*Omega*rc_float32.transform))) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc_float32.noise .+ im*Omega .- im*rc_float32.transform*Omega*rc_float32.transform'))) rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes, 2*nmodes}, nmodes) - @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rc_static.noise) .+ im*Omega .- im*Array(rc_static.transform)'*Omega*Array(rc_static.transform)))) - =# + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rc_static.noise) .+ im*Omega .- im*Array(rc_static.transform)*Omega*Array(rc_static.transform)'))) end end \ No newline at end of file diff --git a/test/test_wigner.jl b/test/test_wigner.jl index 38664a7..aad9729 100644 --- a/test/test_wigner.jl +++ b/test/test_wigner.jl @@ -5,10 +5,10 @@ @testset "symplectic form" begin Omega = symplecticform(2) Omega_static = symplecticform(SMatrix{4, 4}, 2) - test_Omega = [0.0 0.0 1.0 0.0; + test_Omega = [0.0 1.0 0.0 0.0; + -1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0; - -1.0 0.0 0.0 0.0; - 0.0 -1.0 0.0 0.0] + 0.0 0.0 -1.0 0.0] @test isequal(Omega, test_Omega) @test Omega_static isa SMatrix end From bb84d0ab5256ea689ff86e7c3cbd75f31a045e79 Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 17 Nov 2024 20:29:05 -0500 Subject: [PATCH 5/8] update project.toml and CHANGELOG --- CHANGELOG.md | 2 ++ Project.toml | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fc4c007..aea186b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ ## v1.0.3 - dev - Add `prob` function for `Generaldyne` type. +- Add features for random generation of Gaussian states, unitaries, and channels. +- **(breaking)** Switched basis for matrices generated by `symplecticform` from block diagonal to canonical. ## v1.0.2 - 2024-11-03 diff --git a/Project.toml b/Project.toml index 4f05a61..9fb5e66 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gabs" uuid = "0eb812ee-a11f-4f5e-b8d4-bb8a44f06f50" authors = ["Andrew Kille"] -version = "1.0.3-dev" +version = "1.1.0" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" From 7e109cf5dfa6ccb82fbccc35feb33613ef99770d Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 17 Nov 2024 20:35:56 -0500 Subject: [PATCH 6/8] alter Float32 test --- test/test_randoms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_randoms.jl b/test/test_randoms.jl index bf618df..34e1156 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -42,7 +42,7 @@ rspure_float32 = randstate(AbstractArray{Float32}, nmodes, pure = true) @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-3)), real(eigvals(rspure_float32.covar .+ im*Omega))) - @test isapprox(purity(rspure_float32), 1.0, atol = 1e-3) + @test isapprox(Float64(purity(rspure_float32)), 1.0, atol = 1e-3) rs_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) From 10bce54d46af62ebc748bfbbc30dedab48d1c305 Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 17 Nov 2024 20:40:12 -0500 Subject: [PATCH 7/8] add codecov for `randsymplectic` --- src/randoms.jl | 2 +- test/test_randoms.jl | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/randoms.jl b/src/randoms.jl index e13e4eb..91fc2d7 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -100,7 +100,7 @@ function randsymplectic(nmodes::N; passive = false) where {N<:Int} return O * squeezes * O′ end function randsymplectic(::Type{T}, nmodes::N; passive = false) where {T, N<:Int} - symp = randsymplectic(nmodes, passive) + symp = randsymplectic(nmodes, passive = passive) return T(symp) end diff --git a/test/test_randoms.jl b/test/test_randoms.jl index 34e1156..3ac308f 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -19,6 +19,9 @@ S = randsymplectic(nmodes) @test isapprox(S' * Omega * S, Omega, atol = 1e-5) + + S_float32 = randsymplectic(Array{Float32}, nmodes) + @test isapprox(S' * Omega * S, Omega, atol = 1e-5) end @testset "random states" begin From da0d9f376607fedbc7032b280e4363d15111cc47 Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 17 Nov 2024 20:49:04 -0500 Subject: [PATCH 8/8] change floating point tests for rand objects --- test/test_randoms.jl | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/test/test_randoms.jl b/test/test_randoms.jl index 3ac308f..9c21bca 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -20,8 +20,8 @@ S = randsymplectic(nmodes) @test isapprox(S' * Omega * S, Omega, atol = 1e-5) - S_float32 = randsymplectic(Array{Float32}, nmodes) - @test isapprox(S' * Omega * S, Omega, atol = 1e-5) + S_array = randsymplectic(Array, nmodes) + @test isapprox(S_array' * Omega * S_array, Omega, atol = 1e-5) end @testset "random states" begin @@ -37,15 +37,15 @@ @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rspure.covar .+ im*Omega))) @test isapprox(purity(rspure), 1.0, atol = 1e-5) - rs_float32 = randstate(AbstractArray{Float32}, nmodes) - rc_float32 = randchannel(AbstractArray{Float32}, nmodes) - @test rc_float32 isa GaussianChannel - @test rc_float32 * rs_float32 isa GaussianState - @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs_float32.covar .+ im*Omega))) + rs_array = randstate(Array, nmodes) + rc_array = randchannel(Array, nmodes) + @test rc_array isa GaussianChannel + @test rc_array * rs_array isa GaussianState + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs_array.covar .+ im*Omega))) - rspure_float32 = randstate(AbstractArray{Float32}, nmodes, pure = true) - @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-3)), real(eigvals(rspure_float32.covar .+ im*Omega))) - @test isapprox(Float64(purity(rspure_float32)), 1.0, atol = 1e-3) + rspure_array = randstate(Array, nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-3)), real(eigvals(rspure_array.covar .+ im*Omega))) + @test isapprox(purity(rspure_array), 1.0, atol = 1e-3) rs_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) @@ -55,7 +55,7 @@ rspure_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes, pure = true) @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rspure_static.covar) .+ im*Omega))) - @test isapprox(purity(rspure_float32), 1.0, atol = 1e-5) + @test isapprox(purity(rspure_static), 1.0, atol = 1e-5) end @testset "random unitaries" begin @@ -68,12 +68,12 @@ @test isapprox(rupassive.symplectic', inv(rupassive.symplectic), atol = 1e-5) @test isapprox(rupassive.symplectic' * Omega * rupassive.symplectic, Omega, atol = 1e-5) - ru_float32 = randunitary(AbstractArray{Float32}, nmodes) - @test isapprox(ru_float32.symplectic' * Omega * ru_float32.symplectic, Omega, atol = 1e-5) + ru_array = randunitary(Array, nmodes) + @test isapprox(ru_array.symplectic' * Omega * ru_array.symplectic, Omega, atol = 1e-5) - rupassive_float32 = randunitary(nmodes, passive = true) - @test isapprox(rupassive_float32.symplectic', inv(rupassive_float32.symplectic), atol = 1e-5) - @test isapprox(rupassive_float32.symplectic' * Omega * rupassive_float32.symplectic, Omega, atol = 1e-5) + rupassive_array = randunitary(nmodes, passive = true) + @test isapprox(rupassive_array.symplectic', inv(rupassive_array.symplectic), atol = 1e-5) + @test isapprox(rupassive_array.symplectic' * Omega * rupassive_array.symplectic, Omega, atol = 1e-5) ru_static = randunitary(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) @test isapprox(ru_static.symplectic' * Omega * ru_static.symplectic, Omega, atol = 1e-5) @@ -89,8 +89,8 @@ Omega = symplecticform(nmodes); @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc.noise .+ im*Omega .- im*rc.transform*Omega*rc.transform'))) - rc_float32 = randchannel(AbstractArray{Float32}, nmodes) - @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc_float32.noise .+ im*Omega .- im*rc_float32.transform*Omega*rc_float32.transform'))) + rc_array = randchannel(Array, nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc_array.noise .+ im*Omega .- im*rc_array.transform*Omega*rc_array.transform'))) rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes, 2*nmodes}, nmodes) @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rc_static.noise) .+ im*Omega .- im*Array(rc_static.transform)*Omega*Array(rc_static.transform)')))