From c310e27e28965d7721ea8d769f317cf555dcb0ed Mon Sep 17 00:00:00 2001
From: Andrew Kille <68079167+apkille@users.noreply.github.com>
Date: Wed, 25 Dec 2024 03:14:38 -0500
Subject: [PATCH] Add function for computing symplectic spectrum of a Gaussian
 state (#22)

* symplectic spectrum

* update project.toml and changelog.md

* use det in test_states.jl

* adjust tolerances for det test
---
 CHANGELOG.md        |  6 +++++-
 Project.toml        |  2 +-
 src/Gabs.jl         |  2 +-
 src/states.jl       | 13 +++++++++++++
 test/test_states.jl | 19 +++++++++++++++++++
 5 files changed, 39 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4411630..7794b4c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,10 @@
 # News
 
-## v1.2.1 - dev
+## v1.2.2 - 2024-12-25
+
+- Add `sympspectrum` function to compute symplectic spectrum of `GaussianState` object.
+  
+## v1.2.1 - 2024-12-17
 
 - Add `issymplectic` and `isgaussian` checks.
 - **(fix)** Generate correct Williamson matrix in `randstate` generation for `QuadBlockBasis`.
diff --git a/Project.toml b/Project.toml
index 18a3e8f..50d34c4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Gabs"
 uuid = "0eb812ee-a11f-4f5e-b8d4-bb8a44f06f50"
 authors = ["Andrew Kille"]
-version = "1.2.1"
+version = "1.2.2"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
diff --git a/src/Gabs.jl b/src/Gabs.jl
index 6e8a2d6..5205e7b 100644
--- a/src/Gabs.jl
+++ b/src/Gabs.jl
@@ -24,7 +24,7 @@ export
     # wigner functions
     wigner, wignerchar,
     # symplectic form and checks
-    symplecticform, issymplectic, isgaussian,
+    symplecticform, issymplectic, isgaussian, sympspectrum,
     # metrics
     purity
 
diff --git a/src/states.jl b/src/states.jl
index 51879f6..60dfc48 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -712,3 +712,16 @@ function _changebasis(state::GaussianState{B1,M,V}, ::Type{B2}) where {B1<:QuadP
     covar = T * state.covar * transpose(T)
     return GaussianState(B2(nmodes), mean, covar)
 end
+
+"""
+    sympspectrum(state::GaussianState)
+
+Compute the symplectic spectrum of a Gaussian state.
+"""
+function sympspectrum(state::GaussianState)
+    basis = state.basis
+    form = symplecticform(basis)
+    M = form * state.covar
+    spectrum = filter(x -> x > 0, imag.(eigvals(M)))
+    return spectrum
+end
diff --git a/test/test_states.jl b/test/test_states.jl
index 467df4a..3d4006f 100644
--- a/test/test_states.jl
+++ b/test/test_states.jl
@@ -2,6 +2,7 @@
     import Gabs: _changebasis
     using Gabs
     using StaticArrays
+    using LinearAlgebra: det
 
     nmodes = rand(1:5)
     qpairbasis = QuadPairBasis(nmodes)
@@ -114,4 +115,22 @@
         @test ptrace(SVector{2}, SMatrix{2,2}, state_qpair, 1) isa GaussianState
         @test ptrace(SVector{4}, SMatrix{4,4}, state_qpair, [1, 3]) isa GaussianState
     end
+
+    @testset "symplectic spectrum" begin
+        nmodes = rand(1:20)
+        qpairbasis = QuadPairBasis(nmodes)
+        qblockbasis = QuadBlockBasis(nmodes)
+
+        s_qpair = randstate(qpairbasis)
+        s_qblock = randstate(qblockbasis)
+
+        spec_qpair = sympspectrum(s_qpair)
+        spec_qblock = sympspectrum(s_qblock)
+
+        @test all(i > 1 || isapprox(i, 1, atol=1e-5) for i in spec_qpair)
+        @test all(i > 1 || isapprox(i, 1, atol=1e-5) for i in spec_qblock)
+
+        @test isapprox(det(s_qpair.covar), prod(abs2, spec_qpair), atol=1e-3)
+        @test isapprox(det(s_qblock.covar), prod(abs2, spec_qblock), atol=1e-3)
+    end
 end
\ No newline at end of file