Add sprepost to superoperators (#113)

* Add sprepost to superoperators

* Throw ArgumentError in spre and spost for non-square operators

* Bump version number

Author: akirakyle

Project.toml

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

docs/src/api.md

@@ -122,6 +122,10 @@ spre
 spost
 ```
 
+```@docs
+sprepost
+```
+
 ```@docs
 liouvillian
 ```

src/QuantumOpticsBase.jl

@@ -33,7 +33,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
                 current_time, time_shift, time_stretch, time_restrict,
         #superoperators
                 SuperOperator, DenseSuperOperator, DenseSuperOpType,
-                SparseSuperOperator, SparseSuperOpType, spre, spost, liouvillian,
+                SparseSuperOperator, SparseSuperOpType, spre, spost, sprepost, liouvillian,
         #fock
                 FockBasis, number, destroy, create,
                 fockstate, coherentstate, coherentstate!,

src/superoperators.jl

@@ -110,9 +110,16 @@ For operators ``A``, ``B`` the relation
 
 holds. `op` can be a dense or a sparse operator.
 """
-spre(op::AbstractOperator) = SuperOperator((op.basis_l, op.basis_l), (op.basis_r, op.basis_r), tensor(op, identityoperator(op)).data)
+function spre(op::AbstractOperator)
+    if !samebases(op.basis_l, op.basis_r)
+        throw(ArgumentError("It's not clear what spre of a non-square operator should be. See issue #113"))
+    end
+    SuperOperator((op.basis_l, op.basis_l), (op.basis_r, op.basis_r), tensor(op, identityoperator(op)).data)
+end
 
 """
+    spost(op)
+
 Create a super-operator equivalent for left side operator multiplication.
 
 For operators ``A``, ``B`` the relation
@@ -123,8 +130,27 @@ For operators ``A``, ``B`` the relation
 
 holds. `op` can be a dense or a sparse operator.
 """
-spost(op::AbstractOperator) = SuperOperator((op.basis_r, op.basis_r), (op.basis_l, op.basis_l), kron(permutedims(op.data), identityoperator(op).data))
+function spost(op::AbstractOperator)
+    if !samebases(op.basis_l, op.basis_r)
+        throw(ArgumentError("It's not clear what spost of a non-square operator should be. See issue #113"))
+    end
+    SuperOperator((op.basis_r, op.basis_r), (op.basis_l, op.basis_l), kron(permutedims(op.data), identityoperator(op).data))
+end
+
+"""
+    sprepost(op)
 
+Create a super-operator equivalent for left and right side operator multiplication.
+
+For operators ``A``, ``B``, ``C`` the relation
+
+```math
+    \\mathrm{sprepost}(A, B) C = A C B
+```
+
+holds. `A` ond `B` can be dense or a sparse operators.
+"""
+sprepost(A::AbstractOperator, B::AbstractOperator) = SuperOperator((A.basis_l, B.basis_r), (A.basis_r, B.basis_l), kron(permutedims(B.data), A.data))
 
 function _check_input(H::AbstractOperator{B1,B2}, J::Vector, Jdagger::Vector, rates) where {B1,B2}
     for j=J

test/test_apply.jl

@@ -43,12 +43,13 @@ b2 = FockBasis(2)
 k1 = coherentstate(b1, 0.39)
 k2 = coherentstate(b2, 1.4)
 op = projector(k1, k2')
-sOp2 = spre(op)
 
 try
+    sOp2 = spre(op)
     QuantumOpticsBase.apply!(st, 1, sOp2)
 catch e
-    @test typeof(e) <: QuantumInterface.IncompatibleBases
+    @test typeof(e) <: ArgumentError
+    #@test typeof(e) <: QuantumInterface.IncompatibleBases
 end
 
-end #testset
+end #testset

test/test_superoperators.jl

@@ -33,9 +33,14 @@ s_dense = dense(s_)
 
 # Test length
 b1 = GenericBasis(3)
+op = DenseOperator(b1, b1)
+S = spre(op)
+@test length(S) == length(S.data) == (3*3)^2
+S = spost(op)
+@test length(S) == length(S.data) == (3*3)^2
 b2 = GenericBasis(5)
 op = DenseOperator(b1, b2)
-S = spre(op)
+S = sprepost(dagger(op), op)
 @test length(S) == length(S.data) == (3*5)^2
 
 # Test arithmetic
@@ -156,9 +161,19 @@ op2 = DenseOperator(spinbasis, [0.2+0.1im 0.1+2.3im; 0.8+4.0im 0.3+1.4im])
 
 @test spre(sparse(op1))*op2 == op1*op2
 @test spost(sparse(op1))*op2 == op2*op1
+@test sprepost(sparse(op1), op1)*op2 ≈ op1*op2*op1
 
 @test spre(sparse(op1))*sparse(op2) == sparse(op1*op2)
 @test spost(sparse(op1))*sparse(op2) == sparse(op2*op1)
+@test sprepost(sparse(op1), sparse(op1))*sparse(op2) ≈ sparse(op1*op2*op1)
+
+@test sprepost(op1, op2) ≈ spre(op1)*spost(op2)
+b1 = FockBasis(1)
+b2 = FockBasis(5)
+op = fockstate(b1, 0) ⊗ dagger(fockstate(b2, 0))
+@test sprepost(dagger(op), op)*dm(fockstate(b1, 0)) == dm(fockstate(b2, 0))
+@test_throws ArgumentError spre(op)
+@test_throws ArgumentError spost(op)
 
 L = liouvillian(H, J)
 ρ = -1im*(H*ρ₀ - ρ₀*H)