Performance using special structures

[mabuni1998]

@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:

WaveguideQED.jl/src/CavityWaveguideOperator.jl

Lines 50 to 77 in c609dc5

	mutable struct CavityWaveguideAbsorption{B1,B2} <: CavityWaveguideOperator{B1,B2}
	basis_l::B1
	basis_r::B2
	factor::ComplexF64
	op::AbstractOperator
	loc
	indexing::WaveguideIndexing
	function CavityWaveguideAbsorption(basis_l::B1,basis_r::B2,factor::ComplexF64,op::AbstractOperator,loc) where{B1,B2}
	new{B1,B2}(basis_l,basis_r,factor,op,loc,WaveguideIndexing(basis_l,loc))
	end
	end
	"""
	CavityWaveguideEmission{B1,B2} <: CavityWaveguideOperator{B1,B2}
	Structure for fast simultaneous cavity annihilation and waveguide creation operator
	"""
	mutable struct CavityWaveguideEmission{B1,B2} <: CavityWaveguideOperator{B1,B2}
	basis_l::B1
	basis_r::B2
	factor::ComplexF64
	op::AbstractOperator
	loc
	indexing::WaveguideIndexing
	function CavityWaveguideEmission(basis_l::B1,basis_r::B2,factor::ComplexF64,op::AbstractOperator,loc) where{B1,B2}
	new{B1,B2}(basis_l,basis_r,factor,op,loc,WaveguideIndexing(basis_l,loc))
	end
	end

and also for waveguide interactions $w_k^\dagger w_j$:

WaveguideQED.jl/src/WaveguideInteraction.jl

Lines 36 to 43 in c609dc5

	mutable struct WaveguideInteraction{B1,B2} <: AbstractOperator{B1,B2}
	basis_l::B1
	basis_r::B2
	factor::ComplexF64
	op1::AbstractOperator
	op2::AbstractOperator
	loc
	end

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.