Fix eltype propagation in embed and LazySums. (#82)

Author: amilsted

src/operators.jl

@@ -90,7 +90,8 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
         (opsb.basis_r == basis_r.bases[idxsb]) || throw(IncompatibleBases())
     end
 
-    embed_op = tensor([i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...)
+    S = length(operators) > 1 ? mapreduce(eltype, promote_type, operators) : Any
+    embed_op = tensor([i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : identityoperator(T, S, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...)
 
     # Embed all joint-subspace operators.
     idxop_comp = [x for x in zip(indices, operators) if x[1] isa Array]
@@ -116,7 +117,7 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
     reduce(tensor, basis_r.bases[indices]) == op.basis_r || throw(IncompatibleBases())
 
     index_order = [idx for idx in 1:length(basis_l.bases) if idx ∉ indices]
-    all_operators = AbstractOperator[identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in index_order]
+    all_operators = AbstractOperator[identityoperator(T, eltype(op), basis_l.bases[i], basis_r.bases[i]) for i in index_order]
 
     for idx in indices
         pushfirst!(index_order, idx)
@@ -209,21 +210,22 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
     end
     indices, operator_list = zip(operators...)
     operator_list = [operator_list...;]
+    S = mapreduce(eltype, promote_type, operator_list)
     indices_flat = [indices...;]::Vector{Int} # type assertion to help type inference
     start_indices_flat = [i[1] for i in indices]
     complement_indices_flat = Int[i for i=1:N if i ∉ indices_flat]
     operators_flat = AbstractOperator[]
     if all([minimum(I):maximum(I);]==I for I in indices)
         for i in 1:N
             if i in complement_indices_flat
-                push!(operators_flat, identityoperator(T, basis_l.bases[i], basis_r.bases[i]))
+                push!(operators_flat, identityoperator(T, S, basis_l.bases[i], basis_r.bases[i]))
             elseif i in start_indices_flat
                 push!(operators_flat, operator_list[indexin(i, start_indices_flat)[1]])
             end
         end
         return tensor(operators_flat...)
     else
-        complement_operators = [identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in complement_indices_flat]
+        complement_operators = [identityoperator(T, S, basis_l.bases[i], basis_r.bases[i]) for i in complement_indices_flat]
         op = tensor([operator_list; complement_operators]...)
         perm = sortperm([indices_flat; complement_indices_flat])
         return permutesystems(op, perm)

src/operators_lazysum.jl

@@ -55,7 +55,7 @@ function LazySum(factors, operators)
     LazySum(Tf, factors, operators)
 end
 LazySum(::Type{Tf}, operators::AbstractOperator...) where Tf = LazySum(ones(Tf, length(operators)), (operators...,))
-LazySum(operators::AbstractOperator...) = LazySum(ComplexF64, operators...)
+LazySum(operators::AbstractOperator...) = LazySum(mapreduce(eltype, promote_type, operators), operators...)
 LazySum() = throw(ArgumentError("LazySum needs a basis, or at least one operator!"))
 
 Base.copy(x::LazySum) = @samebases LazySum(x.basis_l, x.basis_r, copy(x.factors), copy.(x.operators))

test/test_operators_lazysum.jl

@@ -47,13 +47,19 @@ op2 = sparse(randoperator(b_l, b_r))
 @test 0.1*sparse(op1) + 0.3*op2 == sparse(LazySum([0.1, 0.3], (op1, op2)))
 
 # Test embed
-x1 = randoperator(b1a,b1b)
-y1 = randoperator(b1a,b1b)
-xy1 = LazySum([1., 2.], (x1, y1))
-x = LazySum([1.], (embed(b_l, b_r, 1, x1),))
-y = LazySum([2.], (embed(b_l, b_r, 1, y1),))
-xy = x + y
-@test embed(b_l, b_r, [1], xy1) == xy
+for T in (Float32, Float64, ComplexF32, ComplexF64)
+    x1 = randoperator(T, b1a,b1b)
+    y1 = randoperator(T, b1a,b1b)
+    xy1 = LazySum(T[1., 2.], (x1, y1))
+    x = LazySum(T[1.], (embed(b_l, b_r, 1, x1),))
+    y = LazySum(T[2.], (embed(b_l, b_r, 1, y1),))
+    @test eltype(x) == T
+    @test eltype(y) == T
+    xy = x + y
+    @test embed(b_l, b_r, [1], xy1) == xy
+    @test eltype(xy) == T
+end
+
 
 # Arithmetic operations
 # =====================

test/test_operators_lazytensor.jl

@@ -62,7 +62,15 @@ x = LazyTensor(b_l, b_r, [1, 3], (op1, sparse(op3)), 0.3)
 @test 1e-12 > D(0.3*sparse(op1)⊗I2⊗sparse(op3), sparse(x))
 
 # Test eltype
-@test eltype(x) == ComplexF64
+for T in (Float32, Float64, ComplexF32, ComplexF64)
+    I2_T = identityoperator(T, b2a, b2b)
+    op1_T = randoperator(T, b1a, b1b)
+    op3_T = randoperator(T, b3a, b3b)
+    x_T =  LazyTensor(b_l, b_r, [1, 3], (op1_T, sparse(op3_T)), T(0.3))
+    @test eltype(x_T) == T
+    @test eltype(dense(x_T)) == T
+    @test eltype(sparse(x_T)) == T
+end
 
 # Test suboperators
 @test QuantumOpticsBase.suboperator(x, 1) == op1