From 54875ae2f24d2c863a56a27717298594b709ade0 Mon Sep 17 00:00:00 2001 From: Andrew Kille <68079167+apkille@users.noreply.github.com> Date: Wed, 25 Dec 2024 23:56:57 -0500 Subject: [PATCH] Propagate types in arrays to propagate ForwardDiff duals (#23) * promote types to propagate ForwardDiff duals * add ForwardDiff and FiniteDiff tests --- src/channels.jl | 56 ++++++++++++++----------- src/states.jl | 74 +++++++++++++++++++------------- src/unitaries.jl | 98 ++++++++++++++++++++++++------------------- test/Project.toml | 3 ++ test/test_autodiff.jl | 48 +++++++++++++++++++++ 5 files changed, 183 insertions(+), 96 deletions(-) create mode 100644 test/test_autodiff.jl diff --git a/src/channels.jl b/src/channels.jl index 6ee7460..86baab9 100644 --- a/src/channels.jl +++ b/src/channels.jl @@ -100,16 +100,17 @@ function attenuator(basis::SymplecticBasis{N}, theta::R, n::M) where {N<:Int,R,M end function _attenuator(basis::Union{QuadPairBasis{N},QuadBlockBasis{N}}, theta::R, n::M) where {N<:Int,R<:Real,M<:Int} nmodes = basis.nmodes - disp = zeros(2*nmodes) - transform = Matrix{Float64}(cos(theta) * I, 2*nmodes, 2*nmodes) - noise = Matrix{Float64}((sin(theta))^2 * n * I, 2*nmodes, 2*nmodes) + disp = zeros(R, 2*nmodes) + transform = Matrix{R}(cos(theta) * I, 2*nmodes, 2*nmodes) + noise = Matrix{R}((sin(theta))^2 * n * I, 2*nmodes, 2*nmodes) return disp, transform, noise end function _attenuator(basis::QuadPairBasis{N}, theta::R, n::M) where {N<:Int,R<:Vector,M<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - transform = zeros(2*nmodes, 2*nmodes) - noise = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + transform = zeros(Rt, 2*nmodes, 2*nmodes) + noise = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) ct, st = cos(theta[i]), sin(theta[i]) ni = n[i] @@ -124,9 +125,10 @@ function _attenuator(basis::QuadPairBasis{N}, theta::R, n::M) where {N<:Int,R<:V end function _attenuator(basis::QuadBlockBasis{N}, theta::R, n::M) where {N<:Int,R<:Vector,M<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - transform = zeros(2*nmodes, 2*nmodes) - noise = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + transform = zeros(Rt, 2*nmodes, 2*nmodes) + noise = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) ct, st = cos(theta[i]), sin(theta[i]) ni = n[i] @@ -188,16 +190,17 @@ function amplifier(basis::SymplecticBasis{N}, r::R, n::M) where {N<:Int,R,M} end function _amplifier(basis::Union{QuadPairBasis{N},QuadBlockBasis{N}}, r::R, n::M) where {N<:Int,R<:Real,M<:Int} nmodes = basis.nmodes - disp = zeros(2*nmodes) - transform = Matrix{Float64}(cosh(r) * I, 2*nmodes, 2*nmodes) - noise = Matrix{Float64}((sinh(r))^2 * n * I, 2*nmodes, 2*nmodes) + disp = zeros(R, 2*nmodes) + transform = Matrix{R}(cosh(r) * I, 2*nmodes, 2*nmodes) + noise = Matrix{R}((sinh(r))^2 * n * I, 2*nmodes, 2*nmodes) return disp, transform, noise end function _amplifier(basis::QuadPairBasis{N}, r::R, n::M) where {N<:Int,R<:Vector,M<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - transform = zeros(2*nmodes, 2*nmodes) - noise = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + transform = zeros(Rt, 2*nmodes, 2*nmodes) + noise = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) cr, sr = cos(r[i]), sin(r[i]) ni = n[i] @@ -212,9 +215,10 @@ function _amplifier(basis::QuadPairBasis{N}, r::R, n::M) where {N<:Int,R<:Vector end function _amplifier(basis::QuadBlockBasis{N}, r::R, n::M) where {N<:Int,R<:Vector,M<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - transform = zeros(2*nmodes, 2*nmodes) - noise = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + transform = zeros(Rt, 2*nmodes, 2*nmodes) + noise = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) cr, sr = cos(r[i]), sin(r[i]) ni = n[i] @@ -248,7 +252,8 @@ function _tensor(op1::GaussianChannel{B1,D1,T1}, op2::GaussianChannel{B2,D2,T2}) block1, block2 = Base.OneTo(2*nmodes1), Base.OneTo(2*nmodes2) # initialize direct sum of displacement vectors disp1, disp2 = op1.disp, op2.disp - disp′ = zeros(2*nmodes) + Dt = promote_type(eltype(disp1), eltype(disp2)) + disp′ = zeros(Dt, 2*nmodes) @inbounds for i in block1 disp′[i] = disp1[i] end @@ -257,9 +262,10 @@ function _tensor(op1::GaussianChannel{B1,D1,T1}, op2::GaussianChannel{B2,D2,T2}) end # initialize direct sum of transform and noise matrices trans1, trans2 = op1.transform, op2.transform - transform′ = zeros(2*nmodes, 2*nmodes) + Tt = promote_type(eltype(trans1), eltype(trans2)) + transform′ = zeros(Tt, 2*nmodes, 2*nmodes) noise1, noise2 = op1.noise, op2.noise - noise′ = zeros(2*nmodes, 2*nmodes) + noise′ = zeros(Tt, 2*nmodes, 2*nmodes) @inbounds for i in block1, j in block1 transform′[i,j] = trans1[i,j] noise′[i,j] = noise1[i,j] @@ -281,7 +287,8 @@ function _tensor(op1::GaussianChannel{B1,D1,T1}, op2::GaussianChannel{B2,D2,T2}) block1, block2 = Base.OneTo(nmodes1), Base.OneTo(nmodes2) # initialize direct sum of displacement vectors disp1, disp2 = op1.disp, op2.disp - disp′ = zeros(2*nmodes) + Dt = promote_type(eltype(disp1), eltype(disp2)) + disp′ = zeros(Dt, 2*nmodes) @inbounds for i in block1 disp′[i] = disp1[i] disp′[i+nmodes] = disp1[i+nmodes1] @@ -292,9 +299,10 @@ function _tensor(op1::GaussianChannel{B1,D1,T1}, op2::GaussianChannel{B2,D2,T2}) end # initialize direct sum of transform and noise matrices trans1, trans2 = op1.transform, op2.transform - transform′ = zeros(2*nmodes, 2*nmodes) + Tt = promote_type(eltype(trans1), eltype(trans2)) + transform′ = zeros(Tt, 2*nmodes, 2*nmodes) noise1, noise2 = op1.noise, op2.noise - noise′ = zeros(2*nmodes, 2*nmodes) + noise′ = zeros(Tt, 2*nmodes, 2*nmodes) @inbounds for i in block1, j in block1 transform′[i,j] = trans1[i,j] transform′[i,j+nmodes] = trans1[i,j+nmodes1] diff --git a/src/states.jl b/src/states.jl index 60dfc48..9fe7f59 100644 --- a/src/states.jl +++ b/src/states.jl @@ -158,27 +158,29 @@ end function _coherentstate(basis::QuadPairBasis{N}, alpha::A) where {N<:Int,A<:Number} nmodes = basis.nmodes mean = repeat([sqrt(2)*real(alpha), sqrt(2)*imag(alpha)], nmodes) - covar = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + covar = Matrix{real(A)}(I, 2*nmodes, 2*nmodes) return mean, covar end function _coherentstate(basis::QuadPairBasis{N}, alpha::A) where {N<:Int,A<:Vector} nmodes = basis.nmodes - mean = sqrt(2) * reinterpret(Float64, alpha) - covar = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + Rt = real(eltype(A)) + mean = sqrt(2) * reinterpret(Rt, alpha) + covar = Matrix{Rt}(I, 2*nmodes, 2*nmodes) return mean, covar end function _coherentstate(basis::QuadBlockBasis{N}, alpha::A) where {N<:Int,A<:Number} nmodes = basis.nmodes mean = repeat([sqrt(2)*real(alpha), sqrt(2)*imag(alpha)], inner = nmodes) - covar = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + covar = Matrix{real(A)}(I, 2*nmodes, 2*nmodes) return mean, covar end function _coherentstate(basis::QuadBlockBasis{N}, alpha::A) where {N<:Int,A<:Vector} nmodes = basis.nmodes - re = reinterpret(Float64, alpha) + Rt = real(eltype(A)) + re = reinterpret(Rt, alpha) mean = vcat(@view(re[1:2:end]), @view(re[2:2:end])) mean .*= sqrt(2) - covar = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + covar = Matrix{Rt}(I, 2*nmodes, 2*nmodes) return mean, covar end @@ -228,8 +230,8 @@ function squeezedstate(basis::SymplecticBasis{N}, r::R, theta::R) where {N<:Int, end function _squeezedstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - mean = zeros(2*nmodes) - covar = zeros(2*nmodes, 2*nmodes) + mean = zeros(R, 2*nmodes) + covar = zeros(R, 2*nmodes, 2*nmodes) cr, sr = cosh(2*r), sinh(2*r) ct, st = cos(theta), sin(theta) @inbounds for i in Base.OneTo(nmodes) @@ -242,8 +244,9 @@ function _squeezedstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R end function _squeezedstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - mean = zeros(2*nmodes) - covar = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + mean = zeros(Rt, 2*nmodes) + covar = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) cr, sr = cosh(2*r[i]), sinh(2*r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -256,8 +259,8 @@ function _squeezedstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R end function _squeezedstate(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - mean = zeros(2*nmodes) - covar = zeros(2*nmodes, 2*nmodes) + mean = zeros(R, 2*nmodes) + covar = zeros(R, 2*nmodes, 2*nmodes) cr, sr = cosh(2*r), sinh(2*r) ct, st = cos(theta), sin(theta) @inbounds for i in Base.OneTo(nmodes) @@ -270,8 +273,9 @@ function _squeezedstate(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int, end function _squeezedstate(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - mean = zeros(2*nmodes) - covar = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + mean = zeros(Rt, 2*nmodes) + covar = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) cr, sr = cosh(2*r[i]), sinh(2*r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -335,10 +339,10 @@ function eprstate(basis::SymplecticBasis{N}, r::R, theta::R) where {N<:Int,R} end function _eprstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - mean = zeros(2*nmodes) + mean = zeros(R, 2*nmodes) cr, sr = (1/2)*cosh(2*r), (1/2)*sinh(2*r) ct, st = cos(theta), sin(theta) - covar = zeros(2*nmodes, 2*nmodes) + covar = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) covar[4*i-3, 4*i-3] = cr covar[4*i-3, 4*i-1] = -sr * ct @@ -360,8 +364,9 @@ function _eprstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Rea end function _eprstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - mean = zeros(2*nmodes) - covar = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + mean = zeros(Rt, 2*nmodes) + covar = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) cr, sr = (1/2)*cosh(2*r[i]), (1/2)*sinh(2*r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -386,10 +391,10 @@ function _eprstate(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Vec end function _eprstate(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - mean = zeros(2*nmodes) + mean = zeros(R, 2*nmodes) cr, sr = (1/2)*cosh(2*r), (1/2)*sinh(2*r) ct, st = cos(theta), sin(theta) - covar = zeros(2*nmodes, 2*nmodes) + covar = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) covar[2*i-1, 2*i-1] = cr covar[2*i-1, 2*i] = -sr * ct @@ -411,8 +416,9 @@ function _eprstate(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Re end function _eprstate(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - mean = zeros(2*nmodes) - covar = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + mean = zeros(Rt, 2*nmodes) + covar = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) cr, sr = (1/2)*cosh(2*r[i]), (1/2)*sinh(2*r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -477,12 +483,13 @@ function tensor(state1::GaussianState, state2::GaussianState) end function _tensor(state1::GaussianState{B1,M1,V1}, state2::GaussianState{B2,M2,V2}) where {B1<:QuadPairBasis,B2<:QuadPairBasis,M1,M2,V1,V2} mean1, mean2 = state1.mean, state2.mean + Mt = promote_type(eltype(mean1), eltype(mean2)) basis1, basis2 = state1.basis, state2.basis nmodes1, nmodes2 = basis1.nmodes, basis2.nmodes nmodes = nmodes1 + nmodes2 block1, block2 = Base.OneTo(2*nmodes1), Base.OneTo(2*nmodes2) # initialize direct sum of mean vectors - mean′ = zeros(2*nmodes) + mean′ = zeros(Mt, 2*nmodes) @inbounds for i in block1 mean′[i] = mean1[i] end @@ -491,7 +498,8 @@ function _tensor(state1::GaussianState{B1,M1,V1}, state2::GaussianState{B2,M2,V2 end # initialize direct sum of covariance matrices covar1, covar2 = state1.covar, state2.covar - covar′ = zeros(2*nmodes, 2*nmodes) + Vt = promote_type(eltype(covar1), eltype(covar2)) + covar′ = zeros(Vt, 2*nmodes, 2*nmodes) @inbounds for i in block1, j in block1 covar′[i,j] = covar1[i,j] end @@ -505,13 +513,14 @@ function _tensor(state1::GaussianState{B1,M1,V1}, state2::GaussianState{B2,M2,V2 end function _tensor(state1::GaussianState{B1,M1,V1}, state2::GaussianState{B2,M2,V2}) where {B1<:QuadBlockBasis,B2<:QuadBlockBasis,M1,M2,V1,V2} mean1, mean2 = state1.mean, state2.mean + Mt = promote_type(eltype(mean1), eltype(mean2)) basis1, basis2 = state1.basis, state2.basis nmodes1, nmodes2 = basis1.nmodes, basis2.nmodes nmodes = nmodes1 + nmodes2 block1, block2 = Base.OneTo(nmodes1), Base.OneTo(nmodes2) # initialize direct sum of mean vectors mean1, mean2 = state1.mean, state2.mean - mean′ = zeros(2*nmodes) + mean′ = zeros(Mt, 2*nmodes) @inbounds for i in block1 mean′[i] = mean1[i] mean′[i+nmodes] = mean1[i+nmodes1] @@ -522,7 +531,8 @@ function _tensor(state1::GaussianState{B1,M1,V1}, state2::GaussianState{B2,M2,V2 end # initialize direct sum of covariance matrices covar1, covar2 = state1.covar, state2.covar - covar′ = zeros(2*nmodes, 2*nmodes) + Vt = promote_type(eltype(covar1), eltype(covar2)) + covar′ = zeros(Vt, 2*nmodes, 2*nmodes) @inbounds for i in block1, j in block1 covar′[i,j] = covar1[i,j] covar′[i,j+nmodes] = covar1[i,j+nmodes1] @@ -629,16 +639,18 @@ end function _ptrace(state::GaussianState{B,M,V}, indices::T) where {B<:QuadPairBasis,M,V,T<:AbstractVector} idxlength = length(indices) mean = state.mean + Mt = eltype(mean) covar = state.covar + Vt = eltype(covar) # initialize partial trace of mean vector - mean′ = zeros(2*idxlength) + mean′ = zeros(Mt, 2*idxlength) @inbounds for i in eachindex(indices) idx = indices[i] mean′[2*i-1] = mean[2*idx-1] mean′[2*i] = mean[2*idx] end # initialize partial trace of covariance matrix - covar′ = zeros(2*idxlength, 2*idxlength) + covar′ = zeros(Vt, 2*idxlength, 2*idxlength) @inbounds for i in eachindex(indices) idx = indices[i] covar′[2*i-1, 2*i-1] = covar[2*idx-1, 2*idx-1] @@ -669,16 +681,18 @@ function _ptrace(state::GaussianState{B,M,V}, indices::T) where {B<:QuadBlockBas nmodes = basis.nmodes idxlength = length(indices) mean = state.mean + Mt = eltype(mean) covar = state.covar + Vt = eltype(covar) # initialize partial trace of mean vector - mean′ = zeros(2*idxlength) + mean′ = zeros(Mt, 2*idxlength) @inbounds for i in eachindex(indices) idx = indices[i] mean′[i] = mean[idx] mean′[i+idxlength] = mean[idx+nmodes] end # initialize partial trace of covariance matrix - covar′ = zeros(2*idxlength, 2*idxlength) + covar′ = zeros(Vt, 2*idxlength, 2*idxlength) @inbounds for i in Base.OneTo(idxlength) idx = indices[i] @inbounds for j in i:idxlength diff --git a/src/unitaries.jl b/src/unitaries.jl index 4b9f5d4..ffe3fa8 100644 --- a/src/unitaries.jl +++ b/src/unitaries.jl @@ -49,27 +49,29 @@ end function _displace(basis::QuadPairBasis{N}, alpha::A) where {N<:Int,A<:Number} nmodes = basis.nmodes disp = repeat([sqrt(2)*real(alpha), sqrt(2)*imag(alpha)], nmodes) - symplectic = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + symplectic = Matrix{real(A)}(I, 2*nmodes, 2*nmodes) return disp, symplectic end function _displace(basis::QuadPairBasis{N}, alpha::A) where {N<:Int,A<:Vector} nmodes = basis.nmodes - disp = sqrt(2) * reinterpret(Float64, alpha) - symplectic = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + Rt = real(eltype(A)) + disp = sqrt(2) * reinterpret(Rt, alpha) + symplectic = Matrix{Rt}(I, 2*nmodes, 2*nmodes) return disp, symplectic end function _displace(basis::QuadBlockBasis{N}, alpha::A) where {N<:Int,A<:Number} nmodes = basis.nmodes disp = repeat([sqrt(2)*real(alpha), sqrt(2)*imag(alpha)], inner = nmodes) - symplectic = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + symplectic = Matrix{real(A)}(I, 2*nmodes, 2*nmodes) return disp, symplectic end function _displace(basis::QuadBlockBasis{N}, alpha::A) where {N<:Int,A<:Vector} nmodes = basis.nmodes - re = reinterpret(Float64, alpha) + Rt = real(eltype(A)) + re = reinterpret(Rt, alpha) disp = vcat(@view(re[1:2:end]), @view(re[2:2:end])) disp .*= sqrt(2) - symplectic = Matrix{Float64}(I, 2*nmodes, 2*nmodes) + symplectic = Matrix{Rt}(I, 2*nmodes, 2*nmodes) return disp, symplectic end @@ -123,10 +125,10 @@ function squeeze(basis::SymplecticBasis{N}, r::R, theta::R) where {N<:Int, R} end function _squeeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) cr, sr = cosh(r), sinh(r) ct, st = cos(theta), sin(theta) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) symplectic[2*i-1, 2*i-1] = cr - sr*ct symplectic[2*i-1, 2*i] = -sr * st @@ -137,8 +139,9 @@ function _squeeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Real end function _squeeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) cr, sr = cosh(r[i]), sinh(r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -151,10 +154,10 @@ function _squeeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Vect end function _squeeze(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) cr, sr = cosh(r), sinh(r) ct, st = cos(theta), sin(theta) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) symplectic[i, i] = cr - sr*ct symplectic[i, i+nmodes] = -sr * st @@ -165,8 +168,9 @@ function _squeeze(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Rea end function _squeeze(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) cr, sr = cosh(r[i]), sinh(r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -235,10 +239,10 @@ function twosqueeze(basis::SymplecticBasis{N}, r::R, theta::R) where {N<:Int,R} end function _twosqueeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) cr, sr = cosh(r), sinh(r) ct, st = cos(theta), sin(theta) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) symplectic[4*i-3, 4*i-3] = cr symplectic[4*i-3, 4*i-1] = -sr * ct @@ -260,8 +264,9 @@ function _twosqueeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:R end function _twosqueeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) cr, sr = cosh(r[i]), sinh(r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -286,10 +291,10 @@ function _twosqueeze(basis::QuadPairBasis{N}, r::R, theta::R) where {N<:Int,R<:V end function _twosqueeze(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) cr, sr = cosh(r), sinh(r) ct, st = cos(theta), sin(theta) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) for i in Base.OneTo(Int(nmodes/2)) symplectic[2*i-1, 2*i-1] = cr symplectic[2*i-1, 2*i] = -sr * ct @@ -311,8 +316,9 @@ function _twosqueeze(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<: end function _twosqueeze(basis::QuadBlockBasis{N}, r::R, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) cr, sr = cosh(r[i]), sinh(r[i]) ct, st = cos(theta[i]), sin(theta[i]) @@ -387,9 +393,9 @@ function phaseshift(basis::SymplecticBasis{N}, theta::R) where {N<:Int,R} end function _phaseshift(basis::QuadPairBasis{N}, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) ct, st = cos(theta), sin(theta) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) symplectic[2*i-1, 2*i-1] = ct symplectic[2*i-1, 2*i] = st @@ -400,8 +406,9 @@ function _phaseshift(basis::QuadPairBasis{N}, theta::R) where {N<:Int,R<:Real} end function _phaseshift(basis::QuadPairBasis{N}, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) ct, st = cos(theta[i]), sin(theta[i]) symplectic[2*i-1, 2*i-1] = ct @@ -413,9 +420,9 @@ function _phaseshift(basis::QuadPairBasis{N}, theta::R) where {N<:Int,R<:Vector} end function _phaseshift(basis::QuadBlockBasis{N}, theta::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) ct, st = cos(theta), sin(theta) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) symplectic[i, i] = ct symplectic[i, i+nmodes] = st @@ -426,8 +433,9 @@ function _phaseshift(basis::QuadBlockBasis{N}, theta::R) where {N<:Int,R<:Real} end function _phaseshift(basis::QuadBlockBasis{N}, theta::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) ct, st = cos(theta[i]), sin(theta[i]) symplectic[i, i] = ct @@ -494,9 +502,9 @@ function beamsplitter(basis::SymplecticBasis{N}, transmit::R) where {N<:Int,R} end function _beamsplitter(basis::QuadPairBasis{N}, transmit::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) a1, a2 = sqrt(transmit), sqrt(1 - transmit) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) symplectic[4*i-3, 4*i-3] = a2 symplectic[4*i-3, 4*i-1] = a1 @@ -514,8 +522,9 @@ function _beamsplitter(basis::QuadPairBasis{N}, transmit::R) where {N<:Int,R<:Re end function _beamsplitter(basis::QuadPairBasis{N}, transmit::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) a1, a2 = sqrt(transmit[i]), sqrt(1 - transmit[i]) @@ -535,9 +544,9 @@ function _beamsplitter(basis::QuadPairBasis{N}, transmit::R) where {N<:Int,R<:Ve end function _beamsplitter(basis::QuadBlockBasis{N}, transmit::R) where {N<:Int,R<:Real} nmodes = basis.nmodes - disp = zeros(2*nmodes) + disp = zeros(R, 2*nmodes) a1, a2 = sqrt(transmit), sqrt(1 - transmit) - symplectic = zeros(2*nmodes, 2*nmodes) + symplectic = zeros(R, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(nmodes) symplectic[2*i-1, 2*i-1] = a2 symplectic[2*i-1, 2*i] = a1 @@ -548,8 +557,9 @@ function _beamsplitter(basis::QuadBlockBasis{N}, transmit::R) where {N<:Int,R<:R end function _beamsplitter(basis::QuadBlockBasis{N}, transmit::R) where {N<:Int,R<:Vector} nmodes = basis.nmodes - disp = zeros(2*nmodes) - symplectic = zeros(2*nmodes, 2*nmodes) + Rt = eltype(R) + disp = zeros(Rt, 2*nmodes) + symplectic = zeros(Rt, 2*nmodes, 2*nmodes) @inbounds for i in Base.OneTo(Int(nmodes/2)) a1, a2 = sqrt(transmit[i]), sqrt(1 - transmit[i]) @@ -587,7 +597,8 @@ function _tensor(op1::GaussianUnitary{B1,D1,S1}, op2::GaussianUnitary{B2,D2,S2}) block1, block2 = Base.OneTo(2*nmodes1), Base.OneTo(2*nmodes2) # initialize direct sum of displacement vectors disp1, disp2 = op1.disp, op2.disp - disp′ = zeros(2*nmodes) + Dt = promote_type(eltype(disp1), eltype(disp2)) + disp′ = zeros(Dt, 2*nmodes) @inbounds for i in block1 disp′[i] = disp1[i] end @@ -596,7 +607,8 @@ function _tensor(op1::GaussianUnitary{B1,D1,S1}, op2::GaussianUnitary{B2,D2,S2}) end # initialize direct sum of symplectic matrices symp1, symp2 = op1.symplectic, op2.symplectic - symp′ = zeros(2*nmodes, 2*nmodes) + St = promote_type(eltype(symp1), eltype(symp2)) + symp′ = zeros(St, 2*nmodes, 2*nmodes) @inbounds for i in block1, j in block1 symp′[i,j] = symp1[i,j] end @@ -616,7 +628,8 @@ function _tensor(op1::GaussianUnitary{B1,D1,S1}, op2::GaussianUnitary{B2,D2,S2}) block1, block2 = Base.OneTo(nmodes1), Base.OneTo(nmodes2) # initialize direct sum of displacement vectors disp1, disp2 = op1.disp, op2.disp - disp′ = zeros(2*nmodes) + Dt = promote_type(eltype(disp1), eltype(disp2)) + disp′ = zeros(Dt, 2*nmodes) @inbounds for i in block1 disp′[i] = disp1[i] disp′[i+nmodes] = disp1[i+nmodes1] @@ -627,7 +640,8 @@ function _tensor(op1::GaussianUnitary{B1,D1,S1}, op2::GaussianUnitary{B2,D2,S2}) end # initialize direct sum of symplectic matrices symp1, symp2 = op1.symplectic, op2.symplectic - symp′ = zeros(2*nmodes, 2*nmodes) + St = promote_type(eltype(symp1), eltype(symp2)) + symp′ = zeros(St, 2*nmodes, 2*nmodes) @inbounds for i in block1, j in block1 symp′[i,j] = symp1[i,j] symp′[i,j+nmodes] = symp1[i,j+nmodes1] diff --git a/test/Project.toml b/test/Project.toml index fd0b523..da376ff 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,9 @@ [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/test/test_autodiff.jl b/test/test_autodiff.jl new file mode 100644 index 0000000..6b8da9d --- /dev/null +++ b/test/test_autodiff.jl @@ -0,0 +1,48 @@ +@testitem "Automatic differentiation" begin + using Gabs + using DifferentiationInterface + import ForwardDiff + import FiniteDiff + + + @testset "derivatives" begin + + n = rand(1:20) + xs = rand(Float64) + + function f1(x::R) where {R<:Real} + basis = QuadPairBasis(n) + state = vacuumstate(basis) + op = attenuator(basis, x, 2) + newstate = op * state + return abs(1 - purity(newstate))^2 + end + + findiff_f1 = derivative(f1, AutoFiniteDiff(), xs) + fordiff_f1 = derivative(f1, AutoForwardDiff(), xs) + + @test isapprox(findiff_f1, fordiff_f1, atol = 1e-4) + + end + + @testset "gradients" begin + + n = rand(1:5) + xs = rand(n) + + function f1(x::AbstractVector{<:Real}) + basis = QuadPairBasis(length(x)) + state = squeezedstate(basis, x, x) + op = amplifier(basis, x, ones(length(x))) + newstate = op * state + return abs(1 - purity(newstate))^2 + end + + findiff_f1 = gradient(f1, AutoFiniteDiff(), xs) + fordiff_f1 = gradient(f1, AutoForwardDiff(), xs) + + @test isapprox(findiff_f1, fordiff_f1, atol = 1e-4) + + end + + \ No newline at end of file