From 37af253b749b4f38bfbbccc5edc418b02fe84d5f Mon Sep 17 00:00:00 2001 From: Karol Gietka Date: Thu, 17 Aug 2023 11:10:11 +0200 Subject: [PATCH 01/14] proper squeeze update --- src/spin.jl | 18 ++++++++++++++++++ test/test_spin.jl | 17 +++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/src/spin.jl b/src/spin.jl index 27f268ef..a7660a7d 100644 --- a/src/spin.jl +++ b/src/spin.jl @@ -83,3 +83,21 @@ Spin down state for the given Spin basis. """ spindown(::Type{T}, b::SpinBasis) where T = basisstate(T, b, b.shape[1]) spindown(b::SpinBasis) = spindown(ComplexF64, b) + +""" + squeeze(b::SpinBasis, x::Complex) + +This function generates a spin squeeze operator +along squeezing direction specified by arg(x)/2 + exp(0.5/N*( conj(x)*Jm^2 - x*Jp^2 )) +note that due to the finiteness of the Hilbert space setting a too large + squeezing x will create an oversqueezed state. For nice squeezing + x should be smaller than sqrt(N). +""" +function squeeze(b::SpinBasis,x) + N = length(b)-1 + Jm = sigmam(b)/2 + Jp = sigmap(b)/2 + s = exp(conj(x)*dense(Jm)^2/2/N - x*dense(Jp)^2/2/N) + return s +end \ No newline at end of file diff --git a/test/test_spin.jl b/test/test_spin.jl index 9bfb67d6..40d52f15 100644 --- a/test/test_spin.jl +++ b/test/test_spin.jl @@ -113,4 +113,21 @@ antikommutator(x, y) = x*y + y*x @test 1e-11 > norm(sm*spinup(spinbasis) - spindown(spinbasis)) @test 1e-11 > norm(sp*spindown(spinbasis) - spinup(spinbasis)) +# squeeze oerator test SPIN +Nspins = 500 +b_spin = SpinBasis(Nspins//2) +ss = squeeze(b_spin,rand()*sqrt(Nspins)*exp(1im*rand()*pi)); +st = ss*spindown(b_spin); + +# Heisenberg uncertainty test +@test abs(variance(sigmax(b_spin)/2,st)*variance(sigmay(b_spin)/2,st)) ≥ abs2(expect(sigmaz(b_spin)/2,st))/4 + +# small squeezing test +xx = rand()/10 +ss = squeeze(b_spin,xx*sqrt(Nspins)); +st = ss*spindown(b_spin); + +@test isapprox(2*log(real(variance(sigmax(b_spin)/2,st))/Nspins*4) , -xx*sqrt(Nspins), atol=1e-2) +@test isapprox(2*log(real(variance(sigmay(b_spin)/2,st))/Nspins*4) , xx*sqrt(Nspins), atol=1e-2) + end # testset From dc7edcb80631750ba0ce21b2d1700dabff4beb1c Mon Sep 17 00:00:00 2001 From: christoph Date: Tue, 22 Aug 2023 13:09:29 +0200 Subject: [PATCH 02/14] API fix --- docs/src/api.md | 18 ++++++++++++++++++ src/operators_sparse.jl | 6 ++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index edb392db..2dec73d6 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -594,3 +594,21 @@ DenseChiMatrix ```@docs QuantumOpticsBase.set_printing ``` + +## [LazyTensor functions](@id API: LazyTensor functions) + +```@docs +lazytensor_enable_cache +``` + +```@docs +lazytensor_disable_cache +``` + +```@docs +lazytensor_cachesize +``` + +```@docs +lazytensor_clear_cache +``` \ No newline at end of file diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl index 7765530e..f22fe07b 100644 --- a/src/operators_sparse.jl +++ b/src/operators_sparse.jl @@ -6,14 +6,12 @@ const SparseOpPureType{BL,BR} = Operator{BL,BR,<:SparseMatrixCSC} const SparseOpAdjType{BL,BR} = Operator{BL,BR,<:Adjoint{<:Number,<:SparseMatrixCSC}} const SparseOpType{BL,BR} = Union{SparseOpPureType{BL,BR},SparseOpAdjType{BL,BR}} - """ SparseOperator(b1[, b2, data]) Sparse array implementation of Operator. -The matrix is stored as the julia built-in type `SparseMatrixCSC` -in the `data` field. +The matrix is stored as the julia built-in type `SparseMatrixCSC` in the `data` field. """ SparseOperator(b1::Basis, b2::Basis, data) = Operator(b1, b2, SparseMatrixCSC(data)) SparseOperator(b1::Basis, b2::Basis, data::SparseMatrixCSC) = Operator(b1, b2, data) @@ -89,7 +87,7 @@ identityoperator(::Type{T}, b::Basis) where T<:Number = identityoperator(DataOpe """ diagonaloperator(b::Basis) -Create a diagonal operator of type [`SparseOperator`](@ref). +Create a diagonal operator of type [`SparseOperator`](@ref). """ function diagonaloperator(b::Basis, diag) @assert 1 <= length(diag) <= length(b) From cf51bf053b65db8f7e026fdb396d0615b96c9947 Mon Sep 17 00:00:00 2001 From: Christoph Hotter <52912924+ChristophHotter@users.noreply.github.com> Date: Tue, 22 Aug 2023 13:29:03 +0200 Subject: [PATCH 03/14] Update operators_sparse.jl --- src/operators_sparse.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl index f22fe07b..7c09dd75 100644 --- a/src/operators_sparse.jl +++ b/src/operators_sparse.jl @@ -87,7 +87,7 @@ identityoperator(::Type{T}, b::Basis) where T<:Number = identityoperator(DataOpe """ diagonaloperator(b::Basis) -Create a diagonal operator of type [`SparseOperator`](@ref). +Create a diagonal operator of type [`SparseOperator`](@ref). """ function diagonaloperator(b::Basis, diag) @assert 1 <= length(diag) <= length(b) From 9bb478f2eb535d3805fb21c5982c24677232b901 Mon Sep 17 00:00:00 2001 From: christoph Date: Tue, 22 Aug 2023 17:59:56 +0200 Subject: [PATCH 04/14] squeeze function --- src/spin.jl | 23 +++++++++++------------ test/test_spin.jl | 18 +++++++++--------- 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/src/spin.jl b/src/spin.jl index a7660a7d..f724f723 100644 --- a/src/spin.jl +++ b/src/spin.jl @@ -84,20 +84,19 @@ Spin down state for the given Spin basis. spindown(::Type{T}, b::SpinBasis) where T = basisstate(T, b, b.shape[1]) spindown(b::SpinBasis) = spindown(ComplexF64, b) + """ - squeeze(b::SpinBasis, x::Complex) + squeeze([T=ComplexF64,] b::SpinBasis, z) -This function generates a spin squeeze operator -along squeezing direction specified by arg(x)/2 - exp(0.5/N*( conj(x)*Jm^2 - x*Jp^2 )) -note that due to the finiteness of the Hilbert space setting a too large - squeezing x will create an oversqueezed state. For nice squeezing - x should be smaller than sqrt(N). +Squeezing operator ``S(z)=\\exp{\\left(\\frac{z^*\\hat{J_-}^2 - z\\hat{J}_+}{2 N}\\right)}`` for the +specified Spin-``N/2`` basis with optional data type `T`, computed as the matrix exponential. Too +large squeezing (``|z| > \\sqrt{N}``) will create an oversqueezed state. """ -function squeeze(b::SpinBasis,x) - N = length(b)-1 +function squeeze(::Type{T}, b::SpinBasis, z::Number) where T + N = Int(length(b)-1) + z = T(z)/2N Jm = sigmam(b)/2 Jp = sigmap(b)/2 - s = exp(conj(x)*dense(Jm)^2/2/N - x*dense(Jp)^2/2/N) - return s -end \ No newline at end of file + exp(conj(z)*dense(Jm)^2 - z*dense(Jp)^2) +end +squeeze(b::SpinBasis, z::T) where {T <: Number} = squeeze(ComplexF64, b, z) diff --git a/test/test_spin.jl b/test/test_spin.jl index 40d52f15..97c97925 100644 --- a/test/test_spin.jl +++ b/test/test_spin.jl @@ -113,21 +113,21 @@ antikommutator(x, y) = x*y + y*x @test 1e-11 > norm(sm*spinup(spinbasis) - spindown(spinbasis)) @test 1e-11 > norm(sp*spindown(spinbasis) - spinup(spinbasis)) -# squeeze oerator test SPIN +# squeeze operator test SPIN Nspins = 500 b_spin = SpinBasis(Nspins//2) -ss = squeeze(b_spin,rand()*sqrt(Nspins)*exp(1im*rand()*pi)); -st = ss*spindown(b_spin); +s1 = squeeze(b_spin,0.23*sqrt(Nspins)*exp(1im*0.23*pi)) +s2 = s1*spindown(b_spin); # Heisenberg uncertainty test -@test abs(variance(sigmax(b_spin)/2,st)*variance(sigmay(b_spin)/2,st)) ≥ abs2(expect(sigmaz(b_spin)/2,st))/4 +@test abs(variance(sigmax(b_spin)/2,s2)*variance(sigmay(b_spin)/2,s2)) ≥ abs2(expect(sigmaz(b_spin)/2,s2))/4 # small squeezing test -xx = rand()/10 -ss = squeeze(b_spin,xx*sqrt(Nspins)); -st = ss*spindown(b_spin); +x1 = 0.023 +s1 = squeeze(b_spin,x1*sqrt(Nspins)); +s2 = s1*spindown(b_spin); -@test isapprox(2*log(real(variance(sigmax(b_spin)/2,st))/Nspins*4) , -xx*sqrt(Nspins), atol=1e-2) -@test isapprox(2*log(real(variance(sigmay(b_spin)/2,st))/Nspins*4) , xx*sqrt(Nspins), atol=1e-2) +@test isapprox(2*log(real(variance(sigmax(b_spin)/2,s2))/Nspins*4) , -x1*sqrt(Nspins), atol=1e-2) +@test isapprox(2*log(real(variance(sigmay(b_spin)/2,s2))/Nspins*4) , x1*sqrt(Nspins), atol=1e-2) end # testset From 0b1c86bc01457341b6a0c7c67a8f3c116f7ef3b9 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 11 Sep 2023 10:25:29 -0400 Subject: [PATCH 05/14] Bump actions/checkout from 3 to 4 (#132) Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/breakage.yml | 2 +- .github/workflows/ci-jet.yml | 2 +- .github/workflows/ci.yml | 2 +- .github/workflows/docs.yml | 2 +- .github/workflows/invalidations.yml | 4 ++-- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/breakage.yml b/.github/workflows/breakage.yml index 5d63c92e..01367afb 100644 --- a/.github/workflows/breakage.yml +++ b/.github/workflows/breakage.yml @@ -20,7 +20,7 @@ jobs: pkgversion: [latest] steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 # Install Julia - uses: julia-actions/setup-julia@v1 diff --git a/.github/workflows/ci-jet.yml b/.github/workflows/ci-jet.yml index ffa8d9d7..bb388d79 100644 --- a/.github/workflows/ci-jet.yml +++ b/.github/workflows/ci-jet.yml @@ -17,7 +17,7 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 02b2a3a2..e75f27da 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,7 +23,7 @@ jobs: - os: macOS-latest arch: x86 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index e86cd0ee..56258a2a 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -13,7 +13,7 @@ jobs: build: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@latest with: version: '1' diff --git a/.github/workflows/invalidations.yml b/.github/workflows/invalidations.yml index af4c1580..7e0aefe0 100644 --- a/.github/workflows/invalidations.yml +++ b/.github/workflows/invalidations.yml @@ -19,12 +19,12 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: '1' - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-invalidations@v1 id: invs_pr - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: ref: ${{ github.event.repository.default_branch }} - uses: julia-actions/julia-buildpkg@v1 From 4d99ddf4bc9b0347fa53556125ae53f87abe8628 Mon Sep 17 00:00:00 2001 From: Stefan Krastanov Date: Mon, 11 Sep 2023 11:01:35 -0400 Subject: [PATCH 06/14] Routine static analysis fixes (detected by JET and Aqua) (#135) * misc JET and Aqua fixes * relax jet tests --- Project.toml | 3 ++- src/operators_lazytensor.jl | 7 +++---- src/sparsematrix.jl | 2 +- src/states.jl | 2 +- test/test_jet.jl | 5 +++-- 5 files changed, 10 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 15e41d0a..4ca150fc 100644 --- a/Project.toml +++ b/Project.toml @@ -38,8 +38,9 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" 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", "UnsafeArrays", "FillArrays"] +test = ["LinearAlgebra", "SparseArrays", "Random", "Test", "Aqua", "JET", "Adapt", "Dates", "FFTW", "LRUCache", "Strided", "StridedViews", "UnsafeArrays", "FillArrays"] diff --git a/src/operators_lazytensor.jl b/src/operators_lazytensor.jl index a6802f3b..2b4d1d54 100644 --- a/src/operators_lazytensor.jl +++ b/src/operators_lazytensor.jl @@ -96,7 +96,7 @@ isequal(x::LazyTensor, y::LazyTensor) = samebases(x,y) && isequal(x.indices, y.i # Arithmetic operations -(a::LazyTensor) = LazyTensor(a, -a.factor) -const single_dataoperator{B1,B2} = LazyTensor{B1,B2,F,I,Tuple{T}} where {B1,B2,F,I,T<:DataOperator} +const single_dataoperator{B1,B2} = LazyTensor{B1,B2,F,I,Tuple{T}} where {B1,B2,F,I,T<:DataOperator} function +(a::T1,b::T2) where {T1 <: single_dataoperator{B1,B2},T2 <: single_dataoperator{B1,B2}} where {B1,B2} if a.indices == b.indices op = a.operators[1] * a.factor + b.operators[1] * b.factor @@ -349,7 +349,7 @@ function _tp_matmul_mid!(result, a::AbstractMatrix, loc::Integer, b, α::Number, result_r = Base.ReshapedArray(result, (sz_b_1, size(a, 1), sz_b_3), ()) if a isa FillArrays.Eye - # Square Eyes are skipped higher up. This handles the non-square case. + # Square Eyes are skipped higher up. This handles the non-square case. size(b, loc) == size(a, 2) && size(result, loc) == size(a, 1) || throw(DimensionMismatch("Dimensions of Eye matrix do not match subspace dimensions.")) d = min(size(a)...) @@ -505,7 +505,7 @@ function _explicit_isometries(::Type{T}, used_indices, bl::Basis, br::Basis, shi iso_inds = [i + shift] else push!(isos, Eye{T}(sl, sr)) - push!(iso_inds, i + shift) + push!(iso_inds::Vector{Int}, i + shift) # Help with inference (detected by JET) end end end @@ -752,4 +752,3 @@ _mul_puresparse!(result::DenseOpType{B1,B3},h::LazyTensor{B1,B2,F,I,T},op::Dense _mul_puresparse!(result::DenseOpType{B1,B3},op::DenseOpType{B1,B2},h::LazyTensor{B2,B3,F,I,T},alpha,beta) where {B1,B2,B3,F,I,T} = (_gemm_puresparse(alpha, op.data, h, beta, result.data); result) _mul_puresparse!(result::Ket{B1},a::LazyTensor{B1,B2,F,I,T},b::Ket{B2},alpha,beta) where {B1,B2,F,I,T} = (_gemm_puresparse(alpha, a, b.data, beta, result.data); result) _mul_puresparse!(result::Bra{B2},a::Bra{B1},b::LazyTensor{B1,B2,F,I,T},alpha,beta) where {B1,B2,F,I,T} = (_gemm_puresparse(alpha, a.data, b, beta, result.data); result) - diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 45e08c97..1b4a2b70 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -246,7 +246,7 @@ function _permutedims(x::AbstractSparseMatrix, shape, perm) # TODO upstream as ` shape = (shape...,) shape_perm = ([shape[i] for i in perm]...,) y = spzeros(eltype(x), x.m, x.n) - for I in eachindex(x) + for I in eachindex(x)::CartesianIndices # Help with inference (detected by JET) linear_index = LinearIndices((x.m, x.n))[I.I...] cartesian_index = CartesianIndices(shape)[linear_index] cartesian_index_perm = [cartesian_index[i] for i=perm] diff --git a/src/states.jl b/src/states.jl index 6c6bdb5b..b632c177 100644 --- a/src/states.jl +++ b/src/states.jl @@ -250,4 +250,4 @@ end @inline Base.copyto!(dest::Bra{B1}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {B1,B2,Style<:BraStyle{B2},Axes,F,Args} = throw(IncompatibleBases()) -@inline Base.copyto!(A::T,B::T) where T<:StateVector = (copyto!(A.data,B.data); A) +@inline Base.copyto!(A::T,B::T) where T<:Union{Ket, Bra} = (copyto!(A.data,B.data); A) # Can not use T<:QuantumInterface.StateVector, because StateVector does not imply the existence of a data property diff --git a/test/test_jet.jl b/test/test_jet.jl index 8d33eb57..9dfed8cf 100644 --- a/test/test_jet.jl +++ b/test/test_jet.jl @@ -19,7 +19,7 @@ function (::MayThrowIsOk)(report_type::Type{<:InferenceErrorReport}, @nospeciali end # imported to be declared as modules filtered out from analysis result -using LinearAlgebra, LRUCache, Strided, Dates, SparseArrays +using LinearAlgebra, LRUCache, Strided, StridedViews, Dates, SparseArrays @testset "jet" begin if get(ENV,"JET_TEST","")=="true" @@ -29,11 +29,12 @@ using LinearAlgebra, LRUCache, Strided, Dates, SparseArrays AnyFrameModule(LinearAlgebra), AnyFrameModule(LRUCache), AnyFrameModule(Strided), + AnyFrameModule(StridedViews), AnyFrameModule(Dates), AnyFrameModule(SparseArrays)) ) @show rep - @test length(JET.get_reports(rep)) <= 1 + @test length(JET.get_reports(rep)) <= 6 @test_broken length(JET.get_reports(rep)) == 0 end end # testset From 19a9d53f115991fd9e29204012bf7de3b1e91f31 Mon Sep 17 00:00:00 2001 From: Akira Kyle Date: Tue, 12 Sep 2023 14:38:40 -0600 Subject: [PATCH 07/14] Fix sparse matrix exponential of zero since FastExpm doesn't handle it (#136) --- src/operators_sparse.jl | 6 +++++- src/superoperators.jl | 9 ++++++++- test/test_operators.jl | 1 + test/test_superoperators.jl | 1 + 4 files changed, 15 insertions(+), 2 deletions(-) diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl index 7c09dd75..b52d5ce7 100644 --- a/src/operators_sparse.jl +++ b/src/operators_sparse.jl @@ -56,7 +56,11 @@ If you only need the result of the exponential acting on a vector, consider using much faster implicit methods that do not calculate the entire exponential. """ function exp(op::T; opts...) where {B,T<:SparseOpType{B,B}} - return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...)) + if iszero(op) + return identityoperator(op) + else + return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...)) + end end function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis} diff --git a/src/superoperators.jl b/src/superoperators.jl index 3a137752..8128c3e5 100644 --- a/src/superoperators.jl +++ b/src/superoperators.jl @@ -241,9 +241,16 @@ All optional arguments are passed to `fastExpm` and can be used to specify toler If you only need the result of the exponential acting on an operator, consider using much faster implicit methods that do not calculate the entire exponential. """ -Base.exp(op::SparseSuperOpType; opts...) = SuperOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...)) +function Base.exp(op::SparseSuperOpType; opts...) + if iszero(op) + return identitysuperoperator(op) + else + return SuperOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...)) + end +end # Array-like functions +Base.zero(A::SuperOperator) = SuperOperator(A.basis_l, A.basis_r, zero(A.data)) Base.size(A::SuperOperator) = size(A.data) @inline Base.axes(A::SuperOperator) = axes(A.data) Base.ndims(A::SuperOperator) = 2 diff --git a/test/test_operators.jl b/test/test_operators.jl index 92fc8436..5e180441 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -118,6 +118,7 @@ b = FockBasis(40) alpha = 1+5im H = alpha * create(b) - conj(alpha) * destroy(b) @test exp(sparse(H); threshold=1e-10) ≈ displace(b, alpha) +@test exp(sparse(zero(identityoperator(b)))) ≈ identityoperator(b) @test one(b1).data == Diagonal(ones(b1.shape[1])) @test one(op1).data == Diagonal(ones(b1.shape[1])) diff --git a/test/test_superoperators.jl b/test/test_superoperators.jl index 28043bdb..204264ac 100644 --- a/test/test_superoperators.jl +++ b/test/test_superoperators.jl @@ -221,6 +221,7 @@ Ldense .+= L b = FockBasis(20) L = liouvillian(identityoperator(b), [destroy(b)]) @test exp(sparse(L)).data ≈ exp(dense(L)).data +@test exp(sparse(zero(identitysuperoperator(b)))).data ≈ identitysuperoperator(b).data N = exp(log(2) * sparse(L)) # 50% loss channel ρ = N * dm(coherentstate(b, 1)) @test (1 - real(tr(ρ^2))) < 1e-10 # coherent state remains pure under loss From 032dd824ca25ddf9c82e81b72c919e41952af16b Mon Sep 17 00:00:00 2001 From: Abhishek Bhatt <46929125+Abhishek-1Bhatt@users.noreply.github.com> Date: Wed, 13 Sep 2023 20:46:51 -0400 Subject: [PATCH 08/14] Make `apply!` work when used with any order of subsystem indices (#133) --------- Co-authored-by: Stefan Krastanov --- src/apply.jl | 19 +++++++++++++++++++ test/runtests.jl | 5 +++-- test/test_apply.jl | 25 +++++++++++++++++++++++-- test/test_time_dependent_operators.jl | 10 +++++----- 4 files changed, 50 insertions(+), 9 deletions(-) diff --git a/src/apply.jl b/src/apply.jl index 504578a2..b978cca6 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -1,11 +1,27 @@ function apply!(state::Ket, indices, operation::Operator) op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation) + if (length(indices)>1) + for i in 2:length(indices) + if indices[i]1) + for i in 2:length(indices) + if indices[i]1 + error("Applying SuperOperator to multiple qubits/operators is not supported currently, due to missing tensor product method for SuperOperators") + end op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation) state.data = (op*state).data state diff --git a/test/runtests.jl b/test/runtests.jl index 92456a83..5a719baf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,9 +36,10 @@ names = [ "test_printing.jl", + "test_apply.jl", + "test_aqua.jl", - "test_jet.jl", - "test_apply.jl" + "test_jet.jl" ] detected_tests = filter( diff --git a/test/test_apply.jl b/test/test_apply.jl index 0cf3f61a..a7c4fd12 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -3,13 +3,13 @@ using QuantumInterface using Test @testset "apply" begin - + _b2 = SpinBasis(1//2) _l0 = spinup(_b2) _l1 = spindown(_b2) _x = sigmax(_b2) _y = sigmay(_b2) - +_z = sigmaz(_b2) @test QuantumOpticsBase.apply!(_l0⊗_l1, 1, _x) == _l1⊗_l1 @test QuantumOpticsBase.apply!(_x, 1, _y) == _x @@ -52,4 +52,25 @@ catch e #@test typeof(e) <: QuantumInterface.IncompatibleBases end +# test CNOT₂₋₁ +CNOT = DenseOperator(_b2⊗_b2, Complex.([1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0])) +@test QuantumOpticsBase.apply!(_l0⊗_l1, [2,1], CNOT) == _l1⊗_l1 + +# test operator permutation with 3 qubits/operators for apply! + +# 3-qubit operator permutation +@test QuantumOpticsBase.apply!(_l0⊗_l1⊗_l0, [2,3,1], _x⊗_y⊗_z) ≈ (_y*_l0)⊗(_z*_l1)⊗(_x*_l0) + +# 3-operator operator permutation +@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3,1], _x⊗_y⊗_z) ≈ (_y⊗_z⊗_x)*(_x⊗_y⊗_z)*(_y⊗_z⊗_x)' + +# 3-qubit/operator errors when called for applying superoperator +sOp1 = spre(create(FockBasis(1))) +sOp2 = spost(create(FockBasis(1))) +st1 = coherentstate(FockBasis(1), 1.4) +st2 = coherentstate(FockBasis(1), 0.3) +st3 = coherentstate(FockBasis(1), 0.8) + +@test_throws ErrorException QuantumOpticsBase.apply!(st1⊗st2⊗st3, [2,3,1], sOp1) + end #testset diff --git a/test/test_time_dependent_operators.jl b/test/test_time_dependent_operators.jl index 451c5468..fbc4e664 100644 --- a/test/test_time_dependent_operators.jl +++ b/test/test_time_dependent_operators.jl @@ -90,7 +90,7 @@ using LinearAlgebra, Random o_t1_comp2_ = embed(FockBasis(2) ⊗ GenericBasis(2), [1], o_t1) @test o_t1_comp2 == o_t1_comp2_ - + o_t = TimeDependentSum([f1, f2], [a,n]) @test !QOB.is_const(o_t) @@ -122,7 +122,7 @@ using LinearAlgebra, Random b = randoperator(FockBasis(2)) o2_t = TimeDependentSum([t->cos(t), t->t/3.0], [b, n]) - + # operations o_res = o_t(0.0) + o2_t(0.0) @test isa(o_res, TimeDependentSum) @@ -138,7 +138,7 @@ using LinearAlgebra, Random @test isa(o_res, TimeDependentSum) @test QOB.eval_coefficients(o_res, t) == [-1.0im, -t*3.0] @test all(QOB.suboperators(o_res) .== (a, n)) - + o_res = a + o2_t @test isa(o_res, TimeDependentSum) @test QOB.eval_coefficients(o_res, t) == ComplexF64[1.0, cos(t), t/3.0] @@ -148,7 +148,7 @@ using LinearAlgebra, Random @test isa(o_res, TimeDependentSum) @test QOB.eval_coefficients(o_res, t) == ComplexF64[1.0, -cos(t), -t/3.0] @test all(QOB.suboperators(o_res) .== (a, b, n)) - + o_res = LazySum(a, n) + o2_t @test isa(o_res, TimeDependentSum) @test QOB.eval_coefficients(o_res, t) == ComplexF64[1.0, 1.0, cos(t), t/3.0] @@ -221,7 +221,7 @@ end @test @inferred QOB.eval_coefficients(op_block_shift, 4.9) == (0.0, 0.0) @test @inferred QOB.eval_coefficients(op_block_shift, 5.0) == (1.0, -5.0) @test @inferred QOB.eval_coefficients(op_block_shift, 20.0) == (0.0, 0.0) - + op_stretch_block_shift = time_stretch(op_block_shift, 2.0) @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 9.9) == (0.0, 0.0) @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 10.0) == (1.0, -5.0) From ee6330d46e7a1248defaa9c195b2b6eafe2bd68f Mon Sep 17 00:00:00 2001 From: Stefan Krastanov Date: Thu, 14 Sep 2023 14:07:34 -0400 Subject: [PATCH 09/14] bump julia compat to 1.6 and add 1.6 to CI (#141) * add julia 1.3 to CI * bump julia compat to 1.6 * fix parsing error on julia 1.6 --- .github/workflows/ci.yml | 4 ++++ Project.toml | 2 +- test/test_time_dependent_operators.jl | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e75f27da..f64f7eb2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,6 +22,10 @@ jobs: exclude: - os: macOS-latest arch: x86 + include: + - os: ubuntu-latest + arch: x64 + version: '1.6' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 diff --git a/Project.toml b/Project.toml index 4ca150fc..56f8ab0a 100644 --- a/Project.toml +++ b/Project.toml @@ -24,7 +24,7 @@ LRUCache = "1" QuantumInterface = "0.3.0" Strided = "1, 2" UnsafeArrays = "1" -julia = "1.3" +julia = "1.6" [extras] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/test/test_time_dependent_operators.jl b/test/test_time_dependent_operators.jl index fbc4e664..29b4eaac 100644 --- a/test/test_time_dependent_operators.jl +++ b/test/test_time_dependent_operators.jl @@ -225,7 +225,7 @@ end op_stretch_block_shift = time_stretch(op_block_shift, 2.0) @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 9.9) == (0.0, 0.0) @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 10.0) == (1.0, -5.0) - @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 30.0-√eps()) == (1.0, 5.0-0.5√eps()) + @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 30.0-√eps()) == (1.0, 5.0-0.5*√eps()) @test @inferred QOB.eval_coefficients(op_stretch_block_shift, 30.0) == (0.0, 0.0) op_block = time_restrict(op, 15.0) From 7f0b8a17fe0d1bfca419dfc0dd6b828a1546ba74 Mon Sep 17 00:00:00 2001 From: Stefan Krastanov Date: Thu, 14 Sep 2023 15:41:02 -0400 Subject: [PATCH 10/14] bump version number to 0.4.15 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 56f8ab0a..f9680dd6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "QuantumOpticsBase" uuid = "4f57444f-1401-5e15-980d-4471b28d5678" -version = "0.4.14" +version = "0.4.15" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 97ecd53f420717df31b6b7fe9e8fdc66d106a8ef Mon Sep 17 00:00:00 2001 From: Akira Kyle Date: Sun, 17 Sep 2023 15:49:25 -0600 Subject: [PATCH 11/14] Use FastGaussQuadrature for transform between Position and Fock bases (#142) --- Project.toml | 2 + src/QuantumOpticsBase.jl | 1 - src/polynomials.jl | 96 ---------------------------------------- src/transformations.jl | 41 +++++------------ test/runtests.jl | 1 - test/test_polynomials.jl | 28 ------------ 6 files changed, 13 insertions(+), 156 deletions(-) delete mode 100644 src/polynomials.jl delete mode 100644 test/test_polynomials.jl diff --git a/Project.toml b/Project.toml index f9680dd6..167d41ba 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.4.15" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -19,6 +20,7 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6" Adapt = "1, 2, 3.3" FFTW = "1.2" FastExpm = "1.1.0" +FastGaussQuadrature = "0.5" FillArrays = "0.13, 1" LRUCache = "1" QuantumInterface = "0.3.0" diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl index 1ee5694e..14608f2e 100644 --- a/src/QuantumOpticsBase.jl +++ b/src/QuantumOpticsBase.jl @@ -65,7 +65,6 @@ export Basis, GenericBasis, CompositeBasis, basis, #apply apply! -include("polynomials.jl") include("bases.jl") include("states.jl") include("operators.jl") diff --git a/src/polynomials.jl b/src/polynomials.jl deleted file mode 100644 index ef5615c5..00000000 --- a/src/polynomials.jl +++ /dev/null @@ -1,96 +0,0 @@ -""" - horner(coefficients, x) - -Evaluate the given polynomial at position x using the Horner scheme. - -```math -p(x) = \\sum_{n=0}^N c_n x^n -``` -""" -function horner(coefficients, x) - bn = coefficients[end] - for n=length(coefficients)-1:-1:1 - bn = coefficients[n] + bn*x - end - bn -end - - -module hermite - -""" - a_nk(N) - -Calculate the all coefficients for all Hermite polynomials up to order N. - -```math -H_n(x) = \\sum_{k=0}^n a_{n,k} x^k -``` - -Returns a vector of length N+1 where the n-th entry contains all coefficients -for the n-th Hermite polynomial. -""" -function a(N::T) where T - a = Vector{Vector{T}}(undef, N+1) - a[1] = [1] - a[2] = [0,2] - am = a[2] - for n=2:N - an = zeros(Int, n+1) - a[n+1] = an - if iseven(n) - an[1] = -am[2] - end - an[n+1] = 2*am[n] - if iseven(n) - for k=3:2:n-1 - an[k] = 2*am[k-1] - am[k+1]*k - end - else - for k=2:2:n-1 - an[k] = 2*am[k-1] - am[k+1]*k - end - end - am = an - end - a -end - -""" - A_nk(N) - -Calculate the all scaled coefficients for all Hermite polynomials up to order N. - -The scaled coefficients `A` are connected to the unscaled coefficients `a` by -the relation ``A_{n,k} = \\frac{a_{n,k}}{\\sqrt{2^n n!}}`` - -Returns a vector of length N+1 where the n-th entry contains all scaled -coefficients for the n-th Hermite polynomial. -""" -function A(N::T) where T - A = Vector{Vector{float(T)}}(undef, N+1) - A[1] = [1.] - A[2] = [0., sqrt(2)] - Am = A[2] - for n=2:N - An = zeros(float(T), n+1) - A[n+1] = An - if iseven(n) - An[1] = -Am[2]/sqrt(2*n) - end - An[n+1] = Am[n]*sqrt(2/n) - if iseven(n) - for k=3:2:n-1 - An[k] = Am[k-1]*sqrt(2/n) - Am[k+1]*k/sqrt(2*n) - end - else - for k=2:2:n-1 - An[k] = Am[k-1]*sqrt(2/n) - Am[k+1]*k/sqrt(2*n) - end - end - Am = An - end - A -end - -end # module hermite diff --git a/src/transformations.jl b/src/transformations.jl index 9f0b86c5..bbba8813 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -1,3 +1,5 @@ +import FastGaussQuadrature: hermpoly_rec + """ transform([S=ComplexF64, ]b1::PositionBasis, b2::FockBasis; x0=1) transform([S=ComplexF64, ]b1::FockBasis, b2::PositionBasis; x0=1) @@ -15,38 +17,17 @@ a harmonic trap potential at position ``x``, i.e.: \\frac{1}{\\sqrt{2^n n!}} H_n\\left(\\frac{x}{x_0}\\right) ``` """ -function transform(::Type{S}, b1::PositionBasis, b2::FockBasis; x0=1) where S - T = Matrix{S}(undef, length(b1), length(b2)) - xvec = samplepoints(b1) - A = hermite.A(b2.N) - delta_x = spacing(b1) - c = 1.0/sqrt(x0*sqrt(pi))*sqrt(delta_x) - for i in 1:length(b1) - u = xvec[i]/x0 - C = c*exp(-u^2/2) - for n=b2.offset:b2.N - idx = n-b2.offset+1 - T[i,idx] = C*horner(A[n+1], u) - end +function transform(::Type{S}, bp::PositionBasis, bf::FockBasis; x0=1) where S + T = Matrix{S}(undef, length(bp), length(bf)) + xvec = samplepoints(bp) + C = pi^(-1/4)*sqrt(spacing(bp)/x0) + for i in 1:length(bp) + T[i,:] = C*hermpoly_rec(bf.offset:bf.N, sqrt(2)*xvec[i]/x0) end - DenseOperator(b1, b2, T) + DenseOperator(bp, bf, T) end -function transform(::Type{S}, b1::FockBasis, b2::PositionBasis; x0=1) where S - T = Matrix{S}(undef, length(b1), length(b2)) - xvec = samplepoints(b2) - A = hermite.A(b1.N) - delta_x = spacing(b2) - c = 1.0/sqrt(x0*sqrt(pi))*sqrt(delta_x) - for i in 1:length(b2) - u = xvec[i]/x0 - C = c*exp(-u^2/2) - for n=b1.offset:b1.N - idx = n-b1.offset+1 - T[idx,i] = C*horner(A[n+1], u) - end - end - DenseOperator(b1, b2, T) -end +transform(::Type{S}, bf::FockBasis, bp::PositionBasis; x0=1) where S = + dagger(transform(S, bp, bf; x0=x0)) transform(b1::Basis,b2::Basis;kwargs...) = transform(ComplexF64,b1,b2;kwargs...) diff --git a/test/runtests.jl b/test/runtests.jl index 5a719baf..634a8bc9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,5 @@ names = [ "test_sortedindices.jl", - "test_polynomials.jl", "test_bases.jl", "test_states.jl", diff --git a/test/test_polynomials.jl b/test/test_polynomials.jl deleted file mode 100644 index b9299a37..00000000 --- a/test/test_polynomials.jl +++ /dev/null @@ -1,28 +0,0 @@ -using Test -using QuantumOpticsBase: horner, hermite - -@testset "polynomials" begin - -# Test Horner scheme -c = [0.2, 0.6, 1.7] -x0 = 1.3 -@test horner(c, x0) == c[1] + c[2]*x0 + c[3]*x0^2 - -# Test Hermite polynomials -an = Vector{Vector{Int}}(undef, 8) -an[1] = [1] -an[2] = [0,2] -an[3] = [-2, 0, 4] -an[4] = [0, -12, 0, 8] -an[5] = [12, 0, -48, 0, 16] -an[6] = [0, 120, 0, -160, 0, 32] -an[7] = [-120, 0, 720, 0, -480, 0, 64] -an[8] = [0, -1680, 0, 3360, 0, -1344, 0, 128] -@test hermite.a(7) == an - -A = hermite.A(7) -for n=0:7 - @test A[n+1] ≈ an[n+1]/sqrt(2^n*factorial(n)) -end - -end # testset From 4c6296c57d0ddc2bd8006f0708b6675c0696ef4e Mon Sep 17 00:00:00 2001 From: Stefan Krastanov Date: Mon, 18 Sep 2023 12:16:51 -0400 Subject: [PATCH 12/14] Fixes and tests for `apply!` (#144) --- Project.toml | 2 +- src/apply.jl | 51 +++++++++++++++++++---------------- test/test_apply.jl | 67 +++++++++++++++------------------------------- 3 files changed, 51 insertions(+), 69 deletions(-) diff --git a/Project.toml b/Project.toml index 167d41ba..ce682fa9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "QuantumOpticsBase" uuid = "4f57444f-1401-5e15-980d-4471b28d5678" -version = "0.4.15" +version = "0.4.16" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/apply.jl b/src/apply.jl index b978cca6..24e2d7f7 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -1,27 +1,32 @@ -function apply!(state::Ket, indices, operation::Operator) - op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation) - if (length(indices)>1) - for i in 2:length(indices) - if indices[i]1) - for i in 2:length(indices) - if indices[i]1 - error("Applying SuperOperator to multiple qubits/operators is not supported currently, due to missing tensor product method for SuperOperators") + if is_apply_shortcircuit(state, indices, operation) + state.data = (operation*state).data + return state + else + error("`apply!` does not yet support embedding superoperators acting only on a subsystem of the given state") end - op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation) - state.data = (op*state).data - state end diff --git a/test/test_apply.jl b/test/test_apply.jl index a7c4fd12..36bdd9a9 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,6 +1,7 @@ using QuantumOpticsBase using QuantumInterface using Test +using QuantumOpticsBase: apply! @testset "apply" begin @@ -11,66 +12,42 @@ _x = sigmax(_b2) _y = sigmay(_b2) _z = sigmaz(_b2) -@test QuantumOpticsBase.apply!(_l0⊗_l1, 1, _x) == _l1⊗_l1 -@test QuantumOpticsBase.apply!(_x, 1, _y) == _x +@test apply!(_l0⊗_l1, 1, _x) ≈ _l1⊗_l1 +@test apply!(_x, 1, _y) ≈ _x # Test Operator with IncompatibleBases _l01 = _l0⊗_l1 op = projector(_l0, _l01') - -try - QuantumOpticsBase.apply!(_l0, 1, op) -catch e - @test typeof(e) <: QuantumInterface.IncompatibleBases -end - -try - QuantumOpticsBase.apply!(_x, 1, op) -catch e - @test typeof(e) <: QuantumInterface.IncompatibleBases -end +@test_throws QuantumInterface.IncompatibleBases apply!(_l0, 1, op) +@test_throws QuantumInterface.IncompatibleBases apply!(_x, 1, op) # Test SuperOperator -sOp = spre(create(FockBasis(1))) -st = coherentstate(FockBasis(1), 1.4) - -@test QuantumOpticsBase.apply!(st, 1, sOp).data == (sOp*projector(st)).data +bf1 = FockBasis(1) +sOp = spre(create(bf1)) +st = coherentstate(bf1, 1.4) +@test apply!(st, 1, sOp).data ≈ (sOp*projector(st)).data # Test SuperOperator with IncompatibleBases -b1 = FockBasis(1) -b2 = FockBasis(2) - -k1 = coherentstate(b1, 0.39) -k2 = coherentstate(b2, 1.4) -op = projector(k1, k2') - -try - sOp2 = spre(op) - QuantumOpticsBase.apply!(st, 1, sOp2) -catch e - @test typeof(e) <: ArgumentError - #@test typeof(e) <: QuantumInterface.IncompatibleBases -end +bf2 = FockBasis(2) +op = create(bf2) +sOp2 = spre(op) +@test_throws ArgumentError apply!(st, 1, sOp2) # test CNOT₂₋₁ CNOT = DenseOperator(_b2⊗_b2, Complex.([1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0])) -@test QuantumOpticsBase.apply!(_l0⊗_l1, [2,1], CNOT) == _l1⊗_l1 +@test QuantumOpticsBase.apply!(_l0⊗_l1, [2,1], CNOT) ≈ _l1⊗_l1 # test operator permutation with 3 qubits/operators for apply! +@test QuantumOpticsBase.apply!(_l0⊗_l1⊗_l0, [2,3,1], _x⊗_y⊗_z) ≈ (_z*_l0)⊗(_x*_l1)⊗(_y*_l0) +@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3,1], _x⊗_y⊗_z) ≈ (_z⊗_x⊗_y)*(_x⊗_y⊗_z)*(_z⊗_x⊗_y)' -# 3-qubit operator permutation -@test QuantumOpticsBase.apply!(_l0⊗_l1⊗_l0, [2,3,1], _x⊗_y⊗_z) ≈ (_y*_l0)⊗(_z*_l1)⊗(_x*_l0) - -# 3-operator operator permutation -@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3,1], _x⊗_y⊗_z) ≈ (_y⊗_z⊗_x)*(_x⊗_y⊗_z)*(_y⊗_z⊗_x)' +i = identityoperator(_y) +@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3], _y⊗_z) ≈ (i⊗_y⊗_z)*(_x⊗_y⊗_z)*(i⊗_y⊗_z)' # 3-qubit/operator errors when called for applying superoperator -sOp1 = spre(create(FockBasis(1))) -sOp2 = spost(create(FockBasis(1))) -st1 = coherentstate(FockBasis(1), 1.4) -st2 = coherentstate(FockBasis(1), 0.3) -st3 = coherentstate(FockBasis(1), 0.8) - -@test_throws ErrorException QuantumOpticsBase.apply!(st1⊗st2⊗st3, [2,3,1], sOp1) +sOp1 = spre(create(bf1)) +st1 = coherentstate(bf1, 1.4) +st2 = coherentstate(bf1, 0.3) +@test_broken apply!(st1⊗st2, [1], sOp1) end #testset From 3a637602fbfb5e7777e4b053ea4d6aeab7fdaca0 Mon Sep 17 00:00:00 2001 From: Akira Kyle Date: Thu, 21 Sep 2023 15:16:01 -0600 Subject: [PATCH 13/14] 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 From 06809f719144860f29a9cb91b67006d4ddf96dd1 Mon Sep 17 00:00:00 2001 From: Abhishek Bhatt <46929125+Abhishek-1Bhatt@users.noreply.github.com> Date: Fri, 22 Sep 2023 01:47:29 -0400 Subject: [PATCH 14/14] move nsubsystems to QuantumInterface (#146) Co-authored-by: Stefan Krastanov --- Project.toml | 4 ++-- src/QuantumOpticsBase.jl | 2 +- src/apply.jl | 8 -------- 3 files changed, 3 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index e5a3a61f..6cd6f76b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "QuantumOpticsBase" uuid = "4f57444f-1401-5e15-980d-4471b28d5678" -version = "0.4.16" +version = "0.4.17" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -24,7 +24,7 @@ FastExpm = "1.1.0" FastGaussQuadrature = "0.5" FillArrays = "0.13, 1" LRUCache = "1" -QuantumInterface = "0.3.0" +QuantumInterface = "0.3.3" RandomMatrices = "0.5" Strided = "1, 2" UnsafeArrays = "1" diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl index eea0f2e3..c552186b 100644 --- a/src/QuantumOpticsBase.jl +++ b/src/QuantumOpticsBase.jl @@ -3,7 +3,7 @@ module QuantumOpticsBase using SparseArrays, LinearAlgebra, LRUCache, Strided, UnsafeArrays, FillArrays import LinearAlgebra: mul!, rmul! -import QuantumInterface: dagger, directsum, ⊕, dm, embed, expect, identityoperator, identitysuperoperator, +import QuantumInterface: dagger, directsum, ⊕, dm, embed, nsubsystems, expect, identityoperator, identitysuperoperator, permutesystems, projector, ptrace, reduced, tensor, ⊗, variance, apply!, basis, AbstractSuperOperator # index helpers diff --git a/src/apply.jl b/src/apply.jl index 24e2d7f7..2336cb6d 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -1,11 +1,3 @@ -nsubsystems(s::Ket) = nsubsystems(s.basis) -function nsubsystems(s::Operator) - s.basis_l == s.basis_r || throw(ArgumentError("`nsubsystem(::Operator)` is well defined only if the left and right bases are the same")) - nsubsystems(s.basis_l) -end -nsubsystems(b::CompositeBasis) = length(b.bases) -nsubsystems(b::Basis) = 1 - function is_apply_shortcircuit(state, indices, operation) if nsubsystems(state) == 1 basis(state)==basis(operation) || throw(ArgumentError("`apply!` failed due to incompatible bases of the state and the operation attempted to be applied on it"))