qojulia/waveguideqed.jl

Library for simulating time binned photons in Waveguide Quantum Electrodynamics

waveguideqed.jl's Issues

Performance using special structures

@Krastanov Currently, there are a number of special operator structures that aim to increase performance. For example the operations $a^\dagger w$ and $a w^\dagger$ have the CavitWaveguideOperators: https://github.com/qojulia/WaveguideQED.jl/blob/c609dc544214b994585f05fe3159be518d28391f/src/CavityWaveguideOperator.jl#L50-L77

and also for waveguide interactions $w_k^\dagger w_j$: https://github.com/qojulia/WaveguideQED.jl/blob/c609dc544214b994585f05fe3159be518d28391f/src/WaveguideInteraction.jl#L36-L43

The waveguide interactions have not been optimized much as of now but just using CavityWaveguideOperators instead of a standard LazyTensor gives a speedup of like 1.5.

See the following test:

ξfun2(t1,t2,σ1,σ2,t0) = sqrt(2/σ1)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t1-t0)^2/σ1^2)*sqrt(2/σ2)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t2-t0)^2/σ2^2)

function special()
    times = 0:0.1:10
    be = FockBasis(1)
    bw = WaveguideBasis(2,2,times)
    bw_single = WaveguideBasis(1,2,times)
    b = bw ⊗ be
    a = embed(b,2,destroy(be))
    ad = embed(b,2,create(be))
    wdL = emission(bw,be,1)
    wL = absorption(bw,be,1)
    wdR = emission(bw,be,2) 
    wR = absorption(bw,be,2)

    wdL_in = create(bw,1) ⊗ identityoperator(be)
    wL_in = destroy(bw,1) ⊗ identityoperator(be)
    wdR_in = create(bw,2) ⊗ identityoperator(be)
    wR_in = destroy(bw,2) ⊗ identityoperator(be)

    V = 1
    H = im*sqrt(κ1/dt)*(wL-wdL) + im*sqrt(κ2/dt)*(wR-wdR) + V/dt *(LazyProduct(wdR_in,wL_in)+ LazyProduct(wdL_in,wR_in))
    psi_int  = twophoton(bw,1,ξfun2,1,1,5) ⊗ fockstate(be,1)
    psi_out = waveguide_evolution(times,psi_int,H)
end

function lazy()
    times = 0:0.1:10
    be = FockBasis(1)
    bw = WaveguideBasis(2,2,times)
    bw_single = WaveguideBasis(1,2,times)
    b = bw ⊗ be
    a = embed(b,2,destroy(be))
    ad = embed(b,2,create(be))

    wdL = LazyTensor(b,b,(1,2),(create(bw,1),destroy(be)))
    wL =  LazyTensor(b,b,(1,2),(destroy(bw,1),create(be)))
    wdR =  LazyTensor(b,b,(1,2),(create(bw,2),destroy(be)))
    wR =  LazyTensor(b,b,(1,2),(destroy(bw,2),create(be)))

    wdL_in = create(bw,1) ⊗ identityoperator(be)
    wL_in = destroy(bw,1) ⊗ identityoperator(be)
    wdR_in = create(bw,2) ⊗ identityoperator(be)
    wR_in = destroy(bw,2) ⊗ identityoperator(be)

    V = 1
    H = im*sqrt(κ1/dt)*(wL-wdL) + im*sqrt(κ2/dt)*(wR-wdR) + V/dt *(LazyProduct(wdR_in,wL_in)+ LazyProduct(wdL_in,wR_in))
    psi_int  = twophoton(bw,1,ξfun2,1,1,5) ⊗ fockstate(be,1)
    psi_out = waveguide_evolution(times,psi_int,H)
end

using BenchmarkTools

@btime special()
@btime lazy()

Giving:

930.162 ms (395625 allocations: 31.62 MiB)
1.498 s (219224 allocations: 28.54 MiB)

Right now, we are performing a lot of tests when doing the tensor product between a waveguide operator and any operator see for example:

https://github.com/qojulia/WaveguideQED.jl/blob/c609dc544214b994585f05fe3159be518d28391f/src/CavityWaveguideOperator.jl#L284-L295C4

This is to ensure that if the input is a cavity annihilation or creation operator, then the more efficient CavityWaveguideOperator structure is used. This is, however, very specialized, and maybe could cause problems? Instead, we could just create a LazyTensor per default (which is slower but requires much less code and instead we will be using code already tested in QuantumOpticsBase). Then if one wishes to enhance performance, you have to explicitly create the CavityWaveguideOperator via emission(bc,bw) and absorption(bc,bw). See:

WaveguideQED.jl/src/CavityWaveguideOperator.jl

Lines 131 to 221 in c609dc5

	function absorption(b1::WaveguideBasis{T,Nw},b2::FockBasis,idx) where {T,Nw}
	btotal = tensor(b1,b2)
	return CavityWaveguideAbsorption(btotal,btotal,complex(1.0),destroy(b1,idx),[1,2])
	end
	function absorption(b1::FockBasis,b2::WaveguideBasis{T,Nw},idx) where {T,Nw}
	btotal = tensor(b1,b2)
	return CavityWaveguideAbsorption(btotal,btotal,complex(1.0),destroy(b2,idx),[2,1])
	end
	function absorption(b1::WaveguideBasis{T},b2::FockBasis) where T
	btotal = tensor(b1,b2)
	return CavityWaveguideAbsorption(btotal,btotal,complex(1.0),destroy(b1),[1,2])
	end
	function absorption(b1::FockBasis,b2::WaveguideBasis{T}) where T
	btotal = tensor(b1,b2)
	return CavityWaveguideAbsorption(btotal,btotal,complex(1.0),destroy(b2),[2,1])
	end
	function absorption(a::T,b::Operator) where T<: WaveguideOperatorT
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideAbsorption(btotal,btotal,a.factor,a.operators[1],[a.indices[1],a.indices[1]+1])
	end
	function absorption(a::Operator,b::T) where T<: WaveguideOperatorT
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideAbsorption(btotal,btotal,b.factor,b.operators[1],[b.indices[1]+1,1])
	end
	function absorption(a::T,b::Operator) where T<: WaveguideOperator
	@assert isequal(b,create(basis(b)))
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideAbsorption(btotal,btotal,complex(1.0),a,[1,2])
	end
	function absorption(a::Operator,b::T) where T<: WaveguideOperator
	@assert isequal(a,create(basis(a)))
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideAbsorption(btotal,btotal,complex(1.0),b,[2,1])
	end
	function absorption(a::CompositeBasis,b::T,i::Int) where T<: Union{WaveguideOperator,WaveguideOperatorT}
	btotal = tensor(a,basis(b))
	CavityWaveguideAbsorption(btotal,btotal,complex(1.0),get_waveguide_operators(b)[1],[get_waveguide_location(basis(b))+length(a.shape),i])
	end
	function absorption(a::T,b::CompositeBasis,i::Int) where T<: Union{WaveguideOperator,WaveguideOperatorT}
	btotal = tensor(basis(a),b)
	CavityWaveguideAbsorption(btotal,btotal,complex(1.0),get_waveguide_operators(a)[1],[1,length(basis(a).shape)+i])
	end
	"""
	emission(b1::WaveguideBasis{T},b2::FockBasis) where T
	emission(b1::FockBasis,b2::WaveguideBasis{T}) where T
	Create [`CavityWaveguideEmission`](@ref) that applies `destroy(b::FockBasis)` on `FockBasis` and create(b::WaveguideBasis{T}) on [`WaveguideBasis{T}`](@ref).
	"""
	function emission(b1::WaveguideBasis{T,Nw},b2::FockBasis,idx) where {T,Nw}
	btotal = tensor(b1,b2)
	return CavityWaveguideEmission(btotal,btotal,complex(1.0),create(b1,idx),[1,2])
	end
	function emission(b1::FockBasis,b2::WaveguideBasis{T,Nw},idx) where {T,Nw}
	btotal = tensor(b1,b2)
	return CavityWaveguideEmission(btotal,btotal,complex(1.0),create(b2,idx),[2,1])
	end
	function emission(b1::WaveguideBasis{T,1},b2::FockBasis) where T
	btotal = tensor(b1,b2)
	return CavityWaveguideEmission(btotal,btotal,complex(1.0),create(b1),[1,2])
	end
	function emission(b1::FockBasis,b2::WaveguideBasis{T,1}) where T
	btotal = tensor(b1,b2)
	return CavityWaveguideEmission(btotal,btotal,complex(1.0),create(b2),[2,1])
	end
	function emission(a::T,b::Operator) where T<: WaveguideOperatorT
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideEmission(btotal,btotal,a.factor,a.operators[1],[a.indices[1],a.indices[1]+1])
	end
	function emission(a::Operator,b::T) where T<: WaveguideOperatorT
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideEmission(btotal,btotal,b.factor,b.operators[1],[b.indices[1]+1,1])
	end
	function emission(a::T,b::Operator) where T<: WaveguideOperator
	@assert isequal(b,destroy(basis(b)))
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideEmission(btotal,btotal,complex(1.0),a,[1,2])
	end
	function emission(a::Operator,b::T) where T<: WaveguideOperator
	@assert isequal(a,destroy(basis(a)))
	btotal = tensor(basis(a),basis(b))
	CavityWaveguideEmission(btotal,btotal,complex(1.0),b,[2,1])
	end
	function emission(a::CompositeBasis,b::T,i::Int) where T<: Union{WaveguideOperator,WaveguideOperatorT}
	btotal = tensor(a,basis(b))
	CavityWaveguideEmission(btotal,btotal,complex(1.0),get_waveguide_operators(b)[1],[get_waveguide_location(basis(b))+length(a.shape),i])
	end
	function emission(a::T,b::CompositeBasis,i::Int) where T<: Union{WaveguideOperator,WaveguideOperatorT}
	btotal = tensor(basis(a),b)
	CavityWaveguideEmission(btotal,btotal,complex(1.0),get_waveguide_operators(a)[1],[1,length(basis(a).shape)+i])
	end

What do you think? Should we automatically try to increase performance with the risk that we are reducing flexibility (and possible unexpected errors) or should we just create regular LazyTensor and leave the option for performance enhancement in the functions emission and absorption.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

a few explicit calls to `complex`

There are a few places where you call complex to convert floats to complex, but I am not sure why that would be necessary.

https://github.com/mabuni1998/CavityWaveguide/blob/e8bdfd40788268895f6b489644dd254f84cd8220/Examples/qo_singlephoton.jl#L43

Plus method for LazyTensor and LazySum

I want to implement tensor between QuantumOptics.jl operators and Waveguide operators so that it automatically returns LazyTensor:

function QuantumOpticsBase.:tensor(op1::AbstractOperator,op2::WaveguideOperator) 
    btotal = tensor(op1.basis,op2.basis)
    LazyTensor(btotal,btotal,[1,2],(op1,op2))
end

That is no problem. Then to make these LazyTensors work like normal operators, I want to define a + method that returns LazySum:

function Base.:+(op1::LazyTensor{B1,B2},op2::LazyTensor{B1,B2}) where {B1,B2}
    LazySum(op1,op2)
return

But in QuantumOpticsBase we already have:

function +(a::LazyTensor{B1,B2}, b::LazyTensor{B1,B2}) where {B1,B2}
    if length(a.indices) == 1 && a.indices == b.indices
        op = a.operators[1] * a.factor + b.operators[1] * b.factor
        return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,))
    end
    throw(ArgumentError("Addition of LazyTensor operators is only defined in case both operators act nontrivially on the same, single tensor factor."))
end

Is there a good way to extend this method? I'm thinking of using B1 and B2 to do some dispatch, but they are CompositeBasis, where one has a WaveguideBasis in it, but I don't know how exactly to address this. Finally I also want to define + method for LazyTensor and LazySum, but this is also defined like above, and probably solved the same way:

function Base.:+(op1::LazySum{B1,B2},op2::LazyTensor{B1,B2}) where {B1,B2}
    out = copy(op1)
    push!(out.operators,op2)
    push!(out.factors,1)    
end

Lazy operator

You already have something that works well with the typical dagger and tensor operation.

I would suggest adding a few tests that use it in this fashion and also adding a few tests that uses Lazy Tensor.

When you have checked that this works, it would be a good time to create a brand new type of operator, that is much simpler to work with and that can be used inside of Lazy Tensor. We can call it LazyDagger or something of that type. Let's use this issue to discuss why this might (or might not) be useful and how performant it would be and how to implement it.

Allocation

Hi Stefan,

Can you help me understand why the QuantumOpticsBase._tp_sum_matmul! makes allocations?

At some point it makes a call to: tmp = _tp_sum_get_tmp(op1..., b_data, :_tp_sum_matmul_tmp1)
As I understand it, this should avoid allocations if possible. However, when our code is run, it does not...

It eventually, calls:

function _tp_matmul_get_tmp(::Type{T}, shp::NTuple{N,Int}, sym) where {T,N}
    len = prod(shp)
    use_cache = lazytensor_use_cache()
    key = (sym, taskid(), UInt(len), T)
    if use_cache && Vector{T} <: LazyTensorCacheable
        cached = get!(lazytensor_cache, key) do
            Vector{T}(undef, len)
        end
        # Let's make sure the compiler knows we have the right type
        tmp = Vector{T}(cached)
    else
        tmp = Vector{T}(undef, len)
    end
    Base.ReshapedArray(tmp, shp, ())
end

where I'm guessing cached = get!(lazytensor_cache, key), is where it is supposed to use the cashed tmp vector instead of allocating a new one. However, running the above function I get allocations, and these are the main contributor to the 4GB of allocations during our solver. See the following:

julia>using BenchmarkTools
julia>using QuantumOptics
julia> T = ComplexF64
julia> shp = (2,122)
julia> @btime QuantumOpticsBase._tp_matmul_get_tmp(T,shp,:_tp_sum_matmul_tmp1)
  715.000 ns (3 allocations: 4.02 KiB)
2×122 reshape(::Vector{ComplexF64}, 2, 122) with eltype ComplexF64:
 4.24399e-314+1.35452e-314im       0.0+1.27683e-316im  5.0e-324+2.122e-314im  …  2.122e-314+2.12457e-314im  4.24399e-314+2.57819e-317im  4.4e-323+2.12557e-314im
 1.35448e-314+1.606e-321im    9.0e-323+6.36599e-314im       0.0+0.0im              5.0e-324+2.12457e-314im           0.0+6.0e-323im      7.4e-323+6.36599e-314im

Why does this happen? If we can make it work, then we should have solved the allocation problem.

Basis type

It probably makes sense to create our own basis type. Something like WaveguideBasis instead of FockBasis. When you have that you can also use WaveguideBasis(params) instead of get_cwbasis()

Example of a basis definition: FockBasis

Something like this might make sense (probably with different names for the struct and its properties):

struct WGBasis <: QuantumOptics.Basis
    N::Int # number of steps
    P::Int # maximum number of photons tracked
end

Then the view functions can be somewhat simplified as they will have access to the ket, and the ket itself will have knowledge of the basis, e.g. P=basis(ket).P to find the maximum number of photons.

Views and reshapes of sparsearrays

See the following example where viewing a reshaped sparse matrix doesn't work. Instead, allocating with getindex via b[:,1] allows one to extract the nonzero elements and thus enables a fast loop. This is basically what happens in my code, so if we can solve this problem, then it's done.

a = sparse(zeros(100))
a[1] = 1
b = reshape(a,10,10)
c = view(b,:,1)
d = b[:,1]
println(nonzeros(d))
println(nonzeros(c))

any reason for the `build` directory in the repository?

Most of the code seems to be in src, and build looks like a copy

Performance of mul!

@Krastanov
Update on performance after allocation has been fixed. We still use 1Kb of memory allocation for every multiplication. But this seems to come from indexing etc, and is constant independent of Hilbert size. See below:

using QuantumOptics
using CavityWaveguide
using LinearAlgebra
using BenchmarkTools

times = 0:0.1:10
bw = WaveguideBasis(1,times)
bc = FockBasis(1)
btotal = tensor(bc,bw)
a = sparse(destroy(bc))
ad = sparse(create(bc));
n = ad*a ⊗ identityoperator(bw)
w = destroy(bw)
wd = create(bw);
wda = a ⊗ wd
adw = ad ⊗ w
H = LazySum(param.δ*n,im*sqrt(param.γ/dt)*adw,-im*sqrt(param.γ/dt)*wda,param.x3/4*n*n,-param.x3/4*n)

ξfun(t,σ,t0) = complex(1/(σ*sqrt(2*pi))*exp(-2*log(2)*(t-t0)^2/σ^2))/sqrt(0.2820947917738782)
ξvec = ξfun.(times,1,1)
#Define initial state
ψ_cw = onephoton(bw,ξvec)
psi = fockstate(bc,0) ⊗  ψ_cw 
psi2 = fockstate(bc,0) ⊗  ψ_cw

@btime mul!(psi2,adw,psi,1,1)

Giving as output: 3.962 μs (30 allocations: 1.14 KiB).

With t = 0:0.001:10 we get: 61.100 μs (31 allocations: 1.16 KiB)

Obviously slower, but almost same number of allocations.

The number of allocations is still quite high, though, leading to 69Mb of allocations for a solve (with timestep 0:0.1:10):

julia> @btime ψ = waveguide_evolution(param.times, psi, H)
  224.924 ms (1845632 allocations: 69.59 MiB)

However, running with a smaller timestep (0:0.01:10):

julia> @btime ψ = waveguide_evolution(param.times, psi, H)
  9.946 s (15131058 allocations: 511.66 MiB)

We are now allocating 0.5GB. How come this now scales with the Hilbert space? Is it maybe just the allocation of psi_df used to calculate the differential? Anyways are we happy with the performance now? Or do we want to try and decrease the number of allocations further?

Question: How to extract information from superposition of photon fock states?

As the title indicates, I am curious about the method to extract information (especially phase) from states such as

ξfun(t, σ=1, t0=5) = complex(√(2/σ)*(log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2))
ψ_0 = zerophoton(bw)
ψ_1 = onephoton(bw, ξfun.(times))
ψ_p = ψ_0/√2 + ψ_1/√2
ψ_m = ψ_0/√2 - ψ_1/√2

I have tried expect(wd*w, ψ_p) which indeed gives me 0.5. However, I would like to distinguish between ψ_p and ψ_m. What's the best way to do that?

Thank you very much!

Authorship and license

You have some code comments referring to you authoring the package. A more common/expected method that also permits some automation of bibliography research is to add it in the Project.toml file [here](#authors = ["mabuni1998 [email protected] and contributors"]) and to have it in the license file

overly specific type signatures

In definitions of structs you want concrete types (avoids "boxing" and runtime dispatch or allocations) but in method signatures you do not need to be too specific as types are used only for dispatch there. For instance, this can simply be ::Vector so that people using things like dual numbers for autodifferentiation (instead of complex numbers) can still use it.

https://github.com/mabuni1998/CavityWaveguide/blob/e8bdfd40788268895f6b489644dd254f84cd8220/src/CavityWaveguideObject.jl#L7

Future version of smooth evolution of the coupling

In the solver, the current way to set the time index of the operators is:

WaveguideQED.jl/src/solver.jl

Lines 29 to 33 in 2e4d63c

	function get_hamiltonian(time,psi)
	tidx = min(round(Int,time/dt,RoundUp)+1,nsteps)
	set_waveguidetimeindex!(ops,tidx)
	return H
	end

Since the solver often uses timesteps much smaller than the timestep defined in the basis, the input is often between two-time bins. This is handled by dividing the input (float) with timestep (another float) and rounding. Thus we infer which time bin the input time corresponds to. This is faster(?) than looking into the whole array of times and seeing which two-time bins the input time lies between. However, dividing two floating numbers and turning them into an index seems like bad practice. Any inputs on another implementation @Krastanov ? (the current one works, although I recently had to put min() in because floating point errors made the last index go beyond the maximum time index).

Multiple waveguides

Currently, to model multiple waveguides, we use the specialized basis LeftRightCavityWaveguide:

https://github.com/mabuni1998/CavityWaveguide/blob/cf094c5e8fff16ba9c52084e2134065311f49abc/src/LeftRightWaveguide.jl#L4-L26

However, one could easily imagine that one instead extended the implementation of WaveguideBasis to something like (Np for number of photons, Nw for number of waveguides):

mutable struct WaveguideBasis{Np,Nw} <: QuantumOptics.Basis
    shape::Vector{Int}
    N::Int
    offset::Int
    nsteps::Int
    function WaveguideBasis(Np,Nw,times)
    if Nw==1    
      dim = 0
          for i in 1:Np
              dim = dim + length(times)^i - (length(times)^i-length(times))÷2
          end
          new{Np,Nw}([dim+1], dim, 0,length(times))
      end
    elseif Nw == 2
        dim = 0
        if Np==1
            dim = 1+2*length(times)
        elseif Np==2
            dim = 1+2*length(times)+ length(times)^2 + 2*(length(times)^2 - (length(times)^2-length(times))÷2)
        else
            error("Currently Nw larger than 2 is not supported in WaveguideBasis")
        end
        new{Np,Nw}([dim+1], dim, 0,length(times))
    end
end

Now the waveguide basis would have an additional type, that decides how many waveguides are encoded. This extension would also allow for more than two waveguides, where instead of using :left and :right to index one just uses 1, 2, 3 etc.

Obviously, the above implementation is clunky and not generalized to multiple waveguides, but that could be polished. However, before I make the above implementation, I would like inputs on design options. Should Nw, for example, be an optional argument rather than required since it might be complicated for new users to also specify this? And are there other suggestions?

Solving with matrix exponential

@Krastanov
As we have talked about, we can also solve the differential equation system by continuously applying the unitaryU_n = exp(-i delta t H) to every time-bin. Since our method is only precise up until dt, this should work fairly well and the potential performance gain is enormous. We will be probably sacrificing numerical stability and precision though.

To get a crasp of the performance gain I made the following example where we solve the evolution of a single-photon scattering on a cavity. Due to the low number of excitations, it is sufficient to expand U_n up to first order in t_n. Whether this will always hold is not clear to me.

using QuantumOptics
using WaveguideQED

δ = 0
γ = 1
times = 0:0.1:10
dt = times[2] - times[1]

#Create operators for two photons interacting with cavity
bc = FockBasis(1)
bw = WaveguideBasis(1,times)
a = destroy(bc)
ad = create(bc);
n = ad*a ⊗ identityoperator(bw)
#$w†a and a†w efficient implementation$
wda = emission(bc,bw)
adw = absorption(bc,bw)

#Standard Hamiltonain
H = δ*n + im*sqrt(γ/dt)*(adw-wda)

#Unitary
U = identityoperator(n) - im*dt*(δ*n + im*sqrt(γ/dt)*(adw-wda)) - 1/2*dt^2*( γ/dt*(adw*wda)  )

#Define unitary solver (veeeeeery simple)
function evolve_unitary(psi::Ket,U,times)
    out = copy(psi)
    tmp = copy(psi)
    for i in eachindex(times)
        set_waveguidetimeindex!(U,i)
        mul!(out,U,tmp,1,0)
        tmp.data .= out.data
    end
    out

end

#Define input onephoton state shape
σ = 1
t0 = 5
ξfun(t,σ,t0) = complex(sqrt(2/σ)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2))
ξvec = ξfun.(times,σ,t0)

ψ_cw = onephoton(bw,ξvec)
psi = fockstate(bc,0) ⊗  ψ_cw 

#Solve
@btime ψ = waveguide_evolution(times, psi, H)
@btime ψ = evolve_unitary(psi,U,times)

This gives:

28.266 ms (134309 allocations: 5.53 MiB)
444.100 μs (2832 allocations: 120.38 KiB)

We are gaining like 60 orders of magnitude in performance, and this will be even more apparent for larger systems! To allow for an automatic expansion of the Hamiltonian, you mentioned some tricks we could use. How should we go about constructing the unitary automatically and symbolically?

Activating CavityWaveguide environment

I want to make sure the Project.toml in the hand-in has all the correct dependencies, but since I run the code from the directory above, when I do ] add, I address the environment in that folder. When I cd into the CavityWaveguide directory and activate I get:

julia> ] activate
Activating project at `C:\Users\mabun\.julia\environments\v1.8`
(@v1.8) pkg>

So how do I activate CavityWaveguide as a project?

Typos in README

Thank you very much for the package!

There are some typos in the code demo in README.

First, the second line below should be wd instead of w

# original
w =  destroy(bw)
w =  create(bw)

Second, the first line below should be TwoPhotonView instead of view_twophoton

# original
ψ_double = view_twophoton(ψ);
using PyPlot
fig,ax = subplots(1,1,figsize=(9,4.5))
plot_twophoton!(ax,ψ_double,times)

Finally, would you offer a version of the plot using the Plots package, as the PyPlot package causes some extra trouble, especially in the non-interactive scenario? One of the potential solutions for the README example might be:

using Plots
twophotonstate = ψ_double
cnt1 = contourf(times, times, abs.(twophotonstate).^2, levels=10, color=:turbo, aspect_ratio=:equal, size=(400,400))

I would be very happy to submit a Pull Request for these improvements if a PR is encouraged.

Efficient wda and adw operator

@Krastanov
In the newest version, I implemented a new class of operator called CavityWaveguideOperator inspired by our old approach where the action of $a^\dagger w$ and $w^\dagger a$ was combined. Since this is almost always the only interaction with the waveguide, it's worthwhile to have a special operator which is efficient.

See the following example where I compare the "old" operator $a \otimes wd$ with the "emission" operator which is the new and more efficient implementation. This is a huge performance gain! So large that the limiting factor is no longer the cavity waveguide interaction, leading to a smaller performance gain when solving the whole system.

using QuantumOptics
using CavityWaveguide
using LinearAlgebra
using BenchmarkTools
include("singlecavity_simulations.jl")

param = parameters()

#Set detuning:
param.δ = 0
#Set non linearity:
param.x3 = 0

param.times = 0:0.1:10
dt = param.times[2] -param.times[1]
bw = WaveguideBasis(2,param.times)
bc = FockBasis(1)
a = destroy(bc)
ad = create(bc);
n_l = ad*a ⊗ identityoperator(bw)
w = destroy(bw)
wd = create(bw);

wda_new = emission(bc,bw)
adw_new = absorption(bc,bw)
H_new = param.δ*n_l+im*sqrt(param.γ/dt)*(adw_l-wda_l)+param.x3/4*n_l*n_l-param.x3/4*n_l

wda_old = a ⊗ wd
adw_old = ad ⊗ w
H_old = param.δ*n_f+im*sqrt(param.γ/dt)*(adw_f-wda_f)+param.x3/4*n_f*n_f-param.x3/4*n_f


ξfun(t,σ,t0) = complex(1/(σ*sqrt(2*pi))*exp(-2*log(2)*(t-t0)^2/σ^2))/sqrt(0.2820947917738782)
ξvec = ξfun.(param.times,1,1) * transpose(ξfun.(param.times,1,1))

ψ_cw = twophoton(bw,ξvec)

psi = fockstate(bc,0) ⊗ ψ_cw
out_new = copy(psi_l)
out_old = copy(psi_l)

psi_l.basis.bases[2].timeindex = 100

@btime mul!(out_new,wda_new,psi,1.0,0.0)
@btime mul!(out_old,wda_old,psi,1.0,0.0)
@btime mul!(out_new,adw_new,psi,1.0,0.0)
@btime mul!(out_old,adw_old,psi,1.0,0.0)

@btime ψ = waveguide_evolution(param.times, psi, H_new)
@btime ψ = waveguide_evolution(param.times, psi, H_old)

Giving:

5.214 μs (5 allocations: 304 bytes)
43.000 μs (13 allocations: 496 bytes)
4.586 μs (5 allocations: 304 bytes)
42.200 μs (13 allocations: 496 bytes)
1.431 s (114641 allocations: 12.06 MiB)
2.857 s (458243 allocations: 20.21 MiB)

The implementation is quite simple. I still run the underlying waveguide_mul! which use the efficient view, we implemented yesterday, but the "target" or output is shifted so that that the calculated output is saved in a fashion that mimicks an annihilation. An example is the following. Notice that the two views in waveguide_mul! is shifted, so that the first reads "i" and the second "i-1", which in this case corresponds to absorption or $a^\dagger w$ (creating a photon in cavity and taking it from waveguide):

function waveguide_mul_last!(result, a::CavityWaveguideAbsorption, b, alpha::Number, beta::Number)
    b_data = reshape(b, :, length(basis(a).bases[end]))
    result_data = reshape(result,  :, length(basis(a).bases[end]))
    for i in 2:size(b_data,1)
        waveguide_mul!(view(result_data,i,:),a.op,view(b_data,i-1,:),sqrt(i-1)*alpha*a.factor,beta)
    end
    rmul!(view(result_data,1,:),beta)
end

Corrospondingly for $w^\dagger a$ we have "i+1":

function waveguide_mul_last!(result, a::CavityWaveguideEmission, b, alpha::Number, beta::Number)
    b_data = reshape(b, :, length(basis(a).bases[end]))
    result_data = reshape(result,  :, length(basis(a).bases[end]))
    for i in 1:size(b_data,1)-1
        waveguide_mul!(view(result_data,i,:),a.op,view(b_data,i+1,:),sqrt(i)*alpha*a.factor,beta)
    end
    rmul!(view(result_data,size(b_data,1),:),beta)
end

Handling input/output waveguide in view

I have now implemented a basis for containing twophoton states in input and output waveguides. The new data structure can be viewed via the function:

https://github.com/mabuni1998/CavityWaveguide/blob/36d18ef36dc229c49a820f470cea5ed4aa910097/src/view.jl#L120-L140

Here, the argument type determines which part of the state you view. For example, type=:input returns a view of the two-photon input state. Likewise, type=:output returns a view of the two-photon output state, and type=:inputoutput gives the one-photon in input one photon in output state. type=:singlewaveguide is the standard used when you have only a single waveguide. I know that I could define TwoPhotonView{type} and then have a specialized function for each of the scenarios but for now, I'm just using the if statements here. I'm working on documentation, and the purpose of this issue is not to discuss the rest of the implementation but to discuss specifically how the view is made more friendly to use (if possible). For an example where I use the view, see:
https://github.com/mabuni1998/CavityWaveguide/blob/36d18ef36dc229c49a820f470cea5ed4aa910097/Examples/lodahl_fig2%20IO.jl#L49

Is this the typical way of having options in a function? Sometimes in python you would use string names instead so that the type would be: type="input" and you would dispatch from here. Is there a typical way of doing it in Julia other than the way I did it @Krastanov ?

How to measure the coincidence expectation value

When doing simulations involving HOM dip, I am wondering if we can compute the expectation value for coincidence:

expect(wd1*w1 ⊗ wd2*w2, ψ)

Where ψ is a state with only 2 waveguides, basis defined by WaveguideBasis(2, 2, times) with proper ladder operators.

Currently, Julia throws the error for this expression:

ERROR: LoadError: BoundsError: attempt to access Int64 at index [2]
Stacktrace:
 [1] getindex(x::Int64, i::Int64)
   @ Base ./number.jl:98
 [2] tensor(a::WaveguideInteraction{WaveguideBasis{2, 2}, WaveguideBasis{2, 2}}, b::WaveguideCreate{WaveguideBasis{2, 2}, WaveguideBasis{2, 2}, 2, 2})
   @ WaveguideQED ~/.julia/packages/WaveguideQED/yljRX/src/WaveguideInteraction.jl:104

Is there a way to properly compute the expectation value?

Thank you very much!

Representation of two-photon state

@Krastanov, at some point, there will be performance to gain by only representing the "upper half" of a two-photon state. Right now, we are only changing and interacting with the upper half anyway, but we are using memory to store the lower half. Some indexing will have to be changed, but it's not a big task that will save us N^2/2 of memory, which I guess is somewhat significant. I'm putting it here to have it as a possible task in the near future.

some allocations in `timestep!`

I suspect this one can be made allocation free and faster. It would be a convenient target as it is easy to benchmark independently of the rest.

https://github.com/mabuni1998/CavityWaveguide/blob/e8bdfd40788268895f6b489644dd254f84cd8220/src/CavityWaveguideObject.jl#L101

Order of hilbert spaces impact performance

@Krastanov See the following example where the order of the two hilbert spaces impact the performance of mul!

using QuantumOptics
using CavityWaveguide
using LinearAlgebra
using BenchmarkTools
include("singlecavity_simulations.jl")

param = parameters()
#Set detuning:
param.δ = 0
#Set non linearity:
param.x3 = 0

param.times = 0:0.1:10
dt = param.times[2] -param.times[1]
bw = WaveguideBasis(2,param.times)
bc = FockBasis(1)
a = destroy(bc)
ad = create(bc);
n_l = ad*a ⊗ identityoperator(bw)
w = destroy(bw)
wd = create(bw);

wda_l = a ⊗ wd
adw_l = ad ⊗ w
H_l = param.δ*n_l+im*sqrt(param.γ/dt)*(adw_l-wda_l)+param.x3/4*n_l*n_l-param.x3/4*n_l

n_f = identityoperator(bw) ⊗ (ad*a)
wda_f = wd ⊗ a
adw_f = w ⊗ ad
H_f = param.δ*n_f+im*sqrt(param.γ/dt)*(adw_f-wda_f)+param.x3/4*n_f*n_f-param.x3/4*n_f


ξfun(t,σ,t0) = complex(1/(σ*sqrt(2*pi))*exp(-2*log(2)*(t-t0)^2/σ^2))/sqrt(0.2820947917738782)
ξvec = ξfun.(param.times,1,1) * transpose(ξfun.(param.times,1,1))
#Define initial state
#ψ_cw = onephoton(bw,ξfun,param.times,param.σ,param.t0)
ψ_cw = twophoton(bw,ξvec)

psi_l = fockstate(bc,0) ⊗ ψ_cw
out_l = copy(psi_l)

psi_f = ψ_cw ⊗ fockstate(bc,0)
out_f = copy(psi_f)

psi.basis.bases[2].timeindex = 100

@btime mul!(out_l,wda_l,psi_l,1.0,1.0)
@btime mul!(out_f,wda_f,psi_f,1.0,1.0)
@btime ψ = waveguide_evolution(param.times, psi_l, H_l)
@btime ψ = waveguide_evolution(param.times, psi_f, H_f)
45.300 μs (13 allocations: 496 bytes)
15.400 μs (16 allocations: 624 bytes)
3.090 s (458243 allocations: 20.21 MiB)
1.818 s (572775 allocations: 24.87 MiB)

I guess it has to do with how the array is stored in memory. When the waveguide Hilbert space is the first, then the waveguide array is stored contiguously in memory, while when the cavity Hilbert space is first, the waveguide Hilbert space is split to every second element in memory. Do you have other explanations?