Add finite squeezing keyword argument in the Homodyne interface

CHANGELOG.md

@@ -1,5 +1,11 @@
 # News
 
+## v1.3.4 - 2026-01-05
+
+- Add optional rng arguments to random methods.
+- Implement embed functions for Gaussian unitaries, channels, and linear combinations.
+- Add finite squeezing keyword argument in the `homodyne`, `Homodyne` interface.
+
 ## v1.3.3 - 2025-12-13
 
 - Replace SymplecticFactorizations.jl as a dependency with SymplecticMatrices.jl.

Project.toml

@@ -1,7 +1,7 @@
 name = "Gabs"
 uuid = "0eb812ee-a11f-4f5e-b8d4-bb8a44f06f50"
 authors = ["Andrew Kille"]
-version = "1.3.3"
+version = "1.3.4"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/homodyne.jl

@@ -35,6 +35,10 @@ Iterating the decomposition produces the components `result` and `output`.
 
 Note the measured modes are replaced with vacuum states after the homodyne measurement.
 
+# Keyword arguments
+- `rng::AbstractRNG = Random.default_rng()`: Random number generator that determines a random projection.
+- `squeeze::Real = 1e-12`: Finite squeezing parameter.
+
 # Examples
 ```
 julia> st = squeezedstate(QuadBlockBasis(3), 1.0, pi/4);
@@ -75,15 +79,16 @@ function homodyne(
     state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
     indices::R, 
     angles::G;
-    rng::AbstractRNG = Random.default_rng()
+    rng::AbstractRNG = Random.default_rng(),
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
     basis = state.basis
     nmodes = basis.nmodes
     indlength = length(indices)
     indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR))
     indlength == length(angles) || throw(ArgumentError(Gabs.GENERALDYNE_ERROR))
     # perform conditional mapping of Gaussian quantum state
-    result′, a, A = _homodyne_filter(rng, state, indices, angles)
+    result′, a, A = _homodyne_filter(rng, state, indices, angles; squeeze)
     mean′ = zeros(eltype(Tm), 2*nmodes)
     covar′ = Matrix{eltype(Tc)}((state.ħ/2)*I, 2*nmodes, 2*nmodes)
     # fill in measured modes with vacuum states 
@@ -113,15 +118,16 @@ function homodyne(
     state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
     indices::R, 
     angles::G;
-    rng::AbstractRNG = Random.default_rng()
+    rng::AbstractRNG = Random.default_rng(),
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
     basis = state.basis
     nmodes = basis.nmodes
     indlength = length(indices)
     indlength < nmodes || throw(ArgumentError(Gabs.INDEX_ERROR))
     indlength == length(angles) || throw(ArgumentError(Gabs.GENERALDYNE_ERROR))
     # perform conditional mapping of Gaussian quantum state
-    result′, a, A = _homodyne_filter(rng, state, indices, angles)
+    result′, a, A = _homodyne_filter(rng, state, indices, angles; squeeze)
     mean′ = zeros(eltype(Tm), 2*nmodes)
     covar′ = Matrix{eltype(Tc)}((state.ħ/2)*I, 2*nmodes, 2*nmodes)
     nmodes′ = nmodes - length(indices)
@@ -149,16 +155,16 @@ 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)
+homodyne(rng::AbstractRNG, state::GaussianState{<:QuadPairBasis,Tm,Tc}, indices::R, angles::G; squeeze::Real = 1e-12) where {Tm,Tc,R,G} = homodyne(state, indices, angles; rng, squeeze)
+homodyne(rng::AbstractRNG, state::GaussianState{<:QuadBlockBasis,Tm,Tc}, indices::R, angles::G; squeeze::Real = 1e-12) where {Tm,Tc,R,G} = homodyne(state, indices, angles; rng, squeeze)
 
 """
-    rand([rng::AbstractRNG], ::Type{Homodyne}, state::GaussianState, indices::Vector, angles::Vector; shots = 1) -> Array
+    rand([rng::AbstractRNG], ::Type{Homodyne}, state::GaussianState, indices::Vector, angles::Vector; shots = 1, squeeze = 1e-12) -> 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.
 The number of shots is given by `shots`, which determines how many repeated and random homodyne measurements
-are performed on the quantum system.
+are performed on the quantum system. The `squeeze` parameter determines the finite squeezing performed during the projection.
 
 The output is an `2*length(indices) × shots` array, which contains the measured position and momentum modes columnwise
 for each measurement, the ordering basis of the input Gaussian state `state`.
@@ -181,17 +187,19 @@ function Base.rand(
     indices::R, 
     angles::G; 
     shots::Int = 1,
-    rng::AbstractRNG = Random.default_rng()
+    rng::AbstractRNG = Random.default_rng(),
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
-    return Base.rand(rng, Homodyne, state, indices, angles; shots = shots)
+    return Base.rand(rng, Homodyne, state, indices, angles; shots, squeeze)
 end
 function Base.rand(
     rng::AbstractRNG,
     ::Type{Homodyne}, 
     state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
     indices::R, 
     angles::G; 
-    shots::Int = 1
+    shots::Int = 1,
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
     basis = state.basis
     indlength = length(indices)
@@ -217,12 +225,11 @@ function Base.rand(
     # infinite squeezing along axis defined by `angles`
     @inbounds for i in Base.OneTo(indlength)
         θ = angles[i]
-        sq = 1e-12
         ct, st = cos(θ), sin(θ)
-        B[i,i] += ct^2 * sq + st^2 / sq
-        B[i,i+indlength] += ct * st * (sq - 1 / sq)
-        B[i+indlength,i] += ct * st * (sq - 1 / sq)
-        B[i+indlength,i+indlength] += st^2 * sq + ct^2 / sq
+        B[i,i] += ct^2 * squeeze + st^2 / squeeze
+        B[i,i+indlength] += ct * st * (squeeze - 1 / squeeze)
+        B[i+indlength,i] += ct * st * (squeeze - 1 / squeeze)
+        B[i+indlength,i+indlength] += st^2 * squeeze + ct^2 / squeeze
     end
     # sample from probability distribution by taking the displaced 
     # Cholesky decomposition of the covariance matrix
@@ -242,17 +249,19 @@ function Base.rand(
     indices::R, 
     angles::G;
     shots::Int = 1,
-    rng::AbstractRNG = Random.default_rng()
+    rng::AbstractRNG = Random.default_rng(),
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
-    return Base.rand(rng, Homodyne, state, indices, angles; shots = shots)
+    return Base.rand(rng, Homodyne, state, indices, angles; shots, squeeze)
 end
 function Base.rand(
     rng::AbstractRNG,
     ::Type{Homodyne}, 
     state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
     indices::R, 
     angles::G;
-    shots::Int = 1
+    shots::Int = 1,
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
     basis = state.basis
     nmodes = basis.nmodes
@@ -290,12 +299,11 @@ function Base.rand(
     # infinite squeezing along axis defined by `angles`
     @inbounds for i in Base.OneTo(indlength)
         θ = angles[i]
-        sq = 1e-12
         ct, st = cos(θ), sin(θ)
-        B[i,i] += ct^2 * sq + st^2 / sq
-        B[i,i+indlength] += ct * st * (sq - 1 / sq)
-        B[i+indlength,i] += ct * st * (sq - 1 / sq)
-        B[i+indlength,i+indlength] += st^2 * sq + ct^2 / sq
+        B[i,i] += ct^2 * squeeze + st^2 / squeeze
+        B[i,i+indlength] += ct * st * (squeeze - 1 / squeeze)
+        B[i+indlength,i] += ct * st * (squeeze - 1 / squeeze)
+        B[i+indlength,i+indlength] += st^2 * squeeze + ct^2 / squeeze
     end
     # sample from probability distribution by taking the displaced 
     # Cholesky decomposition of the covariance matrix
@@ -314,7 +322,8 @@ function _homodyne_filter(
     rng::AbstractRNG,
     state::GaussianState{<:QuadPairBasis,Tm,Tc}, 
     indices::R, 
-    angles::G
+    angles::G;
+    squeeze::Real = 1e-12
 ) where {Tm,Tc,R,G}
     basis = state.basis
     indlength = length(indices)
@@ -323,12 +332,11 @@ function _homodyne_filter(
     # infinite squeezing along axis defined by `angles`
     @inbounds for i in Base.OneTo(indlength)
         θ = angles[i]
-        sq = 1e-12
         ct, st = cos(θ), sin(θ)
-        B[2i-1,2i-1] += ct^2 * sq + st^2 / sq
-        B[2i-1,2i] += ct * st * (sq - 1 / sq)
-        B[2i,2i-1] += ct * st * (sq - 1 / sq)
-        B[2i,2i] += st^2 * sq + ct^2 / sq
+        B[2i-1,2i-1] += ct^2 * squeeze + st^2 / squeeze
+        B[2i-1,2i] += ct * st * (squeeze - 1 / squeeze)
+        B[2i,2i-1] += ct * st * (squeeze - 1 / squeeze)
+        B[2i,2i] += st^2 * squeeze + ct^2 / squeeze
     end
     # sample from probability distribution by taking the displaced 
     # Cholesky decomposition of the covariance matrix
@@ -348,7 +356,8 @@ function _homodyne_filter(
     rng::AbstractRNG,
     state::GaussianState{<:QuadBlockBasis,Tm,Tc}, 
     indices::R, 
-    angles::G
+    angles::G;
+    squeeze = 1e-12
 ) where {Tm,Tc,R,G}
     basis = state.basis
     indlength = length(indices)
@@ -357,12 +366,11 @@ function _homodyne_filter(
     # infinite squeezing along axis defined by `angles`
     @inbounds for i in Base.OneTo(indlength)
         θ = angles[i]
-        sq = 1e-12
         ct, st = cos(θ), sin(θ)
-        B[i,i] += ct^2 * sq + st^2 / sq
-        B[i,i+indlength] += ct * st * (sq - 1 / sq)
-        B[i+indlength,i] += ct * st * (sq - 1 / sq)
-        B[i+indlength,i+indlength] += st^2 * sq + ct^2 / sq
+        B[i,i] += ct^2 * squeeze + st^2 / squeeze
+        B[i,i+indlength] += ct * st * (squeeze - 1 / squeeze)
+        B[i+indlength,i] += ct * st * (squeeze - 1 / squeeze)
+        B[i+indlength,i+indlength] += st^2 * squeeze + ct^2 / squeeze
     end
     # sample from probability distribution by taking the displaced 
     # Cholesky decomposition of the covariance matrix

test/test_measurements.jl

@@ -72,17 +72,18 @@
 
         @testset "homodyne" begin
             qpairbasis, qblockbasis = QuadPairBasis(1), QuadBlockBasis(1)
-        
+
             for basis in (qpairbasis, qblockbasis)
                 vac = vacuumstate(basis)
                 @test_throws ArgumentError homodyne(vac, [1, 2], [0.0])
                 @test_throws ArgumentError homodyne(vac, [1], [0.0, π/2])
             end
 
+            squeeze = 1e-14
             # simple test case: measuring 1 mode of 2-mode state
             for basis in (QuadPairBasis(2), QuadBlockBasis(2))
                 st = vacuumstate(basis) ⊗ vacuumstate(basis)
-                M = homodyne(st, [1], [0.0])
+                M = homodyne(st, [1], [0.0]; squeeze)
                 @test M isa Homodyne
                 @test M.result isa Vector
                 @test M.state isa GaussianState
@@ -100,8 +101,8 @@
             st_block = changebasis(QuadBlockBasis, st)
             indices = [2, 4]
             angles = [0.0, π/2]
-            @test size(rand(Homodyne, st, indices, angles, shots=10)) == (4, 10)
-            @test size(rand(Homodyne, st_block, indices, angles, shots=7)) == (4, 7)
+            @test size(rand(Homodyne, st, indices, angles; shots=10, squeeze)) == (4, 10)
+            @test size(rand(Homodyne, st_block, indices, angles; shots=7, squeeze)) == (4, 7)
 
             indices = [1, 2]
             rs_qpair = randstate(QuadPairBasis(4))
@@ -155,18 +156,18 @@
             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])
+            h1 = homodyne(MersenneTwister(seed), base_state, [1], [0.0]; squeeze)
+            h2 = homodyne(MersenneTwister(seed), base_state, [1], [0.0]; squeeze)
             @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])
+            hb1 = homodyne(MersenneTwister(seed), base_state_block, [1], [0.0]; squeeze)
+            hb2 = homodyne(MersenneTwister(seed), base_state_block, [1], [0.0]; squeeze)
             @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)
+            samples_kw = rand(Homodyne, base_state, [1], [π/2]; shots = 3, rng = MersenneTwister(seed), squeeze)
+            samples_pos = rand(MersenneTwister(seed), Homodyne, base_state, [1], [π/2]; shots = 3, squeeze)
             @test samples_kw == samples_pos
 
             @test_throws ArgumentError rand(Homodyne, rs_qpair, collect(1:5), [0.0])