`entanglement_entropy` can be a footgun (it has different meaning for Kets and Operators) -- consider deprecating in favor of separate functions

**Edit (by @Krastanov): This is expected behavior. See first comment below for possible ways forward.** **Original post below:** There seems to be a bug in the entanglement entropy calculation, it seems like it might be a factor of 2 off? Here's the code to reproduce the error: I create the following state and check is normalised: ``` julia> ψ_0 Ket(dim=8) basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)] 0.5773502691896257 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.5773502691896257 + 0.0im 0.0 + 0.0im 0.5773502691896257 + 0.0im julia> ψ_0'*ψ_0 1.0 + 0.0im ``` Now I calculate the entanglement between the A={1,2} and B={3} as follows: ``` julia> entanglement_entropy(dm(ψ_0),[3]) 1.2730283365896256 + 0.0im ``` However, this is wrong. I have done the calculation analytically (two methods, SVD and partial tracing) and it should be $-\frac{2}{3}\log(\frac{2}{3})-\frac{1}{3}\log(\frac{1}{3}) \approx 0.6365471166931516$, which could look like $1.2730283365896256 / 2$ but not exactly? ___ It's interesting cause your ptrace agrees with my calculations since ``` julia> ptrace(dm(ψ_0),[3]).data 4×4 Matrix{ComplexF64}: 0.333333+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.333333+0.0im 0.0+0.0im 0.333333+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.333333+0.0im 0.0+0.0im 0.333333+0.0im julia> eigvals(ptrace(dm(ψ_0),[3]).data) 4-element Vector{Float64}: 0.0 0.0 0.3333333333333333 0.6666666666666666 ``` and by inspection this spectrum does indeed agree with my calculations. I am a bit puzzled cause in the source code you guys have ``` function entropy_vn(rho::DenseOpType{B,B}; tol=1e-15) where B evals::Vector{ComplexF64} = eigvals(rho.data) entr = zero(eltype(rho)) for d ∈ evals if !(abs(d) < tol) entr -= d*log(d) end end return entr end ``` so it should be doing the same provided the spectrum is the same, which it should be since it is being passed the same density matrix. ``` function entanglement_entropy(psi::Ket{B}, partition, entropy_fun=entropy_vn) where B<:CompositeBasis # check that sites are within the range @assert all(partition .<= length(psi.basis.bases)) rho = ptrace(psi, partition) return entropy_fun(rho) end ``` To check I am not crazy: ``` julia> -0.3333333333333333*log(0.3333333333333333)-0.6666666666666666*log(0.6666666666666666) 0.6365141682948128 ``` I guess it could be that I am miss using the `partition` argument, I didn't find the documentation very clear.