From f9c3a8c2bdb4a397977028f33fa0b97c89455271 Mon Sep 17 00:00:00 2001
From: Matias Bundgaard-Nielsen <56378038+mabuni1998@users.noreply.github.com>
Date: Fri, 7 Apr 2023 17:26:39 -0400
Subject: [PATCH 1/5] Arithemtic between Operators and LazyOperators

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

Note that tensor(a::Operator,b::LazyTensor) only behaves well if Operator does not contain a CompositeBasis. This could probably be fixed by specifying the type more explicitly and returning error if it does contain CompositeBasis. Other solutions are welcome.
---
 src/LazyOperatorAritmetic.jl | 107 +++++++++++++++++++++++++++++++++++
 1 file changed, 107 insertions(+)
 create mode 100644 src/LazyOperatorAritmetic.jl

diff --git a/src/LazyOperatorAritmetic.jl b/src/LazyOperatorAritmetic.jl
new file mode 100644
index 00000000..f9ae77c7
--- /dev/null
+++ b/src/LazyOperatorAritmetic.jl
@@ -0,0 +1,107 @@
+
+function QuantumOpticsBase.:+(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    LazySum(a,b)
+end
+
+function QuantumOpticsBase.:-(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    LazySum([1,-1],[a,b])
+end
+
+function QuantumOpticsBase.:+(a::LazyTensor{B1,B2},b::Operator{B1,B2}) where {B1,B2}
+    LazySum(a) + b
+end
+function QuantumOpticsBase.:+(a::Operator{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    +(b,a)
+end
+function QuantumOpticsBase.:+(a::LazyProduct{B1,B2},b::Operator{B1,B2}) where {B1,B2}
+    LazySum(a) + b
+end
+function QuantumOpticsBase.:+(a::Operator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    +(b,a)
+end
+function QuantumOpticsBase.:-(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    LazySum(a) - b
+end
+function QuantumOpticsBase.:-(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    a-LazySum(b)
+end
+function QuantumOpticsBase.:+(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    LazySum(a) + LazySum(b)
+end
+function QuantumOpticsBase.:+(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    +(b,a)
+end
+
+function QuantumOpticsBase.:⊗(a::LazyTensor,b::Operator)
+    if isequal(b,identityoperator(basis(b)))
+        btotal = basis(a) ⊗ basis(b)
+        LazyTensor(btotal,btotal,[a.indices...],(a.operators...,),a.factor)
+    else
+        a ⊗ LazyTensor(b.basis_l,b.basis_r,[1],(b,),1)
+    end
+end
+
+function QuantumOpticsBase.:⊗(a::Operator,b::LazyTensor)
+    if isequal(a,identityoperator(basis(a)))
+        btotal = basis(a) ⊗ basis(b)
+        LazyTensor(btotal,btotal,[b.indices...].+1 ,(b.operators...,),b.factor)
+    else
+        LazyTensor(a.basis_l,a.basis_r,[1],(a,),1) ⊗ b
+    end
+end
+
+
+
+function QuantumOpticsBase.:⊗(a::Operator,b::LazySum)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    ops = Vector{AbstractOperator}(undef,length(b.operators))
+    for i in eachindex(ops)
+        ops[i] = a ⊗ b.operators[i]
+    end
+    LazySum(btotal_l,btotal_r,b.factors,ops)
+end
+function QuantumOpticsBase.:⊗(a::LazySum,b::Operator)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    ops = Vector{AbstractOperator}(undef,length(a.operators))
+    for i in eachindex(ops)
+        ops[i] = a.operators[i] ⊗ b
+    end
+    LazySum(btotal_l,btotal_r,a.factors,ops)
+end
+
+
+function QuantumOpticsBase.:⊗(a::Operator,b::LazyProduct)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    ops = Vector{AbstractOperator}(undef,length(b.operators))
+    for i in eachindex(ops)
+        ops[i] = a ⊗ b.operators[i]
+    end
+    LazyProduct(btotal_l,btotal_r,b.factor,ops)
+end
+function QuantumOpticsBase.:⊗(a::LazyProduct,b::Operator)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    ops = Vector{AbstractOperator}(undef,length(a.operators))
+    for i in eachindex(ops)
+        ops[i] = a.operators[i] ⊗ b
+    end
+    LazyProduct(btotal_l,btotal_r,a.factor,ops)
+end
+
+function QuantumOpticsBase.:⊗(a::AbstractOperator,b::LazyTensor)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    indices = (1,(b.indices .+ 1)...)
+    ops = (a,b.operators...)
+    LazyTensor(btotal_l,btotal_r,indices,ops,b.factor)
+end
+function QuantumOpticsBase.:⊗(a::LazyTensor,b::AbstractOperator)
+    btotal_l = a.basis_l ⊗ b.basis_l
+    btotal_r = a.basis_r ⊗ b.basis_r
+    indices = (a.indices...,length(a.indices)+1)
+    ops = (a.operators...,b)
+    LazyTensor(btotal_l,btotal_r,indices,ops,a.factor)
+end

From f4b872eaaa7b014afa9ad6d502ac81ca756dd805 Mon Sep 17 00:00:00 2001
From: mabuni1998 <mabuni1998@gmail.com>
Date: Sat, 8 Apr 2023 13:36:27 -0400
Subject: [PATCH 2/5] Tuples now used in lazyarithmetic

---
 src/QuantumOpticsBase.jl                      |  1 +
 ...ritmetic.jl => operators_lazyaritmetic.jl} | 26 +++++++------------
 2 files changed, 10 insertions(+), 17 deletions(-)
 rename src/{LazyOperatorAritmetic.jl => operators_lazyaritmetic.jl} (84%)

diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl
index 70d7bc5b..b964a12c 100644
--- a/src/QuantumOpticsBase.jl
+++ b/src/QuantumOpticsBase.jl
@@ -80,5 +80,6 @@ include("metrics.jl")
 include("spinors.jl")
 include("phasespace.jl")
 include("printing.jl")
+include("operators_lazyaritmetic.jl")
 
 end # module
diff --git a/src/LazyOperatorAritmetic.jl b/src/operators_lazyaritmetic.jl
similarity index 84%
rename from src/LazyOperatorAritmetic.jl
rename to src/operators_lazyaritmetic.jl
index f9ae77c7..f84b79ec 100644
--- a/src/LazyOperatorAritmetic.jl
+++ b/src/operators_lazyaritmetic.jl
@@ -51,23 +51,16 @@ function QuantumOpticsBase.:⊗(a::Operator,b::LazyTensor)
 end
 
 
-
 function QuantumOpticsBase.:⊗(a::Operator,b::LazySum)
     btotal_l = a.basis_l ⊗ b.basis_l
     btotal_r = a.basis_r ⊗ b.basis_r
-    ops = Vector{AbstractOperator}(undef,length(b.operators))
-    for i in eachindex(ops)
-        ops[i] = a ⊗ b.operators[i]
-    end
+    ops = ([a ⊗ op op in b.operators]...,)
     LazySum(btotal_l,btotal_r,b.factors,ops)
 end
 function QuantumOpticsBase.:⊗(a::LazySum,b::Operator)
     btotal_l = a.basis_l ⊗ b.basis_l
     btotal_r = a.basis_r ⊗ b.basis_r
-    ops = Vector{AbstractOperator}(undef,length(a.operators))
-    for i in eachindex(ops)
-        ops[i] = a.operators[i] ⊗ b
-    end
+    ops = ([op ⊗ b op in a.operators]...,)
     LazySum(btotal_l,btotal_r,a.factors,ops)
 end
 
@@ -75,19 +68,13 @@ end
 function QuantumOpticsBase.:⊗(a::Operator,b::LazyProduct)
     btotal_l = a.basis_l ⊗ b.basis_l
     btotal_r = a.basis_r ⊗ b.basis_r
-    ops = Vector{AbstractOperator}(undef,length(b.operators))
-    for i in eachindex(ops)
-        ops[i] = a ⊗ b.operators[i]
-    end
+    ops = ([a ⊗ op op in b.operators]...,)
     LazyProduct(btotal_l,btotal_r,b.factor,ops)
 end
 function QuantumOpticsBase.:⊗(a::LazyProduct,b::Operator)
     btotal_l = a.basis_l ⊗ b.basis_l
     btotal_r = a.basis_r ⊗ b.basis_r
-    ops = Vector{AbstractOperator}(undef,length(a.operators))
-    for i in eachindex(ops)
-        ops[i] = a.operators[i] ⊗ b
-    end
+    ops = ([op ⊗ b op in a.operators]...,)
     LazyProduct(btotal_l,btotal_r,a.factor,ops)
 end
 
@@ -105,3 +92,8 @@ function QuantumOpticsBase.:⊗(a::LazyTensor,b::AbstractOperator)
     ops = (a.operators...,b)
     LazyTensor(btotal_l,btotal_r,indices,ops,a.factor)
 end
+
+
+
+
+

From 690fff2714b7c550adba1b45722dc27aa9335cca Mon Sep 17 00:00:00 2001
From: mabuni1998 <mabuni1998@gmail.com>
Date: Mon, 10 Apr 2023 14:22:27 -0400
Subject: [PATCH 3/5] Tests for arithmetic of Lazy operators

Added test for addition and substraction of LazyTensor and LazyProduct.
Also tensor of LazyTensor, LazySum, and LazyProduct.
---
 src/QuantumOpticsBase.jl           |  2 +-
 src/operators_lazyaritmetic.jl     | 99 ------------------------------
 src/operators_lazyproduct.jl       | 30 +++++++++
 src/operators_lazysum.jl           | 13 ++++
 src/operators_lazytensor.jl        | 48 +++++++++++----
 test/test_abstractdata.jl          | 10 +--
 test/test_operators_lazyproduct.jl | 28 ++++++++-
 test/test_operators_lazysum.jl     |  9 +++
 test/test_operators_lazytensor.jl  | 54 +++++++++++++++-
 9 files changed, 173 insertions(+), 120 deletions(-)
 delete mode 100644 src/operators_lazyaritmetic.jl

diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl
index b964a12c..ed1da309 100644
--- a/src/QuantumOpticsBase.jl
+++ b/src/QuantumOpticsBase.jl
@@ -80,6 +80,6 @@ include("metrics.jl")
 include("spinors.jl")
 include("phasespace.jl")
 include("printing.jl")
-include("operators_lazyaritmetic.jl")
+#include("operators_lazyaritmetic.jl")
 
 end # module
diff --git a/src/operators_lazyaritmetic.jl b/src/operators_lazyaritmetic.jl
deleted file mode 100644
index f84b79ec..00000000
--- a/src/operators_lazyaritmetic.jl
+++ /dev/null
@@ -1,99 +0,0 @@
-
-function QuantumOpticsBase.:+(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
-    LazySum(a,b)
-end
-
-function QuantumOpticsBase.:-(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
-    LazySum([1,-1],[a,b])
-end
-
-function QuantumOpticsBase.:+(a::LazyTensor{B1,B2},b::Operator{B1,B2}) where {B1,B2}
-    LazySum(a) + b
-end
-function QuantumOpticsBase.:+(a::Operator{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
-    +(b,a)
-end
-function QuantumOpticsBase.:+(a::LazyProduct{B1,B2},b::Operator{B1,B2}) where {B1,B2}
-    LazySum(a) + b
-end
-function QuantumOpticsBase.:+(a::Operator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    +(b,a)
-end
-function QuantumOpticsBase.:-(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    LazySum(a) - b
-end
-function QuantumOpticsBase.:-(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    a-LazySum(b)
-end
-function QuantumOpticsBase.:+(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    LazySum(a) + LazySum(b)
-end
-function QuantumOpticsBase.:+(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    +(b,a)
-end
-
-function QuantumOpticsBase.:⊗(a::LazyTensor,b::Operator)
-    if isequal(b,identityoperator(basis(b)))
-        btotal = basis(a) ⊗ basis(b)
-        LazyTensor(btotal,btotal,[a.indices...],(a.operators...,),a.factor)
-    else
-        a ⊗ LazyTensor(b.basis_l,b.basis_r,[1],(b,),1)
-    end
-end
-
-function QuantumOpticsBase.:⊗(a::Operator,b::LazyTensor)
-    if isequal(a,identityoperator(basis(a)))
-        btotal = basis(a) ⊗ basis(b)
-        LazyTensor(btotal,btotal,[b.indices...].+1 ,(b.operators...,),b.factor)
-    else
-        LazyTensor(a.basis_l,a.basis_r,[1],(a,),1) ⊗ b
-    end
-end
-
-
-function QuantumOpticsBase.:⊗(a::Operator,b::LazySum)
-    btotal_l = a.basis_l ⊗ b.basis_l
-    btotal_r = a.basis_r ⊗ b.basis_r
-    ops = ([a ⊗ op op in b.operators]...,)
-    LazySum(btotal_l,btotal_r,b.factors,ops)
-end
-function QuantumOpticsBase.:⊗(a::LazySum,b::Operator)
-    btotal_l = a.basis_l ⊗ b.basis_l
-    btotal_r = a.basis_r ⊗ b.basis_r
-    ops = ([op ⊗ b op in a.operators]...,)
-    LazySum(btotal_l,btotal_r,a.factors,ops)
-end
-
-
-function QuantumOpticsBase.:⊗(a::Operator,b::LazyProduct)
-    btotal_l = a.basis_l ⊗ b.basis_l
-    btotal_r = a.basis_r ⊗ b.basis_r
-    ops = ([a ⊗ op op in b.operators]...,)
-    LazyProduct(btotal_l,btotal_r,b.factor,ops)
-end
-function QuantumOpticsBase.:⊗(a::LazyProduct,b::Operator)
-    btotal_l = a.basis_l ⊗ b.basis_l
-    btotal_r = a.basis_r ⊗ b.basis_r
-    ops = ([op ⊗ b op in a.operators]...,)
-    LazyProduct(btotal_l,btotal_r,a.factor,ops)
-end
-
-function QuantumOpticsBase.:⊗(a::AbstractOperator,b::LazyTensor)
-    btotal_l = a.basis_l ⊗ b.basis_l
-    btotal_r = a.basis_r ⊗ b.basis_r
-    indices = (1,(b.indices .+ 1)...)
-    ops = (a,b.operators...)
-    LazyTensor(btotal_l,btotal_r,indices,ops,b.factor)
-end
-function QuantumOpticsBase.:⊗(a::LazyTensor,b::AbstractOperator)
-    btotal_l = a.basis_l ⊗ b.basis_l
-    btotal_r = a.basis_r ⊗ b.basis_r
-    indices = (a.indices...,length(a.indices)+1)
-    ops = (a.operators...,b)
-    LazyTensor(btotal_l,btotal_r,indices,ops,a.factor)
-end
-
-
-
-
-
diff --git a/src/operators_lazyproduct.jl b/src/operators_lazyproduct.jl
index dd373282..039fe29a 100644
--- a/src/operators_lazyproduct.jl
+++ b/src/operators_lazyproduct.jl
@@ -56,6 +56,26 @@ isequal(x::LazyProduct{B1,B2}, y::LazyProduct{B1,B2}) where {B1,B2} = (samebases
 # Arithmetic operations
 -(a::T) where T<:LazyProduct = T(a.operators,a.ket_l,a.bra_r, -a.factor)
 
+
+function +(a::LazyProduct{B1,B2},b::Operator{B1,B2}) where {B1,B2}
+    LazySum(a) + b
+end
+function +(a::Operator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    +(b,a)
+end
+function -(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    LazySum(a) - b
+end
+function +(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    LazySum(a) + LazySum(b)
+end
+function -(a::LazyProduct{B1,B2},b::Operator{B1,B2}) where {B1,B2}
+    LazySum(a) - b
+end
+function -(a::Operator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
+    a - LazySum(b)
+end
+
 *(a::LazyProduct{B1,B2}, b::LazyProduct{B2,B3}) where {B1,B2,B3} = LazyProduct((a.operators..., b.operators...), a.factor*b.factor)
 *(a::LazyProduct, b::Number) = LazyProduct(a.operators, a.factor*b)
 *(a::Number, b::LazyProduct) = LazyProduct(b.operators, a*b.factor)
@@ -74,6 +94,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))
 
 
+#Assume same basis
+function tensor(a::Operator{B1,B1},b::LazyProduct{B, B, F, T, KTL, BTR}) where {B1,B, F, T, KTL, BTR}
+    ops = ([(i == 1 ? a : identityoperator(a)) ⊗ op for (i,op) in enumerate(b.operators)]...,)
+    LazyProduct(ops,b.factor)
+end
+function tensor(a::LazyProduct{B, B, F, T, KTL, BTR},b::Operator{B1,B1}) where {B1,B, F, T, KTL, BTR}
+    ops = ([op ⊗ (i == 1 ? b : identityoperator(b)) 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)
diff --git a/src/operators_lazysum.jl b/src/operators_lazysum.jl
index 5ba01731..bcc49898 100644
--- a/src/operators_lazysum.jl
+++ b/src/operators_lazysum.jl
@@ -106,6 +106,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)
diff --git a/src/operators_lazytensor.jl b/src/operators_lazytensor.jl
index 1e980f9f..36425b08 100644
--- a/src/operators_lazytensor.jl
+++ b/src/operators_lazytensor.jl
@@ -96,22 +96,48 @@ 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
-        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."))
+function +(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    LazySum(a,b)
+end
+function -(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    LazySum((1,-1),(a,b))
+end
+function +(a::LazyTensor{B1,B2},b::Operator{B1,B2}) where {B1,B2}
+    LazySum(a) + b
+end
+function +(a::Operator{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    +(b,a)
+end
+function -(a::LazyTensor{B1,B2},b::Operator{B1,B2}) where {B1,B2}
+    LazySum(a) - b
+end
+function -(a::Operator{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
+    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
-        op = a.operators[1] * a.factor - b.operators[1] * b.factor
-        return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,))
+function tensor(a::LazyTensor{B1,B1},b::Operator{B2,B2}) where {B1,B2}
+    if isequal(b,identityoperator(basis(b)))
+        btotal = basis(a) ⊗ basis(b)
+        LazyTensor(btotal,btotal,a.indices,(a.operators...,),a.factor)
+    elseif B2 <: CompositeBasis
+        throw(ArgumentError("tensor(a::LazyTensor{B1,B1},b::Operator{B2,B2}) is not implemented for B2 being CompositeBasis "))
+    else
+        a ⊗ LazyTensor(b.basis_l,b.basis_r,[1],(b,),1)
+    end
+end
+function tensor(a::Operator{B1,B1},b::LazyTensor{B2,B2})  where {B1,B2}
+    if isequal(a,identityoperator(basis(a)))
+        btotal = basis(a) ⊗ basis(b)
+        LazyTensor(btotal,btotal,b.indices.+length(basis(a).shape) ,(b.operators...,),b.factor)
+    elseif B1 <: CompositeBasis
+        throw(ArgumentError("tensor(a::Operator{B1,B1},b::LazyTensor{B2,B2}) is not implemented for B1 being CompositeBasis "))
+    else
+        LazyTensor(a.basis_l,a.basis_r,[1],(a,),1) ⊗ b
     end
-    throw(ArgumentError("Subtraction of LazyTensor operators is only defined in case both operators act nontrivially on the same, single tensor factor."))
 end
 
+
+
 function *(a::LazyTensor{B1,B2}, b::LazyTensor{B2,B3}) where {B1,B2,B3}
     indices = sort(union(a.indices, b.indices))
     # ops = Vector{AbstractOperator}(undef, length(indices))
diff --git a/test/test_abstractdata.jl b/test/test_abstractdata.jl
index b4a9bc5d..cbf0033a 100644
--- a/test/test_abstractdata.jl
+++ b/test/test_abstractdata.jl
@@ -412,9 +412,9 @@ 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)))
 
-# Addition
-@test_throws ArgumentError op1 + op2
-@test_throws ArgumentError op1 - op2
+# Addition Addition of LazyTensor now returns LazySum
+#@test_throws ArgumentError op1 + op2
+#@test_throws ArgumentError op1 - op2
 @test D(-op1_, -op1, 1e-12)
 
 # Test multiplication
@@ -448,7 +448,9 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 
 # Addition
-@test_throws ArgumentError op1 + op2
+#Commented following line since addition of LazyProduct returns LazySum and is allowed.
+#@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
diff --git a/test/test_operators_lazyproduct.jl b/test/test_operators_lazyproduct.jl
index f2ee3d07..347677bf 100644
--- a/test/test_operators_lazyproduct.jl
+++ b/test/test_operators_lazyproduct.jl
@@ -62,9 +62,22 @@ 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
+#Following is commented since addition between LazyProduct now returns LazySum
+#@test_throws ArgumentError op1 + op2
+
+# Test Addition
 @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 Subtraction
+@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 multiplication
 @test_throws DimensionMismatch op1a*op1a
@@ -78,6 +91,17 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 # Test division
 @test 1e-14 > D(op1/7, op1_/7)
 
+#Test Tensor product (NOTE: it is assumed that BL == BR of both involved operators.)
+op1_l = randoperator(b_l, b_l)
+op2_l = randoperator(b_l, b_l)
+op_r = randoperator(b_r, b_r)
+op = op_r ⊗ LazyProduct([op1_l, sparse(op2_l)])*0.1
+op_ = op_r ⊗ (0.1*op1_l*op2_l)
+@test 1e-11 > D(op,op_)
+op = LazyProduct([op1_l, sparse(op2_l)])*0.1 ⊗ op_r
+op_ = (0.1*op1_l*op2_l) ⊗ op_r 
+@test 1e-11 > D(op,op_)
+
 # Test identityoperator
 Idense = identityoperator(DenseOpType, b_l)
 id = identityoperator(LazyProduct, b_l)
diff --git a/test/test_operators_lazysum.jl b/test/test_operators_lazysum.jl
index 508b47b8..7e17ef99 100644
--- a/test/test_operators_lazysum.jl
+++ b/test/test_operators_lazysum.jl
@@ -105,6 +105,15 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 # 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 tuples vs. vectors
 @test (op1+op1).operators isa Tuple
 @test (op1+op2).operators isa Tuple
diff --git a/test/test_operators_lazytensor.jl b/test/test_operators_lazytensor.jl
index bbaef1b9..e660f1a4 100644
--- a/test/test_operators_lazytensor.jl
+++ b/test/test_operators_lazytensor.jl
@@ -113,11 +113,25 @@ fac = randn()
 @test dense(op3 + fac * op3) ≈ dense(op3) * (1 + fac)
 @test dense(op3 - fac * op3) ≈ dense(op3) * (1 - fac)
 
-# Forbidden addition
-@test_throws ArgumentError op1 + op2
-@test_throws ArgumentError op1 - op2
+#Commented since addition between lazytensors now return LazySum and its therefore always allowed.
+# Forbidden addition. 
+#@test_throws ArgumentError op1 + op2
+#@test_throws ArgumentError op1 - op2
 @test 1e-14 > D(-op1_, -op1)
 
+#Test addition
+@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 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 multiplication
 @test_throws DimensionMismatch op1*op2
 @test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@@ -413,4 +427,38 @@ Lop1 = LazyTensor(b1^2, b2^2, 2, sparse(randoperator(b1, b2)))
 @test_throws DimensionMismatch Lop1*dense(Lop1)
 @test_throws DimensionMismatch Lop1*sparse(Lop1)
 
+
+#Test Tensor (only bl == br allowed)
+subop1 = randoperator(b1a, b1a)
+subop2 = randoperator(b2a, b2a)
+subop3 = randoperator(b3a, b3a)
+I1 = dense(identityoperator(b1a, b1a))
+I2 = dense(identityoperator(b2a, b2a))
+I3 = dense(identityoperator(b3a, b3a))
+op1 = LazyTensor(b_l, b_l, [1, 3], (subop1, sparse(subop3)), 0.1)
+op1_ = 0.1*subop1 ⊗ I2 ⊗ subop3
+op2 = LazyTensor(b_l, b_l, [2, 3], (sparse(subop2), subop3), 0.7)
+op2_ = 0.7*I1 ⊗ subop2 ⊗ subop3
+op3 = 0.3*LazyTensor(b_l, b_l, 3, subop3)
+op3_ = 0.3*I1 ⊗ I2 ⊗ subop3
+
+#Cannot tensor CompositeBasis with LazyTensor if its not identity:
+@test_throws ArgumentError op1 ⊗ op1_
+@test_throws ArgumentError op2_ ⊗ op2
+
+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)
+
+op1_tensor =  (I1 ⊗ I3) ⊗ op1
+op1_tensor_ =  (I1 ⊗ I3) ⊗ op1_
+op2_tensor =   op1 ⊗ (I1 ⊗ I3)
+op2_tensor_ =   op1_ ⊗ (I1 ⊗ I3)
+@test 1e-14 > D(op1_tensor,op1_tensor_)
+@test 1e-14 > D(op1_tensor_,op1_tensor)
+
+
 end # testset

From 64da84f45beb5fc813b4cc5e4afb5ed7a6bf0fc8 Mon Sep 17 00:00:00 2001
From: mabuni1998 <mabuni1998@gmail.com>
Date: Wed, 12 Apr 2023 20:54:50 -0400
Subject: [PATCH 4/5] Lazy arithmetics with BL!=BR and test added.

Abstract type LazyOperator added as supertype for LazySum, LazyTensor, and LazyProduct.  +(a::LazyOperator,b::LazyOperator) implemented to make LazyOperators contagious.

added identityoperator for composite basis since previously the following would return false:
```
b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(6)
b3b = GenericBasis(5)
I1 = dense(identityoperator(b1a, b1b))
I2 = dense(identityoperator(b2a, b2b))
I3 = dense(identityoperator(b3a, b3b))

a = tensor(I1,I3)
a == identityoperator(a)
```
---
 src/QuantumOpticsBase.jl           |  1 -
 src/operators_dense.jl             |  3 ++
 src/operators_lazyproduct.jl       | 32 +++----------
 src/operators_lazysum.jl           | 19 ++++++--
 src/operators_lazytensor.jl        | 58 +++++++++++------------
 test/test_abstractdata.jl          |  6 +--
 test/test_operators_lazyproduct.jl | 26 ++++++-----
 test/test_operators_lazysum.jl     | 14 ++++++
 test/test_operators_lazytensor.jl  | 74 ++++++++++++------------------
 9 files changed, 111 insertions(+), 122 deletions(-)

diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl
index ed1da309..70d7bc5b 100644
--- a/src/QuantumOpticsBase.jl
+++ b/src/QuantumOpticsBase.jl
@@ -80,6 +80,5 @@ include("metrics.jl")
 include("spinors.jl")
 include("phasespace.jl")
 include("printing.jl")
-#include("operators_lazyaritmetic.jl")
 
 end # module
diff --git a/src/operators_dense.jl b/src/operators_dense.jl
index c3a2dbb2..831f5f1a 100644
--- a/src/operators_dense.jl
+++ b/src/operators_dense.jl
@@ -234,6 +234,9 @@ permutesystems(a::AdjointOperator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:Co
 
 identityoperator(::Type{S}, ::Type{T}, b1::Basis, b2::Basis) where {S<:DenseOpType,T<:Number} =
     Operator(b1, b2, Matrix{T}(I, length(b1), length(b2)))
+identityoperator(::Type{S}, ::Type{T}, b1::CompositeBasis, b2::CompositeBasis) where {S<:DenseOpType,T<:Number} =
+    Operator(b1, b2, kron([Matrix{T}(I, length(b1.bases[i]), length(b2.bases[i])) for i in reverse(1:length(b1.bases))]...))
+
 
 """
     projector(a::Ket, b::Bra)
diff --git a/src/operators_lazyproduct.jl b/src/operators_lazyproduct.jl
index 039fe29a..3027f761 100644
--- a/src/operators_lazyproduct.jl
+++ b/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
@@ -56,26 +56,6 @@ isequal(x::LazyProduct{B1,B2}, y::LazyProduct{B1,B2}) where {B1,B2} = (samebases
 # Arithmetic operations
 -(a::T) where T<:LazyProduct = T(a.operators,a.ket_l,a.bra_r, -a.factor)
 
-
-function +(a::LazyProduct{B1,B2},b::Operator{B1,B2}) where {B1,B2}
-    LazySum(a) + b
-end
-function +(a::Operator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    +(b,a)
-end
-function -(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    LazySum(a) - b
-end
-function +(a::LazyProduct{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    LazySum(a) + LazySum(b)
-end
-function -(a::LazyProduct{B1,B2},b::Operator{B1,B2}) where {B1,B2}
-    LazySum(a) - b
-end
-function -(a::Operator{B1,B2},b::LazyProduct{B1,B2}) where {B1,B2}
-    a - LazySum(b)
-end
-
 *(a::LazyProduct{B1,B2}, b::LazyProduct{B2,B3}) where {B1,B2,B3} = LazyProduct((a.operators..., b.operators...), a.factor*b.factor)
 *(a::LazyProduct, b::Number) = LazyProduct(a.operators, a.factor*b)
 *(a::Number, b::LazyProduct) = LazyProduct(b.operators, a*b.factor)
@@ -94,13 +74,13 @@ 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))
 
 
-#Assume same basis
-function tensor(a::Operator{B1,B1},b::LazyProduct{B, B, F, T, KTL, BTR}) where {B1,B, F, T, KTL, BTR}
-    ops = ([(i == 1 ? a : identityoperator(a)) ⊗ op for (i,op) in enumerate(b.operators)]...,)
+
+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{B, B, F, T, KTL, BTR},b::Operator{B1,B1}) where {B1,B, F, T, KTL, BTR}
-    ops = ([op ⊗ (i == 1 ? b : identityoperator(b)) for (i,op) in enumerate(a.operators)]...,)
+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
 
diff --git a/src/operators_lazysum.jl b/src/operators_lazysum.jl
index bcc49898..dcb52584 100644
--- a/src/operators_lazysum.jl
+++ b/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) = promote_type(eltype(x.factors), mapreduce(eltype, promote_type, x.operators))
@@ -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
diff --git a/src/operators_lazytensor.jl b/src/operators_lazytensor.jl
index 36425b08..98b49787 100644
--- a/src/operators_lazytensor.jl
+++ b/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,48 +96,46 @@ 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}
-    LazySum(a,b)
-end
-function -(a::LazyTensor{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
-    LazySum((1,-1),(a,b))
-end
-function +(a::LazyTensor{B1,B2},b::Operator{B1,B2}) where {B1,B2}
-    LazySum(a) + b
-end
-function +(a::Operator{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
-    +(b,a)
-end
-function -(a::LazyTensor{B1,B2},b::Operator{B1,B2}) where {B1,B2}
-    LazySum(a) - b
+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
+    LazySum(a) + LazySum(b)
 end
-function -(a::Operator{B1,B2},b::LazyTensor{B1,B2}) where {B1,B2}
-    a - LazySum(b)
+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
+    LazySum(a) - LazySum(b)
 end
 
-function tensor(a::LazyTensor{B1,B1},b::Operator{B2,B2}) where {B1,B2}
-    if isequal(b,identityoperator(basis(b)))
-        btotal = basis(a) ⊗ basis(b)
-        LazyTensor(btotal,btotal,a.indices,(a.operators...,),a.factor)
-    elseif B2 <: CompositeBasis
-        throw(ArgumentError("tensor(a::LazyTensor{B1,B1},b::Operator{B2,B2}) is not implemented for B2 being CompositeBasis "))
+function tensor(a::LazyTensor{B1,B2},b::AbstractOperator{B3,B4}) where {B1,B2,B3,B4}
+    if isequal(b,identityoperator(b))
+        btotal_l = a.basis_l ⊗ b.basis_l
+        btotal_r = a.basis_r ⊗ b.basis_r
+        LazyTensor(btotal_l,btotal_r,a.indices,(a.operators...,),a.factor)
+    elseif 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::Operator{B1,B1},b::LazyTensor{B2,B2})  where {B1,B2}
-    if isequal(a,identityoperator(basis(a)))
-        btotal = basis(a) ⊗ basis(b)
-        LazyTensor(btotal,btotal,b.indices.+length(basis(a).shape) ,(b.operators...,),b.factor)
-    elseif B1 <: CompositeBasis
-        throw(ArgumentError("tensor(a::Operator{B1,B1},b::LazyTensor{B2,B2}) is not implemented for B1 being CompositeBasis "))
+function tensor(a::AbstractOperator{B1,B2},b::LazyTensor{B3,B4})  where {B1,B2,B3,B4}
+    if isequal(a,identityoperator(a))
+        btotal_l = a.basis_l ⊗ b.basis_l
+        btotal_r = a.basis_r ⊗ b.basis_r
+        LazyTensor(btotal_l,btotal_r,b.indices.+length(a.basis_l.shape) ,(b.operators...,),b.factor)
+    elseif 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))
     # ops = Vector{AbstractOperator}(undef, length(indices))
diff --git a/test/test_abstractdata.jl b/test/test_abstractdata.jl
index cbf0033a..659c1fca 100644
--- a/test/test_abstractdata.jl
+++ b/test/test_abstractdata.jl
@@ -412,9 +412,7 @@ 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)))
 
-# Addition Addition of LazyTensor now returns LazySum
-#@test_throws ArgumentError op1 + op2
-#@test_throws ArgumentError op1 - op2
+# Addition
 @test D(-op1_, -op1, 1e-12)
 
 # Test multiplication
@@ -448,8 +446,6 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 
 # Addition
-#Commented following line since addition of LazyProduct returns LazySum and is allowed.
-#@test_throws ArgumentError op1 + op2
 @test D(2.1*op1 + 0.3*op2, 2.1*op1_+0.3*op2_)
 @test D(-op1_, -op1)
 
diff --git a/test/test_operators_lazyproduct.jl b/test/test_operators_lazyproduct.jl
index 347677bf..2b26f11c 100644
--- a/test/test_operators_lazyproduct.jl
+++ b/test/test_operators_lazyproduct.jl
@@ -50,33 +50,40 @@ 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)))
 
-#Following is commented since addition between LazyProduct now returns LazySum
-#@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
@@ -91,15 +98,12 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))
 # Test division
 @test 1e-14 > D(op1/7, op1_/7)
 
-#Test Tensor product (NOTE: it is assumed that BL == BR of both involved operators.)
-op1_l = randoperator(b_l, b_l)
-op2_l = randoperator(b_l, b_l)
-op_r = randoperator(b_r, b_r)
-op = op_r ⊗ LazyProduct([op1_l, sparse(op2_l)])*0.1
-op_ = op_r ⊗ (0.1*op1_l*op2_l)
+#Test Tensor product
+op = op1a ⊗ op1
+op_ = op1a ⊗ op1_
 @test 1e-11 > D(op,op_)
-op = LazyProduct([op1_l, sparse(op2_l)])*0.1 ⊗ op_r
-op_ = (0.1*op1_l*op2_l) ⊗ op_r 
+op = op1 ⊗ op2a
+op_ = op1_ ⊗ op2a 
 @test 1e-11 > D(op,op_)
 
 # Test identityoperator
diff --git a/test/test_operators_lazysum.jl b/test/test_operators_lazysum.jl
index 7e17ef99..8b146e1f 100644
--- a/test/test_operators_lazysum.jl
+++ b/test/test_operators_lazysum.jl
@@ -113,6 +113,20 @@ 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
diff --git a/test/test_operators_lazytensor.jl b/test/test_operators_lazytensor.jl
index e660f1a4..e033af2f 100644
--- a/test/test_operators_lazytensor.jl
+++ b/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,29 +101,31 @@ 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)
-
-#Commented since addition between lazytensors now return LazySum and its therefore always allowed.
-# Forbidden addition. 
-#@test_throws ArgumentError op1 + op2
-#@test_throws ArgumentError op1 - op2
-@test 1e-14 > D(-op1_, -op1)
+# Addition one operator at different index
+@test dense(op3 + fac * op4) ≈ dense(op3_ + fac*op4_)
+@test dense(op3 - fac * op4) ≈ dense(op3_ - fac*op4_)
 
 #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_)
@@ -132,6 +133,27 @@ fac = randn()
 @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)
+
+op1_tensor =  (I1 ⊗ I3) ⊗ op1
+op1_tensor_ =  (I1 ⊗ I3) ⊗ op1_
+op2_tensor =   op1 ⊗ (I1 ⊗ I3)
+op2_tensor_ =   op1_ ⊗ (I1 ⊗ I3)
+@test 1e-14 > D(op1_tensor,op1_tensor_)
+@test 1e-14 > D(op1_tensor_,op1_tensor)
+
+#Cannot tensor CompositeBasis with LazyTensor if its not identity:
+@test_throws ArgumentError op1 ⊗ op1_
+@test_throws ArgumentError op2_ ⊗ op2
+
+
+
 # Test multiplication
 @test_throws DimensionMismatch op1*op2
 @test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@@ -426,39 +448,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)
-
-
-#Test Tensor (only bl == br allowed)
-subop1 = randoperator(b1a, b1a)
-subop2 = randoperator(b2a, b2a)
-subop3 = randoperator(b3a, b3a)
-I1 = dense(identityoperator(b1a, b1a))
-I2 = dense(identityoperator(b2a, b2a))
-I3 = dense(identityoperator(b3a, b3a))
-op1 = LazyTensor(b_l, b_l, [1, 3], (subop1, sparse(subop3)), 0.1)
-op1_ = 0.1*subop1 ⊗ I2 ⊗ subop3
-op2 = LazyTensor(b_l, b_l, [2, 3], (sparse(subop2), subop3), 0.7)
-op2_ = 0.7*I1 ⊗ subop2 ⊗ subop3
-op3 = 0.3*LazyTensor(b_l, b_l, 3, subop3)
-op3_ = 0.3*I1 ⊗ I2 ⊗ subop3
-
-#Cannot tensor CompositeBasis with LazyTensor if its not identity:
-@test_throws ArgumentError op1 ⊗ op1_
-@test_throws ArgumentError op2_ ⊗ op2
-
-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)
-
-op1_tensor =  (I1 ⊗ I3) ⊗ op1
-op1_tensor_ =  (I1 ⊗ I3) ⊗ op1_
-op2_tensor =   op1 ⊗ (I1 ⊗ I3)
-op2_tensor_ =   op1_ ⊗ (I1 ⊗ I3)
-@test 1e-14 > D(op1_tensor,op1_tensor_)
-@test 1e-14 > D(op1_tensor_,op1_tensor)
-
-
-end # testset

From d34baa8336411ce0d4679b1e7a35c763efd6f8a1 Mon Sep 17 00:00:00 2001
From: mabuni1998 <mabuni1998@gmail.com>
Date: Sat, 15 Apr 2023 09:29:12 -0400
Subject: [PATCH 5/5] LazyTensor identityoperator check removed.

Also removed identityopreator method for CompositeBasis
---
 src/operators_dense.jl            |  3 ---
 src/operators_lazytensor.jl       | 12 ++----------
 test/test_operators_lazytensor.jl |  9 +--------
 3 files changed, 3 insertions(+), 21 deletions(-)

diff --git a/src/operators_dense.jl b/src/operators_dense.jl
index 831f5f1a..c3a2dbb2 100644
--- a/src/operators_dense.jl
+++ b/src/operators_dense.jl
@@ -234,9 +234,6 @@ permutesystems(a::AdjointOperator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:Co
 
 identityoperator(::Type{S}, ::Type{T}, b1::Basis, b2::Basis) where {S<:DenseOpType,T<:Number} =
     Operator(b1, b2, Matrix{T}(I, length(b1), length(b2)))
-identityoperator(::Type{S}, ::Type{T}, b1::CompositeBasis, b2::CompositeBasis) where {S<:DenseOpType,T<:Number} =
-    Operator(b1, b2, kron([Matrix{T}(I, length(b1.bases[i]), length(b2.bases[i])) for i in reverse(1:length(b1.bases))]...))
-
 
 """
     projector(a::Ket, b::Bra)
diff --git a/src/operators_lazytensor.jl b/src/operators_lazytensor.jl
index 98b49787..34626111 100644
--- a/src/operators_lazytensor.jl
+++ b/src/operators_lazytensor.jl
@@ -113,22 +113,14 @@ function -(a::T1,b::T2) where {T1 <: single_dataoperator{B1,B2},T2 <: single_dat
 end
 
 function tensor(a::LazyTensor{B1,B2},b::AbstractOperator{B3,B4}) where {B1,B2,B3,B4}
-    if isequal(b,identityoperator(b))
-        btotal_l = a.basis_l ⊗ b.basis_l
-        btotal_r = a.basis_r ⊗ b.basis_r
-        LazyTensor(btotal_l,btotal_r,a.indices,(a.operators...,),a.factor)
-    elseif B3 <: CompositeBasis || B4 <: CompositeBasis
+    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 isequal(a,identityoperator(a))
-        btotal_l = a.basis_l ⊗ b.basis_l
-        btotal_r = a.basis_r ⊗ b.basis_r
-        LazyTensor(btotal_l,btotal_r,b.indices.+length(a.basis_l.shape) ,(b.operators...,),b.factor)
-    elseif B1 <: CompositeBasis || B2 <: CompositeBasis
+    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
diff --git a/test/test_operators_lazytensor.jl b/test/test_operators_lazytensor.jl
index e033af2f..c7cef8b9 100644
--- a/test/test_operators_lazytensor.jl
+++ b/test/test_operators_lazytensor.jl
@@ -141,14 +141,7 @@ op2_tensor_ =   op1_ ⊗ subop1
 @test 1e-14 > D(op1_tensor,op1_tensor_)
 @test 1e-14 > D(op1_tensor_,op1_tensor)
 
-op1_tensor =  (I1 ⊗ I3) ⊗ op1
-op1_tensor_ =  (I1 ⊗ I3) ⊗ op1_
-op2_tensor =   op1 ⊗ (I1 ⊗ I3)
-op2_tensor_ =   op1_ ⊗ (I1 ⊗ I3)
-@test 1e-14 > D(op1_tensor,op1_tensor_)
-@test 1e-14 > D(op1_tensor_,op1_tensor)
-
-#Cannot tensor CompositeBasis with LazyTensor if its not identity:
+#Cannot tensor CompositeBasis with LazyTensor:
 @test_throws ArgumentError op1 ⊗ op1_
 @test_throws ArgumentError op2_ ⊗ op2