Functionality of expect -> expect_waveguide

expect(O,psi) previously returned a sum over all time indeces. But this behaviour was unstable when combined with other basises. Now instead expect_waveguide(op,psi) gives this over all waveguide indeces instead.
Updated docs accordingly

Project.toml

@@ -1,7 +1,7 @@
 name = "WaveguideQED"
 uuid = "c4555495-0e8d-488d-8e6a-965ecefe9055"
 authors = ["Matias Bundgaard-Nielsen <[email protected]>, Mikkel Heuck <[email protected]>, and Stefan Krastanov <[email protected]>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

docs/Manifest.toml

@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.9.0"
+julia_version = "1.10.0"
 manifest_format = "2.0"
 project_hash = "d6a5993505edb88a10fbee1b0a4d47aa45a4e462"
 
@@ -190,7 +190,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.2+0"
+version = "1.0.5+1"
 
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
@@ -633,21 +633,26 @@ version = "1.0.0"
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.3"
+version = "0.6.4"
 
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "7.84.0+0"
+version = "8.4.0+0"
 
 [[deps.LibGit2]]
-deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[deps.LibGit2_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
+uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
+version = "1.6.4+0"
+
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.10.2+0"
+version = "1.11.0+1"
 
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -772,7 +777,7 @@ version = "1.1.7"
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.2+0"
+version = "2.28.2+1"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -785,7 +790,7 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2022.10.11"
+version = "2023.1.10"
 
 [[deps.MuladdMacro]]
 git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
@@ -829,12 +834,12 @@ version = "1.12.9"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.21+4"
+version = "0.3.23+2"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+0"
+version = "0.8.1+2"
 
 [[deps.OpenSSL]]
 deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"]
@@ -892,7 +897,7 @@ version = "2.7.1"
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.PoissonRandom]]
 deps = ["Random"]
@@ -993,7 +998,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[deps.Random]]
-deps = ["SHA", "Serialization"]
+deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.Random123]]
@@ -1174,6 +1179,7 @@ version = "1.1.1"
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+version = "1.10.0"
 
 [[deps.SparseDiffTools]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "Reexport", "Requires", "SciMLOperators", "Setfield", "SparseArrays", "StaticArrayInterface", "StaticArrays", "Tricks", "VertexSafeGraphs"]
@@ -1240,7 +1246,7 @@ version = "1.4.0"
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.StatsAPI]]
 deps = ["LinearAlgebra"]
@@ -1303,9 +1309,9 @@ deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[deps.SuiteSparse_jll]]
-deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.10.1+6"
+version = "7.2.1+1"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["DocStringExtensions"]
@@ -1422,7 +1428,7 @@ version = "0.2.0"
 deps = ["FFTW", "LinearAlgebra", "Parameters", "PrecompileTools", "QuantumOptics", "QuantumOpticsBase", "SparseArrays", "Strided", "UnsafeArrays"]
 path = ".."
 uuid = "c4555495-0e8d-488d-8e6a-965ecefe9055"
-version = "0.2.1"
+version = "0.2.4"
 
 [[deps.WignerSymbols]]
 deps = ["HalfIntegers", "LRUCache", "Primes", "RationalRoots"]
@@ -1439,7 +1445,7 @@ version = "0.4.9"
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.13+0"
+version = "1.2.13+1"
 
 [[deps.ZygoteRules]]
 deps = ["ChainRulesCore", "MacroTools"]
@@ -1450,14 +1456,14 @@ version = "0.2.3"
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.7.0+0"
+version = "5.8.0+1"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.48.0+0"
+version = "1.52.0+1"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+0"
+version = "17.4.0+2"

docs/src/tutorial.md

@@ -98,13 +98,13 @@ If we also want to know the number of photons in the waveguide state as a functi
 n = (ad*a) ⊗ identityoperator(bw)
 n_w = identityoperator(bc) ⊗ (create(bw)*destroy(bw))
 function exp_na_and_nw(time,psi)
-    (expect(n,psi),expect(n_w,psi))
+    (expect(n,psi),expect_waveguide(n_w,psi))
 end
 ψ_out,na,nw = waveguide_evolution(times,ψ_in,H,fout=exp_na_and_nw)
 nothing #hide
 ```
 
-Where `expect(n_w,psi)` calculates the expectation value of all times of the pulse at each timestep: $\mathrm{expect(n_w,psi)} = \bra{\psi} \sum_k  I \otimes w_k^\dagger w_k  \ket{\psi}$
+Where `expect_waveguide(n_w,psi)` calculates the expectation value of all times of the pulse at each timestep: $\mathrm{expect_waveguide(n_w,psi)} = \bra{\psi} \sum_k  I \otimes w_k^\dagger w_k  \ket{\psi}$
 
 This can be plotted as:
 

src/WaveguideInteraction.jl

@@ -55,6 +55,45 @@ end
 
 Base.:*(x::WaveguideInteraction{BL,BR},y::WaveguideInteraction{BL,BR}) where {BL,BR} = LazyProduct((x,y),x.factor*y.factor)
 
+
+
+
+function Base.:*(a::O1,b::O2) where {O1 <: WaveguideOperator,O2 <: WaveguideOperator}
+    WaveguideInteraction(basis(a),basis(a),complex(1.0),a,b,1)
+end
+
+function Base.:*(a::O1,b::O2) where {O1 <: WaveguideInteraction,O2 <: WaveguideOperator}
+    LazyProduct((a,b),a.factor*b.factor)
+end
+function Base.:*(a::O1,b::O2) where {O1 <: WaveguideOperator,O2 <: WaveguideInteraction}
+    LazyProduct((a,b),a.factor*b.factor)
+end
+
+function Base.:*(a::T1,b::T2) where {T1<: WaveguideOperatorT,T2<: WaveguideOperatorT}
+    @assert a.indices[1] == b.indices[1]
+    WaveguideInteraction(basis(a),basis(a),a.factor*b.factor,a.operators[1],b.operators[1],a.indices[1])
+end
+
+const TensorWaveguideInteraction = LazyTensor{B,B,F,V,T} where {B,F,V,T<:Tuple{Vararg{Union{WaveguideInteraction,WaveguideOperator,CavityWaveguideOperator}}}}
+const WaveguideSum = LazySum{B,B,F,T} where {B,F,T<:Tuple{Vararg{Union{WaveguideInteraction,WaveguideOperator,CavityWaveguideOperator}}}}
+
+function zerophoton_projector(psi::Ket)
+    zerophoton_projector(basis(psi))
+end
+
+function zerophoton_projector(basis::WaveguideBasis)
+    dm(zerophoton(basis))
+end
+
+function zerophoton_projector(basis::Basis)
+    identityoperator(basis)
+end
+
+function zerophoton_projector(b::CompositeBasis)
+    tensor([zerophoton_projector(b.bases[i]) for i in 1:length(b.bases)]...)
+end
+
+
 identityoperator(x::WaveguideInteraction) = identityoperator(x.basis_l)
 function identityoperator(::Type{T},::Type{ComplexF64}, b1::Basis, b2::Basis) where {T<:WaveguideInteraction}
     @assert b1==b2
@@ -297,55 +336,3 @@ function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,Np,idx1},bo
     @inbounds result[1+(idx1-1)*nsteps+timeindex] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex] 
     result
 end
-
-
-function Base.:*(a::O1,b::O2) where {O1 <: WaveguideOperator,O2 <: WaveguideOperator}
-    WaveguideInteraction(basis(a),basis(a),complex(1.0),a,b,1)
-end
-function Base.:*(a::T1,b::T2) where {T1<: WaveguideOperatorT,T2<: WaveguideOperatorT}
-    @assert a.indices[1] == b.indices[1]
-    WaveguideInteraction(basis(a),basis(a),a.factor*b.factor,a.operators[1],b.operators[1],a.indices[1])
-end
-
-const TensorWaveguideInteraction = LazyTensor{B,B,F,V,T} where {B,F,V,T<:Tuple{Vararg{Union{WaveguideInteraction,WaveguideOperator,CavityWaveguideOperator}}}}
-const WaveguideSum = LazySum{B,B,F,T} where {B,F,T<:Tuple{Vararg{Union{WaveguideInteraction,WaveguideOperator,CavityWaveguideOperator}}}}
-
-function zerophoton_projector(psi::Ket)
-    zerophoton_projector(basis(psi))
-end
-
-function zerophoton_projector(basis::WaveguideBasis)
-    dm(zerophoton(basis))
-end
-
-function zerophoton_projector(basis::Basis)
-    identityoperator(basis)
-end
-
-function zerophoton_projector(b::CompositeBasis)
-    tensor([zerophoton_projector(b.bases[i]) for i in 1:length(b.bases)]...)
-end
-
-function expect(a::T,psi::Ket) where T<:Union{WaveguideOperator,WaveguideInteraction,TensorWaveguideInteraction,CavityWaveguideOperator}
-    out = 0
-    tmp_ket = Ket(psi.basis)
-    projector = zerophoton_projector(psi)
-    set_waveguidetimeindex!(a,1)
-    mul!(tmp_ket,a,psi,a.factor,0)
-    zero_part = dot(psi.data,(projector*tmp_ket).data)
-    out += dot(psi.data,tmp_ket.data)
-    for i in 2:get_nsteps(basis(a))
-        set_waveguidetimeindex!(a,i)
-        mul!(tmp_ket,a,psi,a.factor,0)
-        out += (dot(psi.data,tmp_ket.data) - zero_part)
-    end
-    out
-end
-
-function expect(a::T,psi::Ket) where T<:Union{WaveguideSum}
-    out = 0
-    for (i,O) in enumerate(a.operators)
-        out +=  expect(O,psi)*a.factors[i]  
-    end
-    out
-end

src/WaveguideOperator.jl

@@ -71,6 +71,8 @@ Base.:+(a::LazySum{B1,B2},b::WaveguideOperator{B1,B2}) where {B1,B2} = LazySum(a
 Base.:-(a::WaveguideOperator{B1,B2},b::LazySum{B1,B2}) where {B1,B2}= LazySum(a) - LazySum(b)
 Base.:-(a::LazySum{B1,B2},b::WaveguideOperator{B1,B2}) where {B1,B2}= LazySum(a) - LazySum(b)
 
+Base.:*(a::WaveguideOperator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2} = LazyProduct(a) * b
+Base.:*(a::LazyProduct{B1,B2},b::WaveguideOperator{B1,B2}) where {B1,B2} = a * LazyProduct(b)
 """
     destroy(basis::WaveguideBasis{Np,1}) where {Np}
     destroy(basis::WaveguideBasis{Np,Nw},i::Int) where {Np,Nw}
@@ -178,6 +180,27 @@ function identityoperator(::Type{T},::Type{ComplexF64}, b1::Basis, b2::Basis) wh
     identityoperator(b1)
 end
 
+"""
+    expect_waveguide(O, psi,times)
+
+Returns the sum of expectation values of O over all times in times. This is usefull to extract information from waveguide states.
+For example the number of photons can be computed as: `expect_waveguide(create(bw)*destroy(bw), psi, times)`, where BW is the waveguidebasis. Mathematically this is equivalent to ``sum_n \\langle \\psi | O(t_n) | \\psi \\rangle``.
+For ``O(t_k) = w_k^\\dagger w_k`` acting on a single photon state ``|\\psi \\rangle = \\sum_k \\xi(t_k) w_k^\\dagger |0\\rangle``, one would therefore have: ``\\langle O \\rangle = \\sum_k |\\xi(t_k)|^2 = 1``, where the last follows from ``\\xi(t)`` being normalized.
+"""
+function expect_waveguide(O,psi)
+    psi_c = copy(psi)
+    expval = 0
+    ops = get_waveguide_operators(O)
+    for i in 1:get_nsteps(basis(O))
+        set_waveguidetimeindex!(ops,i)
+        mul!(psi_c, O, psi)
+        expval += dot(psi_c.data,psi.data)
+    end
+    return expval
+end
+
+
+
 """
     mul!(result::Ket{B1}, a::LazyTensor{B1,B2,F,I,T}, b::Ket{B2}, alpha, beta)
     

src/WaveguideQED.jl

@@ -8,7 +8,7 @@ import LinearAlgebra: axpy!, dot, mul!, rmul!,I
 import QuantumOpticsBase: create, dagger, destroy, expect, identityoperator, tensor,set_time!
 
 export TwoPhotonTimestepView,TwoWaveguideTimestepView,OnePhotonView,TwoPhotonView,TwoWaveguideView,
-    WaveguideBasis,zerophoton,onephoton,twophoton,view_waveguide,get_waveguidetimeindex,set_waveguidetimeindex!,get_dt,get_nsteps,get_waveguide_location,get_waveguide_basis,get_number_of_waveguides,get_waveguide_operators,
+    WaveguideBasis,zerophoton,onephoton,twophoton,view_waveguide,get_waveguidetimeindex,set_waveguidetimeindex!,get_dt,get_nsteps,get_waveguide_location,get_waveguide_basis,get_number_of_waveguides,get_waveguide_operators,expect_waveguide,
     WaveguideOperator,WaveguideDestroy,WaveguideCreate,
     CavityWaveguideAbsorption,CavityWaveguideEmission,emission,absorption,
     WaveguideInteraction,
@@ -17,7 +17,7 @@ export TwoPhotonTimestepView,TwoWaveguideTimestepView,OnePhotonView,TwoPhotonVie
     detect_single_click,detect_single_click!,LazyTensorKet,LazyTensorBra,LazySumKet,get_all_projectors,detect_double_click,detect_double_click!,Detector,
     plot_twophoton!,
     WaveguideTransform,effective_hamiltonian,fftket,
-    destroy,create,tensor,⊗,dagger,identityoperator,expect,set_time!
+    destroy,create,tensor,⊗,dagger,identityoperator,set_time!
 
 include("view.jl")
 include("basis.jl")

src/basis.jl

@@ -830,5 +830,3 @@ Returns [`WaveguideBasis`](@ref) from `CompositeBasis.bases`
 function get_waveguide_basis(basis::CompositeBasis)
     basis.bases[findall(x->typeof(x)<:WaveguideBasis,basis.bases)]
 end
-
-