diff --git a/CHANGELOG.md b/CHANGELOG.md index d48259f..810b368 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # News +## Unreleased + +- Add `heterodyne` measurement interface for Gaussian states. Heterodyne measurement simultaneously measures both quadratures of one or more modes. + ## v1.3.4 - 2026-01-05 - Add optional rng arguments to random methods. diff --git a/src/Gabs.jl b/src/Gabs.jl index 5dd5f1c..89f2498 100644 --- a/src/Gabs.jl +++ b/src/Gabs.jl @@ -15,7 +15,7 @@ export # types GaussianState, GaussianUnitary, GaussianChannel, GaussianLinearCombination, # Gaussian measurements - generaldyne, Generaldyne, homodyne, Homodyne, + generaldyne, Generaldyne, homodyne, Homodyne, heterodyne, Heterodyne, # symplectic representations QuadPairBasis, QuadBlockBasis, changebasis, # operations @@ -66,6 +66,8 @@ include("generaldyne.jl") include("homodyne.jl") +include("heterodyne.jl") + include("wigner.jl") include("metrics.jl") diff --git a/src/heterodyne.jl b/src/heterodyne.jl new file mode 100644 index 0000000..cc6fbc9 --- /dev/null +++ b/src/heterodyne.jl @@ -0,0 +1,337 @@ +struct Heterodyne{R,S<:GaussianState} <: Gabs.AbstractGaussianMeasurement + result::R + state::S + function Heterodyne(r::R, s::S) where {R,S<:GaussianState} + return new{R,S}(r, s) + end +end + +# iteration for destructuring into components +Base.iterate(F::Heterodyne) = (F.result, Val(:state)) +Base.iterate(F::Heterodyne, ::Val{:state}) = (F.state, Val(:done)) +Base.iterate(F::Heterodyne, ::Val{:done}) = nothing + +# printing method +function Base.show( + io::IO, + mime::MIME{Symbol("text/plain")}, + H::Heterodyne{<:Any,<:GaussianState} +) + Base.summary(io, H); println(io) + println(io, "result:") + Base.show(io, mime, H.result) + println(io, "\noutput state:") + Base.show(io, mime, H.state) +end + +""" + heterodyne(state::GaussianState, indices::Vector) -> Heterodyne + heterodyne(state::GaussianState, index::Int) -> Heterodyne + +Compute the heterodyne measurement of the subsystem of a Gaussian state `state` +indicated by `indices` and return a `Heterodyne` object. +The `result` (phase-space quadratures) and mapped state `output` can be obtained +from the Heterodyne object `M` via `M.result` and `M.output`. +Iterating the decomposition produces the components `result` and `output`. + +Heterodyne measurement simultaneously measures both quadratures of the mode(s) with unit gain. +The measured modes are replaced with vacuum states after the heterodyne measurement. + +# Keyword arguments +- `rng::AbstractRNG = Random.default_rng()`: Random number generator that determines a random projection. + +# Examples +``` +julia> st = squeezedstate(QuadBlockBasis(3), 1.0, pi/4); + +julia> M = heterodyne(st, [1, 3]) +Heterodyne{Matrix{ComplexF64}, GaussianState{QuadBlockBasis{Int64}, Vector{Float64}, Matrix{Float64}}} +result: +2×2 Matrix{ComplexF64}: + 0.123456+0.234567im 0.345678+0.456789im + 0.567890+0.678901im 0.789012+0.890123im +output state: +GaussianState for 3 modes. + symplectic basis: QuadBlockBasis +mean: 6-element Vector{Float64}: + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 +covariance: 6×6 Matrix{Float64}: + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 0.0 0.0 1.0 + +julia> result, state = M; # destructuring via iteration + +julia> result == M.result && state == M.state +true +``` +""" +function heterodyne( + state::GaussianState{<:QuadPairBasis,Tm,Tc}, + indices::R; + rng::AbstractRNG = Random.default_rng() +) where {Tm,Tc,R} + basis = state.basis + nmodes = basis.nmodes + indlength = length(indices) + indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR)) + # perform conditional mapping of Gaussian quantum state + result′, a, A = _heterodyne_filter(rng, state, indices) + mean′ = zeros(eltype(Tm), 2*nmodes) + covar′ = Matrix{eltype(Tc)}((state.ħ/2)*I, 2*nmodes, 2*nmodes) + # fill in measured modes with vacuum states + notindices = setdiff(1:nmodes, indices) + @inbounds for i in eachindex(notindices) + idx = notindices[i] + copyto!(@view(mean′[2idx-1:2idx]), @view(a[2i-1:2i])) + @inbounds for j in i:length(notindices) + otheridx = notindices[j] + covar′[2*idx-1, 2*otheridx-1] = A[2*i-1, 2*j-1] + covar′[2*idx-1, 2*otheridx] = A[2*i-1, 2*j] + covar′[2*idx, 2*otheridx-1] = A[2*i, 2*j-1] + covar′[2*idx, 2*otheridx] = A[2*i, 2*j] + covar′[2*otheridx-1, 2*idx-1] = A[2*j-1, 2*i-1] + covar′[2*otheridx-1, 2*idx] = A[2*j-1, 2*i] + covar′[2*otheridx, 2*idx-1] = A[2*j, 2*i-1] + covar′[2*otheridx, 2*idx] = A[2*j, 2*i] + end + end + # promote output array type to ensure it matches the input array type + mean′′ = Gabs._promote_output_vector(Tm, mean′, 2*nmodes) + covar′′ = Gabs._promote_output_matrix(Tc, covar′, 2*nmodes) + state′ = GaussianState(basis, mean′′, covar′′, ħ = state.ħ) + return Heterodyne(result′, state′) +end +function heterodyne( + state::GaussianState{<:QuadBlockBasis,Tm,Tc}, + indices::R; + rng::AbstractRNG = Random.default_rng() +) where {Tm,Tc,R} + basis = state.basis + nmodes = basis.nmodes + indlength = length(indices) + indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR)) + # perform conditional mapping of Gaussian quantum state + result′, a, A = _heterodyne_filter(rng, state, indices) + mean′ = zeros(eltype(Tm), 2*nmodes) + covar′ = Matrix{eltype(Tc)}((state.ħ/2)*I, 2*nmodes, 2*nmodes) + nmodes′ = nmodes - length(indices) + # fill in measured modes with vacuum states + notindices = setdiff(1:nmodes, indices) + @inbounds for i in eachindex(notindices) + idx = notindices[i] + mean′[idx] = a[i] + mean′[idx+nmodes] = a[i+nmodes′] + @inbounds for j in i:length(notindices) + otheridx = notindices[j] + covar′[idx,otheridx] = A[i,j] + covar′[otheridx,idx] = A[j,i] + covar′[idx+nmodes,otheridx] = A[i+nmodes′,j] + covar′[idx,otheridx+nmodes] = A[i,j+nmodes′] + covar′[otheridx,idx+nmodes] = A[j,i+nmodes′] + covar′[otheridx+nmodes,idx] = A[j+nmodes′,i] + covar′[idx+nmodes,otheridx+nmodes] = A[i+nmodes′,j+nmodes′] + covar′[otheridx+nmodes,idx+nmodes] = A[j+nmodes′,i+nmodes′] + end + end + # promote output array type to ensure it matches the input array type + mean′′ = Gabs._promote_output_vector(Tm, mean′, 2*nmodes) + covar′′ = Gabs._promote_output_matrix(Tc, covar′, 2*nmodes) + state′ = GaussianState(basis, mean′′, covar′′, ħ = state.ħ) + return Heterodyne(result′, state′) +end +heterodyne(rng::AbstractRNG, state::GaussianState{<:QuadPairBasis,Tm,Tc}, indices::R) where {Tm,Tc,R} = heterodyne(state, indices; rng) +heterodyne(rng::AbstractRNG, state::GaussianState{<:QuadBlockBasis,Tm,Tc}, indices::R) where {Tm,Tc,R} = heterodyne(state, indices; rng) + +""" + Base.rand(::Type{Heterodyne}, state::GaussianState, indices; shots=1, rng=Random.default_rng()) + +Generate random heterodyne measurement outcomes for a Gaussian state. +Returns a matrix with shape (nmodes_measured, shots) where each column is a complex outcome +representing the two-quadrature measurement. +""" +function Base.rand( + rng::AbstractRNG, + ::Type{Heterodyne}, + state::GaussianState{<:QuadPairBasis,Tm,Tc}, + indices::R; + shots::Int = 1 +) where {Tm,Tc,R} + basis = state.basis + indlength = length(indices) + indlength < basis.nmodes || throw(ArgumentError(Gabs.INDEX_ERROR)) + mean, covar = state.mean, state.covar + # write mean and covariance matrix of measured modes to vector `b` and matrix `B`, respectively + b, B = zeros(2*indlength), zeros(2*indlength, 2*indlength) + @inbounds for i in eachindex(indices) + idx = indices[i] + b[2i-1:2i] .= @view(mean[2idx-1:2idx]) + @inbounds for j in eachindex(indices) + otheridx = indices[j] + if idx == otheridx + B[2i-1:2i, 2i-1:2i] .= @view(covar[2idx-1:2idx, 2idx-1:2idx]) + else + B[2i-1:2i, 2j-1:2j] .= @view(covar[2idx-1:2idx, 2otheridx-1:2otheridx]) + B[2j-1:2j, 2i-1:2i] .= @view(covar[2otheridx-1:2otheridx, 2idx-1:2idx]) + end + end + end + # add vacuum noise (unit gain heterodyne measurement) + @inbounds for i in Base.OneTo(indlength) + B[2i-1, 2i-1] += 1.0 + B[2i, 2i] += 1.0 + end + # sample from probability distribution by taking the displaced + # Cholesky decomposition of the covariance matrix + symB = Symmetric(B) + L = cholesky(symB).L + buf = zeros(2*indlength) + results = zeros(ComplexF64, indlength, shots) + @inbounds for shot in Base.OneTo(shots) + mul!(buf, L, randn!(rng, buf)) + buf .+= b + @inbounds for i in Base.OneTo(indlength) + results[i, shot] = complex(buf[2i-1], buf[2i]) + end + end + return results +end +function Base.rand( + rng::AbstractRNG, + ::Type{Heterodyne}, + state::GaussianState{<:QuadBlockBasis,Tm,Tc}, + indices::R; + shots::Int = 1 +) where {Tm,Tc,R} + basis = state.basis + nmodes = basis.nmodes + indlength = length(indices) + indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR)) + nmodes′ = nmodes - indlength + mean, covar = state.mean, state.covar + # write mean and covariance matrix of measured modes to vector `b` and matrix `B`, respectively + b, B = zeros(2*indlength), zeros(2*indlength, 2*indlength) + @inbounds for i in eachindex(indices) + idx = indices[i] + b[i] = mean[idx] + b[i+indlength] = mean[idx+nmodes] + @inbounds for j in eachindex(indices) + otheridx = indices[j] + if idx == otheridx + B[i, i] = covar[idx, idx] + B[i+indlength, i] = covar[idx+nmodes, idx] + B[i, i+indlength] = covar[idx, idx+nmodes] + B[i+indlength, i+indlength] = covar[idx+nmodes, idx+nmodes] + else + B[i, j] = covar[idx, otheridx] + B[i+indlength, j] = covar[idx+nmodes, otheridx] + B[i, j+indlength] = covar[idx, otheridx+nmodes] + B[i+indlength, j+indlength] = covar[idx+nmodes, otheridx+nmodes] + + B[j, i] = covar[otheridx, idx] + B[j+indlength, i] = covar[otheridx+nmodes, idx] + B[j, i+indlength] = covar[otheridx, idx+nmodes] + B[j+indlength, i+indlength] = covar[otheridx+nmodes, idx+nmodes] + end + end + end + # add vacuum noise (unit gain heterodyne measurement) + @inbounds for i in Base.OneTo(indlength) + B[i, i] += 1.0 + B[i+indlength, i+indlength] += 1.0 + end + # sample from probability distribution by taking the displaced + # Cholesky decomposition of the covariance matrix + symB = Symmetric(B) + L = cholesky(symB).L + buf = zeros(2*indlength) + results = zeros(ComplexF64, indlength, shots) + @inbounds for shot in Base.OneTo(shots) + mul!(buf, L, randn!(rng, buf)) + buf .+= b + @inbounds for i in Base.OneTo(indlength) + results[i, shot] = complex(buf[i], buf[i+indlength]) + end + end + return results +end +function Base.rand( + ::Type{Heterodyne}, + state::GaussianState, + indices::R; + shots::Int = 1, + rng::AbstractRNG = Random.default_rng() +) where {R} + return Base.rand(rng, Heterodyne, state, indices; shots) +end + +function _heterodyne_filter( + rng::AbstractRNG, + state::GaussianState{<:QuadPairBasis,Tm,Tc}, + indices::R +) where {Tm,Tc,R} + basis = state.basis + indlength = length(indices) + nmodes′ = basis.nmodes - indlength + a, b, A, B, C = _part_state(state, indices) + # add unit-gain vacuum noise for heterodyne measurement + @inbounds for i in Base.OneTo(indlength) + B[2i-1,2i-1] += 1.0 + B[2i,2i] += 1.0 + end + # sample from probability distribution by taking the displaced + # Cholesky decomposition of the covariance matrix + symB = Symmetric(B) + L = cholesky(symB).L + resultmean = L * randn(rng, 2*indlength) + b + meandiff = resultmean - b + # conditional mapping + buf = C * inv(symB) + a .+= buf * meandiff + A .-= buf * C' + # convert to complex representation and promote output array type + result_complex = zeros(ComplexF64, indlength) + @inbounds for i in Base.OneTo(indlength) + result_complex[i] = complex(resultmean[2i-1], resultmean[2i]) + end + return result_complex, a, A +end +function _heterodyne_filter( + rng::AbstractRNG, + state::GaussianState{<:QuadBlockBasis,Tm,Tc}, + indices::R +) where {Tm,Tc,R} + basis = state.basis + indlength = length(indices) + nmodes′ = basis.nmodes - indlength + a, b, A, B, C = _part_state(state, indices) + # add unit-gain vacuum noise for heterodyne measurement + @inbounds for i in Base.OneTo(indlength) + B[i,i] += 1.0 + B[i+indlength,i+indlength] += 1.0 + end + # sample from probability distribution by taking the displaced + # Cholesky decomposition of the covariance matrix + symB = Symmetric(B) + L = cholesky(symB).L + resultmean = L * randn(rng, 2*indlength) + b + meandiff = resultmean - b + # conditional mapping + buf = C * inv(symB) + a .+= buf * meandiff + A .-= buf * C' + # convert to complex representation + result_complex = zeros(ComplexF64, indlength) + @inbounds for i in Base.OneTo(indlength) + result_complex[i] = complex(resultmean[i], resultmean[i+indlength]) + end + return result_complex, a, A +end diff --git a/test/test_measurements.jl b/test/test_measurements.jl index eb51569..a728eeb 100644 --- a/test/test_measurements.jl +++ b/test/test_measurements.jl @@ -174,4 +174,102 @@ @test_throws ArgumentError rand(Homodyne, rs_qblock, collect(1:5), [π/2]) end end + + @testset "heterodyne" begin + + @testset "heterodyne" begin + qpairbasis, qblockbasis = QuadPairBasis(1), QuadBlockBasis(1) + + for basis in (qpairbasis, qblockbasis) + vac = vacuumstate(basis) + @test_throws ArgumentError heterodyne(vac, [1, 2]) + end + + # simple test case: measuring 1 mode of 2-mode state + for basis in (QuadPairBasis(2), QuadBlockBasis(2)) + st = vacuumstate(basis) ⊗ vacuumstate(basis) + M = heterodyne(st, [1]) + @test M isa Heterodyne + @test M.result isa Vector{ComplexF64} + @test M.state isa GaussianState + + result, state = M + @test result == M.result + @test state == M.state + + # check if measured mode is replaced with vacuum + @test isapprox(state.mean[1:2], zeros(2), atol=1e-12) + @test isapprox(state.covar[1:2, 1:2], Matrix{Float64}(I,2,2), atol=1e-12) + end + + st = squeezedstate(QuadPairBasis(4), 0.5, π/2) + st_block = changebasis(QuadBlockBasis, st) + indices = [2, 4] + @test size(rand(Heterodyne, st, indices; shots=10)) == (2, 10) + @test size(rand(Heterodyne, st_block, indices; shots=7)) == (2, 7) + + indices = [1, 2] + rs_qpair = randstate(QuadPairBasis(4)) + rs_qblock = changebasis(QuadBlockBasis, rs_qpair) + M_qpair = heterodyne(rs_qpair, indices) + M_qblock = heterodyne(rs_qblock, indices) + + # extract analytical conditional update + xA, xB, VA, VB, VAB = Gabs._part_state(rs_qpair, indices) + B = copy(VB) + # add unit-gain vacuum noise for heterodyne (no squeezing) + for i in 1:2 + B[2i-1,2i-1] += 1.0 + B[2i,2i] += 1.0 + end + L = cholesky(Symmetric(B)).L + resultmean = L * randn(4) + xB + xA′ = xA .+ VAB * (inv(B) * (resultmean - xB)) + VA′ = VA .- VAB * (inv(B) * transpose(VAB)) + + nm = 4 + out_mean = zeros(2nm) + notindices = setdiff(1:nm, indices) + for i in eachindex(notindices) + idx = notindices[i] + out_mean[2idx-1:2idx] .= @view(xA′[2i-1:2i]) + end + out_covar = Matrix{Float64}(I, 2nm, 2nm) + for i in eachindex(notindices), j in i:length(notindices) + idx, jdx = notindices[i], notindices[j] + out_covar[2idx-1:2idx, 2jdx-1:2jdx] .= @view(VA′[2i-1:2i, 2j-1:2j]) + out_covar[2jdx-1:2jdx, 2idx-1:2idx] .= @view(VA′[2j-1:2j, 2i-1:2i]) + end + expected_state = GaussianState(QuadPairBasis(nm), out_mean, out_covar) + @test isapprox(M_qpair.state.covar, expected_state.covar, atol=1e-8) + @test isapprox(M_qblock.state.covar, changebasis(QuadBlockBasis, expected_state).covar, atol=1e-8) + + sstatic = vacuumstate(SVector{2}, SMatrix{2,2}, QuadPairBasis(1)) + statestatic = sstatic ⊗ sstatic ⊗ sstatic ⊗ sstatic + hstatic = heterodyne(statestatic, [2]) + @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 = heterodyne(MersenneTwister(seed), base_state, [1]) + h2 = heterodyne(MersenneTwister(seed), base_state, [1]) + @test h1.result == h2.result + @test h1.state == h2.state + + hb1 = heterodyne(MersenneTwister(seed), base_state_block, [1]) + hb2 = heterodyne(MersenneTwister(seed), base_state_block, [1]) + @test hb1.result == hb2.result + @test hb1.state == hb2.state + + samples_kw = rand(Heterodyne, base_state, [1]; shots = 3, rng = MersenneTwister(seed)) + samples_pos = rand(MersenneTwister(seed), Heterodyne, base_state, [1]; shots = 3) + @test samples_kw == samples_pos + + @test_throws ArgumentError rand(Heterodyne, rs_qpair, collect(1:5)) + @test_throws ArgumentError rand(Heterodyne, rs_qblock, collect(1:5)) + end + end end \ No newline at end of file