From 3a637602fbfb5e7777e4b053ea4d6aeab7fdaca0 Mon Sep 17 00:00:00 2001
From: Akira Kyle <akira@akirakyle.com>
Date: Thu, 21 Sep 2023 15:16:01 -0600
Subject: [PATCH] Add haar random states and unitaries (#143)

* Add randstate_haar and randunitary_haar

* Silence Aqua ambiguities test for StatsBase pulled in by RandomMatrices
---
 Project.toml                   |  5 ++++-
 src/QuantumOpticsBase.jl       |  1 +
 src/state_definitions.jl       | 15 +++++++++++++++
 test/test_aqua.jl              |  4 +++-
 test/test_state_definitions.jl | 10 ++++++++++
 5 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index ce682fa9..e5a3a61f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,7 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RandomMatrices = "2576dda1-a324-5b11-aa66-c48ed7e3c618"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
 UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
@@ -24,6 +25,7 @@ FastGaussQuadrature = "0.5"
 FillArrays = "0.13, 1"
 LRUCache = "1"
 QuantumInterface = "0.3.0"
+RandomMatrices = "0.5"
 Strided = "1, 2"
 UnsafeArrays = "1"
 julia = "1.6"
@@ -39,10 +41,11 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
 StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 
 [targets]
-test = ["LinearAlgebra", "SparseArrays", "Random", "Test", "Aqua", "JET", "Adapt", "Dates", "FFTW", "LRUCache", "Strided", "StridedViews", "UnsafeArrays", "FillArrays"]
+test = ["LinearAlgebra", "SparseArrays", "Random", "Test", "Aqua", "JET", "Adapt", "Dates", "FFTW", "LRUCache", "Strided", "StridedViews", "UnsafeArrays", "FillArrays", "StatsBase"]
diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl
index 14608f2e..eea0f2e3 100644
--- a/src/QuantumOpticsBase.jl
+++ b/src/QuantumOpticsBase.jl
@@ -41,6 +41,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
                 displace, displace_analytical, displace_analytical!,
                 squeeze,
         randstate, randoperator, thermalstate, coherentthermalstate, phase_average, passive_state,
+        randstate_haar, randunitary_haar,
         #spin
                 SpinBasis, sigmax, sigmay, sigmaz, sigmap, sigmam, spinup, spindown,
         #subspace
diff --git a/src/state_definitions.jl b/src/state_definitions.jl
index 34086a41..657e186e 100644
--- a/src/state_definitions.jl
+++ b/src/state_definitions.jl
@@ -1,3 +1,4 @@
+using RandomMatrices
 """
     randstate([T=ComplexF64,] basis)
 
@@ -10,6 +11,13 @@ function randstate(::Type{T}, b::Basis) where T
 end
 randstate(b) = randstate(ComplexF64, b)
 
+"""
+    randstate_haar(b::Basis)
+
+Returns a Haar random pure state of dimension `length(b)` by applying a Haar random unitary to a fixed pure state.
+"""
+randstate_haar(b::Basis) = Ket(b, rand(Haar(2), length(b), 2)[:,1])
+
 """
     randoperator([T=ComplexF64,] b1[, b2])
 
@@ -20,6 +28,13 @@ randoperator(b1::Basis, b2::Basis) = randoperator(ComplexF64, b1, b2)
 randoperator(::Type{T}, b::Basis) where T = randoperator(T, b, b)
 randoperator(b) = randoperator(ComplexF64, b)
 
+"""
+    randunitary_haar(b::Basis)
+
+Returns a Haar random unitary matrix of dimension `length(b)`.
+"""
+randunitary_haar(b::Basis) = DenseOperator(b, b, rand(Haar(2), length(b), 2))
+
 """
     thermalstate(H,T)
 
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index b3ecf9ad..bd5289af 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -2,10 +2,12 @@ using Test
 using QuantumOpticsBase
 using Aqua
 using FillArrays
+using StatsBase
 
 @testset "aqua" begin
     Aqua.test_all(QuantumOpticsBase;
-                  ambiguities=(exclude=[FillArrays.reshape],),  # Due to https://github.com/JuliaArrays/FillArrays.jl/issues/105#issuecomment-1518406018
+                  ambiguities=(exclude=[FillArrays.reshape, # Due to https://github.com/JuliaArrays/FillArrays.jl/issues/105#issuecomment-1518406018
+                                        StatsBase.TestStat, StatsBase.:(==) , StatsBase.sort!],),  # Due to https://github.com/JuliaStats/StatsBase.jl/issues/861
                   piracy=(broken=true,)
                   )
     # manual piracy check to exclude identityoperator
diff --git a/test/test_state_definitions.jl b/test/test_state_definitions.jl
index 678d97b5..6cd3eeca 100644
--- a/test/test_state_definitions.jl
+++ b/test/test_state_definitions.jl
@@ -1,5 +1,6 @@
 using Test
 using QuantumOpticsBase
+using Random
 
 @testset "state_definitions" begin
 
@@ -37,4 +38,13 @@ end
 @test isapprox(expect(destroy(b),rpas),0, atol=1e-14)
 @test isapprox(entropy_vn(rpas),S, atol=1e-14)
 
+Random.seed!(0)
+
+for n in 1:5:50
+    b = FockBasis(n)
+    @test isapprox(norm(randstate_haar(b)), 1, atol=1e-8)
+    U = randunitary_haar(b)
+    @test isapprox(dagger(U)*U, identityoperator(b), atol=1e-8)
+end
+
 end # testset