Fix ptrace and operators

Author: albertomercurio

src/qobj/arithmetic_and_attributes.jl

@@ -515,7 +515,7 @@ function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{
     length(QO.dims) == 1 && return QO
 
     ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, SVector(sel))
-    return QuantumObject(ρtr, dims = dkeep)
+    return QuantumObject(ρtr, type=Operator, dims = dkeep)
 end
 
 ptrace(QO::QuantumObject{<:AbstractArray,BraQuantumObject}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
@@ -524,7 +524,7 @@ function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::U
     length(QO.dims) == 1 && return QO
 
     ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, SVector(sel))
-    return QuantumObject(ρtr, dims = dkeep)
+    return QuantumObject(ρtr, type=Operator, dims = dkeep)
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 
@@ -547,7 +547,7 @@ function _ptrace_ket(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
 
     vmat = reshape(QO, reverse(dims)...)
     topermute = nd + 1 .- sel_qtrace
-    vmat = PermutedDimsArray(vmat, topermute)
+    vmat = permutedims(vmat, topermute) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
     vmat = reshape(vmat, prod(dkeep), prod(dtrace))
 
     return vmat * vmat', dkeep
@@ -576,14 +576,14 @@ function _ptrace_oper(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
 
     ρmat = reshape(QO, reverse(vcat(dims, dims))...)
     topermute = 2 * nd + 1 .- qtrace_sel
-    ρmat = PermutedDimsArray(ρmat, reverse(topermute))
+    ρmat = permutedims(ρmat, reverse(topermute)) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
 
     ## TODO: Check if it works always
 
     # ρmat = row_major_reshape(ρmat, prod(dtrace), prod(dtrace), prod(dkeep), prod(dkeep))
     # res = dropdims(mapslices(tr, ρmat, dims=(1,2)), dims=(1,2))
     ρmat = reshape(ρmat, prod(dkeep), prod(dkeep), prod(dtrace), prod(dtrace))
-    res = dropdims(mapslices(tr, ρmat, dims = (3, 4)), dims = (3, 4))
+    res = map(tr, eachslice(ρmat, dims = (1, 2)))
 
     return res, dkeep
 end

src/qobj/operators.jl

@@ -13,7 +13,7 @@ export tunneling
 export qft
 
 @doc raw"""
-    rand_unitary(dimensions, distribution=:haar)
+    rand_unitary(dimensions, distribution=Val(:haar))
 
 Returns a random unitary [`QuantumObject`](@ref).
 
@@ -27,10 +27,13 @@ The `distribution` specifies which of the method used to obtain the unitary matr
 
 # References
 1. [F. Mezzadri, How to generate random matrices from the classical compact groups, arXiv:math-ph/0609050 (2007)](https://arxiv.org/abs/math-ph/0609050)
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `rand_unitary(dimensions, Val(distribution))` instead of `rand_unitary(dimensions, distribution)`. Also, put `dimensions` as `Tuple` or `SVector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-rand_unitary(dimensions::Int, distribution::Symbol = :haar) = rand_unitary(SVector(dimensions), Val(distribution))
-rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, distribution::Symbol = :haar) =
-    rand_unitary(dimensions, Val(distribution))
+rand_unitary(dimensions::Int, distribution::Union{Symbol,Val} = Val(:haar)) = rand_unitary(SVector(dimensions), makeVal(distribution))
+rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, distribution::Union{Symbol,Val} = Val(:haar)) =
+    rand_unitary(dimensions, makeVal(distribution))
 function rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{:haar})
     N = prod(dimensions)
 
@@ -224,7 +227,7 @@ function phase(N::Int, ϕ0::Real = 0)
 end
 
 @doc raw"""
-    jmat(j::Real, which::Symbol)
+    jmat(j::Real, which::Union{Symbol,Val})
 
 Generate higher-order Spin-`j` operators, where `j` is the spin quantum number and can be a non-negative integer or half-integer
 
@@ -250,7 +253,18 @@ Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=false
 2×2 SparseMatrixCSC{ComplexF64, Int64} with 1 stored entry:
      ⋅          ⋅    
  1.0+0.0im      ⋅
+
+julia> jmat(1.5, Val(:z))
+Quantum Object:   type=Operator   dims=[4]   size=(4, 4)   ishermitian=true
+4×4 SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
+ 1.5+0.0im      ⋅           ⋅           ⋅    
+     ⋅      0.5+0.0im       ⋅           ⋅    
+     ⋅          ⋅      -0.5+0.0im       ⋅    
+     ⋅          ⋅           ⋅      -1.5+0.0im
 ```
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `jmat(j, Val(which))` instead of `jmat(j, which)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 jmat(j::Real, which::Symbol) = jmat(j, Val(which))
 jmat(j::Real) = (jmat(j, Val(:x)), jmat(j, Val(:y)), jmat(j, Val(:z)))
@@ -427,8 +441,8 @@ d_j = \sigma_z^{\otimes j} \otimes \sigma_{-} \otimes I^{\otimes N-j-1}
 
 Note that the site index `j` should satisfy: `0 ≤ j ≤ N - 1`.
 
-> [!IMPORTANT]
-> If you want to keep type stability, it is recommended to use `fdestroy(Val(N), j)` instead of `fdestroy(N, j)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `fdestroy(Val(N), j)` instead of `fdestroy(N, j)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 fdestroy(N::Union{Int,Val}, j::Int) = _Jordan_Wigner(N, j, sigmam())
 
@@ -444,8 +458,8 @@ d_j^\dagger = \sigma_z^{\otimes j} \otimes \sigma_{+} \otimes I^{\otimes N-j-1}
 
 Note that the site index `j` should satisfy: `0 ≤ j ≤ N - 1`.
 
-> [!IMPORTANT]
-> If you want to keep type stability, it is recommended to use `fcreate(Val(N), j)` instead of `fcreate(N, j)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `fcreate(Val(N), j)` instead of `fcreate(N, j)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 fcreate(N::Union{Int,Val}, j::Int) = _Jordan_Wigner(N, j, sigmap())
 
@@ -470,7 +484,7 @@ end
 
 Generates the projection operator ``\hat{O} = \dyad{i}{j}`` with Hilbert space dimension `N`.
 """
-projection(N::Int, i::Int, j::Int) = QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N))
+projection(N::Int, i::Int, j::Int) = QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type=Operator)
 
 @doc raw"""
     tunneling(N::Int, m::Int=1; sparse::Union{Bool,Val{<:Bool}}=Val(false))
@@ -485,8 +499,8 @@ where ``N`` is the number of basis states in the Hilbert space, and ``m`` is the
 
 If `sparse=true`, the operator is returned as a sparse matrix, otherwise a dense matrix is returned.
 
-> [!IMPORTANT]
-> If you want to keep type stability, it is recommended to use `tunneling(N, m, Val(sparse))` instead of `tunneling(N, m, sparse)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `tunneling(N, m, Val(sparse))` instead of `tunneling(N, m, sparse)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function tunneling(N::Int, m::Int = 1; sparse::Union{Bool,Val} = Val(false))
     (m < 1) && throw(ArgumentError("The number of excitations (m) cannot be less than 1"))
@@ -522,6 +536,9 @@ The `dimensions` can be either the following types:
 ```
 
 where ``\omega = \exp(\frac{2 \pi i}{N})``.
+
+!!! warning "Beware of type-stability!"
+    It is highly recommended to use `qft(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 qft(dimensions::Int) = QuantumObject(_qft_op(dimensions), Operator, dimensions)
 qft(dimensions::Union{AbstractVector{T},Tuple}) where {T} =

src/qobj/states.jl

@@ -15,41 +15,50 @@ Returns a zero [`Ket`](@ref) vector with given argument `dimensions`.
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
 - `dimensions::Union{AbstractVector{Int}, Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+
+!!! warning "Beware of type-stability!"
+    It is highly recommended to use `zero_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 zero_ket(dimensions::Int) = QuantumObject(zeros(ComplexF64, dimensions), Ket, dimensions)
 zero_ket(dimensions::Union{AbstractVector{Int},Tuple}) =
     QuantumObject(zeros(ComplexF64, prod(dimensions)), Ket, dimensions)
 
 @doc raw"""
-    fock(N::Int, pos::Int=0; dims::Union{Int,AbstractVector{Int},Tuple}=N, sparse::Union{Bool,Val}=Val(false))
+    fock(N::Int, j::Int=0; dims::Union{Int,AbstractVector{Int},Tuple}=N, sparse::Union{Bool,Val}=Val(false))
 
 Generates a fock state ``\ket{\psi}`` of dimension `N`. 
 
 It is also possible to specify the list of dimensions `dims` if different subsystems are present.
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `fock(N, j, dims=dims, sparse=Val(sparse))` instead of `fock(N, j, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function fock(
     N::Int,
-    pos::Int = 0;
+    j::Int = 0;
     dims::Union{Int,AbstractVector{Int},Tuple} = N,
     sparse::Union{Bool,Val} = Val(false),
 )
     if getVal(makeVal(sparse))
-        array = sparsevec([pos + 1], [1.0 + 0im], N)
+        array = sparsevec([j + 1], [1.0 + 0im], N)
     else
         array = zeros(ComplexF64, N)
-        array[pos+1] = 1
+        array[j+1] = 1
     end
     return QuantumObject(array; type = Ket, dims = dims)
 end
 
 @doc raw"""
-    basis(N::Int, pos::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple}=N)
+    basis(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple}=N)
 
 Generates a fock state like [`fock`](@ref).
 
 It is also possible to specify the list of dimensions `dims` if different subsystems are present.
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `basis(N, j, dims=dims)` with `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-basis(N::Int, pos::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple} = N) = fock(N, pos, dims = dims)
+basis(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple} = N) = fock(N, j, dims = dims)
 
 @doc raw"""
     coherent(N::Int, α::Number)
@@ -68,6 +77,9 @@ Generate a random normalized [`Ket`](@ref) vector with given argument `dimension
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
 - `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `rand_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 rand_ket(dimensions::Int) = rand_ket(SVector(dimensions))
 function rand_ket(dimensions::Union{AbstractVector{Int},Tuple})
@@ -82,6 +94,9 @@ end
 Density matrix representation of a Fock state.
 
 Constructed via outer product of [`fock`](@ref).
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `fock_dm(N, pos, dims=dims, sparse=Val(sparse))` instead of `fock_dm(N, pos, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function fock_dm(
     N::Int,
@@ -113,8 +128,8 @@ Density matrix for a thermal state (generating thermal state probabilities) with
 - `n::Real`: Expectation value for number of particles in the thermal state.
 - `sparse::Union{Bool,Val}`: If `true`, return a sparse matrix representation.
 
-> [!IMPORTANT]
-> If you want to keep type stability, it is recommended to use `thermal_dm(N, n, Val(sparse))` instead of `thermal_dm(N, n, sparse)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `thermal_dm(N, n, sparse=Val(sparse))` instead of `thermal_dm(N, n, sparse=sparse)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function thermal_dm(N::Int, n::Real; sparse::Union{Bool,Val} = Val(false))
     β = log(1.0 / n + 1.0)
@@ -135,6 +150,9 @@ Returns the maximally mixed density matrix with given argument `dimensions`.
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
 - `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `maximally_mixed_dm(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 maximally_mixed_dm(dimensions::Int) = QuantumObject(I(dimensions) / complex(dimensions), Operator, SVector(dimensions))
 function maximally_mixed_dm(dimensions::Union{AbstractVector{Int},Tuple})
@@ -153,6 +171,9 @@ The `dimensions` can be either the following types:
 
 The default keyword argument `rank = prod(dimensions)` (full rank).
 
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `rand_dm(dimensions; rank=rank)` with `dimensions` as `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
+
 # References
 - [J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6.3 (1965): 440-449](https://doi.org/10.1063/1.1704292)
 - [K. Życzkowski, et al., Generating random density matrices, Journal of Mathematical Physics 52, 062201 (2011)](http://dx.doi.org/10.1063/1.3595693)
@@ -284,8 +305,8 @@ Returns the `n`-qubit [W-state](https://en.wikipedia.org/wiki/W_state):
 \frac{1}{\sqrt{n}} \left( |100...0\rangle + |010...0\rangle + \cdots + |00...01\rangle \right)
 ```
 
-> [!IMPORTANT]
-> If you want to keep type stability, it is recommended to use `w_state(Val(n))` instead of `w_state(n)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `w_state(Val(n))` instead of `w_state(n)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
 """
 function w_state(::Val{n}) where {n}
     nzind = 2 .^ (0:(n-1)) .+ 1
@@ -305,8 +326,8 @@ Returns the generalized `n`-qudit [Greenberger–Horne–Zeilinger (GHZ) state](
 
 Here, `d` specifies the dimension of each qudit. Default to `d=2` (qubit).
 
-> [!IMPORTANT]
-> If you want to keep type stability, it is recommended to use `ghz_state(Val(n); kwargs...)` instead of `ghz_state(n; kwargs...)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `ghz_state(Val(n))` instead of `ghz_state(n)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) for more details.
 """
 function ghz_state(::Val{n}; d::Int = 2) where {n}
     nzind = collect((0:(d-1)) .* Int((d^n - 1) / (d - 1)) .+ 1)