Fixes and tests for apply!

Project.toml

@@ -1,6 +1,6 @@
 name = "QuantumOpticsBase"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "0.4.15"
+version = "0.4.16"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/apply.jl

@@ -1,27 +1,32 @@
-function apply!(state::Ket, indices, operation::Operator)
-    op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation)
-    if (length(indices)>1)
-        for i in 2:length(indices)
-            if indices[i]<indices[i-1]
-                op = permutesystems(op, indices)
-                break
-            end
-        end
+nsubsystems(s::Ket) = nsubsystems(s.basis)
+function nsubsystems(s::Operator)
+    s.basis_l == s.basis_r || throw(ArgumentError("`nsubsystem(::Operator)` is well defined only if the left and right bases are the same"))
+    nsubsystems(s.basis_l)
+end
+nsubsystems(b::CompositeBasis) = length(b.bases)
+nsubsystems(b::Basis) = 1
+
+function is_apply_shortcircuit(state, indices, operation)
+    if nsubsystems(state) == 1
+        basis(state)==basis(operation) || throw(ArgumentError("`apply!` failed due to incompatible bases of the state and the operation attempted to be applied on it"))
     end
+    basis(state)==basis(operation) || return false
+    j = 1
+    for i in indices
+        i == j || return false
+        j+=1
+    end
+    return j-1 == length(indices)
+end
+
+function apply!(state::Ket, indices, operation::Operator)
+    op = is_apply_shortcircuit(state, indices, operation) ? operation : embed(basis(state), indices, operation)
     state.data = (op*state).data
     state
 end
 
 function apply!(state::Operator, indices, operation::Operator)
-    op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation)
-    if (length(indices)>1)
-        for i in 2:length(indices)
-            if indices[i]<indices[i-1]
-                op = permutesystems(op, indices)
-                break
-            end
-        end
-    end
+    op = is_apply_shortcircuit(state, indices, operation) ? operation : embed(basis(state), indices, operation)
     state.data = (op*state*op').data
     state
 end
@@ -31,10 +36,10 @@ function apply!(state::Ket, indices, operation::T) where {T<:AbstractSuperOperat
 end
 
 function apply!(state::Operator, indices, operation::T) where {T<:AbstractSuperOperator}
-    if length(indices)>1
-        error("Applying SuperOperator to multiple qubits/operators is not supported currently, due to missing tensor product method for SuperOperators")
+    if is_apply_shortcircuit(state, indices, operation)
+        state.data = (operation*state).data
+        return state
+    else
+        error("`apply!` does not yet support embedding superoperators acting only on a subsystem of the given state")
     end
-    op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation)
-    state.data = (op*state).data
-    state
 end

test/test_apply.jl

@@ -1,6 +1,7 @@
 using QuantumOpticsBase
 using QuantumInterface
 using Test
+using QuantumOpticsBase: apply!
 
 @testset "apply" begin
 
@@ -11,66 +12,42 @@ _x = sigmax(_b2)
 _y = sigmay(_b2)
 _z = sigmaz(_b2)
 
-@test QuantumOpticsBase.apply!(_l0⊗_l1, 1, _x) == _l1⊗_l1
-@test QuantumOpticsBase.apply!(_x, 1, _y) == _x
+@test apply!(_l0⊗_l1, 1, _x) ≈ _l1⊗_l1
+@test apply!(_x, 1, _y) ≈ _x
 
 # Test Operator with IncompatibleBases
 _l01 = _l0⊗_l1
 op = projector(_l0, _l01')
-
-try
-    QuantumOpticsBase.apply!(_l0, 1, op)
-catch e
-    @test typeof(e) <: QuantumInterface.IncompatibleBases
-end
-
-try
-    QuantumOpticsBase.apply!(_x, 1, op)
-catch e
-    @test typeof(e) <: QuantumInterface.IncompatibleBases
-end
+@test_throws QuantumInterface.IncompatibleBases apply!(_l0, 1, op)
+@test_throws QuantumInterface.IncompatibleBases apply!(_x, 1, op)
 
 # Test SuperOperator
-sOp = spre(create(FockBasis(1)))
-st = coherentstate(FockBasis(1), 1.4)
-
-@test QuantumOpticsBase.apply!(st, 1, sOp).data == (sOp*projector(st)).data
+bf1 = FockBasis(1)
+sOp = spre(create(bf1))
+st = coherentstate(bf1, 1.4)
+@test apply!(st, 1, sOp).data ≈ (sOp*projector(st)).data
 
 # Test SuperOperator with IncompatibleBases
-b1 = FockBasis(1)
-b2 = FockBasis(2)
-
-k1 = coherentstate(b1, 0.39)
-k2 = coherentstate(b2, 1.4)
-op = projector(k1, k2')
-
-try
-    sOp2 = spre(op)
-    QuantumOpticsBase.apply!(st, 1, sOp2)
-catch e
-    @test typeof(e) <: ArgumentError
-    #@test typeof(e) <: QuantumInterface.IncompatibleBases
-end
+bf2 = FockBasis(2)
+op = create(bf2)
+sOp2 = spre(op)
+@test_throws ArgumentError apply!(st, 1, sOp2)
 
 # test CNOT₂₋₁
 CNOT = DenseOperator(_b2⊗_b2, Complex.([1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0]))
-@test QuantumOpticsBase.apply!(_l0⊗_l1, [2,1], CNOT) == _l1⊗_l1
+@test QuantumOpticsBase.apply!(_l0⊗_l1, [2,1], CNOT) ≈ _l1⊗_l1
 
 # test operator permutation with 3 qubits/operators for apply!
+@test QuantumOpticsBase.apply!(_l0⊗_l1⊗_l0, [2,3,1], _x⊗_y⊗_z) ≈ (_z*_l0)⊗(_x*_l1)⊗(_y*_l0)
+@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3,1], _x⊗_y⊗_z) ≈ (_z⊗_x⊗_y)*(_x⊗_y⊗_z)*(_z⊗_x⊗_y)'
 
-# 3-qubit operator permutation
-@test QuantumOpticsBase.apply!(_l0⊗_l1⊗_l0, [2,3,1], _x⊗_y⊗_z) ≈ (_y*_l0)⊗(_z*_l1)⊗(_x*_l0)
-
-# 3-operator operator permutation
-@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3,1], _x⊗_y⊗_z) ≈ (_y⊗_z⊗_x)*(_x⊗_y⊗_z)*(_y⊗_z⊗_x)'
+i = identityoperator(_y)
+@test QuantumOpticsBase.apply!(_x⊗_y⊗_z, [2,3], _y⊗_z) ≈ (i⊗_y⊗_z)*(_x⊗_y⊗_z)*(i⊗_y⊗_z)'
 
 # 3-qubit/operator errors when called for applying superoperator
-sOp1 = spre(create(FockBasis(1)))
-sOp2 = spost(create(FockBasis(1)))
-st1 = coherentstate(FockBasis(1), 1.4)
-st2 = coherentstate(FockBasis(1), 0.3)
-st3 = coherentstate(FockBasis(1), 0.8)
-
-@test_throws ErrorException QuantumOpticsBase.apply!(st1⊗st2⊗st3, [2,3,1], sOp1)
+sOp1 = spre(create(bf1))
+st1 = coherentstate(bf1, 1.4)
+st2 = coherentstate(bf1, 0.3)
+@test_broken apply!(st1⊗st2, [1], sOp1)
 
 end #testset