Request For Comment: Arithemtic between Operators and LazyOperators (#86)

Extend some of the current operation with LazyOperators to allow for smoother interface with e.g. LazyTensor. As an example +(a::LazyTensor,LazyTensor) now returns LazySum.

Author: mabuni1998

src/operators_lazyproduct.jl

@@ -11,7 +11,7 @@ complex factor is stored in the `factor` field which allows for fast
 multiplication with numbers.
 """
 
-mutable struct LazyProduct{BL,BR,F,T,KTL,BTR} <: AbstractOperator{BL,BR}
+mutable struct LazyProduct{BL,BR,F,T,KTL,BTR} <: LazyOperator{BL,BR}
     basis_l::BL
     basis_r::BR
     factor::F
@@ -74,6 +74,16 @@ permutesystems(op::LazyProduct, perm::Vector{Int}) = LazyProduct(([permutesystem
 identityoperator(::Type{LazyProduct}, ::Type{S}, b1::Basis, b2::Basis) where S<:Number = LazyProduct(identityoperator(S, b1, b2))
 
 
+
+function tensor(a::Operator{B1,B2},b::LazyProduct{B3, B4, F, T, KTL, BTR}) where {B1,B2,B3,B4, F, T, KTL, BTR}
+    ops = ([(i == 1 ? a : identityoperator(a.basis_r) ) ⊗ op for (i,op) in enumerate(b.operators)]...,)
+    LazyProduct(ops,b.factor)
+end
+function tensor(a::LazyProduct{B1, B2, F, T, KTL, BTR},b::Operator{B3,B4}) where {B1,B2,B3,B4, F, T, KTL, BTR}
+    ops = ([op ⊗ (i == length(a.operators) ? b : identityoperator(a.basis_l) ) for (i,op) in enumerate(a.operators)]...,)
+    LazyProduct(ops,a.factor)
+end
+
 function mul!(result::Ket{B1},a::LazyProduct{B1,B2},b::Ket{B2},alpha,beta) where {B1,B2}
     if length(a.operators)==1
         mul!(result,a.operators[1],b,a.factor*alpha,beta)

src/operators_lazysum.jl

@@ -9,6 +9,11 @@ function _check_bases(basis_l, basis_r, operators)
     end
 end
 
+"""
+Abstract class for all Lazy type operators ([`LazySum`](@ref), [`LazyProduct`](@ref), and [`LazyTensor`](@ref))
+"""
+abstract type LazyOperator{BL,BR} <: AbstractOperator{BL,BR} end
+
 """
     LazySum([Tf,] [factors,] operators)
     LazySum([Tf,] basis_l, basis_r, [factors,] [operators])
@@ -28,7 +33,7 @@ of operator-state operations, such as simulating time evolution. A `Vector` can
 reduce compile-time overhead when doing arithmetic on `LazySum`s, such as
 summing many `LazySum`s together.
 """
-mutable struct LazySum{BL,BR,F,T} <: AbstractOperator{BL,BR}
+mutable struct LazySum{BL,BR,F,T} <: LazyOperator{BL,BR}
     basis_l::BL
     basis_r::BR
     factors::F
@@ -57,6 +62,7 @@ end
 LazySum(::Type{Tf}, operators::AbstractOperator...) where Tf = LazySum(ones(Tf, length(operators)), (operators...,))
 LazySum(operators::AbstractOperator...) = LazySum(mapreduce(eltype, promote_type, operators), operators...)
 LazySum() = throw(ArgumentError("LazySum needs a basis, or at least one operator!"))
+LazySum(x::LazySum) = x
 
 Base.copy(x::LazySum) = @samebases LazySum(x.basis_l, x.basis_r, copy(x.factors), copy.(x.operators))
 Base.eltype(x::LazySum) = mapreduce(eltype, promote_type, x.operators; init=eltype(x.factors))
@@ -81,8 +87,9 @@ function +(a::LazySum{B1,B2}, b::LazySum{B1,B2}) where {B1,B2}
     @samebases LazySum(a.basis_l, a.basis_r, factors, ops)
 end
 +(a::LazySum{B1,B2}, b::LazySum{B3,B4}) where {B1,B2,B3,B4} = throw(IncompatibleBases())
-+(a::LazySum, b::AbstractOperator) = a + LazySum(b)
-+(a::AbstractOperator, b::LazySum) = LazySum(a) + b
++(a::LazyOperator, b::AbstractOperator) = LazySum(a) + LazySum(b)
++(a::AbstractOperator, b::LazyOperator) = LazySum(a) + LazySum(b)
++(a::O1, b::O2) where {O1<:LazyOperator,O2<:LazyOperator} = LazySum(a) + LazySum(b)
 
 -(a::LazySum) = @samebases LazySum(a.basis_l, a.basis_r, -a.factors, a.operators)
 function -(a::LazySum{B1,B2}, b::LazySum{B1,B2}) where {B1,B2}
@@ -92,8 +99,10 @@ function -(a::LazySum{B1,B2}, b::LazySum{B1,B2}) where {B1,B2}
     @samebases LazySum(a.basis_l, a.basis_r, factors, ops)
 end
 -(a::LazySum{B1,B2}, b::LazySum{B3,B4}) where {B1,B2,B3,B4} = throw(IncompatibleBases())
--(a::LazySum, b::AbstractOperator) = a - LazySum(b)
--(a::AbstractOperator, b::LazySum) = LazySum(a) - b
+-(a::LazyOperator, b::AbstractOperator) = LazySum(a) - LazySum(b)
+-(a::AbstractOperator, b::LazyOperator) = LazySum(a) - LazySum(b)
+-(a::O1, b::O2) where {O1<:LazyOperator,O2<:LazyOperator} = LazySum(a) - LazySum(b)
+
 
 function *(a::LazySum, b::Number)
     factors = b*a.factors
@@ -106,6 +115,19 @@ function /(a::LazySum, b::Number)
     @samebases LazySum(a.basis_l, a.basis_r, factors, a.operators)
 end
 
+function tensor(a::Operator,b::LazySum)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    ops = ([a ⊗ op for op in b.operators]...,)
+    LazySum(btotal_l,btotal_r,b.factors,ops)
+end
+function tensor(a::LazySum,b::Operator)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    ops = ([op ⊗ b for op in a.operators]...,)
+    LazySum(btotal_l,btotal_r,a.factors,ops)
+end
+
 function dagger(op::LazySum)
     ops = dagger.(op.operators)
     @samebases LazySum(op.basis_r, op.basis_l, conj.(op.factors), ops)

src/operators_lazytensor.jl

@@ -11,7 +11,7 @@ must be sorted.
 Additionally, a factor is stored in the `factor` field which allows for fast
 multiplication with numbers.
 """
-mutable struct LazyTensor{BL,BR,F,I,T} <: AbstractOperator{BL,BR}
+mutable struct LazyTensor{BL,BR,F,I,T} <: LazyOperator{BL,BR}
     basis_l::BL
     basis_r::BR
     factor::F
@@ -96,21 +96,37 @@ isequal(x::LazyTensor, y::LazyTensor) = samebases(x,y) && isequal(x.indices, y.i
 # Arithmetic operations
 -(a::LazyTensor) = LazyTensor(a, -a.factor)
 
-function +(a::LazyTensor{B1,B2}, b::LazyTensor{B1,B2}) where {B1,B2}
-    if length(a.indices) == 1 && a.indices == b.indices
+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
         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."))
+    LazySum(a) + LazySum(b)
 end
-
-function -(a::LazyTensor{B1,B2}, b::LazyTensor{B1,B2}) where {B1,B2}
-    if length(a.indices) == 1 && a.indices == b.indices
+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
         return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,))
     end
-    throw(ArgumentError("Subtraction of LazyTensor operators is only defined in case both operators act nontrivially on the same, single tensor factor."))
+    LazySum(a) - LazySum(b)
+end
+
+function tensor(a::LazyTensor{B1,B2},b::AbstractOperator{B3,B4}) where {B1,B2,B3,B4}
+    if B3 <: CompositeBasis || B4 <: CompositeBasis
+        throw(ArgumentError("tensor(a::LazyTensor{B1,B2},b::AbstractOperator{B3,B4}) is not implemented for B3 or B4 being CompositeBasis unless b is identityoperator "))
+    else
+        a ⊗ LazyTensor(b.basis_l,b.basis_r,[1],(b,),1)
+    end
 end
+function tensor(a::AbstractOperator{B1,B2},b::LazyTensor{B3,B4})  where {B1,B2,B3,B4}
+    if B1 <: CompositeBasis || B2 <: CompositeBasis
+        throw(ArgumentError("tensor(a::AbstractOperator{B1,B2},b::LazyTensor{B3,B4}) is not implemented for B1 or B2 being CompositeBasis unless b is identityoperator "))
+    else
+        LazyTensor(a.basis_l,a.basis_r,[1],(a,),1) ⊗ b
+    end
+end
+
 
 function *(a::LazyTensor{B1,B2}, b::LazyTensor{B2,B3}) where {B1,B2,B3}
     indices = sort(union(a.indices, b.indices))

test/test_abstractdata.jl

@@ -413,8 +413,6 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 
 # Addition
-@test_throws ArgumentError op1 + op2
-@test_throws ArgumentError op1 - op2
 @test D(-op1_, -op1, 1e-12)
 
 # Test multiplication
@@ -448,7 +446,7 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 
 # Addition
-@test_throws ArgumentError op1 + op2
+@test D(2.1*op1 + 0.3*op2, 2.1*op1_+0.3*op2_)
 @test D(-op1_, -op1)
 
 # Test multiplication

test/test_operators_lazyproduct.jl

@@ -50,21 +50,41 @@ op1b = randoperator(b_r, b_l)
 op2a = randoperator(b_l, b_r)
 op2b = randoperator(b_r, b_l)
 op3a = randoperator(b_l, b_l)
+op3b = randoperator(b_r, b_r)
+
 op1 = LazyProduct([op1a, sparse(op1b)])*0.1
 op1_ = 0.1*(op1a*op1b)
 op2 = LazyProduct([sparse(op2a), op2b], 0.3)
 op2_ = 0.3*(op2a*op2b)
 op3 = LazyProduct(op3a)
 op3_ = op3a
+op4 = LazyProduct(op3b)
+op4_ = op3b
 
 x1 = Ket(b_l, rand(ComplexF64, length(b_l)))
 x2 = Ket(b_l, rand(ComplexF64, length(b_l)))
 xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 
-# Addition
-@test_throws ArgumentError op1 + op2
+# Test Addition
+@test_throws QuantumOpticsBase.IncompatibleBases op1 + dagger(op4)
 @test 1e-14 > D(-op1_, -op1)
+@test 1e-14 > D(op1+op2, op1_+op2_)
+@test 1e-14 > D(op1+op2_, op1_+op2_)
+@test 1e-14 > D(op1_+op2, op1_+op2_)
+#Check for unallowed addition:
+@test_throws QuantumOpticsBase.IncompatibleBases LazyProduct([op1a, sparse(op1a)])+LazyProduct([sparse(op2b), op2b], 0.3)
+
+# Test Subtraction
+@test_throws QuantumOpticsBase.IncompatibleBases op1 - dagger(op4)
+@test 1e-14 > D(op1 - op2, op1_ - op2_)
+@test 1e-14 > D(op1 - op2_, op1_ - op2_)
+@test 1e-14 > D(op1_ - op2, op1_ - op2_)
+@test 1e-14 > D(op1 + (-op2), op1_ - op2_)
+@test 1e-14 > D(op1 + (-1*op2), op1_ - op2_)
+#Check for unallowed subtraction:
+@test_throws QuantumOpticsBase.IncompatibleBases LazyProduct([op1a, sparse(op1a)])-LazyProduct([sparse(op2b), op2b], 0.3)
+
 
 # Test multiplication
 @test_throws DimensionMismatch op1a*op1a
@@ -78,6 +98,14 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 # Test division
 @test 1e-14 > D(op1/7, op1_/7)
 
+#Test Tensor product
+op = op1a ⊗ op1
+op_ = op1a ⊗ op1_
+@test 1e-11 > D(op,op_)
+op = op1 ⊗ op2a
+op_ = op1_ ⊗ op2a 
+@test 1e-11 > D(op,op_)
+
 # Test identityoperator
 Idense = identityoperator(DenseOpType, b_l)
 id = identityoperator(LazyProduct, b_l)

test/test_operators_lazysum.jl

@@ -134,6 +134,29 @@ ops_tuple = (symmetric_op,hermitian_op,adjoint_op)
 # Test division
 @test 1e-14 > D(op1/7, op1_/7)
 
+#Test Tensor
+op1_tensor = op1a ⊗ op1
+op2_tensor = op1 ⊗ op1a
+op1_tensor_ = op1a ⊗ op1_
+op2_tensor_ = op1_ ⊗ op1a
+@test 1e-14 > D(op1_tensor,op1_tensor_)
+@test 1e-14 > D(op1_tensor_,op1_tensor)
+
+#Test addition of with and of LazyTensor and LazyProduct
+subop1 = randoperator(b1a, b1b)
+subop2 = randoperator(b2a, b2b)
+subop3 = randoperator(b3a, b3b)
+I1 = dense(identityoperator(b1a, b1b))
+I2 = dense(identityoperator(b2a, b2b))
+I3 = dense(identityoperator(b3a, b3b))
+op1_LT = LazyTensor(b_l, b_r, [1, 3], (subop1, sparse(subop3)), 0.1)
+op1_LT_ = 0.1*subop1 ⊗ I2 ⊗ subop3
+op1_LP  = LazyProduct([op1a])*0.1
+op1_LP_ = op1a*0.1
+@test 1e-14 > D(op1_LT+op1_LP,op1_LT_+op1_LP_)
+@test 1e-14 > D(op1_LT+op1,op1_LT_+op1_)
+@test 1e-14 > D(op1_LP+op1,op1_LP_+op1_)
+
 # Test tuples vs. vectors
 @test (op1+op1).operators isa Tuple
 @test (op1+op2).operators isa Tuple

test/test_operators_lazytensor.jl

@@ -10,7 +10,6 @@ mutable struct test_lazytensor{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR}
 end
 Base.eltype(::test_lazytensor) = ComplexF64
 
-@testset "operators-lazytensor" begin
 
 Random.seed!(0)
 
@@ -102,21 +101,51 @@ op2 = LazyTensor(b_l, b_r, [2, 3], (sparse(subop2), subop3), 0.7)
 op2_ = 0.7*I1 ⊗ subop2 ⊗ subop3
 op3 = 0.3*LazyTensor(b_l, b_r, 3, subop3)
 op3_ = 0.3*I1 ⊗ I2 ⊗ subop3
+op4 = 0.4*LazyTensor(b_l, b_r, 2, subop2)
+op4_ = 0.4*I1 ⊗ subop2 ⊗ I3
 
 x1 = Ket(b_r, rand(ComplexF64, length(b_r)))
 x2 = Ket(b_r, rand(ComplexF64, length(b_r)))
 xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 
-# Allowed addition
+# Addition one operator at same index
 fac = randn()
 @test dense(op3 + fac * op3) ≈ dense(op3) * (1 + fac)
 @test dense(op3 - fac * op3) ≈ dense(op3) * (1 - fac)
+# Addition one operator at different index
+@test dense(op3 + fac * op4) ≈ dense(op3_ + fac*op4_)
+@test dense(op3 - fac * op4) ≈ dense(op3_ - fac*op4_)
 
-# Forbidden addition
-@test_throws ArgumentError op1 + op2
-@test_throws ArgumentError op1 - op2
+#Test addition
+@test_throws QuantumOpticsBase.IncompatibleBases op1 + dagger(op2)
 @test 1e-14 > D(-op1_, -op1)
+@test 1e-14 > D(op1+op2, op1_+op2_)
+@test 1e-14 > D(op1+op2_, op1_+op2_)
+@test 1e-14 > D(op1_+op2, op1_+op2_)
+
+#Test substraction
+@test_throws QuantumOpticsBase.IncompatibleBases op1 - dagger(op2)
+@test 1e-14 > D(op1 - op2, op1_ - op2_)
+@test 1e-14 > D(op1 - op2_, op1_ - op2_)
+@test 1e-14 > D(op1_ - op2, op1_ - op2_)
+@test 1e-14 > D(op1 + (-op2), op1_ - op2_)
+@test 1e-14 > D(op1 + (-1*op2), op1_ - op2_)
+
+
+#Test tensor
+op1_tensor =  subop1 ⊗ op1
+op1_tensor_ =  subop1 ⊗ op1_
+op2_tensor =   op1 ⊗ subop1
+op2_tensor_ =   op1_ ⊗ subop1
+@test 1e-14 > D(op1_tensor,op1_tensor_)
+@test 1e-14 > D(op1_tensor_,op1_tensor)
+
+#Cannot tensor CompositeBasis with LazyTensor:
+@test_throws ArgumentError op1 ⊗ op1_
+@test_throws ArgumentError op2_ ⊗ op2
+
+
 
 # Test multiplication
 @test_throws DimensionMismatch op1*op2
@@ -412,5 +441,3 @@ Lop1 = LazyTensor(b1^2, b2^2, 2, sparse(randoperator(b1, b2)))
 @test_throws DimensionMismatch sparse(Lop1)*Lop1
 @test_throws DimensionMismatch Lop1*dense(Lop1)
 @test_throws DimensionMismatch Lop1*sparse(Lop1)
-
-end # testset