Lazy identities via FillArrays

Project.toml

@@ -5,6 +5,7 @@ version = "0.3.11"
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
@@ -16,6 +17,7 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 [compat]
 Adapt = "1, 2, 3.3"
 FFTW = "1.2"
+FillArrays = "0.13, 1"
 LRUCache = "1"
 QuantumInterface = "0.1.0"
 Strided = "1, 2"

src/QuantumOpticsBase.jl

@@ -1,6 +1,6 @@
 module QuantumOpticsBase
 
-using SparseArrays, LinearAlgebra, LRUCache, Strided, UnsafeArrays
+using SparseArrays, LinearAlgebra, LRUCache, Strided, UnsafeArrays, FillArrays
 import LinearAlgebra: mul!, rmul!
 
 import QuantumInterface: dagger, directsum, ⊕, dm, embed, expect, permutesystems,
@@ -17,7 +17,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
         #operators_dense
                 Operator, DenseOperator, DenseOpType, projector, dm,
         #operators_sparse
-                SparseOperator, diagonaloperator, SparseOpType,
+                SparseOperator, diagonaloperator, SparseOpType, EyeOpType,
         #operators_lazysum
                 LazySum,
         #operators_lazyproduct

src/operators.jl

@@ -350,7 +350,6 @@ to be used in the identity matrix.
 identityoperator(::Type{T}, ::Type{S}, b1::Basis, b2::Basis) where {T<:AbstractOperator,S} = throw(ArgumentError("Identity operator not defined for operator type $T."))
 identityoperator(::Type{T}, ::Type{S}, b::Basis) where {T<:AbstractOperator,S} = identityoperator(T,S,b,b)
 identityoperator(::Type{T}, bases::Basis...) where T<:AbstractOperator = identityoperator(T,eltype(T),bases...)
-identityoperator(b::Basis) = identityoperator(b,b)
 identityoperator(op::T) where {T<:AbstractOperator} = identityoperator(T, op.basis_l, op.basis_r)
 
 # Catch case where eltype cannot be inferred from type; this is a bit hacky

src/operators_sparse.jl

@@ -60,11 +60,24 @@ function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:Composite
     SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
 end
 
-identityoperator(::Type{T}, ::Type{S}, b1::Basis, b2::Basis) where {T<:DataOperator,S<:Number} =
+identityoperator(::Type{T}, ::Type{S}, b1::Basis, b2::Basis) where {T<:SparseOpType,S<:Number} =
     SparseOperator(b1, b2, sparse(one(S)*I, length(b1), length(b2)))
+identityoperator(::Type{T}, ::Type{S}, b::Basis) where {T<:SparseOpType,S<:Number} =
+    SparseOperator(b, b, sparse(one(S)*I, length(b), length(b)))
+
+const EyeOpPureType{BL,BR} = Operator{BL,BR,<:Eye}
+const EyeOpAdjType{BL,BR} = Operator{BL,BR,<:Adjoint{<:Number,<:Eye}}
+const EyeOpType{BL,BR} = Union{EyeOpPureType{BL,BR},EyeOpAdjType{BL,BR}}
+
+identityoperator(::Type{T}, ::Type{S}, b1::Basis, b2::Basis) where {T<:DataOperator,S<:Number} =
+    Operator(b1, b2, Eye{S}(length(b1), length(b2)))
+identityoperator(::Type{T}, ::Type{S}, b::Basis) where {T<:DataOperator,S<:Number} =
+    Operator(b, b, Eye{S}(length(b)))
+
 identityoperator(::Type{T}, b1::Basis, b2::Basis) where T<:Number = identityoperator(DataOperator, T, b1, b2)
-identityoperator(::Type{T}, b::Basis) where T<:Number = identityoperator(T, b, b)
+identityoperator(::Type{T}, b::Basis) where T<:Number = identityoperator(DataOperator, T, b)
 identityoperator(b1::Basis, b2::Basis) = identityoperator(ComplexF64, b1, b2)
+identityoperator(b::Basis) = identityoperator(ComplexF64, b)
 
 """
     diagonaloperator(b::Basis)

test/test_aqua.jl

@@ -3,7 +3,8 @@ using QuantumOpticsBase
 using Aqua
 
 @testset "aqua" begin
-    Aqua.test_all(QuantumOpticsBase,
-                  piracy=false  # TODO: Due to Base methods in QuantumOpticsBase, for types defined in QuantumInterface
+    Aqua.test_all(QuantumOpticsBase;
+                  piracy=false,  # TODO: Due to Base methods in QuantumOpticsBase, for types defined in QuantumInterface
+                  ambiguities=(exclude=[Base.reshape],)  # FIXME: Temporarily work-around ambiguities from FillArrays
                   )
 end # testset

test/test_operators_lazytensor.jl

@@ -378,7 +378,6 @@ result = deepcopy(result_)
 QuantumOpticsBase.mul!(result,op,state,alpha,beta)
 @test 1e-12 > D(result, alpha*op_*state + beta*result_)
 
-
 # Test gemm errors
 test_op = test_lazytensor(b1a, b1a, rand(2, 2))
 test_lazy = LazyTensor(tensor(b1a, b1a), [1, 2], (test_op, test_op))
@@ -414,3 +413,40 @@ Lop1 = LazyTensor(b1^2, b2^2, 2, sparse(randoperator(b1, b2)))
 @test_throws DimensionMismatch Lop1*sparse(Lop1)
 
 end # testset
+
+@testset "LazyTensor: explicit isometries" begin
+
+D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))
+D(x1::StateVector, x2::StateVector) = norm(x2-x1)
+
+# Test explicit identities and isometries
+bl = FockBasis(2) ⊗ GenericBasis(2) ⊗ SpinBasis(1//2) ⊗ GenericBasis(1) ⊗ GenericBasis(2)
+br = FockBasis(2) ⊗ GenericBasis(2) ⊗ SpinBasis(1//2) ⊗ GenericBasis(2) ⊗ GenericBasis(1)
+
+iso = identityoperator(bl.bases[5], br.bases[5])
+
+n1 = LazyTensor(bl, br, (1,3), (number(bl.bases[1]), sigmax(bl.bases[3])))
+n1_sp = LazyTensor(bl, br, (1,2,3,5), (number(bl.bases[1]), identityoperator(bl.bases[2]), sigmax(bl.bases[3]), iso))
+n1_de = LazyTensor(bl, br, (1,2,3,5), (dense(number(bl.bases[1])), identityoperator(bl.bases[2]), sigmax(bl.bases[3]), iso))
+
+@test dense(n1) == dense(n1_sp)
+@test dense(n1) == dense(n1_de)
+
+state = randoperator(br,br)
+
+@test 1e-12 > D(n1 * state, n1_sp * state)
+@test 1e-12 > D(n1 * state, n1_de * state)
+
+out = randoperator(bl,br)
+alpha = randn()
+beta = randn()
+out_ref = mul!(copy(out), n1, state, alpha, beta)
+@test 1e-12 > D(out_ref, mul!(copy(out), n1_sp, state, alpha, beta))
+@test 1e-12 > D(out_ref, mul!(copy(out), n1_de, state, alpha, beta))
+
+state = randoperator(bl,bl)
+out_ref = mul!(copy(out), state, n1, alpha, beta)
+@test 1e-12 > D(out_ref, mul!(copy(out), state, n1_sp, alpha, beta))
+@test 1e-12 > D(out_ref, mul!(copy(out), state, n1_de, alpha, beta))
+
+end

test/test_operators_sparse.jl

@@ -112,6 +112,11 @@ I = identityoperator(SparseOpType, b_l)
 @test 1e-11 > D(xbra1*I, xbra1)
 @test I == identityoperator(SparseOpType, b1a) ⊗ identityoperator(SparseOpType, b2a) ⊗ identityoperator(SparseOpType, b3a)
 
+IEye = identityoperator(b_l)
+@test isa(IEye, EyeOpType)
+@test sparse(IEye) == I
+Icomp = identityoperator(b1a) ⊗ identityoperator(b2a) ⊗ identityoperator(b3a)
+@test IEye == Icomp
 
 # Test tr and normalize
 op = sparse(DenseOperator(GenericBasis(3), [1 3 2;5 2 2;-1 2 5]))