Removed "should_upstream.jl" and added QUBase 0.4 as dependency

Author: mabuni1998

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.0"
 manifest_format = "2.0"
-project_hash = "223688a101f8a9d7186d2b2c1de1b4a7a04b14df"
+project_hash = "4dde56137bf92cdd4b71288bb369d001941db81e"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -780,15 +780,15 @@ version = "0.1.0"
 
 [[deps.QuantumOptics]]
 deps = ["Arpack", "DiffEqBase", "DiffEqCallbacks", "FFTW", "ForwardDiff", "IterativeSolvers", "LinearAlgebra", "LinearMaps", "OrdinaryDiffEq", "QuantumOpticsBase", "Random", "RecursiveArrayTools", "Reexport", "SparseArrays", "StochasticDiffEq", "WignerSymbols"]
-git-tree-sha1 = "81d7cf56eddb8b7c6b65c9da638d10bf2f10d7c2"
+git-tree-sha1 = "ff110fa403b1cd68b97527e42fb1333ab0b67969"
 uuid = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
-version = "1.0.9"
+version = "1.0.10"
 
 [[deps.QuantumOpticsBase]]
-deps = ["Adapt", "FFTW", "LRUCache", "LinearAlgebra", "QuantumInterface", "Random", "SparseArrays", "Strided", "UnsafeArrays"]
-git-tree-sha1 = "c0fe25d57ffe78ec8a56f831e1631d448e5973d0"
+deps = ["Adapt", "FFTW", "FillArrays", "LRUCache", "LinearAlgebra", "QuantumInterface", "Random", "SparseArrays", "Strided", "UnsafeArrays"]
+git-tree-sha1 = "0d916919aa637558cc8e146dfaffa60130219db7"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "0.3.12"
+version = "0.4.0"
 
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]

Project.toml

@@ -17,7 +17,7 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 Parameters = "0.12"
 PrecompileTools = "1"
 QuantumOptics = "1"
-QuantumOpticsBase = "0.3.9"
+QuantumOpticsBase = "0.4"
 Strided = "1,2"
 UnsafeArrays = "1"
-julia = "1.8"
+julia = "1.9"

README.md

@@ -1,12 +1,12 @@
 # WaveguideQED.jl
 <a href="https://qojulia.github.io/WaveguideQED.jl/dev/"><img src="https://img.shields.io/badge/docs-stable-blue.svg" alt="Documentation of latest stable version"></a> 
-<a href="https://codecov.io/gh/mabuni1998/WaveguideQED.jl"><img src="https://img.shields.io/codecov/c/gh/mabuni1998/WaveguideQED.jl?label=codecov" alt="Test coverage from codecov"></a>
+<a href="https://codecov.io/gh/qojulia/WaveguideQED.jl"><img src="https://img.shields.io/codecov/c/gh/qojulia/WaveguideQED.jl?label=codecov" alt="Test coverage from codecov"></a>
 
 
-A julia package for simulating quantum states of photon wavepackets using a discrete-time formalism [Phys. Rev. A 101, 042322](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.042322). Package works as an extension to [QuantumOptics.jl](https://qojulia.org/) where basises and operators from WaveguideQED.jl can be used together with operators and basises from QuantumOpics.jl. 
+A Julia package for simulating quantum states of photon wavepackets using a discrete-time formalism [Phys. Rev. A 101, 042322](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.042322). The package works as an extension to [QuantumOptics.jl](https://qojulia.org/) where bases and operators from WaveguideQED.jl can be used together with operators and bases from QuantumOpics.jl. 
 
 ### Example of usage:
-Define a waveguide basis, containing a two photon wavepacket for a time interval 0 to 20 with 0.2 timesteps:
+Define a waveguide basis, containing a two-photon wavepacket for a time interval from 0 to 20 with timesteps of 0.2:
 
 
 ```julia
@@ -32,7 +32,7 @@ wda = a ⊗ wd
 adw = ad ⊗ w
 ```
 
-Finally, we can define an initial twophoton gaussian wavepacket state with view_twophoton and zero photons in the cavity, an Hamiltonian, and simulate the evolution:
+Finally, we can define an initial two-photon Gaussian wavepacket state with view_twophoton and zero photons in the cavity, a Hamiltonian, and simulate the evolution:
 
 
 ```julia
@@ -44,7 +44,7 @@ H = im*sqrt(1/dt)*(adw-wda)
 ψ = waveguide_evolution(times, psi, H)
 ```
 
-Plotting the twophoton state is also simple:
+Plotting the two-photon state is also simple:
 
 
 ```julia

docs/src/detection.md

@@ -26,7 +26,7 @@ $$\begin{align*}
 &= \frac{1}{2} \left ( W_c^\dagger(\xi_a) W_c^\dagger(\xi_b) \ket{0}_c - W_d^\dagger(\xi_a) W_d^\dagger(\xi_b) \ket{0}_d + \int_{t_0}^{t_{end}} \mathrm{d}t \int_{t_0}^{t_{end}} \mathrm{d}t' \left [ \xi_a(t)\xi_b(t') - \xi_a(t')\xi_b(t) \right] \ket{1}_c\ket{1}_d \right)
 \end{align*}$$
 
-where we introduced $$W_{c/d}^\dagger(\xi_a) W_{c/d}^\dagger(\xi_b) \ket{0}_{c/d} = int_{t_0}^{t_{end}} \mathrm{d}t \int_{t_0}^{t_{end}} \mathrm{d}t' \xi_a(t)\xi_b(t') w_{c/d}^\dagger(t) w_{c/d}^\dagger(t') \ket{0}_{c/d}$$. $$W_{c/d}^\dagger(\xi_a) W_{c/d}^\dagger(\xi_b) \ket{0}_{c/d}$$ thus corresponds to both photons going into the same direction. It is also evident that if $$\xi_a(t)\xi_b(t') - \xi_a(t')\xi_b(t) \right = 0$$ then we will have no photons in waveguide c and d simultaneously. This condition is exactly fulfilled if the photon in waveguide a is indistinguishable from the photon in waveguide b. This also means that if the photons ARE distinguishable, we will start to see photon occurring in waveguide c and d simultaneously. All this and more can be simulated in the code, and in the next section, we walk through how to set the above example up in the code.
+where we introduced $$W_{c/d}^\dagger(\xi_a) W_{c/d}^\dagger(\xi_b) \ket{0}_{c/d} = int_{t_0}^{t_{end}} \mathrm{d}t \int_{t_0}^{t_{end}} \mathrm{d}t' \xi_a(t)\xi_b(t') w_{c/d}^\dagger(t) w_{c/d}^\dagger(t') \ket{0}_{c/d}$$. $$W_{c/d}^\dagger(\xi_a) W_{c/d}^\dagger(\xi_b) \ket{0}_{c/d}$$ thus corresponds to both photons going into the same direction. It is also evident that if $$\xi_a(t)\xi_b(t') - \xi_a(t')\xi_b(t) = 0$$ then we will have no photons in waveguide c and d simultaneously. This condition is exactly fulfilled if the photon in waveguide a is indistinguishable from the photon in waveguide b. This also means that if the photons ARE distinguishable, we will start to see photons occurring in waveguides c and d simultaneously. All this and more can be simulated in the code, and in the next section, we walk through how to set the above example up in the code.
 
 ## Beamsplitter and detection in WaveguideQED.jl
 
@@ -60,7 +60,7 @@ Dminus = Detector(wa/sqrt(2),-wb/sqrt(2))
 nothing #hide
 ``` 
 
-The [`Detector`](@ref) applies the first operator (`wa/sqrt(2)`) to the first `Ket` in LazyTensorKet (`waveguide_a`) and the second operator (L`$\pm$ wb/sqrt(2)`) to the second `Ket` in `LazyTensorKet` (waveguide_b). The probability of detecting a photon in the detectors can then be calculated by:
+The [`Detector`](@ref) applies the first operator (`wa/sqrt(2)`) to the first `Ket` in LazyTensorKet (`waveguide_a`) and the second operator ($$\pm $$ `wb/sqrt(2)`) to the second `Ket` in `LazyTensorKet` (waveguide_b). The probability of detecting a photon in the detectors can then be calculated by:
 
 ```@repl detection
 p_plus = Dplus * ψ_total

docs/src/example_lodahl.md

@@ -93,7 +93,7 @@ end
 nothing #hide
 ```
 
-Finally, we can plot the scattered wavefunctions, and we note that this matches Fig. 3 in Ref. ^[1]:
+Finally, we can plot the scattered wavefunctions, and we note that this matches Fig. 3 in Ref.[^1]:
 
 ```@example lodahl
 using PyPlot; #hide

docs/src/index.md

@@ -4,8 +4,8 @@ WaveguideQED.jl is a package for simulating continuous fockstates in waveguides.
 
 ## Dev docs
 Added functionalities:
-* [`WaveguideBasis`](@ref) for representing the waveguide space and the related generator functions: [`zerophoton`](@ref), [`onephoton`](@ref), and [`twophoton`](@ref). Also see [`OnePhotonView`](@ref), [`TwoPhotonView`](@ref), and [`plot_twophoton!`](@ref) for viewing the waveguide data for plotting. Note that [`WaveguideBasis`](@ref) can contain multiple waveguides.
-* [`WaveguideOperator`](@ref) which are specialized operators allowing efficient annihilation and creation operators at each time-bin in the waveguide. They are created by giving a basis to [`WaveguideQED.destroy`](@ref) and [`WaveguideQED.create`](@ref)
+* [`WaveguideBasis`](@ref) for representing the waveguide Hilbert space and the related functions for generating states in this Hilbert space: [`zerophoton`](@ref), [`onephoton`](@ref), and [`twophoton`](@ref). Also see [`OnePhotonView`](@ref), [`TwoPhotonView`](@ref), and [`plot_twophoton!`](@ref) for viewing the waveguide states and plotting them. Note that [`WaveguideBasis`](@ref) can contain multiple waveguides.
+* [`WaveguideOperator`](@ref) are specialized operators allowing efficient annihilation and creation operators at each time-bin in the waveguide. They are created by giving a basis to [`WaveguideQED.destroy`](@ref) and [`WaveguideQED.create`](@ref)
 * Since the interaction between the waveguide time-bin mode $k$ and cavity/emitter is given as: $a^\dagger w_k - a w_k^\dagger$ there are specially optimized functions for doing these operations called [`CavityWaveguideOperator`](@ref) which are created using a fockbasis and a waveguide basis and the functions [`emission`](@ref) and [`absorption`](@ref).
 * [`Detector`](@ref), [`LazyTensorKet`](@ref), and [`LazySumKet`](@ref), together with [`detect_single_click`](@ref) and [`detect_double_click`](@ref) allow one to do a beamsplitter interference and subsequent detection of photons coming from two waveguides. 
 

docs/src/references.md

@@ -1,16 +1,16 @@
 # Suggested readings
 
-## Theory and background
+### Theory and background
 The theory of time-bin continuous fockstates is introduced in [Heuck2020Photon-photonCavities](@cite) where it is used to derive the equations of motion for Waveguide QED systems. The numerical method in this library is heavily based on this approach, where instead of deriving the equations, the systems are solved by applying the operators themselves. The time-bin method is also used in [Heuck2020](@cite) and [Krastanov2022](@cite).
 
-## Similar or other approaches
+### Similar approaches
 
-The following is a list of approaches or methods that are trying to solve similar problems that can be treated with this library.
+The following is a list of approaches that are trying to solve problems that can also be treated with this library.
 * [HughesWQEDMPP2021](@cite) considers feedback in waveguide systems and uses a space-discretized waveguide picture with Monte Carlo trajectories
-* [Fischer2018](@cite) relates many approaches to solving WaveguideQED problems with each other and also introduces a framework that aims to deal with similar problems through master equations, photon counting and tracing.
+* [Fischer2018](@cite) relates many approaches to solving WaveguideQED problems with each other and also introduces a framework that aims to deal with similar problems through master equations, photon counting and tracing out the waveguide.
 * The SLH formalism introduced in [Kiilerich2019](@cite) and [Kiilerich2019](@cite) uses cascaded cavities to simulate quantum pulses. Further work also includes: [Yang2022](@cite) [Christiansen2023](@cite)
 
-## Papers where we reproduce results from
+### Papers where we reproduce results from
 * The theoretical results in [DynamicalPhotonLodahl2022](@cite) are reproduced in [Scattering on a two-level system](https://qojulia.github.io/WaveguideQED.jl/dev/example_lodahl/)
 
 # References

src/WaveguideOperator.jl

@@ -174,13 +174,28 @@ function mul!(result::Ket{B1}, a::LazyTensor{B1,B2,F,I,T}, b::Ket{B2}, alpha, be
     b_data = Base.ReshapedArray(b.data, QuantumOpticsBase._comp_size(basis(b)), ())
     result_data = Base.ReshapedArray(result.data, QuantumOpticsBase._comp_size(basis(result)), ())
 
-    tp_ops = (tuple(( (isa(op,DataOperator) ? op.data : op) for op in a.operators)...), a.indices)
-    iso_ops = QuantumOpticsBase._explicit_isometries(a.indices, a.basis_l, a.basis_r)
+    tp_ops = QuantumOpticsBase._tpops_tuple(a)
+    iso_ops = QuantumOpticsBase._explicit_isometries(eltype(a), a.indices, a.basis_l, a.basis_r)
 
     QuantumOpticsBase._tp_sum_matmul!(result_data, tp_ops, iso_ops, b_data, alpha * a.factor, beta)
     result
 end
 
+function QuantumOpticsBase._tpops_tuple(a::LazyTensor{B1,B2,F,I,T}; shift=0, op_transform=identity)  where {B1<:Basis,B2<:Basis, F,I,T<:Tuple{Vararg{AbstractOperator}}}
+    length(a.operators) == 0 == length(a.indices) && return ()
+    op_pairs = tuple((((isa(op,DataOperator) ? op_transform(op.data) : op_transform(op)), i + shift) for (op, i) in zip(a.operators, a.indices))...)
+
+    # Filter out identities:
+    # This induces a non-trivial cost only if _is_square_eye is not inferrable.
+    # This happens if we have Eyes that are not SquareEyes.
+    # This can happen if the user constructs LazyTensor operators including
+    # explicit identityoperator(b,b).
+    filtered = filter(p->!QuantumOpticsBase._is_square_eye(p[1]), op_pairs)
+    return filtered
+end
+QuantumOpticsBase._is_square_eye(a::AbstractOperator) = false
+
+
 #Called from _tp_sum_matmul!
 #Makes sure operator works on correct part of tensor.
 function QuantumOpticsBase.:_tp_matmul!(result, a::WaveguideOperator, loc::Integer, b, α::Number, β::Number)

src/WaveguideQED.jl

@@ -22,7 +22,7 @@ include("WaveguideOperator.jl")
 include("CavityWaveguideOperator.jl")
 include("WaveguideInteraction.jl")
 include("solver.jl")
-include("should_upstream.jl")
+#include("should_upstream.jl")
 include("detection.jl")
 include("plotting.jl")
 include("precompile.jl")

src/basis.jl

@@ -724,6 +724,7 @@ end
 - `norm::Bool=true`: normalize the resulting wavepacket.
 """
 twophoton(b::WaveguideBasis{T,Nw},idx::I,ξ::Matrix;norm=true) where {T,Nw,I<:Union{Vector{Int64},Tuple{Vararg{Int64}}}} = twophoton(b,idx[1],idx[2],ξ,norm=norm)
+twophoton(b::WaveguideBasis{1,Nw},args...;norm=true) where {Nw} = throw(ArgumentError("The input basis b only contains onephoton and I can't create twophotons. Create a new basis containing two-photons using WaveguideBasis(2,times)"))
 
 """
     twophoton(b::WaveguideBasis{T,Nw},idx::I,ξ::Function,times,args...;norm=true) where {T,Nw,I<:Union{Vector{Int64},Tuple{Vararg{Int64}}}} = twophoton(b,idx[1],idx[2],ξ,times,args...,norm=norm)