Remove lazy wrapper.

src/QuantumOptics.jl

@@ -9,8 +9,9 @@ export bases, Basis, GenericBasis, CompositeBasis,
         operators, Operator, expect, identityoperator, ptrace, embed,
         operators_dense, DenseOperator, projector,
         operators_sparse, SparseOperator, diagonaloperator,
-        operators_lazy, lazy, LazyWrapper,
-                LazyTensor, LazySum, LazyProduct,
+        operators_lazysum, LazySum,
+        operators_lazyproduct, LazyProduct,
+        operators_lazytensor, LazyTensor,
         super, DenseSuperOperator, SparseSuperOperator,
                 spre, spost, liouvillian,
         fock, FockBasis, number, destroy, create,
@@ -40,7 +41,9 @@ include("operators.jl")
 include("operators_dense.jl")
 include("sparsematrix.jl")
 include("operators_sparse.jl")
-include("operators_lazy.jl")
+include("operators_lazysum.jl")
+include("operators_lazyproduct.jl")
+include("operators_lazytensor.jl")
 include("superoperators.jl")
 include("spin.jl")
 include("fock.jl")
@@ -71,7 +74,9 @@ using .states
 using .operators
 using .operators_dense
 using .operators_sparse
-using .operators_lazy
+using .operators_lazysum
+using .operators_lazyproduct
+using .operators_lazytensor
 using .superoperators
 using .spin
 using .fock

src/correlationexpansion.jl

@@ -10,7 +10,8 @@ using ..sortedindices
 using ..bases
 using ..operators
 using ..operators_dense
-using ..operators_lazy
+using ..operators_lazysum
+using ..operators_lazytensor
 using ..ode_dopri
 
 
@@ -300,7 +301,7 @@ function gemm!(alpha, a::LazyTensor, b::CorrelationExpansion, beta, result::Corr
         if isempty(I_)
             copy!(op.data, b.correlations[mask].data)
         else
-            operators = operators_lazy.suboperators(a, I_)
+            operators = operators_lazytensor.suboperators(a, I_)
             sortedindices.reducedindices!(I_, I)
             a_ = LazyTensor(op.basis_l, op.basis_r, I_, operators)
             gemm!(alpha, a_, b.correlations[mask], beta, op)
@@ -330,7 +331,7 @@ function gemm!(alpha, a::CorrelationExpansion, b::LazyTensor, beta, result::Corr
         if isempty(I_)
             copy!(op.data, a.correlations[mask].data)
         else
-            operators = operators_lazy.suboperators(b, I_)
+            operators = operators_lazytensor.suboperators(b, I_)
             sortedindices.reducedindices!(I_, I)
             b_ = LazyTensor(op.basis_l, op.basis_r, I_, operators)
             gemm!(alpha, a.correlations[mask], b_, beta, op)

src/cumulantexpansion.jl

@@ -4,7 +4,8 @@ using ..bases
 using ..states
 using ..operators
 using ..operators_dense
-using ..operators_lazy
+using ..operators_lazysum
+using ..operators_lazytensor
 using ..ode_dopri
 
 import Base: *, full

src/operators_lazyproduct.jl

@@ -1,3 +1,13 @@
+module operators_lazyproduct
+
+import Base: ==, *, /, +, -
+import ..operators: dagger, identityoperator,
+                    trace, ptrace, normalize!, tensor, permutesystems
+
+using ..bases, ..states, ..operators, ..operators_dense
+
+export LazyProduct
+
 """
 Lazy evaluation of product of operators.
 
@@ -30,9 +40,6 @@ Base.sparse(op::LazyProduct) = op.factor*prod(sparse(op_i) for op_i in op.operat
 *(a::LazyProduct, b::LazyProduct) = 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)
-*(a::LazyWrapper, b::LazyWrapper) = LazyProduct([a.operator, b.operator], a.factor*b.factor)
-*(a::LazyWrapper, b::LazyProduct) = LazyProduct([a.operator; b.operators], a.factor*b.factor)
-*(a::LazyProduct, b::LazyWrapper) = LazyProduct([a.operators; b.operator], a.factor*b.factor)
 
 /(a::LazyProduct, b::Number) = LazyProduct(a.operators, a.factor/b)
 
@@ -69,3 +76,5 @@ function operators.gemv!(alpha, a::Bra, b::LazyProduct, beta, result::Bra)
     end
     operators.gemv!(alpha, tmp1, b.operators[end], beta, result)
 end
+
+end # module

src/operators_lazysum.jl

@@ -1,3 +1,13 @@
+module operators_lazysum
+
+import Base: ==, *, /, +, -
+import ..operators: dagger, identityoperator,
+                    trace, ptrace, normalize!, tensor, permutesystems
+
+using ..bases, ..states, ..operators, ..operators_dense
+
+export LazySum
+
 """
 Lazy evaluation of sum of operators.
 
@@ -34,15 +44,9 @@ Base.sparse(op::LazySum) = sum(a*sparse(op_i) for (a, op_i) in zip(op.factors, o
 /(a::LazySum, b::Number) = LazySum(a.factors/b, a.operators)
 
 +(a::LazySum, b::LazySum) = (operators.check_samebases(a,b); LazySum([a.factors; b.factors], [a.operators; b.operators]))
-+(a::LazyWrapper, b::LazyWrapper) = (operators.check_samebases(a,b); LazySum([a.factor, b.factor], [a.operator, b.operator]))
-+(a::LazySum, b::LazyWrapper) = (operators.check_samebases(a,b); LazySum([a.factors; b.factor], [a.operators; b.operator]))
-+(a::LazyWrapper, b::LazySum) = (operators.check_samebases(a,b); LazySum([a.factor; b.factors], [a.operator; b.operators]))
 
 -(a::LazySum) = LazySum(-a.factors, a.operators)
 -(a::LazySum, b::LazySum) = (operators.check_samebases(a,b); LazySum([a.factors; -b.factors], [a.operators; b.operators]))
--(a::LazyWrapper, b::LazyWrapper) = (operators.check_samebases(a,b); LazySum([a.factor, -b.factor], [a.operator, b.operator]))
--(a::LazySum, b::LazyWrapper) = (operators.check_samebases(a,b); LazySum([a.factors; -b.factor], [a.operators; b.operator]))
--(a::LazyWrapper, b::LazySum) = (operators.check_samebases(a,b); LazySum([a.factor; -b.factors], [a.operator; b.operators]))
 
 
 identityoperator(::Type{LazySum}, b1::Basis, b2::Basis) = LazySum(identityoperator(b1, b2))
@@ -78,4 +82,6 @@ function operators.gemv!(alpha, a::Bra, b::LazySum, beta, result::Bra)
     for i=2:length(b.operators)
         operators.gemv!(alpha*b.factors[i], a, b.operators[i], Complex(1.), result)
     end
-end
+end
+
+end # module

src/operators_lazytensor.jl

@@ -1,3 +1,14 @@
+module operators_lazytensor
+
+import Base: ==, *, /, +, -
+import ..operators: dagger, identityoperator,
+                    trace, ptrace, normalize!, tensor, permutesystems
+
+using ..sortedindices
+using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
+
+export LazyTensor
+
 """
 Lazy implementation of a tensor product of operators.
 
@@ -147,9 +158,6 @@ function ptrace(op::LazyTensor, indices::Vector{Int})
 end
 
 tensor(a::LazyTensor, b::LazyTensor) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [a.indices; b.indices+length(a.basis_l.bases)], Operator[a.operators; b.operators], a.factor*b.factor)
-tensor(a::LazyWrapper, b::LazyWrapper) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [1, 2], Operator[a.operator, b.operator], a.factor*b.factor)
-tensor(a::LazyTensor, b::LazyWrapper) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [a.indices; length(a.basis_l.bases)+1], [a.operators; b.operator], a.factor*b.factor)
-tensor(a::LazyWrapper, b::LazyTensor) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [1; b.indices], [a.operator; b.operators], a.factor*b.factor)
 
 function permutesystems(op::LazyTensor, perm::Vector{Int})
     b_l = permutesystems(op.basis_l, perm)
@@ -173,7 +181,7 @@ function _gemm_recursive_dense_lazy(i_k::Int, N_k::Int, K::Int, J::Int, val::Com
         return nothing
     end
     if i_k in indices
-        h_i = operators_lazy.suboperator(h, i_k)
+        h_i = operators_lazytensor.suboperator(h, i_k)
         if isa(h_i, SparseOperator)
             h_i_data = h_i.data::SparseMatrixCSC{Complex128,Int}
             @inbounds for k=1:h_i_data.n
@@ -223,7 +231,7 @@ function _gemm_recursive_lazy_dense(i_k::Int, N_k::Int, K::Int, J::Int, val::Com
         return nothing
     end
     if i_k in indices
-        h_i = operators_lazy.suboperator(h, i_k)
+        h_i = suboperator(h, i_k)
         if isa(h_i, SparseOperator)
             h_i_data = h_i.data::SparseMatrixCSC{Complex128,Int}
             @inbounds for k=1:h_i_data.n
@@ -299,3 +307,5 @@ function operators.gemv!(alpha, a::Bra, b::LazyTensor, beta, result::Bra)
     result_data = reshape(result.data, 1, length(result.data))
     gemm(alpha, a_data, b, beta, result_data)
 end
+
+end # module

src/particle.jl

@@ -3,7 +3,7 @@ module particle
 import Base.==
 import ..operators
 
-using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse, ..operators_lazy
+using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
 
 export PositionBasis, MomentumBasis,
         gaussianstate,

test/test_correlationexpansion.jl

@@ -42,7 +42,7 @@ op1_ = ce.approximate(op1, S2 ∪ S3 ∪ S4)
 @test 1e-13 > D(op1, op1_)
 
 # Test multiplication
-h = 0.5*lazy(randop(b1)) ⊗ lazy(randop(b2)) ⊗ lazy(randop(b3)) ⊗ lazy(randop(b4))
+h = 0.5*LazyTensor(b, [1, 2, 3, 4], [randop(b1), randop(b2), randop(b3), randop(b4)])
 
 @test 1e-13 > D(full(h)*0.3*full(op1_), h*(0.3*op1_))
 @test 1e-13 > D(full(op1_)*0.3*full(h), (op1_*0.3)*h)
@@ -70,8 +70,8 @@ x = full(h)*op1
 rho = randdo(b1) ⊗ randdo(b2) ⊗ randdo(b3) ⊗ randdo(b4)
 rho_ce = ce.approximate(rho, S2 ∪ S3 ∪ S4)
 
-j1 = lazy(randop(b1)) ⊗ lazy(randop(b2)) ⊗ lazy(randop(b3)) ⊗ lazy(randop(b4))
-j2 = lazy(randop(b1)) ⊗ lazy(randop(b2)) ⊗ lazy(randop(b3)) ⊗ lazy(randop(b4))
+j1 = LazyTensor(b, [1, 2, 3, 4], [randop(b1), randop(b2), randop(b3), randop(b4)])
+j2 = LazyTensor(b, [1, 2, 3, 4], [randop(b1), randop(b2), randop(b3), randop(b4)])
 J = LazyTensor[j1, j2]
 v = rand(Float64, length(J))
 Γ = v * transpose(v)