declare identitysuperoperator (#16)

qojulia · Jul 16, 2023 · cd247d8 · cd247d8 · Krastanov · Jul 16, 2023
1 parent e73f39a
commit cd247d8
Show file tree

Hide file tree

Showing 5 changed files with 15 additions and 6 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # News
 
+## v0.3.1
+
+- `identitysuperoperator` declared.
+
 ## v0.3.0
 
 - Redo some of the `identityoperator` methods.

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumInterface"
 uuid = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 authors = ["QuantumInterface.jl contributors"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

diff --git a/src/identityoperator.jl b/src/identityoperator.jl
@@ -20,3 +20,6 @@ identityoperator(op::T) where {T<:AbstractOperator} = identityoperator(T, op.bas
 identityoperator(::Type{T}, ::Type{Any}, b1::Basis, b2::Basis) where T<:AbstractOperator = identityoperator(T, ComplexF64, b1, b2)
 
 identityoperator(b1::Basis, b2::Basis) = identityoperator(ComplexF64, b1, b2)
+
+"""Prepare the identity superoperator over a given space."""
+function identitysuperoperator end
diff --git a/src/sparse.jl b/src/sparse.jl
@@ -9,26 +9,27 @@ Convert an arbitrary operator into a [`SparseOperator`](@ref).
 """
 sparse(a::AbstractOperator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) from QuantumOptics instead."))
 
-function ptrace(x::AbstractSparseMatrix, shape_nd, indices)
+function ptrace(x::AbstractSparseMatrix, shape_nd, indices) # A fairly QuantumOptics-specific method, left here to avoid piracy
     shape_nd = (shape_nd...,)
     N = div(length(shape_nd), 2)
     shape_2d = (x.m, x.n)
     shape_nd_after = ([i ∈ indices || i-N ∈ indices ? 1 : shape_nd[i] for i=1:2*N]...,)
     shape_2d_after = (prod(shape_nd_after[1:N]), prod(shape_nd_after[N+1:end]))
     I_nd_after_max = CartesianIndex(shape_nd_after...)
     y = spzeros(eltype(x), shape_2d_after...)
-    for I in eachindex(x)
-        I_nd = sub2sub(shape_2d, shape_nd, I)
+    for I in eachindex(x)::CartesianIndices # Manual type assertions to help JET
+        println(I.I)
+        I_nd = sub2sub(shape_2d, shape_nd, I)::CartesianIndex
         if I_nd.I[indices] != I_nd.I[indices .+ N]
             continue
         end
-        I_after = sub2sub(shape_nd_after, shape_2d_after, min(I_nd, I_nd_after_max))
+        I_after = sub2sub(shape_nd_after, shape_2d_after, min(I_nd, I_nd_after_max)::CartesianIndex)
         y[I_after] += x[I]
     end
     y
 end
 
-function sub2sub(shape1, shape2, I)
+function sub2sub(shape1, shape2, I::CartesianIndex)::CartesianIndex # Manual type assertions to help JET
     linearindex = LinearIndices(shape1)[I.I...]
     CartesianIndices(shape2)[linearindex]
 end
diff --git a/test/test_jet.jl b/test/test_jet.jl
@@ -1,4 +1,5 @@
 using QuantumInterface
+using Test
 using JET
 
 using JET: ReportPass, BasicPass, InferenceErrorReport, UncaughtExceptionReport, MethodErrorReport