From 851abe0e642c34e475942548eb64eb6969a5f35a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Designolle?= Date: Mon, 25 Nov 2024 15:41:29 +0100 Subject: [PATCH] Apply runic in src/ and test/ --- src/basic.jl | 102 +++++++-------- src/entanglement.jl | 96 +++++++------- src/entropy.jl | 30 ++--- src/games.jl | 28 ++-- src/incompatibility.jl | 28 ++-- src/measurements.jl | 278 ++++++++++++++++++++-------------------- src/multilinear.jl | 80 ++++++------ src/nonlocal.jl | 72 +++++------ src/random.jl | 20 +-- src/seesaw.jl | 100 +++++++-------- src/states.jl | 29 ++--- src/supermaps.jl | 10 +- test/entanglement.jl | 26 ++-- test/entropy.jl | 2 +- test/incompatibility.jl | 14 +- test/measurements.jl | 4 +- test/nonlocal.jl | 12 +- test/norms.jl | 4 +- test/random.jl | 8 +- test/supermaps.jl | 8 +- 20 files changed, 475 insertions(+), 476 deletions(-) diff --git a/src/basic.jl b/src/basic.jl index 451566d..98466b8 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -3,7 +3,7 @@ Produces a ket of dimension `d` with nonzero element `i`. """ -function ket(::Type{T}, i::Integer, d::Integer = 2) where {T<:Number} +function ket(::Type{T}, i::Integer, d::Integer = 2) where {T <: Number} psi = zeros(T, d) psi[i] = 1 return psi @@ -26,7 +26,7 @@ export ketbra Produces a projector onto the basis state `i` in dimension `d`. """ -function proj(::Type{T}, i::Integer, d::Integer = 2) where {T<:Number} +function proj(::Type{T}, i::Integer, d::Integer = 2) where {T <: Number} p = Hermitian(zeros(T, d, d)) p[i, i] = 1 return p @@ -34,7 +34,7 @@ end proj(i::Integer, d::Integer = 2) = proj(Bool, i, d) export proj -const Measurement{T} = Vector{Hermitian{T,Matrix{T}}} +const Measurement{T} = Vector{Hermitian{T, Matrix{T}}} export Measurement """ @@ -54,7 +54,7 @@ Converts a set of measurements in the common tensor format into a matrix of (her By default, the second argument is fixed by the size of `A`. It can also contain custom number of outcomes if there are measurements with less outcomes. """ -function tensor_to_povm(Aax::Array{T,4}, o::Vector{Int64} = fill(size(Aax, 3), size(Aax, 4))) where {T} +function tensor_to_povm(Aax::Array{T, 4}, o::Vector{Int64} = fill(size(Aax, 3), size(Aax, 4))) where {T} return [[Hermitian(Aax[:, :, a, x]) for a in 1:o[x]] for x in axes(Aax, 4)] end export tensor_to_povm @@ -64,7 +64,7 @@ export tensor_to_povm Converts a matrix of (hermitian) matrices into a set of measurements in the common tensor format. """ -function povm_to_tensor(Axa::Vector{Measurement{T}}) where {T<:Number} +function povm_to_tensor(Axa::Vector{Measurement{T}}) where {T <: Number} d, o, m = _measurements_parameters(Axa) Aax = zeros(T, d, d, maximum(o), m) for x in eachindex(Axa) @@ -76,7 +76,7 @@ function povm_to_tensor(Axa::Vector{Measurement{T}}) where {T<:Number} end export povm_to_tensor -function _measurements_parameters(Axa::Vector{Measurement{T}}) where {T<:Number} +function _measurements_parameters(Axa::Vector{Measurement{T}}) where {T <: Number} @assert !isempty(Axa) # dimension on which the measurements act d = size(Axa[1][1], 1) @@ -93,11 +93,11 @@ _measurements_parameters(Aa::Measurement) = _measurements_parameters([Aa]) Checks if the measurement defined by A is valid (hermitian, semi-definite positive, and normalized). """ -function test_povm(E::Vector{<:AbstractMatrix{T}}) where {T<:Number} +function test_povm(E::Vector{<:AbstractMatrix{T}}) where {T <: Number} !all(ishermitian.(E)) && return false d = size(E[1], 1) !(sum(E) ≈ I(d)) && return false - for i = 1:length(E) + for i in 1:length(E) minimum(eigvals(E[i])) < -_rtol(T) && return false end return true @@ -111,10 +111,10 @@ Constructs the shift operator X of dimension `d` to the power `p`. Reference: [Generalized Clifford algebra](https://en.wikipedia.org/wiki/Generalized_Clifford_algebra) """ -function shift(::Type{T}, d::Integer, p::Integer = 1) where {T<:Number} +function shift(::Type{T}, d::Integer, p::Integer = 1) where {T <: Number} X = zeros(T, d, d) - for i in 0:d-1 - X[mod(i + p, d)+1, i+1] = 1 + for i in 0:(d - 1) + X[mod(i + p, d) + 1, i + 1] = 1 end return X end @@ -128,21 +128,21 @@ Constructs the clock operator Z of dimension `d` to the power `p`. Reference: [Generalized Clifford algebra](https://en.wikipedia.org/wiki/Generalized_Clifford_algebra) """ -function clock(::Type{T}, d::Integer, p::Integer = 1) where {T<:Number} +function clock(::Type{T}, d::Integer, p::Integer = 1) where {T <: Number} z = zeros(T, d) ω = _root_unity(T, d) - for i in 0:d-1 + for i in 0:(d - 1) exponent = mod(i * p, d) if exponent == 0 - z[i+1] = 1 + z[i + 1] = 1 elseif 4exponent == d - z[i+1] = im + z[i + 1] = im elseif 2exponent == d - z[i+1] = -1 + z[i + 1] = -1 elseif 4exponent == 3d - z[i+1] = -im + z[i + 1] = -im else - z[i+1] = ω^exponent + z[i + 1] = ω^exponent end end return Diagonal(z) @@ -160,17 +160,17 @@ Vectors of integers between 0 and 3 or strings of I, X, Y, Z automatically generate Kronecker products of the corresponding operators. """ -function pauli(::Type{T}, i::Integer) where {T<:Number} +function pauli(::Type{T}, i::Integer) where {T <: Number} return gell_mann(T, i ÷ 2 + 1, i % 2 + 1, 2) end -function pauli(::Type{T}, ind::Vector{<:Integer}) where {T<:Number} +function pauli(::Type{T}, ind::Vector{<:Integer}) where {T <: Number} if length(ind) == 1 return pauli(T, ind[1]) else return kron([pauli(T, i) for i in ind]...) end end -function pauli(::Type{T}, str::String) where {T<:Number} +function pauli(::Type{T}, str::String) where {T <: Number} ind = Int64[] for c in str if c in ['I', 'i', '1'] @@ -187,7 +187,7 @@ function pauli(::Type{T}, str::String) where {T<:Number} end return pauli(T, ind) end -pauli(::Type{T}, c::Char) where {T<:Number} = pauli(T, string(c)) +pauli(::Type{T}, c::Char) where {T <: Number} = pauli(T, string(c)) pauli(i::Integer) = pauli(ComplexF64, i) pauli(ind::Vector{<:Integer}) = pauli(ComplexF64, ind) pauli(str::String) = pauli(ComplexF64, str) @@ -202,7 +202,7 @@ Constructs the set `G` of generalized Gell-Mann matrices in dimension `d` such t Reference: [Generalizations of Pauli matrices](https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices) """ -function gell_mann(::Type{T}, d::Integer = 3) where {T<:Number} +function gell_mann(::Type{T}, d::Integer = 3) where {T <: Number} return [gell_mann(T, i, j, d) for j in 1:d, i in 1:d][:] # d=2, ρ = 1/2(σ0 + n*σ) # d=3, ρ = 1/3(I + √3 n*λ) @@ -221,7 +221,7 @@ Constructs the set `i`,`j`th Gell-Mann matrix of dimension `d`. Reference: [Generalizations of Pauli matrices](https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices) """ -function gell_mann(::Type{T}, i::Integer, j::Integer, d::Integer = 3) where {T<:Number} +function gell_mann(::Type{T}, i::Integer, j::Integer, d::Integer = 3) where {T <: Number} return gell_mann!(zeros(T, d, d), i, j, d) end gell_mann(i::Integer, j::Integer, d::Integer = 3) = gell_mann(ComplexF64, i, j, d) @@ -231,7 +231,7 @@ gell_mann(i::Integer, j::Integer, d::Integer = 3) = gell_mann(ComplexF64, i, j, In-place version of `gell_mann`. """ -function gell_mann!(res::AbstractMatrix{T}, i::Integer, j::Integer, d::Integer = 3) where {T<:Number} +function gell_mann!(res::AbstractMatrix{T}, i::Integer, j::Integer, d::Integer = 3) where {T <: Number} if i < j res[i, j] = 1 res[j, i] = 1 @@ -244,21 +244,21 @@ function gell_mann!(res::AbstractMatrix{T}, i::Integer, j::Integer, d::Integer = end elseif i == d tmp = _sqrt(T, 2) / _sqrt(T, d * (d - 1)) - for k in 1:d-1 + for k in 1:(d - 1) res[k, k] = tmp end res[d, d] = -(d - 1) * tmp else - gell_mann!(view(res, 1:d-1, 1:d-1), i, j, d - 1) + gell_mann!(view(res, 1:(d - 1), 1:(d - 1)), i, j, d - 1) end return res end export gell_mann! -_rtol(::Type{T}) where {T<:Number} = Base.rtoldefault(real(T)) +_rtol(::Type{T}) where {T <: Number} = Base.rtoldefault(real(T)) -_eps(::Type{T}) where {T<:Number} = _realeps(real(T)) -_realeps(::Type{T}) where {T<:AbstractFloat} = eps(T) +_eps(::Type{T}) where {T <: Number} = _realeps(real(T)) +_realeps(::Type{T}) where {T <: AbstractFloat} = eps(T) _realeps(::Type{<:Real}) = 0 """ @@ -266,13 +266,13 @@ _realeps(::Type{<:Real}) = 0 Zeroes out real or imaginary parts of `M` that are smaller than `tol`. """ -function cleanup!(M::AbstractArray{T}; tol = _eps(T)) where {T<:Number} +function cleanup!(M::AbstractArray{T}; tol = _eps(T)) where {T <: Number} wrapper = Base.typename(typeof(M)).wrapper cleanup!(parent(M); tol) return wrapper(M) end -function cleanup!(M::Array{T}; tol = _eps(T)) where {T<:Number} +function cleanup!(M::Array{T}; tol = _eps(T)) where {T <: Number} if isbitstype(T) M2 = reinterpret(real(T), M) #this is a no-op when T<:Real _cleanup!(M2; tol) @@ -288,7 +288,7 @@ end export cleanup! function _cleanup!(M; tol) - return M[abs.(M). tol) dec.U[:, 1:rank] end -_orthonormal_range_svd(A::AbstractMatrix; tol::Union{Real,Nothing} = nothing) = +_orthonormal_range_svd(A::AbstractMatrix; tol::Union{Real, Nothing} = nothing) = _orthonormal_range_svd!(deepcopy(A); tol = tol) -function _orthonormal_range_qr(A::SA.AbstractSparseMatrix{T,M}; tol::Union{Real,Nothing} = nothing) where {T<:Number,M} +function _orthonormal_range_qr(A::SA.AbstractSparseMatrix{T, M}; tol::Union{Real, Nothing} = nothing) where {T <: Number, M} dec = qr(A) tol = isnothing(tol) ? maximum(abs.(dec.R)) * _eps(T) : tol rank = sum(abs.(Diagonal(dec.R)) .> tol) @@ -325,10 +325,10 @@ otherwise uses an SVD and returns a dense matrix (`mode = 1`). Input `A` will be Tolerance `tol` is used to compute the rank and is automatically set if not provided. """ function orthonormal_range( - A::SA.AbstractMatrix{T}; - mode::Integer = -1, - tol::Union{Real,Nothing} = nothing -) where {T<:Number} + A::SA.AbstractMatrix{T}; + mode::Integer = -1, + tol::Union{Real, Nothing} = nothing + ) where {T <: Number} mode == 1 && SA.issparse(A) && throw(ArgumentError("SVD does not work with sparse matrices, use a dense matrix.")) mode == -1 && (mode = SA.issparse(A) ? 0 : 1) @@ -383,12 +383,12 @@ This function returns a generator, which can then be used e.g. in for loops with entire basis at once. If you need a vector, call `collect` on it. """ function n_body_basis( - n::Integer, - n_parties::Integer; - sb::AbstractVector{<:AbstractMatrix} = [pauli(1), pauli(2), pauli(3)], - sparse::Bool = true, - eye::AbstractMatrix = I(size(sb[1], 1)) -) + n::Integer, + n_parties::Integer; + sb::AbstractVector{<:AbstractMatrix} = [pauli(1), pauli(2), pauli(3)], + sparse::Bool = true, + eye::AbstractMatrix = I(size(sb[1], 1)) + ) (n >= 0 && n_parties >= 2) || throw(ArgumentError("Number of parties must be ≥ 2 and n ≥ 0.")) n <= n_parties || throw(ArgumentError("Number of parties cannot be larger than n.")) @@ -397,8 +397,8 @@ function n_body_basis( idx_eye = length(sb) ( Base.kron(sb[t]...) - for p in Combinatorics.with_replacement_combinations(1:nb, n) - for t in Combinatorics.multiset_permutations([p; repeat([idx_eye], n_parties - n)], n_parties) + for p in Combinatorics.with_replacement_combinations(1:nb, n) + for t in Combinatorics.multiset_permutations([p; repeat([idx_eye], n_parties - n)], n_parties) ) end export n_body_basis diff --git a/src/entanglement.jl b/src/entanglement.jl index 0325506..0b9dc1d 100644 --- a/src/entanglement.jl +++ b/src/entanglement.jl @@ -7,7 +7,7 @@ end """ schmidt_decomposition(ψ::AbstractVector, dims::AbstractVector{<:Integer} = _equal_sizes(ψ)) - + Produces the Schmidt decomposition of `ψ` with subsystem dimensions `dims`. If the argument `dims` is omitted equally-sized subsystems are assumed. Returns the (sorted) Schmidt coefficients λ and isometries U, V such that kron(U', V')*`ψ` is of Schmidt form. Reference: [Schmidt decomposition](https://en.wikipedia.org/wiki/Schmidt_decomposition). @@ -62,7 +62,7 @@ function entanglement_entropy(ρ::AbstractMatrix{T}, dims::AbstractVector = _equ JuMP.@variable(model, h) JuMP.@objective(model, Min, h / log(Rs(2))) - JuMP.@constraint(model, [h; σvec; ρvec] in Hypatia.EpiTrRelEntropyTriCone{Rs,Ts}(1 + 2 * vec_dim)) + JuMP.@constraint(model, [h; σvec; ρvec] in Hypatia.EpiTrRelEntropyTriCone{Rs, Ts}(1 + 2 * vec_dim)) JuMP.set_optimizer(model, Hypatia.Optimizer{Rs}) JuMP.set_silent(model) JuMP.optimize!(model) @@ -98,8 +98,8 @@ function _test_entanglement_entropy_qubit(h, ρ, σ) R = typeof(h) λ, U = eigen(σ) g = zeros(R, 4, 4) - for j = 1:4 - for i = 1:j-1 + for j in 1:4 + for i in 1:(j - 1) g[i, j] = (λ[i] - λ[j]) / log(λ[i] / λ[j]) end g[j, j] = λ[j] @@ -109,8 +109,8 @@ function _test_entanglement_entropy_qubit(h, ρ, σ) λ2, U2 = eigen(σT) phi = partial_transpose(ketbra(U2[:, 1]), 2, [2, 2]) G = zero(U) - for i = 1:4 - for j = 1:4 + for i in 1:4 + for j in 1:4 G += g[i, j] * ketbra(U[:, i]) * phi * ketbra(U[:, j]) end end @@ -138,21 +138,21 @@ If a state ``ρ`` with local dimensions ``d_A`` and ``d_B`` has Schmidt number ` a PSD matrix ``ω`` in the extended space ``AA′B′B``, where ``A′`` and ``B^′`` have dimension ``s``, such that ``ω / s`` is separable against ``AA′|B′B`` and ``Π† ω Π = ρ``, where ``Π = 1_A ⊗ s ψ^+ ⊗ 1_B``, and ``ψ^+`` is a non-normalized maximally entangled state. Separabiity is tested with the DPS hierarchy, -with `n` controlling the how many copies of the ``B′B`` subsystem are used. +with `n` controlling the how many copies of the ``B′B`` subsystem are used. References: Hulpke, Bruss, Lewenstein, Sanpera [arXiv:quant-ph/0401118](https://arxiv.org/abs/quant-ph/0401118)\ - Weilenmann, Dive, Trillo, Aguilar, Navascués [arXiv:1912.10056](https://arxiv.org/abs/1912.10056) +Weilenmann, Dive, Trillo, Aguilar, Navascués [arXiv:1912.10056](https://arxiv.org/abs/1912.10056) """ function schmidt_number( - ρ::AbstractMatrix{T}, - s::Integer = 2, - dims::AbstractVector{<:Integer} = _equal_sizes(ρ), - n::Integer = 1; - ppt::Bool = true, - verbose::Bool = false, - solver = Hypatia.Optimizer{_solver_type(T)} -) where {T<:Number} + ρ::AbstractMatrix{T}, + s::Integer = 2, + dims::AbstractVector{<:Integer} = _equal_sizes(ρ), + n::Integer = 1; + ppt::Bool = true, + verbose::Bool = false, + solver = Hypatia.Optimizer{_solver_type(T)} + ) where {T <: Number} ishermitian(ρ) || throw(ArgumentError("State must be Hermitian")) s >= 1 || throw(ArgumentError("Schmidt number must be ≥ 1")) if s == 1 @@ -198,13 +198,13 @@ export schmidt_number Lower bounds the random robustness of state `ρ` with subsystem dimensions `dims` using level `n` of the DPS hierarchy. Argument `ppt` indicates whether to include the partial transposition constraints. """ function random_robustness( - ρ::AbstractMatrix{T}, - dims::AbstractVector{<:Integer} = _equal_sizes(ρ), - n::Integer = 1; - ppt::Bool = true, - verbose::Bool = false, - solver = Hypatia.Optimizer{_solver_type(T)} -) where {T<:Number} + ρ::AbstractMatrix{T}, + dims::AbstractVector{<:Integer} = _equal_sizes(ρ), + n::Integer = 1; + ppt::Bool = true, + verbose::Bool = false, + solver = Hypatia.Optimizer{_solver_type(T)} + ) where {T <: Number} ishermitian(ρ) || throw(ArgumentError("State must be Hermitian")) is_complex = (T <: Complex) @@ -240,14 +240,14 @@ References: Doherty, Parrilo, Spedalieri [arXiv:quant-ph/0308032](https://arxiv.org/abs/quant-ph/0308032) """ function _dps_constraints!( - model::JuMP.GenericModel{T}, - ρ::AbstractMatrix, - dims::AbstractVector{<:Integer}, - n::Integer; - ppt::Bool = true, - is_complex::Bool = true, - isometry::AbstractMatrix = I(size(ρ, 1)) -) where {T} + model::JuMP.GenericModel{T}, + ρ::AbstractMatrix, + dims::AbstractVector{<:Integer}, + n::Integer; + ppt::Bool = true, + is_complex::Bool = true, + isometry::AbstractMatrix = I(size(ρ, 1)) + ) where {T} ishermitian(ρ) || throw(ArgumentError("State must be Hermitian")) dA, dB = dims @@ -267,11 +267,11 @@ function _dps_constraints!( JuMP.@variable(model, symmetric_meat[1:d, 1:d] in psd_cone) lifted = wrapper(V * symmetric_meat * V') - JuMP.@expression(model, reduced, partial_trace(lifted, 3:n+1, ext_dims)) + JuMP.@expression(model, reduced, partial_trace(lifted, 3:(n + 1), ext_dims)) JuMP.@constraint(model, witness_constraint, ρ == wrapper(isometry' * reduced * isometry)) if ppt - for i in 2:n+1 + for i in 2:(n + 1) JuMP.@constraint(model, partial_transpose(lifted, 2:i, ext_dims) in psd_cone) end end @@ -282,7 +282,7 @@ function _fully_decomposable_witness_constraints!(model, dims, W) biparts = Combinatorics.partitions(1:nparties, 2) dim = prod(dims) - Ps = [JuMP.@variable(model, [1:dim,1:dim] in JuMP.HermitianPSDCone()) for _ in 1:length(biparts)] + Ps = [JuMP.@variable(model, [1:dim, 1:dim] in JuMP.HermitianPSDCone()) for _ in 1:length(biparts)] JuMP.@constraint(model, tr(W) == 1) # this can be used instead of tr(W) = 1 if we want a GME entanglement quantifier (see ref.) @@ -319,17 +319,17 @@ with the PPT criterion. If the state is a PPT mixture, returns a 0 matrix instea Reference: Jungnitsch, Moroder, Guehne [arXiv:quant-ph/0401118](https://arxiv.org/abs/quant-ph/0401118) """ function ppt_mixture( - ρ::AbstractMatrix{T}, - dims::AbstractVector{<:Integer}; - verbose::Bool = false, - solver = Hypatia.Optimizer{_solver_type(T)} -) where {T <: Number} + ρ::AbstractMatrix{T}, + dims::AbstractVector{<:Integer}; + verbose::Bool = false, + solver = Hypatia.Optimizer{_solver_type(T)} + ) where {T <: Number} dim = size(ρ, 1) prod(dims) == size(ρ, 1) || throw(ArgumentError("State dimension does not agree with local dimensions.")) model = JuMP.GenericModel{_solver_type(T)}() - JuMP.@variable(model, W[1:dim,1:dim], Hermitian) + JuMP.@variable(model, W[1:dim, 1:dim], Hermitian) _fully_decomposable_witness_constraints!(model, dims, W) _min_dotprod!(model, ρ, W, solver, verbose) @@ -356,7 +356,7 @@ defining such a witness. More precisely, if a list of observables ``O_i`` is provided in the parameter `obs`, the witness will be of the form ``∑_i α_i O_i`` and detects ρ only using these observables. For example, using only two-body operators (and lower order) -one can call +one can call ```julia-repl julia> two_body_basis = collect(Iterators.flatten(n_body_basis(i, 3) for i in 0:2)) @@ -366,18 +366,18 @@ julia> ppt_mixture(state_ghz(2, 3), [2, 2, 2], two_body_basis) Reference: Jungnitsch, Moroder, Guehne [arXiv:quant-ph/0401118](https://arxiv.org/abs/quant-ph/0401118) """ function ppt_mixture( - ρ::AbstractMatrix{T}, - dims::AbstractVector{<:Integer}, - obs::AbstractVector{<:AbstractMatrix}; - verbose::Bool = false, - solver = Hypatia.Optimizer{_solver_type(T)} -) where {T <: Number} + ρ::AbstractMatrix{T}, + dims::AbstractVector{<:Integer}, + obs::AbstractVector{<:AbstractMatrix}; + verbose::Bool = false, + solver = Hypatia.Optimizer{_solver_type(T)} + ) where {T <: Number} dim = size(ρ, 1) prod(dims) == size(ρ, 1) || throw(ArgumentError("State dimension does not agree with local dimensions.")) model = JuMP.GenericModel{_solver_type(T)}() - JuMP.@variable(model, w_coeffs[1:length(obs)]) + JuMP.@variable(model, w_coeffs[1:length(obs)]) W = sum(w_coeffs[i] * obs[i] for i in eachindex(w_coeffs)) _fully_decomposable_witness_constraints!(model, dims, W) diff --git a/src/entropy.jl b/src/entropy.jl index 1c33385..daa7d58 100644 --- a/src/entropy.jl +++ b/src/entropy.jl @@ -7,7 +7,7 @@ Computes the (quantum) relative entropy tr(`ρ` (log `ρ` - log `σ`)) between p Reference: [Quantum relative entropy](https://en.wikipedia.org/wiki/Quantum_relative_entropy). """ -function relative_entropy(base::Real, ρ::AbstractMatrix{T}, σ::AbstractMatrix{T}) where {T<:Number} +function relative_entropy(base::Real, ρ::AbstractMatrix{T}, σ::AbstractMatrix{T}) where {T <: Number} if size(ρ) != size(σ) throw(ArgumentError("ρ and σ have the same size.")) end @@ -24,7 +24,7 @@ function relative_entropy(base::Real, ρ::AbstractMatrix{T}, σ::AbstractMatrix{ logσ_λ = _log.(Ref(base), σ_λ) d = size(ρ, 1) h = zero(real(T)) - @inbounds for j = 1:d, i = 1:d + @inbounds for j in 1:d, i in 1:d h += ρ_λ[i] * (logρ_λ[i] - logσ_λ[j]) * m[i, j] end return h @@ -39,7 +39,7 @@ Computes the relative entropy D(`p`||`q`) = Σᵢpᵢlog(pᵢ/qᵢ) between two Reference: [Relative entropy](https://en.wikipedia.org/wiki/Relative_entropy). """ -function relative_entropy(base::Real, p::AbstractVector{T}, q::AbstractVector{T}) where {T<:Real} +function relative_entropy(base::Real, p::AbstractVector{T}, q::AbstractVector{T}) where {T <: Real} if length(p) != length(q) throw(ArgumentError("`p` and q must have the same length.")) end @@ -48,7 +48,7 @@ function relative_entropy(base::Real, p::AbstractVector{T}, q::AbstractVector{T} end logp = _log.(Ref(base), p) logq = _log.(Ref(base), q) - h = sum(p[i] * (logp[i] - logq[i]) for i = 1:length(p)) + h = sum(p[i] * (logp[i] - logq[i]) for i in 1:length(p)) return h end relative_entropy(p::AbstractVector, q::AbstractVector) = relative_entropy(2, p, q) @@ -79,7 +79,7 @@ function entropy(base::Real, ρ::AbstractMatrix{T}) where {T <: Number} if any(λ .< -_rtol(T)) throw(ArgumentError("ρ must be positive semidefinite.")) end - h = -sum(λ[i] * _log(base, λ[i]) for i = 1:size(ρ, 1)) + h = -sum(λ[i] * _log(base, λ[i]) for i in 1:size(ρ, 1)) return h end entropy(ρ::AbstractMatrix) = entropy(2, ρ) @@ -92,11 +92,11 @@ Computes the Shannon entropy -Σᵢpᵢlog(pᵢ) of a non-negative vector `p` us Reference: [Entropy (information theory)](https://en.wikipedia.org/wiki/Entropy_(information_theory)). """ -function entropy(base::Real, p::AbstractVector{T}) where {T<:Real} +function entropy(base::Real, p::AbstractVector{T}) where {T <: Real} if any(p .< T(0)) throw(ArgumentError("p must be non-negative.")) end - h = -sum(p[i] * _log(base, p[i]) for i = 1:length(p)) + h = -sum(p[i] * _log(base, p[i]) for i in 1:length(p)) return h end entropy(p::AbstractVector) = entropy(2, p) @@ -120,14 +120,14 @@ Computes the conditional Shannon entropy H(A|B) of the joint probability distrib Reference: [Conditional entropy](https://en.wikipedia.org/wiki/Conditional_entropy). """ -function conditional_entropy(base::Real, pAB::AbstractMatrix{T}) where {T<:Real} +function conditional_entropy(base::Real, pAB::AbstractMatrix{T}) where {T <: Real} nA, nB = size(pAB) if any(pAB .< T(0)) throw(ArgumentError("pAB must be non-negative.")) end h = T(0) pB = sum(pAB; dims = 1) - for a = 1:nA, b = 1:nB + for a in 1:nA, b in 1:nB h -= pAB[a, b] * _log(base, pAB[a, b] / pB[b]) end return h @@ -142,17 +142,17 @@ Computes the conditional von Neumann entropy of `rho` with subsystem dimensions Reference: [Conditional quantum entropy](https://en.wikipedia.org/wiki/Conditional_quantum_entropy). """ function conditional_entropy( - base::Real, - rho::AbstractMatrix, - csys::AbstractVector{<:Integer}, - dims::AbstractVector{<:Integer} -) + base::Real, + rho::AbstractMatrix, + csys::AbstractVector{<:Integer}, + dims::AbstractVector{<:Integer} + ) isempty(csys) && return entropy(base, rho) length(csys) == length(dims) && return zero(real(eltype(rho))) remove = Vector{eltype(csys)}(undef, length(dims) - length(csys)) # To condition on csys we trace out the rest counter = 0 - for i = 1:length(dims) + for i in 1:length(dims) if !(i in csys) counter += 1 remove[counter] = i diff --git a/src/games.jl b/src/games.jl index 718f1b2..681ef31 100644 --- a/src/games.jl +++ b/src/games.jl @@ -14,9 +14,9 @@ function chsh(::Type{T}, d::Integer = 2) where {T} element = inv(T(d^2)) end - for a = 0:d-1, b = 0:d-1, x = 0:d-1, y = 0:d-1 + for a in 0:(d - 1), b in 0:(d - 1), x in 0:(d - 1), y in 0:(d - 1) if mod(a + b + x * y, d) == 0 - G[a+1, b+1, x+1, y+1] = element + G[a + 1, b + 1, x + 1, y + 1] = element end end @@ -41,9 +41,9 @@ function cglmp(::Type{T}, d::Integer = 3) where {T} normalization = inv(T(4 * (d - 1))) end - for a = 0:d-1, b = 0:d-1, x = 0:1, y = 0:1, k = 0:d-2 + for a in 0:(d - 1), b in 0:(d - 1), x in 0:1, y in 0:1, k in 0:(d - 2) if mod(a - b, d) == mod((-1)^mod(x + y, 2) * k + x * y, d) - G[a+1, b+1, x+1, y+1] = normalization * (d - 1 - k) + G[a + 1, b + 1, x + 1, y + 1] = normalization * (d - 1 - k) end end @@ -60,19 +60,19 @@ inn22 Bell functional in Collins-Gisin notation. Local bound 1. Reference: Śliwa, [arXiv:quant-ph/0305190](https://arxiv.org/abs/quant-ph/0305190) """ function inn22(::Type{T}, n) where {T} - C = zeros(T, n+1,n+1) - for x = 1:n - for y = 1:n - if x+y <= n + 1 - C[x+1,y+1] = -1 - elseif x+y == n + 2 - C[x+1,y+1] = 1 + C = zeros(T, n + 1, n + 1) + for x in 1:n + for y in 1:n + if x + y <= n + 1 + C[x + 1, y + 1] = -1 + elseif x + y == n + 2 + C[x + 1, y + 1] = 1 end end end - C[1,2] = 1 - C[2,1] = 1 - return C + C[1, 2] = 1 + C[2, 1] = 1 + return C end inn22(n::Integer = 3) = inn22(Int, n) export inn22 diff --git a/src/incompatibility.jl b/src/incompatibility.jl index 99f3f41..482604d 100644 --- a/src/incompatibility.jl +++ b/src/incompatibility.jl @@ -14,12 +14,12 @@ Returns the parent POVM if `return_parent = true`. Reference: Designolle, Farkas, Kaniewski, [arXiv:1906.00448](https://arxiv.org/abs/1906.00448) """ function incompatibility_robustness( - A::Vector{Measurement{T}}, - measure::String = "g"; - verbose = false, - return_parent = false, - solver = Hypatia.Optimizer{_solver_type(T)} -) where {T<:Number} + A::Vector{Measurement{T}}, + measure::String = "g"; + verbose = false, + return_parent = false, + solver = Hypatia.Optimizer{_solver_type(T)} + ) where {T <: Number} @assert measure in ["d", "r", "p", "jm", "g"] d, o, m = _measurements_parameters(A) is_complex = T <: Complex @@ -44,8 +44,8 @@ function incompatibility_robustness( end # constraints - lhs = zero(JuMP.GenericAffExpr{stT,JuMP.GenericVariableRef{stT}}) - rhs = zero(JuMP.GenericAffExpr{stT,JuMP.GenericVariableRef{stT}}) + lhs = zero(JuMP.GenericAffExpr{stT, JuMP.GenericVariableRef{stT}}) + rhs = zero(JuMP.GenericAffExpr{stT, JuMP.GenericVariableRef{stT}}) if measure in ["d", "r", "p"] con = JuMP.@constraint(model, [j in CartesianIndices(o)], sum(X[x][j.I[x]] for x in 1:m) in cone) JuMP.add_to_expression!(lhs, 1) @@ -69,7 +69,7 @@ function incompatibility_robustness( JuMP.@constraint(model, X[x][a] in cone) end end - if measure == "p" + if measure == "p" JuMP.add_to_expression!(rhs, ξ[x]) end end @@ -112,7 +112,7 @@ Computes the incompatibility depolarizing robustness of the measurements in the Reference: Designolle, Farkas, Kaniewski, [arXiv:1906.00448](https://arxiv.org/abs/1906.00448) """ -function incompatibility_robustness_depolarizing(A::Vector{Measurement{T}}; kwargs...) where {T<:Number} +function incompatibility_robustness_depolarizing(A::Vector{Measurement{T}}; kwargs...) where {T <: Number} return incompatibility_robustness(A, "d"; kwargs...) end export incompatibility_robustness_depolarizing @@ -124,7 +124,7 @@ Computes the incompatibility random robustness of the measurements in the vector Reference: Designolle, Farkas, Kaniewski, [arXiv:1906.00448](https://arxiv.org/abs/1906.00448) """ -function incompatibility_robustness_random(A::Vector{Measurement{T}}; kwargs...) where {T<:Number} +function incompatibility_robustness_random(A::Vector{Measurement{T}}; kwargs...) where {T <: Number} return incompatibility_robustness(A, "r"; kwargs...) end export incompatibility_robustness_random @@ -136,7 +136,7 @@ Computes the incompatibility probabilistic robustness of the measurements in the Reference: Designolle, Farkas, Kaniewski, [arXiv:1906.00448](https://arxiv.org/abs/1906.00448) """ -function incompatibility_robustness_probabilistic(A::Vector{Measurement{T}}; kwargs...) where {T<:Number} +function incompatibility_robustness_probabilistic(A::Vector{Measurement{T}}; kwargs...) where {T <: Number} return incompatibility_robustness(A, "p"; kwargs...) end export incompatibility_robustness_probabilistic @@ -148,7 +148,7 @@ Computes the incompatibility jointly measurable robustness of the measurements i Reference: Designolle, Farkas, Kaniewski, [arXiv:1906.00448](https://arxiv.org/abs/1906.00448) """ -function incompatibility_robustness_jointly_measurable(A::Vector{Measurement{T}}; kwargs...) where {T<:Number} +function incompatibility_robustness_jointly_measurable(A::Vector{Measurement{T}}; kwargs...) where {T <: Number} return incompatibility_robustness(A, "jm"; kwargs...) end export incompatibility_robustness_jointly_measurable @@ -160,7 +160,7 @@ Computes the incompatibility generalized robustness of the measurements in the v Reference: Designolle, Farkas, Kaniewski, [arXiv:1906.00448](https://arxiv.org/abs/1906.00448) """ -function incompatibility_robustness_generalized(A::Vector{Measurement{T}}; kwargs...) where {T<:Number} +function incompatibility_robustness_generalized(A::Vector{Measurement{T}}; kwargs...) where {T <: Number} return incompatibility_robustness(A, "g"; kwargs...) end export incompatibility_robustness_generalized diff --git a/src/measurements.jl b/src/measurements.jl index cd0295f..d0d26e4 100644 --- a/src/measurements.jl +++ b/src/measurements.jl @@ -1,30 +1,30 @@ -_root_unity(::Type{T}, n::Integer) where {T<:Number} = exp(2 * im * real(T)(π) / n) -_sqrt(::Type{T}, n::Integer) where {T<:Number} = sqrt(real(T)(n)) +_root_unity(::Type{T}, n::Integer) where {T <: Number} = exp(2 * im * real(T)(π) / n) +_sqrt(::Type{T}, n::Integer) where {T <: Number} = sqrt(real(T)(n)) # MUBs -function mub_prime(::Type{T}, p::Integer) where {T<:Number} +function mub_prime(::Type{T}, p::Integer) where {T <: Number} γ = _root_unity(T, p) inv_sqrt_p = inv(_sqrt(T, p)) - B = [Matrix{T}(undef, p, p) for _ in 1:p+1] + B = [Matrix{T}(undef, p, p) for _ in 1:(p + 1)] B[1] .= I(p) if p == 2 B[2] .= [1 1; 1 -1] .* inv_sqrt_p B[3] .= [1 1; im -im] .* inv_sqrt_p else - for k in 0:p-1 - fill!(B[k+2], inv_sqrt_p) - for t in 0:p-1, j in 0:p-1 + for k in 0:(p - 1) + fill!(B[k + 2], inv_sqrt_p) + for t in 0:(p - 1), j in 0:(p - 1) exponent = mod(j * (t + k * j), p) if exponent == 0 continue elseif 4exponent == p - B[k+2][j+1, t+1] *= im + B[k + 2][j + 1, t + 1] *= im elseif 2exponent == p - B[k+2][j+1, t+1] *= -1 + B[k + 2][j + 1, t + 1] *= -1 elseif 4exponent == 3p - B[k+2][j+1, t+1] *= -im + B[k + 2][j + 1, t + 1] *= -im else - B[k+2][j+1, t+1] *= γ^exponent + B[k + 2][j + 1, t + 1] *= γ^exponent end end end @@ -33,29 +33,29 @@ function mub_prime(::Type{T}, p::Integer) where {T<:Number} end mub_prime(p::Integer) = mub_prime(ComplexF64, p) -function mub_prime_power(::Type{T}, p::Integer, r::Integer) where {T<:Number} +function mub_prime_power(::Type{T}, p::Integer, r::Integer) where {T <: Number} d = Int64(p^r) γ = _root_unity(T, p) inv_sqrt_d = inv(_sqrt(T, d)) - B = [zeros(T, d, d) for _ in 1:d+1] + B = [zeros(T, d, d) for _ in 1:(d + 1)] B[1] .= I(d) f, x = Nemo.finite_field(p, r, "x") - pow = [x^i for i in 0:r-1] - el = [sum(digits(i; base = p, pad = r) .* pow) for i in 0:d-1] + pow = [x^i for i in 0:(r - 1)] + el = [sum(digits(i; base = p, pad = r) .* pow) for i in 0:(d - 1)] if p == 2 - for i in 1:d, k in 0:d-1, q in 0:d-1 + for i in 1:d, k in 0:(d - 1), q in 0:(d - 1) aux = one(T) q_bin = digits(q; base = 2, pad = r) - for m in 0:r-1, n in 0:r-1 - aux *= conj(im^_tr_ff(el[i] * el[q_bin[m+1]*2^m+1] * el[q_bin[n+1]*2^n+1])) + for m in 0:(r - 1), n in 0:(r - 1) + aux *= conj(im^_tr_ff(el[i] * el[q_bin[m + 1] * 2^m + 1] * el[q_bin[n + 1] * 2^n + 1])) end - B[i+1][:, k+1] += (-1)^_tr_ff(el[q+1] * el[k+1]) * aux * B[1][:, q+1] * inv_sqrt_d + B[i + 1][:, k + 1] += (-1)^_tr_ff(el[q + 1] * el[k + 1]) * aux * B[1][:, q + 1] * inv_sqrt_d end else inv_two = inv(2 * one(f)) - for i in 1:d, k in 0:d-1, q in 0:d-1 - B[i+1][:, k+1] += - γ^_tr_ff(-el[q+1] * el[k+1]) * γ^_tr_ff(el[i] * el[q+1] * el[q+1] * inv_two) * B[1][:, q+1] * inv_sqrt_d + for i in 1:d, k in 0:(d - 1), q in 0:(d - 1) + B[i + 1][:, k + 1] += + γ^_tr_ff(-el[q + 1] * el[k + 1]) * γ^_tr_ff(el[i] * el[q + 1] * el[q + 1] * inv_two) * B[1][:, q + 1] * inv_sqrt_d end end return B @@ -75,7 +75,7 @@ The output contains 1+minᵢ pᵢ^rᵢ bases, where `d` = ∏ᵢ pᵢ^rᵢ. Reference: Durt, Englert, Bengtsson, Życzkowski, [arXiv:1004.3348](https://arxiv.org/abs/1004.3348). """ -function mub(::Type{T}, d::Integer) where {T<:Number} +function mub(::Type{T}, d::Integer) where {T <: Number} # the dimension d can be any integer greater than two @assert d ≥ 2 f = collect(Nemo.factor(Int64(d))) # Nemo.factor requires d to be an Int64 (or UInt64) @@ -100,7 +100,7 @@ mub(d::Integer) = mub(ComplexF64, d) export mub # Select a specific subset with k bases -function mub(::Type{T}, d::Integer, k::Integer, s::Integer = 1) where {T<:Number} +function mub(::Type{T}, d::Integer, k::Integer, s::Integer = 1) where {T <: Number} B = mub(T, d) subs = collect(Iterators.take(Combinatorics.combinations(1:length(B), k), s)) sub = subs[end] @@ -113,7 +113,7 @@ mub(d::Integer, k::Integer, s::Integer = 1) = mub(ComplexF64, d, k, s) Checks if the input bases are mutually unbiased. """ -function test_mub(B::Vector{Matrix{T}}) where {T<:Number} +function test_mub(B::Vector{Matrix{T}}) where {T <: Number} d = size(B[1], 1) k = length(B) inv_d = inv(T(d)) @@ -145,11 +145,11 @@ Reference: Appleby, Yadsan-Appleby, Zauner, [arXiv:1209.1813](http://arxiv.org/a function sic_povm(::Type{T}, d::Integer) where {T} fiducial = _fiducial_WH(real(T), d) vecs = Vector{Vector{T}}(undef, d^2) - for p in 0:d-1 + for p in 0:(d - 1) Xp = shift(T, d, p) - for q in 0:d-1 + for q in 0:(d - 1) Zq = clock(T, d, q) - vecs[d*p+q+1] = Xp * Zq * fiducial + vecs[d * p + q + 1] = Xp * Zq * fiducial end end sqrt_d = _sqrt(T, d) @@ -166,12 +166,12 @@ export sic_povm Checks if `vecs` is a vector of `d²` vectors |vᵢ⟩ such that |vᵢ⟩⟨vᵢ| forms a SIC-POVM of dimension `d`. """ -function test_sic(vecs::Vector{Vector{T}}) where {T<:Number} +function test_sic(vecs::Vector{Vector{T}}) where {T <: Number} d = length(vecs[1]) length(vecs) == d^2 || throw(ArgumentError("Number of vectors must be d² = $(d^2), got $(length(vecs)).")) normalization = inv(T(d^2)) symmetry = inv(T(d^2 * (d + 1))) - for j in 1:d^2, i in 1:j + for j in 1:(d^2), i in 1:j inner_product = abs2(dot(vecs[i], vecs[j])) if i == j deviation = abs2(inner_product - normalization) @@ -191,11 +191,11 @@ export test_sic Does the Naimark dilation of a rank-1 POVM given as a vector of vectors. This is the minimal dilation. """ -function dilate_povm(vecs::Vector{Vector{T}}) where {T<:Number} +function dilate_povm(vecs::Vector{Vector{T}}) where {T <: Number} d = length(vecs[1]) n = length(vecs) V = Matrix{T}(undef, n, d) - for j = 1:d + for j in 1:d for i in 1:n V[i, j] = conj(vecs[i][j]) end @@ -236,12 +236,12 @@ function _fiducial_WH(::Type{T}, d::Integer) where {T} return [ T(1) / 40 * (-a + 5) * r + T(1) / 20 * (-a + 5), ((-T(1) / 40 * a * r - T(1) / 40 * a) * b + (T(1) / 80 * (a - 5) * r + T(1) / 40 * (a - 5))) * im + - -T(1) / 40 * a * b + - T(1) / 80 * (-a + 5) * r, + -T(1) / 40 * a * b + + T(1) / 80 * (-a + 5) * r, T(1) / 40 * (a - 5) * r * im, ((T(1) / 40 * a * r + T(1) / 40 * a) * b + (T(1) / 80 * (a - 5) * r + T(1) / 40 * (a - 5))) * im + - T(1) / 40 * a * b + - T(1) / 80 * (-a + 5) * r + T(1) / 40 * a * b + + T(1) / 80 * (-a + 5) * r, ] elseif d == 5 a = sqrt(T(3)) @@ -252,46 +252,46 @@ function _fiducial_WH(::Type{T}, d::Integer) where {T} T(1) / 60 * (-a + 3) * r + T(1) / 12 * (-a + 3), ( ((T(1) / 120 * a * r + T(1) / 40 * a) * t + (-T(1) / 240 * a * r + T(1) / 240 * a)) * b + - (T(1) / 120 * (a - 3) * r * t + (T(1) / 32 * (-a + 1) * r + T(1) / 32 * (-a + 1))) + (T(1) / 120 * (a - 3) * r * t + (T(1) / 32 * (-a + 1) * r + T(1) / 32 * (-a + 1))) ) * im + - ((-T(1) / 120 * a * r - T(1) / 120 * a) * t - T(1) / 120 * a) * b + - T(1) / 40 * (a - 1) * r * t + - T(1) / 160 * (-a + 3) * r + - T(1) / 96 * (a - 3), + ((-T(1) / 120 * a * r - T(1) / 120 * a) * t - T(1) / 120 * a) * b + + T(1) / 40 * (a - 1) * r * t + + T(1) / 160 * (-a + 3) * r + + T(1) / 96 * (a - 3), ( ( (-T(1) / 80 * a * r + T(1) / 240 * (-7 * a + 6)) * t + - (T(1) / 480 * (-a + 9) * r + T(1) / 480 * (-a + 15)) + (T(1) / 480 * (-a + 9) * r + T(1) / 480 * (-a + 15)) ) * b + (T(1) / 120 * (-a + 3) * r * t + (T(1) / 32 * (-a + 1) * r + T(1) / 32 * (-a + 1))) ) * im + - ( + ( (T(1) / 240 * (2 * a + 3) * r + T(1) / 240 * (4 * a + 3)) * t + - (T(1) / 480 * (-a - 3) * r + T(1) / 160 * (-a - 5)) + (T(1) / 480 * (-a - 3) * r + T(1) / 160 * (-a - 5)) ) * b + - T(1) / 40 * (-a + 1) * r * t + - T(1) / 160 * (-a + 3) * r + - T(1) / 96 * (a - 3), + T(1) / 40 * (-a + 1) * r * t + + T(1) / 160 * (-a + 3) * r + + T(1) / 96 * (a - 3), ( ( (T(1) / 80 * a * r + T(1) / 240 * (7 * a - 6)) * t + - (T(1) / 480 * (a - 9) * r + T(1) / 480 * (a - 15)) + (T(1) / 480 * (a - 9) * r + T(1) / 480 * (a - 15)) ) * b + (T(1) / 120 * (-a + 3) * r * t + (T(1) / 32 * (-a + 1) * r + T(1) / 32 * (-a + 1))) ) * im + - ( + ( (T(1) / 240 * (-2 * a - 3) * r + T(1) / 240 * (-4 * a - 3)) * t + - (T(1) / 480 * (a + 3) * r + T(1) / 160 * (a + 5)) + (T(1) / 480 * (a + 3) * r + T(1) / 160 * (a + 5)) ) * b + - T(1) / 40 * (-a + 1) * r * t + - T(1) / 160 * (-a + 3) * r + - T(1) / 96 * (a - 3), + T(1) / 40 * (-a + 1) * r * t + + T(1) / 160 * (-a + 3) * r + + T(1) / 96 * (a - 3), ( ((-T(1) / 120 * a * r - T(1) / 40 * a) * t + (T(1) / 240 * a * r - T(1) / 240 * a)) * b + - (T(1) / 120 * (a - 3) * r * t + (T(1) / 32 * (-a + 1) * r + T(1) / 32 * (-a + 1))) + (T(1) / 120 * (a - 3) * r * t + (T(1) / 32 * (-a + 1) * r + T(1) / 32 * (-a + 1))) ) * im + - ((T(1) / 120 * a * r + T(1) / 120 * a) * t + T(1) / 120 * a) * b + - T(1) / 40 * (a - 1) * r * t + - T(1) / 160 * (-a + 3) * r + - T(1) / 96 * (a - 3) + ((T(1) / 120 * a * r + T(1) / 120 * a) * t + T(1) / 120 * a) * b + + T(1) / 40 * (a - 1) * r * t + + T(1) / 160 * (-a + 3) * r + + T(1) / 96 * (a - 3), ] elseif d == 6 a = sqrt(T(21)) @@ -300,89 +300,89 @@ function _fiducial_WH(::Type{T}, d::Integer) where {T} return [ ( T(1) / 1008 * (-a + 3) * r[2] * b[1] * b[2]^2 + - T(1) / 504 * (a - 6) * r[2] * b[1] * b[2] + - (T(1) / 168 * (a - 2) * r[2] - T(1) / 168 * a) * b[1] + T(1) / 504 * (a - 6) * r[2] * b[1] * b[2] + + (T(1) / 168 * (a - 2) * r[2] - T(1) / 168 * a) * b[1] ) * im + - T(1) / 504 * (-a + 5) * r[2] * b[2]^2 + - T(1) / 504 * (a + 1) * r[2] * b[2] + - T(1) / 504 * (-a - 13) * r[2] + - T(1) / 168 * (a + 21), + T(1) / 504 * (-a + 5) * r[2] * b[2]^2 + + T(1) / 504 * (a + 1) * r[2] * b[2] + + T(1) / 504 * (-a - 13) * r[2] + + T(1) / 168 * (a + 21), ( (T(1) / 2016 * (-a + 3) * r[2] * b[1] + (T(1) / 504 * (a - 7) * r[2] + T(1) / 336 * (a - 5))) * b[2]^2 + - ( + ( (T(1) / 1008 * (a - 6) * r[2] + T(1) / 672 * (-a + 7)) * b[1] + - (T(1) / 1008 * (-a - 7) * r[2] + T(1) / 336 * (-a - 1)) + (T(1) / 1008 * (-a - 7) * r[2] + T(1) / 336 * (-a - 1)) ) * b[2] + - ( + ( (T(1) / 672 * (a + 3) * r[2] + T(1) / 672 * (-3 * a + 7)) * b[1] + - (T(1) / 504 * (-3 * a + 14) * r[2] + T(1) / 168 * (-a - 4)) + (T(1) / 504 * (-3 * a + 14) * r[2] + T(1) / 168 * (-a - 4)) ) ) * im + - ( + ( (T(1) / 1008 * (-a + 3) * r[2] + T(1) / 672 * (a - 3)) * b[1] + - (T(1) / 1008 * (-3 * a + 7) * r[2] - T(1) / 84) + (T(1) / 1008 * (-3 * a + 7) * r[2] - T(1) / 84) ) * b[2]^2 + - ( + ( (T(1) / 2016 * (a - 3) * r[2] - T(1) / 336) * b[1] + - (T(1) / 1008 * (5 * a - 7) * r[2] + T(1) / 336 * (a - 5)) + (T(1) / 1008 * (5 * a - 7) * r[2] + T(1) / 336 * (a - 5)) ) * b[2] + - (T(1) / 224 * (a - 5) * r[2] + T(1) / 672 * (-7 * a + 19)) * b[1] + - T(1) / 72 * (a - 4) * r[2] + - T(1) / 168 * (-a + 22), + (T(1) / 224 * (a - 5) * r[2] + T(1) / 672 * (-7 * a + 19)) * b[1] + + T(1) / 72 * (a - 4) * r[2] + + T(1) / 168 * (-a + 22), ( (T(1) / 672 * (-a + 3) * r[2] * b[1] + T(1) / 336 * (a - 5)) * b[2]^2 + - (T(1) / 336 * r[2] * b[1] + T(1) / 336 * (-a - 1)) * b[2] + - ( + (T(1) / 336 * r[2] * b[1] + T(1) / 336 * (-a - 1)) * b[2] + + ( (T(1) / 672 * (5 * a - 19) * r[2] + T(1) / 224 * (-a + 7)) * b[1] + - (T(1) / 336 * (5 * a - 21) * r[2] + T(1) / 336 * (-5 * a + 13)) + (T(1) / 336 * (5 * a - 21) * r[2] + T(1) / 336 * (-5 * a + 13)) ) ) * im + - (T(1) / 672 * (-a + 3) * b[1] + T(1) / 1008 * (a - 5) * r[2]) * b[2]^2 + - (T(1) / 336 * b[1] + T(1) / 1008 * (-a - 1) * r[2]) * b[2] + - (T(1) / 672 * (-a + 7) * r[2] + T(1) / 672 * (5 * a - 19)) * b[1] + - T(1) / 1008 * (-5 * a + 13) * r[2] + - T(1) / 336 * (5 * a - 21), + (T(1) / 672 * (-a + 3) * b[1] + T(1) / 1008 * (a - 5) * r[2]) * b[2]^2 + + (T(1) / 336 * b[1] + T(1) / 1008 * (-a - 1) * r[2]) * b[2] + + (T(1) / 672 * (-a + 7) * r[2] + T(1) / 672 * (5 * a - 19)) * b[1] + + T(1) / 1008 * (-5 * a + 13) * r[2] + + T(1) / 336 * (5 * a - 21), ( (T(1) / 504 * (a - 3) * r[2] * b[1] + T(1) / 504 * (a - 1) * r[2]) * b[2]^2 + - (T(1) / 1008 * (-a + 3) * r[2] * b[1] + T(1) / 252 * (-a + 2) * r[2]) * b[2] + - (T(1) / 168 * (-a + 4) * r[2] * b[1] + T(1) / 504 * (-9 * a + 11) * r[2]) + (T(1) / 1008 * (-a + 3) * r[2] * b[1] + T(1) / 252 * (-a + 2) * r[2]) * b[2] + + (T(1) / 168 * (-a + 4) * r[2] * b[1] + T(1) / 504 * (-9 * a + 11) * r[2]) ) * im + - (T(1) / 1008 * (a - 3) * r[2] * b[1] - T(1) / 126 * r[2]) * b[2]^2 + - (T(1) / 1008 * (a - 9) * r[2] * b[1] + T(1) / 504 * (a - 5) * r[2]) * b[2] + - T(1) / 168 * (-a + 2) * r[2] * b[1] + - T(1) / 504 * (-5 * a + 23) * r[2], + (T(1) / 1008 * (a - 3) * r[2] * b[1] - T(1) / 126 * r[2]) * b[2]^2 + + (T(1) / 1008 * (a - 9) * r[2] * b[1] + T(1) / 504 * (a - 5) * r[2]) * b[2] + + T(1) / 168 * (-a + 2) * r[2] * b[1] + + T(1) / 504 * (-5 * a + 23) * r[2], ( (T(1) / 2016 * (-a + 3) * r[2] * b[1] + T(1) / 336 * (-a + 1)) * b[2]^2 + - (T(1) / 2016 * (-a + 9) * r[2] * b[1] + T(1) / 168 * (a - 2)) * b[2] + - ( + (T(1) / 2016 * (-a + 9) * r[2] * b[1] + T(1) / 168 * (a - 2)) * b[2] + + ( (T(1) / 672 * (a + 3) * r[2] + T(1) / 672 * (a - 21)) * b[1] + - (T(1) / 168 * a * r[2] + T(1) / 168 * (3 * a - 16)) + (T(1) / 168 * a * r[2] + T(1) / 168 * (3 * a - 16)) ) ) * im + - (T(1) / 672 * (-a + 3) * b[1] + T(1) / 1008 * (a + 7) * r[2]) * b[2]^2 + - (T(1) / 672 * (a - 5) * b[1] + T(1) / 504 * (-2 * a + 7) * r[2]) * b[2] + - (T(1) / 672 * (3 * a - 7) * r[2] + T(1) / 672 * (a - 5)) * b[1] + - T(1) / 504 * (-a - 28) * r[2] + - T(1) / 168 * a, + (T(1) / 672 * (-a + 3) * b[1] + T(1) / 1008 * (a + 7) * r[2]) * b[2]^2 + + (T(1) / 672 * (a - 5) * b[1] + T(1) / 504 * (-2 * a + 7) * r[2]) * b[2] + + (T(1) / 672 * (3 * a - 7) * r[2] + T(1) / 672 * (a - 5)) * b[1] + + T(1) / 504 * (-a - 28) * r[2] + + T(1) / 168 * a, ( (T(1) / 672 * (a - 3) * b[1] + (T(1) / 1008 * (-a + 1) * r[2] + T(1) / 84)) * b[2]^2 + - ( + ( (T(1) / 672 * (a - 7) * r[2] + T(1) / 672 * (-a + 5)) * b[1] + - (T(1) / 504 * (a - 2) * r[2] + T(1) / 336 * (-a + 5)) + (T(1) / 504 * (a - 2) * r[2] + T(1) / 336 * (-a + 5)) ) * b[2] + - ( + ( (T(1) / 672 * (a - 7) * r[2] + T(1) / 672 * (-3 * a + 5)) * b[1] + - (T(1) / 1008 * (3 * a - 11) * r[2] + T(1) / 336 * (-a - 23)) + (T(1) / 1008 * (3 * a - 11) * r[2] + T(1) / 336 * (-a - 23)) ) ) * im + - (T(1) / 672 * (-a + 3) * r[2] * b[1] + (T(1) / 252 * r[2] + T(1) / 336 * (a - 1))) * b[2]^2 + - ( + (T(1) / 672 * (-a + 3) * r[2] * b[1] + (T(1) / 252 * r[2] + T(1) / 336 * (a - 1))) * b[2]^2 + + ( (T(1) / 672 * (a - 5) * r[2] + T(1) / 672 * (a - 7)) * b[1] + - (T(1) / 1008 * (-a + 5) * r[2] + T(1) / 168 * (-a + 2)) + (T(1) / 1008 * (-a + 5) * r[2] + T(1) / 168 * (-a + 2)) ) * b[2] + - (T(1) / 672 * (3 * a - 5) * r[2] + T(1) / 672 * (a - 7)) * b[1] + - T(1) / 1008 * (-a - 23) * r[2] + - T(1) / 336 * (-3 * a + 11) + (T(1) / 672 * (3 * a - 5) * r[2] + T(1) / 672 * (a - 7)) * b[1] + + T(1) / 1008 * (-a - 23) * r[2] + + T(1) / 336 * (-3 * a + 11), ] elseif d == 7 a = sqrt(T(2)) @@ -393,54 +393,54 @@ function _fiducial_WH(::Type{T}, d::Integer) where {T} T(1) / 14 * (a + 1), ( (T(1) / 196 * (a - 4) * r * t^2 + T(1) / 392 * (3 * a + 2) * r * t + T(1) / 392 * (a + 3) * r) * b + - (T(1) / 196 * (a + 6) * r * t^2 + T(1) / 392 * (-3 * a - 4) * r * t + T(1) / 392 * (-5 * a - 9) * r) + (T(1) / 196 * (a + 6) * r * t^2 + T(1) / 392 * (-3 * a - 4) * r * t + T(1) / 392 * (-5 * a - 9) * r) ) * im + - (-T(1) / 28 * a * t^2 + T(1) / 56 * (a - 2) * t + T(1) / 56 * (a - 1)) * b + - T(1) / 28 * (3 * a + 6) * t^2 + - T(1) / 56 * (-a - 4) * t + - T(1) / 56 * (-3 * a - 5), + (-T(1) / 28 * a * t^2 + T(1) / 56 * (a - 2) * t + T(1) / 56 * (a - 1)) * b + + T(1) / 28 * (3 * a + 6) * t^2 + + T(1) / 56 * (-a - 4) * t + + T(1) / 56 * (-3 * a - 5), ( (T(1) / 196 * (3 * a + 2) * r * t^2 + T(1) / 196 * (-2 * a + 1) * r * t + T(1) / 784 * (a - 4) * r) * - b + - (T(1) / 196 * (-3 * a - 4) * r * t^2 + T(1) / 196 * (a - 1) * r * t + T(1) / 784 * (-5 * a - 2) * r) + b + + (T(1) / 196 * (-3 * a - 4) * r * t^2 + T(1) / 196 * (a - 1) * r * t + T(1) / 784 * (-5 * a - 2) * r) ) * im + - (T(1) / 28 * (a - 2) * t^2 + T(1) / 28 * t - T(1) / 112 * a) * b + - T(1) / 28 * (-a - 4) * t^2 + - T(1) / 28 * (-a - 1) * t + - T(1) / 112 * (a + 6), + (T(1) / 28 * (a - 2) * t^2 + T(1) / 28 * t - T(1) / 112 * a) * b + + T(1) / 28 * (-a - 4) * t^2 + + T(1) / 28 * (-a - 1) * t + + T(1) / 112 * (a + 6), ( (T(1) / 98 * (2 * a - 1) * r * t^2 + T(1) / 392 * (-a + 4) * r * t + T(1) / 784 * (-11 * a + 2) * r) * - b + (T(1) / 98 * (a - 1) * r * t^2 + T(1) / 392 * (a + 6) * r * t + T(1) / 784 * (-13 * a - 8) * r) + b + (T(1) / 98 * (a - 1) * r * t^2 + T(1) / 392 * (a + 6) * r * t + T(1) / 784 * (-13 * a - 8) * r) ) * im + - (-T(1) / 14 * t^2 + T(1) / 56 * a * t + T(1) / 112 * (-a + 6)) * b + - T(1) / 14 * (-a - 1) * t^2 + - T(1) / 56 * (3 * a + 6) * t + - T(1) / 112 * a, + (-T(1) / 14 * t^2 + T(1) / 56 * a * t + T(1) / 112 * (-a + 6)) * b + + T(1) / 14 * (-a - 1) * t^2 + + T(1) / 56 * (3 * a + 6) * t + + T(1) / 112 * a, ( (T(1) / 98 * (-2 * a + 1) * r * t^2 + T(1) / 392 * (a - 4) * r * t + T(1) / 784 * (11 * a - 2) * r) * - b + (T(1) / 98 * (a - 1) * r * t^2 + T(1) / 392 * (a + 6) * r * t + T(1) / 784 * (-13 * a - 8) * r) + b + (T(1) / 98 * (a - 1) * r * t^2 + T(1) / 392 * (a + 6) * r * t + T(1) / 784 * (-13 * a - 8) * r) ) * im + - (T(1) / 14 * t^2 - T(1) / 56 * a * t + T(1) / 112 * (a - 6)) * b + - T(1) / 14 * (-a - 1) * t^2 + - T(1) / 56 * (3 * a + 6) * t + - T(1) / 112 * a, + (T(1) / 14 * t^2 - T(1) / 56 * a * t + T(1) / 112 * (a - 6)) * b + + T(1) / 14 * (-a - 1) * t^2 + + T(1) / 56 * (3 * a + 6) * t + + T(1) / 112 * a, ( (T(1) / 196 * (-3 * a - 2) * r * t^2 + T(1) / 196 * (2 * a - 1) * r * t + T(1) / 784 * (-a + 4) * r) * - b + - (T(1) / 196 * (-3 * a - 4) * r * t^2 + T(1) / 196 * (a - 1) * r * t + T(1) / 784 * (-5 * a - 2) * r) + b + + (T(1) / 196 * (-3 * a - 4) * r * t^2 + T(1) / 196 * (a - 1) * r * t + T(1) / 784 * (-5 * a - 2) * r) ) * im + - (T(1) / 28 * (-a + 2) * t^2 - T(1) / 28 * t + T(1) / 112 * a) * b + - T(1) / 28 * (-a - 4) * t^2 + - T(1) / 28 * (-a - 1) * t + - T(1) / 112 * (a + 6), + (T(1) / 28 * (-a + 2) * t^2 - T(1) / 28 * t + T(1) / 112 * a) * b + + T(1) / 28 * (-a - 4) * t^2 + + T(1) / 28 * (-a - 1) * t + + T(1) / 112 * (a + 6), ( (T(1) / 196 * (-a + 4) * r * t^2 + T(1) / 392 * (-3 * a - 2) * r * t + T(1) / 392 * (-a - 3) * r) * b + - (T(1) / 196 * (a + 6) * r * t^2 + T(1) / 392 * (-3 * a - 4) * r * t + T(1) / 392 * (-5 * a - 9) * r) + (T(1) / 196 * (a + 6) * r * t^2 + T(1) / 392 * (-3 * a - 4) * r * t + T(1) / 392 * (-5 * a - 9) * r) ) * im + - (T(1) / 28 * a * t^2 + T(1) / 56 * (-a + 2) * t + T(1) / 56 * (-a + 1)) * b + - T(1) / 28 * (3 * a + 6) * t^2 + - T(1) / 56 * (-a - 4) * t + - T(1) / 56 * (-3 * a - 5) + (T(1) / 28 * a * t^2 + T(1) / 56 * (-a + 2) * t + T(1) / 56 * (-a + 1)) * b + + T(1) / 28 * (3 * a + 6) * t^2 + + T(1) / 56 * (-a - 4) * t + + T(1) / 56 * (-3 * a - 5), ] else throw(ArgumentError(string("Invalid input dimension d = ", d))) diff --git a/src/multilinear.jl b/src/multilinear.jl index b03ef62..0fb3df2 100644 --- a/src/multilinear.jl +++ b/src/multilinear.jl @@ -11,9 +11,9 @@ end function _tidx!(tidx::AbstractVector{<:Integer}, idx::Integer, dims::Vector{<:Integer}) nsys = length(dims) - cidx = idx - 1 # Current index + cidx = idx - 1 # Current index dr = prod(dims) - for k = 1:nsys + for k in 1:nsys # Everytime you increase a tensor index you shift by the product of remaining dimensions dr ÷= dims[k] tidx[k] = (cidx ÷ dr) + 1 @@ -48,28 +48,28 @@ for (T, limit, wrapper) in [(:AbstractMatrix, :dY, :identity), (:(Hermitian), :j, :(Hermitian)), (:(Symmetric), :j, :(Symmetric))] @eval begin function partial_trace( - X::$T, - remove::AbstractVector{<:Integer}, - dims::AbstractVector{<:Integer} = _equal_sizes(X) - ) + X::$T, + remove::AbstractVector{<:Integer}, + dims::AbstractVector{<:Integer} = _equal_sizes(X) + ) isempty(remove) && return X length(remove) == length(dims) && return fill(tr(X), 1, 1) - keep = Vector{eltype(remove)}(undef, length(dims) - length(remove)) # Systems kept + keep = Vector{eltype(remove)}(undef, length(dims) - length(remove)) # Systems kept counter = 0 - for i = 1:length(dims) + for i in 1:length(dims) if !(i in remove) counter += 1 keep[counter] = i end end - dimsY = dims[keep] # The tensor dimensions of Y + dimsY = dims[keep] # The tensor dimensions of Y dimsR = dims[remove] # The tensor dimensions of the traced out systems dY = prod(dimsY) # Dimension of Y dR = prod(dimsR) # Dimension of system traced out Y = similar(X, (dY, dY)) # Final output Y - tXi = Vector{Int64}(undef, length(dims)) # Tensor indexing of X for column + tXi = Vector{Int64}(undef, length(dims)) # Tensor indexing of X for column tXj = Vector{Int64}(undef, length(dims)) # Tensor indexing of X for row @views tXikeep = tXi[keep] @@ -78,16 +78,16 @@ for (T, limit, wrapper) in @views tXjremove = tXj[remove] # We loop through Y and find the corresponding element - @inbounds for j = 1:dY + @inbounds for j in 1:dY # Find current column tensor index for Y _tidx!(tXjkeep, j, dimsY) - for i = 1:$limit + for i in 1:$limit # Find current row tensor index for Y _tidx!(tXikeep, i, dimsY) - # Now loop through the diagonal of the traced out systems + # Now loop through the diagonal of the traced out systems Y[i, j] = 0 - for k = 1:dR + for k in 1:dR _tidx!(tXiremove, k, dimsR) _tidx!(tXjremove, k, dimsR) @@ -128,16 +128,16 @@ for (T, wrapper) in [(:AbstractMatrix, :identity), (:(Hermitian), :(Hermitian)), (:(Symmetric), :(Symmetric))] @eval begin function partial_transpose( - X::$T, - transp::AbstractVector{<:Integer}, - dims::AbstractVector{<:Integer} = _equal_sizes(X) - ) + X::$T, + transp::AbstractVector{<:Integer}, + dims::AbstractVector{<:Integer} = _equal_sizes(X) + ) isempty(transp) && return X length(transp) == length(dims) && return transpose(X) - keep = Vector{eltype(transp)}(undef, length(dims) - length(transp)) # Systems kept + keep = Vector{eltype(transp)}(undef, length(dims) - length(transp)) # Systems kept counter = 0 - for i = 1:length(dims) + for i in 1:length(dims) if !(i in transp) counter += 1 keep[counter] = i @@ -147,15 +147,15 @@ for (T, wrapper) in d = size(X, 1) # Dimension of the final output Y Y = similar(X, (d, d)) # Final output Y - tXi = Vector{Int64}(undef, length(dims)) # Tensor indexing of X for row + tXi = Vector{Int64}(undef, length(dims)) # Tensor indexing of X for row tXj = Vector{Int64}(undef, length(dims)) # Tensor indexing of X for column - tYi = Vector{Int64}(undef, length(dims)) # Tensor indexing of Y for row + tYi = Vector{Int64}(undef, length(dims)) # Tensor indexing of Y for row tYj = Vector{Int64}(undef, length(dims)) # Tensor indexing of Y for column @inbounds for j in 1:d _tidx!(tYj, j, dims) - for i in 1:j-1 + for i in 1:(j - 1) _tidx!(tYi, i, dims) for k in keep @@ -228,10 +228,10 @@ end Permutes the order of the subsystems of vector `X` with subsystem dimensions `dims` in-place according to the permutation `perm`. If the argument `dims` is omitted two equally-sized subsystems are assumed. """ function permute_systems!( - X::AbstractVector{T}, - perm::AbstractVector{<:Integer}, - dims::AbstractVector{<:Integer} = _equal_sizes(X) -) where {T} + X::AbstractVector{T}, + perm::AbstractVector{<:Integer}, + dims::AbstractVector{<:Integer} = _equal_sizes(X) + ) where {T} perm == 1:length(perm) && return X X == 1:length(X) && return _idxperm!(X, perm, dims) @@ -242,10 +242,10 @@ end export permute_systems! @doc """ - permute_systems(X::AbstractMatrix, perm::AbstractVector, dims::AbstractVector = _equal_sizes(X)) + permute_systems(X::AbstractMatrix, perm::AbstractVector, dims::AbstractVector = _equal_sizes(X)) - Permutes the order of the subsystems of the square matrix `X`, which is composed by square subsystems of dimensions `dims`, according to the permutation `perm`. If the argument `dims` is omitted two equally-sized subsystems are assumed. - """ permute_systems( +Permutes the order of the subsystems of the square matrix `X`, which is composed by square subsystems of dimensions `dims`, according to the permutation `perm`. If the argument `dims` is omitted two equally-sized subsystems are assumed. +""" permute_systems( X::AbstractMatrix, perm::AbstractVector, dims::AbstractVector = _equal_sizes(X); @@ -255,11 +255,11 @@ for (T, wrapper) in [(:AbstractMatrix, :identity), (:(Hermitian), :(Hermitian)), (:(Symmetric), :(Symmetric))] @eval begin function permute_systems( - X::$T, - perm::AbstractVector{<:Integer}, - dims::AbstractVector{<:Integer} = _equal_sizes(X); - rows_only::Bool = false - ) + X::$T, + perm::AbstractVector{<:Integer}, + dims::AbstractVector{<:Integer} = _equal_sizes(X); + rows_only::Bool = false + ) perm == 1:length(perm) && return X p = _idxperm(perm, dims) @@ -272,7 +272,7 @@ end permute_systems(X::AbstractMatrix, perm::Vector, dims::Matrix) Permutes the order of the subsystems of the matrix `X`, which is composed by subsystems of dimensions `dims`, according to the permutation `perm`. -`dims` should be a n x 2 matrix where `dims[i, 1]` is the number of rows of subsystem i, and `dims[i,2]` is its number of columns. +`dims` should be a n x 2 matrix where `dims[i, 1]` is the number of rows of subsystem i, and `dims[i,2]` is its number of columns. """ function permute_systems(X::AbstractMatrix, perm::AbstractVector{<:Integer}, dims::Matrix{<:Integer}) perm == 1:length(perm) && return X @@ -290,10 +290,10 @@ Unitary that permutes subsystems of dimension `dims` according to the permutatio If `dims` is an Integer, assumes there are `length(perm)` subsystems of equal dimensions `dims`. """ function permutation_matrix( - dims::Union{Integer,AbstractVector{<:Integer}}, - perm::AbstractVector{<:Integer}; - is_sparse::Bool = true -) + dims::Union{Integer, AbstractVector{<:Integer}}, + perm::AbstractVector{<:Integer}; + is_sparse::Bool = true + ) dims = dims isa Integer ? fill(dims, length(perm)) : dims d = prod(dims) id = is_sparse ? SA.sparse(I, (d, d)) : I(d) diff --git a/src/nonlocal.jl b/src/nonlocal.jl index fdc1ebf..dea4dc6 100644 --- a/src/nonlocal.jl +++ b/src/nonlocal.jl @@ -6,7 +6,7 @@ as a 4-dimensional array. Reference: Araújo, Hirsch, and Quintino, [arXiv:2005.13418](https://arxiv.org/abs/2005.13418). """ -function local_bound(G::Array{T,4}) where {T<:Real} +function local_bound(G::Array{T, 4}) where {T <: Real} oa, ob, ia, ib = size(G) if oa^ia < ob^ib @@ -17,7 +17,7 @@ function local_bound(G::Array{T,4}) where {T<:Real} G = permutedims(G, (1, 3, 2, 4)) squareG = reshape(G, oa * ia, ob * ib) - offset = Vector(1 .+ ob * (0:ib-1)) + offset = Vector(1 .+ ob * (0:(ib - 1))) @views initial_score = sum(maximum(reshape(sum(squareG[:, offset]; dims = 2), oa, ia); dims = 1)) #compute initial_score for the all-zeros strategy to serve as a reference point chunks = _partition(ob^ib - 1, Threads.nthreads()) @@ -30,14 +30,14 @@ function local_bound(G::Array{T,4}) where {T<:Real} end export local_bound -function _local_bound_single(initial_score::T, chunk, sizeG, offset, squareG::Array{T,2}) where {T} +function _local_bound_single(initial_score::T, chunk, sizeG, offset, squareG::Array{T, 2}) where {T} oa, ob, ia, ib = sizeG score = initial_score ind = digits(chunk[1]; base = ob, pad = ib) offset_ind = zeros(Int64, ib) Galice = zeros(T, oa * ia, 1) maxvec = zeros(T, 1, ia) - for b = chunk[1]:chunk[2] + for b in chunk[1]:chunk[2] offset_ind .= ind .+ offset @views sum!(Galice, squareG[:, offset_ind]) squareGalice = Base.ReshapedArray(Galice, (oa, ia), ()) @@ -55,9 +55,9 @@ end If `n ≥ k` partitions the set `1:n` into `k` parts as equally sized as possible. Otherwise partitions it into `n` parts of size 1. """ -function _partition(n::T, k::T) where {T<:Integer} +function _partition(n::T, k::T) where {T <: Integer} num_parts = min(k, n) - parts = Vector{Tuple{T,T}}(undef, num_parts) + parts = Vector{Tuple{T, T}}(undef, num_parts) base_size = div(n, k) num_larger = rem(n, k) if num_larger > 0 @@ -67,17 +67,17 @@ function _partition(n::T, k::T) where {T<:Integer} end i = 2 while i ≤ num_larger - parts[i] = (1, base_size + 1) .+ parts[i-1][2] + parts[i] = (1, base_size + 1) .+ parts[i - 1][2] i += 1 end while i ≤ num_parts - parts[i] = (1, base_size) .+ parts[i-1][2] + parts[i] = (1, base_size) .+ parts[i - 1][2] i += 1 end return parts end -#copied from QETLAB +# copied from QETLAB function _update_odometer!(ind::AbstractVector{<:Integer}, upper_lim::Integer) # Start by increasing the last index by 1. ind_len = length(ind) @@ -86,13 +86,13 @@ function _update_odometer!(ind::AbstractVector{<:Integer}, upper_lim::Integer) # Now we work the "odometer": repeatedly set each digit to 0 if it # is too high and carry the addition to the left until we hit a # digit that *isn't* too high. - for j = ind_len:-1:1 + for j in ind_len:-1:1 # If we've hit the upper limit in this entry, move onto the next # entry. if ind[j] ≥ upper_lim ind[j] = 0 if j ≥ 2 - ind[j-1] += 1 + ind[j - 1] += 1 else # we're at the left end of the vector; just stop return end @@ -122,7 +122,7 @@ export tsirelson_bound Takes a bipartite Bell functional `V` in full probability notation and transforms it to Collins-Gisin notation. If `behaviour` is `true` do instead the transformation for behaviours. Doesn't assume normalization. """ -function fp2cg(V::AbstractArray{T,4}, behaviour::Bool = false) where {T} +function fp2cg(V::AbstractArray{T, 4}, behaviour::Bool = false) where {T} oa, ob, ia, ib = size(V) alice_pars = ia * (oa - 1) + 1 bob_pars = ib * (ob - 1) + 1 @@ -133,25 +133,25 @@ function fp2cg(V::AbstractArray{T,4}, behaviour::Bool = false) where {T} if ~behaviour CG[1, 1] = sum(V[oa, ob, :, :]) - for a = 1:oa-1, x = 1:ia + for a in 1:(oa - 1), x in 1:ia CG[aindex(a, x), 1] = sum(V[a, ob, x, :] - V[oa, ob, x, :]) end - for b = 1:ob-1, y = 1:ib + for b in 1:(ob - 1), y in 1:ib CG[1, bindex(b, y)] = sum(V[oa, b, :, y] - V[oa, ob, :, y]) end - for a = 1:oa-1, b = 1:ob-1, x = 1:ia, y = 1:ib + for a in 1:(oa - 1), b in 1:(ob - 1), x in 1:ia, y in 1:ib CG[aindex(a, x), bindex(b, y)] = V[a, b, x, y] - V[a, ob, x, y] - V[oa, b, x, y] + V[oa, ob, x, y] end else CG[1, 1] = sum(V) / (ia * ib) - for x = 1:ia, a = 1:oa-1 - CG[aindex(a, x), 1] = sum(V[a, b, x, y] for b = 1:ob, y = 1:ib) / ib + for x in 1:ia, a in 1:(oa - 1) + CG[aindex(a, x), 1] = sum(V[a, b, x, y] for b in 1:ob, y in 1:ib) / ib end - for y = 1:ib, b = 1:ob-1 - CG[1, bindex(b, y)] = sum(V[a, b, x, y] for a = 1:oa, x = 1:ia) / ia + for y in 1:ib, b in 1:(ob - 1) + CG[1, bindex(b, y)] = sum(V[a, b, x, y] for a in 1:oa, x in 1:ia) / ia end - for x = 1:ia, y = 1:ib - CG[aindex(1, x):aindex(oa - 1, x), bindex(1, y):bindex(ob - 1, y)] = V[1:oa-1, 1:ob-1, x, y] + for x in 1:ia, y in 1:ib + CG[aindex(1, x):aindex(oa - 1, x), bindex(1, y):bindex(ob - 1, y)] = V[1:(oa - 1), 1:(ob - 1), x, y] end end return CG @@ -169,19 +169,19 @@ function cg2fp(CG::Matrix{T}, scenario::Vector{<:Integer}, behaviour::Bool = fal aindex(a, x) = 1 + a + (x - 1) * (oa - 1) bindex(b, y) = 1 + b + (y - 1) * (ob - 1) - V = Array{T,4}(undef, (oa, ob, ia, ib)) + V = Array{T, 4}(undef, (oa, ob, ia, ib)) if ~behaviour - for x = 1:ia, y = 1:ib + for x in 1:ia, y in 1:ib V[oa, ob, x, y] = CG[1, 1] / (ia * ib) end - for x = 1:ia, y = 1:ib, b = 1:ob-1 + for x in 1:ia, y in 1:ib, b in 1:(ob - 1) V[oa, b, x, y] = CG[1, 1] / (ia * ib) + CG[1, bindex(b, y)] / ia end - for x = 1:ia, y = 1:ib, a = 1:oa-1 + for x in 1:ia, y in 1:ib, a in 1:(oa - 1) V[a, ob, x, y] = CG[1, 1] / (ia * ib) + CG[aindex(a, x), 1] / ib end - for x = 1:ia, y = 1:ib, a = 1:oa-1, b = 1:ob-1 + for x in 1:ia, y in 1:ib, a in 1:(oa - 1), b in 1:(ob - 1) V[a, b, x, y] = CG[1, 1] / (ia * ib) + CG[aindex(a, x), 1] / ib + @@ -189,12 +189,12 @@ function cg2fp(CG::Matrix{T}, scenario::Vector{<:Integer}, behaviour::Bool = fal CG[aindex(a, x), bindex(b, y)] end else - for x = 1:ia, y = 1:ib - V[1:oa-1, 1:ob-1, x, y] = CG[aindex(1, x):aindex(oa - 1, x), bindex(1, y):bindex(ob - 1, y)] - V[1:oa-1, ob, x, y] = + for x in 1:ia, y in 1:ib + V[1:(oa - 1), 1:(ob - 1), x, y] = CG[aindex(1, x):aindex(oa - 1, x), bindex(1, y):bindex(ob - 1, y)] + V[1:(oa - 1), ob, x, y] = CG[aindex(1, x):aindex(oa - 1, x), 1] - sum(CG[aindex(1, x):aindex(oa - 1, x), bindex(1, y):bindex(ob - 1, y)]; dims = 2) - V[oa, 1:ob-1, x, y] = + V[oa, 1:(ob - 1), x, y] = CG[1, bindex(1, y):bindex(ob - 1, y)] - vec(sum(CG[aindex(1, x):aindex(oa - 1, x), bindex(1, y):bindex(ob - 1, y)]; dims = 1)) V[oa, ob, x, y] = @@ -213,10 +213,10 @@ export cg2fp Applies N sets of measurements onto a state `rho` to form a probability array. """ function probability_tensor( - rho::Hermitian{T1,Matrix{T1}}, - first_Aax::Vector{Measurement{T2}}, # needed so that T2 is not unbounded - other_Aax::Vector{Measurement{T2}}... -) where {T1,T2} + rho::Hermitian{T1, Matrix{T1}}, + first_Aax::Vector{Measurement{T2}}, # needed so that T2 is not unbounded + other_Aax::Vector{Measurement{T2}}... + ) where {T1, T2} T = real(promote_type(T1, T2)) all_Aax = (first_Aax, other_Aax...) N = length(all_Aax) @@ -267,10 +267,10 @@ Convert a 2x...x2xmx...xm probability array into - a mx...xm correlation array (no marginals) - a (m+1)x...x(m+1) correlation array (marginals). """ -function correlation_tensor(p::AbstractArray{T,N2}; marg::Bool = true) where {T} where {N2} +function correlation_tensor(p::AbstractArray{T, N2}; marg::Bool = true) where {T} where {N2} @assert iseven(N2) N = N2 ÷ 2 - m = size(p)[N+1:end] # numbers of inputs per party + m = size(p)[(N + 1):end] # numbers of inputs per party o = size(p)[1:N] # numbers of outputs per party @assert collect(o) == 2ones(Int64, N) res = zeros(T, (marg ? m .+ 1 : m)...) diff --git a/src/random.jl b/src/random.jl index 0c18790..5310510 100644 --- a/src/random.jl +++ b/src/random.jl @@ -1,5 +1,5 @@ #often we don't actually need the variance of the normal variables to be 1, so we don't need to waste time diving everything by sqrt(2) -_randn(::Type{Complex{T}}) where {T<:AbstractFloat} = Complex{T}(randn(T), randn(T)) +_randn(::Type{Complex{T}}) where {T <: AbstractFloat} = Complex{T}(randn(T), randn(T)) _randn(::Type{T}) where {T} = randn(T) _randn(::Type{T}, dim1::Integer, dims::Integer...) where {T} = _randn!(Array{T}(undef, dim1, dims...)) function _randn!(A::AbstractArray{T}) where {T} @@ -48,16 +48,16 @@ If `T` is a real type the output is instead a Haar-random (real) orthogonal matr Reference: Stewart, [doi:10.1137/0717034](https://doi.org/10.1137/0717034). """ -function random_unitary(::Type{T}, d::Integer) where {T<:Number} +function random_unitary(::Type{T}, d::Integer) where {T <: Number} z = Matrix{T}(undef, d, d) - @inbounds for j = 1:d - for i = j:d + @inbounds for j in 1:d + for i in j:d z[i, j] = _randn(T) end end τ = Vector{T}(undef, d) s = Vector{T}(undef, d) - @inbounds for k = 1:d #this is a partial QR decomposition where we don't apply the reflection to the rest of the matrix + @inbounds for k in 1:d #this is a partial QR decomposition where we don't apply the reflection to the rest of the matrix @views x = z[k:d, k] τ[k] = LinearAlgebra.reflector!(x) s[k] = sign(real(x[1])) @@ -75,17 +75,17 @@ Produces a random POVM of dimension `d` with `n` outcomes and rank `min(k, d)`. Reference: Heinosaari et al., [arXiv:1902.04751](https://arxiv.org/abs/1902.04751). """ -function random_povm(::Type{T}, d::Integer, n::Integer, k::Integer = d) where {T<:Number} - d ≤ n * k || throw(ArgumentError("We need d ≤ n*k, but got d = $(d) and n*k = $(n*k)")) - E = [Matrix{T}(undef, (d, d)) for _ = 1:n] - for i = 1:n +function random_povm(::Type{T}, d::Integer, n::Integer, k::Integer = d) where {T <: Number} + d ≤ n * k || throw(ArgumentError("We need d ≤ n*k, but got d = $(d) and n*k = $(n * k)")) + E = [Matrix{T}(undef, (d, d)) for _ in 1:n] + for i in 1:n G = randn(T, (d, k)) mul!(E[i], G, G') end S = sum(E) rootinvS = Hermitian(S)^-0.5 #don't worry, the probability of getting a singular S is zero mat = Matrix{T}(undef, (d, d)) - for i = 1:n + for i in 1:n mul!(mat, rootinvS, E[i]) mul!(E[i], mat, rootinvS) end diff --git a/src/seesaw.jl b/src/seesaw.jl index 7568b4f..6346be9 100644 --- a/src/seesaw.jl +++ b/src/seesaw.jl @@ -6,10 +6,10 @@ Maximizes bipartite Bell functional in Collins-Gisin notation `CG` using the see If `oa` == `ob` == 2 the heuristic reduces to a bunch of eigenvalue problems. Otherwise semidefinite programming is needed and we use the assemblage version of seesaw. -References: Pál and Vértesi, [arXiv:1006.3032](https://arxiv.org/abs/1006.3032), +References: Pál and Vértesi, [arXiv:1006.3032](https://arxiv.org/abs/1006.3032), section II.B.1 of Tavakoli et al., [arXiv:2307.02551](https://arxiv.org/abs/2307.02551) """ -function seesaw(CG::Matrix{T}, scenario::Vector{<:Integer}, d::Integer) where {T<:Real} +function seesaw(CG::Matrix{T}, scenario::Vector{<:Integer}, d::Integer) where {T <: Real} R = _solver_type(T) CG = R.(CG) T2 = Complex{R} @@ -32,8 +32,8 @@ function seesaw(CG::Matrix{T}, scenario::Vector{<:Integer}, d::Integer) where {T if new_ω - ω <= minimumincrease || i > maxiter ω = new_ω ψ = state_phiplus_ket(T2, d; coeff = λ) - A = [[A[x]] for x = 1:ia] #rather inconvenient format - B = [[B[y]] for y = 1:ib] #but consistent with the general case + A = [[A[x]] for x in 1:ia] #rather inconvenient format + B = [[B[y]] for y in 1:ib] #but consistent with the general case break end ω = new_ω @@ -41,8 +41,8 @@ function seesaw(CG::Matrix{T}, scenario::Vector{<:Integer}, d::Integer) where {T end else B = Vector{Measurement{T2}}(undef, ib) - for y = 1:ib - B[y] = random_povm(T2, d, ob)[1:ob-1] + for y in 1:ib + B[y] = random_povm(T2, d, ob)[1:(ob - 1)] end local ψ, A ω = -R(Inf) @@ -63,17 +63,17 @@ function seesaw(CG::Matrix{T}, scenario::Vector{<:Integer}, d::Integer) where {T end export seesaw -function _optimise_alice_assemblage(CG::Matrix{R}, scenario, B; solver = Hypatia.Optimizer{R}) where {R<:AbstractFloat} +function _optimise_alice_assemblage(CG::Matrix{R}, scenario, B; solver = Hypatia.Optimizer{R}) where {R <: AbstractFloat} oa, ob, ia, ib = scenario d = size(B[1][1], 1) model = JuMP.GenericModel{R}() - ρxa = [[JuMP.@variable(model, [1:d, 1:d] in JuMP.HermitianPSDCone()) for _ in 1:oa-1] for _ in 1:ia] #assemblage + ρxa = [[JuMP.@variable(model, [1:d, 1:d] in JuMP.HermitianPSDCone()) for _ in 1:(oa - 1)] for _ in 1:ia] #assemblage ρ_B = JuMP.@variable(model, [1:d, 1:d], Hermitian) #auxiliary quantum state JuMP.@constraint(model, tr(ρ_B) == 1) - for x = 1:ia - JuMP.@constraint(model, ρ_B - sum(ρxa[x][a] for a in 1:oa-1) in JuMP.HermitianPSDCone()) + for x in 1:ia + JuMP.@constraint(model, ρ_B - sum(ρxa[x][a] for a in 1:(oa - 1)) in JuMP.HermitianPSDCone()) end ω = _compute_value_assemblage(CG, scenario, ρxa, ρ_B, B) @@ -82,18 +82,18 @@ function _optimise_alice_assemblage(CG::Matrix{R}, scenario, B; solver = Hypatia JuMP.set_optimizer(model, solver) JuMP.set_silent(model) JuMP.optimize!(model) - value_ρxa = [[Hermitian(JuMP.value.(ρxa[x][a])) for a = 1:oa-1] for x = 1:ia] + value_ρxa = [[Hermitian(JuMP.value.(ρxa[x][a])) for a in 1:(oa - 1)] for x in 1:ia] return JuMP.value(ω), value_ρxa, Hermitian(JuMP.value.(ρ_B)) end -function _optimise_bob_povm(CG::Matrix{R}, scenario, ρxa, ρ_B; solver = Hypatia.Optimizer{R}) where {R<:AbstractFloat} +function _optimise_bob_povm(CG::Matrix{R}, scenario, ρxa, ρ_B; solver = Hypatia.Optimizer{R}) where {R <: AbstractFloat} oa, ob, ia, ib = scenario d = size(ρ_B, 1) model = JuMP.GenericModel{R}() - B = [[JuMP.@variable(model, [1:d, 1:d] in JuMP.HermitianPSDCone()) for _ in 1:ob-1] for _ in 1:ib] #povm - for y = 1:ib - JuMP.@constraint(model, I - sum(B[y][b] for b in 1:ob-1) in JuMP.HermitianPSDCone()) + B = [[JuMP.@variable(model, [1:d, 1:d] in JuMP.HermitianPSDCone()) for _ in 1:(ob - 1)] for _ in 1:ib] #povm + for y in 1:ib + JuMP.@constraint(model, I - sum(B[y][b] for b in 1:(ob - 1)) in JuMP.HermitianPSDCone()) end ω = _compute_value_assemblage(CG, scenario, ρxa, ρ_B, B) @@ -102,29 +102,29 @@ function _optimise_bob_povm(CG::Matrix{R}, scenario, ρxa, ρ_B; solver = Hypati JuMP.set_optimizer(model, solver) JuMP.set_silent(model) JuMP.optimize!(model) - B = [[Hermitian(JuMP.value.(B[y][b])) for b = 1:ob-1] for y = 1:ib] + B = [[Hermitian(JuMP.value.(B[y][b])) for b in 1:(ob - 1)] for y in 1:ib] return JuMP.value(ω), B end -function _compute_value_assemblage(CG::Matrix{R}, scenario, ρxa, ρ_B, B) where {R<:AbstractFloat} +function _compute_value_assemblage(CG::Matrix{R}, scenario, ρxa, ρ_B, B) where {R <: AbstractFloat} oa, ob, ia, ib = scenario aind(a, x) = 1 + a + (x - 1) * (oa - 1) bind(b, y) = 1 + b + (y - 1) * (ob - 1) - ω = CG[1, 1] * one(JuMP.GenericAffExpr{R,JuMP.GenericVariableRef{R}}) - for a = 1:oa-1 - for x = 1:ia - tempB = sum(CG[aind(a, x), bind(b, y)] * B[y][b] for b = 1:ob-1 for y = 1:ib) + ω = CG[1, 1] * one(JuMP.GenericAffExpr{R, JuMP.GenericVariableRef{R}}) + for a in 1:(oa - 1) + for x in 1:ia + tempB = sum(CG[aind(a, x), bind(b, y)] * B[y][b] for b in 1:(ob - 1) for y in 1:ib) ω += real(dot(tempB, ρxa[x][a])) end end - for a = 1:oa-1 - for x = 1:ia + for a in 1:(oa - 1) + for x in 1:ia ω += CG[aind(a, x), 1] * real(tr(ρxa[x][a])) end end - for b = 1:ob-1 - for y = 1:ib + for b in 1:(ob - 1) + for y in 1:ib ω += CG[1, bind(b, y)] * real(dot(B[y][b], ρ_B)) end end @@ -137,28 +137,28 @@ function _decompose_assemblage(scenario, ρxa, ρ_B::AbstractMatrix{T}) where {T d = size(ρ_B, 1) λ, U = eigen(ρ_B) ψ = zeros(T, d^2) - for i = 1:d + for i in 1:d @views ψ .+= sqrt(λ[i]) * kron(conj(U[:, i]), U[:, i]) end W = zeros(T, d, d) - for i = 1:d + for i in 1:d if λ[i] >= _eps(T) @views W .+= ketbra(U[:, i]) / sqrt(λ[i]) end end - A = [[Hermitian(conj(W * ρxa[x][a] * W)) for a = 1:oa-1] for x = 1:ia] + A = [[Hermitian(conj(W * ρxa[x][a] * W)) for a in 1:(oa - 1)] for x in 1:ia] return ψ, A end function _optimise_alice_projectors!(CG::Matrix, λ::Vector, A, B) ia, ib = size(CG) .- 1 d = length(λ) - for x = 1:ia - for j = 1:d - for i = 1:j - A[x].data[i, j] = CG[x+1, 1] * (i == j) * abs2(λ[i]) - for y = 1:ib - A[x].data[i, j] += CG[x+1, y+1] * λ[i] * conj(λ[j] * B[y][i, j]) + for x in 1:ia + for j in 1:d + for i in 1:j + A[x].data[i, j] = CG[x + 1, 1] * (i == j) * abs2(λ[i]) + for y in 1:ib + A[x].data[i, j] += CG[x + 1, y + 1] * λ[i] * conj(λ[j] * B[y][i, j]) end end end @@ -169,12 +169,12 @@ end function _optimise_bob_projectors!(CG::Matrix, λ::Vector, A, B) ia, ib = size(CG) .- 1 d = length(λ) - for y = 1:ib - for j = 1:d - for i = 1:j - B[y].data[i, j] = CG[1, y+1] * (i == j) * abs2(λ[i]) - for x = 1:ia - B[y].data[i, j] += CG[x+1, y+1] * λ[i] * conj(λ[j] * A[x][i, j]) + for y in 1:ib + for j in 1:d + for i in 1:j + B[y].data[i, j] = CG[1, y + 1] * (i == j) * abs2(λ[i]) + for x in 1:ia + B[y].data[i, j] += CG[x + 1, y + 1] * λ[i] * conj(λ[j] * A[x][i, j]) end end end @@ -185,7 +185,7 @@ end function _positive_projection!(M::AbstractMatrix{T}) where {T} λ, U = eigen!(M) fill!(M, 0) - for i = 1:length(λ) + for i in 1:length(λ) if λ[i] > _rtol(T) @views M.data .+= ketbra(U[:, i]) end @@ -197,17 +197,17 @@ function _optimise_state!(CG::Matrix, λ::Vector{T}, A, B) where {T} ia, ib = size(CG) .- 1 d = length(λ) M = Hermitian(zeros(T, d, d)) - for x = 1:ia - M += CG[x+1, 1] * real(Diagonal(A[x])) + for x in 1:ia + M += CG[x + 1, 1] * real(Diagonal(A[x])) end - for y = 1:ib - M += CG[1, y+1] * real(Diagonal(B[y])) + for y in 1:ib + M += CG[1, y + 1] * real(Diagonal(B[y])) end - for y = 1:ib - for x = 1:ia - for j = 1:d - for i = 1:j - M.data[i, j] += CG[x+1, y+1] * A[x].data[i, j] * B[y].data[i, j] + for y in 1:ib + for x in 1:ia + for j in 1:d + for i in 1:j + M.data[i, j] += CG[x + 1, y + 1] * A[x].data[i, j] * B[y].data[i, j] end end end diff --git a/src/states.jl b/src/states.jl index dd5d08a..ebc7778 100644 --- a/src/states.jl +++ b/src/states.jl @@ -31,7 +31,7 @@ export white_noise! Produces the vector of the maximally entangled state Φ⁺ of local dimension `d`. """ -function state_phiplus_ket(::Type{T}, d::Integer = 2; kwargs...) where {T<:Number} +function state_phiplus_ket(::Type{T}, d::Integer = 2; kwargs...) where {T <: Number} return state_ghz_ket(T, d, 2; kwargs...) end state_phiplus_ket(d::Integer = 2; kwargs...) = state_phiplus_ket(ComplexF64, d; kwargs...) @@ -42,7 +42,7 @@ export state_phiplus_ket Produces the maximally entangled state Φ⁺ of local dimension `d` with visibility `v`. """ -function state_phiplus(::Type{T}, d::Integer = 2; v::Real = 1) where {T<:Number} +function state_phiplus(::Type{T}, d::Integer = 2; v::Real = 1) where {T <: Number} rho = ketbra(state_phiplus_ket(T, d; coeff = one(T))) parent(rho) ./= d return white_noise!(rho, v) @@ -55,9 +55,9 @@ export state_phiplus Produces the vector of the maximally entangled state ψ⁻ of local dimension `d`. """ -function state_psiminus_ket(::Type{T}, d::Integer = 2; coeff = inv(_sqrt(T, d))) where {T<:Number} +function state_psiminus_ket(::Type{T}, d::Integer = 2; coeff = inv(_sqrt(T, d))) where {T <: Number} psi = zeros(T, d^2) - psi[d.+(d-1)*(0:d-1)] .= (-1) .^ (0:d-1) .* coeff + psi[d .+ (d - 1) * (0:(d - 1))] .= (-1) .^ (0:(d - 1)) .* coeff return psi end state_psiminus_ket(d::Integer = 2; kwargs...) = state_psiminus_ket(ComplexF64, d; kwargs...) @@ -68,7 +68,7 @@ export state_psiminus_ket Produces the maximally entangled state ψ⁻ of local dimension `d` with visibility `v`. """ -function state_psiminus(::Type{T}, d::Integer = 2; v::Real = 1) where {T<:Number} +function state_psiminus(::Type{T}, d::Integer = 2; v::Real = 1) where {T <: Number} rho = ketbra(state_psiminus_ket(T, d; coeff = one(T))) parent(rho) ./= d return white_noise!(rho, v) @@ -81,10 +81,10 @@ export state_psiminus Produces the vector of the GHZ state local dimension `d`. """ -function state_ghz_ket(::Type{T}, d::Integer = 2, N::Integer = 3; coeff = inv(_sqrt(T, d))) where {T<:Number} +function state_ghz_ket(::Type{T}, d::Integer = 2, N::Integer = 3; coeff = inv(_sqrt(T, d))) where {T <: Number} psi = zeros(T, d^N) spacing = (1 - d^N) ÷ (1 - d) - psi[1:spacing:d^N] .= coeff + psi[1:spacing:(d^N)] .= coeff return psi end state_ghz_ket(d::Integer = 2, N::Integer = 3; kwargs...) = state_ghz_ket(ComplexF64, d, N; kwargs...) @@ -95,7 +95,7 @@ export state_ghz_ket Produces the GHZ state of local dimension `d` with visibility `v`. """ -function state_ghz(::Type{T}, d::Integer = 2, N::Integer = 3; v::Real = 1, kwargs...) where {T<:Number} +function state_ghz(::Type{T}, d::Integer = 2, N::Integer = 3; v::Real = 1, kwargs...) where {T <: Number} return white_noise!(ketbra(state_ghz_ket(T, d, N; kwargs...)), v) end state_ghz(d::Integer = 2, N::Integer = 3; kwargs...) = state_ghz(ComplexF64, d, N; kwargs...) @@ -106,9 +106,9 @@ export state_ghz Produces the vector of the `N`-partite W state. """ -function state_w_ket(::Type{T}, N::Integer = 3; coeff = inv(_sqrt(T, N))) where {T<:Number} +function state_w_ket(::Type{T}, N::Integer = 3; coeff = inv(_sqrt(T, N))) where {T <: Number} psi = zeros(T, 2^N) - psi[2 .^ (0:N-1).+1] .= coeff + psi[2 .^ (0:(N - 1)) .+ 1] .= coeff return psi end state_w_ket(N::Integer = 3; kwargs...) = state_w_ket(ComplexF64, N; kwargs...) @@ -119,7 +119,7 @@ export state_w_ket Produces the `N`-partite W state with visibility `v`. """ -function state_w(::Type{T}, N::Integer = 3; v::Real = 1, kwargs...) where {T<:Number} +function state_w(::Type{T}, N::Integer = 3; v::Real = 1, kwargs...) where {T <: Number} return white_noise!(ketbra(state_w_ket(T, N; kwargs...)), v) end state_w(N::Integer = 3; kwargs...) = state_w(ComplexF64, N; kwargs...) @@ -130,7 +130,7 @@ export state_w Produces the isotropic state of local dimension `d` with visibility `v`. """ -function isotropic(v::T, d::Integer = 2) where {T<:Real} +function isotropic(v::T, d::Integer = 2) where {T <: Real} return state_phiplus(T, d; v) end export isotropic @@ -142,7 +142,7 @@ Produces the vector of the `N`-partite `N`-level singlet state. Reference: Adán Cabello, [arXiv:quant-ph/0203119](https://arxiv.org/abs/quant-ph/0203119) """ -function state_super_singlet_ket(::Type{T}, N::Integer = 3; coeff = inv(_sqrt(T, factorial(N)))) where {T<:Number} +function state_super_singlet_ket(::Type{T}, N::Integer = 3; coeff = inv(_sqrt(T, factorial(N)))) where {T <: Number} psi = zeros(T, N^N) for per in Combinatorics.permutations(1:N) tmp = kron(ket.((T,), per, (N,))...) @@ -169,11 +169,10 @@ This state is invariant under simultaneous rotations on all parties: `(U⊗… Reference: Adán Cabello, [arXiv:quant-ph/0203119](https://arxiv.org/abs/quant-ph/0203119) """ -function state_super_singlet(::Type{T}, N::Integer = 3; v::Real = 1) where {T<:Number} +function state_super_singlet(::Type{T}, N::Integer = 3; v::Real = 1) where {T <: Number} rho = ketbra(state_super_singlet_ket(T, N; coeff = one(T))) parent(rho) ./= factorial(N) return white_noise!(rho, v) end state_super_singlet(N::Integer = 3; kwargs...) = state_super_singlet(ComplexF64, N; kwargs...) export state_super_singlet - diff --git a/src/supermaps.jl b/src/supermaps.jl index 0b22745..3bb429c 100644 --- a/src/supermaps.jl +++ b/src/supermaps.jl @@ -1,7 +1,7 @@ #extract from T the kind of float to be used in the conic solver -_solver_type(::Type{T}) where {T<:AbstractFloat} = T -_solver_type(::Type{Complex{T}}) where {T<:AbstractFloat} = T -_solver_type(::Type{T}) where {T<:Number} = Float64 +_solver_type(::Type{T}) where {T <: AbstractFloat} = T +_solver_type(::Type{Complex{T}}) where {T <: AbstractFloat} = T +_solver_type(::Type{T}) where {T <: Number} = Float64 """ choi(K::Vector{<:AbstractMatrix}) @@ -26,13 +26,13 @@ function diamond_norm(J::AbstractMatrix{T}, dims::AbstractVector; solver = Hypat din, dout = dims model = JuMP.GenericModel{_solver_type(T)}() if is_complex - JuMP.@variable(model, Y[1:din*dout, 1:din*dout], Hermitian) + JuMP.@variable(model, Y[1:(din * dout), 1:(din * dout)], Hermitian) JuMP.@variable(model, σ[1:din, 1:din], Hermitian) bigσ = Hermitian(kron(σ, I(dout))) JuMP.@constraint(model, bigσ - Y in JuMP.HermitianPSDCone()) JuMP.@constraint(model, bigσ + Y in JuMP.HermitianPSDCone()) else - JuMP.@variable(model, Y[1:din*dout, 1:din*dout], Symmetric) + JuMP.@variable(model, Y[1:(din * dout), 1:(din * dout)], Symmetric) JuMP.@variable(model, σ[1:din, 1:din], Symmetric) bigσ = Symmetric(kron(σ, I(dout))) JuMP.@constraint(model, bigσ - Y in JuMP.PSDCone()) diff --git a/test/entanglement.jl b/test/entanglement.jl index 01b5ca4..d444313 100644 --- a/test/entanglement.jl +++ b/test/entanglement.jl @@ -14,13 +14,13 @@ for R in [Float64, Double64] #Float128 and BigFloat take too long Random.seed!(8) #makes all states entangled ψ = random_state_ket(R, 6) - @test entanglement_entropy(ψ, [2, 3]) ≈ entanglement_entropy(ketbra(ψ), [2, 3])[1] atol = 1e-3 rtol = 1e-3 + @test entanglement_entropy(ψ, [2, 3]) ≈ entanglement_entropy(ketbra(ψ), [2, 3])[1] atol = 1.0e-3 rtol = 1.0e-3 ρ = random_state(R, 4) h, σ = entanglement_entropy(ρ) @test Ket._test_entanglement_entropy_qubit(h, ρ, σ) T = Complex{R} ψ = random_state_ket(T, 6) - @test entanglement_entropy(ψ, [2, 3]) ≈ entanglement_entropy(ketbra(ψ), [2, 3])[1] atol = 1e-3 rtol = 1e-3 + @test entanglement_entropy(ψ, [2, 3]) ≈ entanglement_entropy(ketbra(ψ), [2, 3])[1] atol = 1.0e-3 rtol = 1.0e-3 ρ = random_state(T, 4) h, σ = entanglement_entropy(ρ) @test Ket._test_entanglement_entropy_qubit(h, ρ, σ) @@ -31,22 +31,22 @@ ρ = state_ghz(R, 2, 2) s, W = random_robustness(ρ) @test eltype(W) == R - @test s ≈ 0.5 atol = 1e-5 rtol = 1e-5 - @test dot(ρ, W) ≈ -s atol = 1e-5 rtol = 1e-5 + @test s ≈ 0.5 atol = 1.0e-5 rtol = 1.0e-5 + @test dot(ρ, W) ≈ -s atol = 1.0e-5 rtol = 1.0e-5 T = Complex{R} ρ = state_ghz(T, 2, 2) s, W = random_robustness(ρ) @test eltype(W) == T - @test s ≈ 0.5 atol = 1e-5 rtol = 1e-5 - @test dot(ρ, W) ≈ -s atol = 1e-5 rtol = 1e-5 + @test s ≈ 0.5 atol = 1.0e-5 rtol = 1.0e-5 + @test dot(ρ, W) ≈ -s atol = 1.0e-5 rtol = 1.0e-5 end d = 3 - @test isapprox(schmidt_number(state_ghz(ComplexF64, d, 2), 2), 1 / 15, atol = 1e-3, rtol = 1e-3) + @test isapprox(schmidt_number(state_ghz(ComplexF64, d, 2), 2), 1 / 15, atol = 1.0e-3, rtol = 1.0e-3) @test isapprox( schmidt_number(state_ghz(Float64, d, 2), 2, [d, d], 2; solver = SCS.Optimizer), 1 / 15, - atol = 1e-3, - rtol = 1e-3 + atol = 1.0e-3, + rtol = 1.0e-3 ) @test schmidt_number(random_state(Float64, 6), 2, [2, 3], 1; solver = SCS.Optimizer) <= 0 end @@ -55,15 +55,15 @@ ρ = state_ghz(R, 2, 3) v, W = ppt_mixture(ρ, [2, 2, 2]) - @test isapprox(v, 0.4285, atol = 1e-3) + @test isapprox(v, 0.4285, atol = 1.0e-3) full_body_basis = collect(Iterators.flatten(n_body_basis(i, 3) for i in 0:3)) v, w = ppt_mixture(ρ, [2, 2, 2], full_body_basis) - @test isapprox(v, 0.4285, atol = 1e-3) - @test isapprox(sum(w[i] * full_body_basis[i] for i in eachindex(w)), W, atol = 1e-3) + @test isapprox(v, 0.4285, atol = 1.0e-3) + @test isapprox(sum(w[i] * full_body_basis[i] for i in eachindex(w)), W, atol = 1.0e-3) two_body_basis = collect(Iterators.flatten(n_body_basis(i, 3) for i in 0:2)) v, w = ppt_mixture(state_w(3), [2, 2, 2], two_body_basis) - @test isapprox(v, 0.6960, atol = 1e-3) + @test isapprox(v, 0.696, atol = 1.0e-3) end end end diff --git a/test/entropy.jl b/test/entropy.jl index 67444c4..1a64781 100644 --- a/test/entropy.jl +++ b/test/entropy.jl @@ -29,7 +29,7 @@ U = random_unitary(Complex{R}, 4) ρ2 = Hermitian(U * [ρ zeros(3, 1); zeros(1, 3) 0] * U') σ2 = Hermitian(U * [σ zeros(3, 1); zeros(1, 3) 0] * U') - @test relative_entropy(ρ, σ) ≈ relative_entropy(ρ2, σ2) atol = 1e-8 rtol = sqrt(Base.rtoldefault(R)) + @test relative_entropy(ρ, σ) ≈ relative_entropy(ρ2, σ2) atol = 1.0e-8 rtol = sqrt(Base.rtoldefault(R)) p = rand(R) q = rand(R) @test binary_relative_entropy(p, q) ≈ binary_relative_entropy(ℯ, p, q) / log(R(2)) diff --git a/test/incompatibility.jl b/test/incompatibility.jl index 1600553..bcfcc4e 100644 --- a/test/incompatibility.jl +++ b/test/incompatibility.jl @@ -1,11 +1,11 @@ @testset "Incompatibility " begin A = povm(mub(2)) - @test abs(incompatibility_robustness_depolarizing(A) - 1/√3) < 1e-6 - @test abs(incompatibility_robustness_random(A) - 1/√3) < 1e-6 - @test abs(incompatibility_robustness_probabilistic(A) - 1/√3) < 1e-6 - @test abs(incompatibility_robustness_jointly_measurable(A) - (√3 - 1)) < 1e-6 - @test abs(incompatibility_robustness_generalized(A) - (1 + 1/√3) / 2) < 1e-6 + @test abs(incompatibility_robustness_depolarizing(A) - 1 / √3) < 1.0e-6 + @test abs(incompatibility_robustness_random(A) - 1 / √3) < 1.0e-6 + @test abs(incompatibility_robustness_probabilistic(A) - 1 / √3) < 1.0e-6 + @test abs(incompatibility_robustness_jointly_measurable(A) - (√3 - 1)) < 1.0e-6 + @test abs(incompatibility_robustness_generalized(A) - (1 + 1 / √3) / 2) < 1.0e-6 # other types - @test abs(incompatibility_robustness_depolarizing(povm(broadcast.(Float64, mub(2, 2)))) - 1/√2) < 1e-6 - @test abs(incompatibility_robustness_depolarizing(povm((mub(Complex{Double64}, 2)))) - 1/√3) < 1e-10 + @test abs(incompatibility_robustness_depolarizing(povm(broadcast.(Float64, mub(2, 2)))) - 1 / √2) < 1.0e-6 + @test abs(incompatibility_robustness_depolarizing(povm((mub(Complex{Double64}, 2)))) - 1 / √3) < 1.0e-10 end diff --git a/test/measurements.jl b/test/measurements.jl index e006499..0d44087 100644 --- a/test/measurements.jl +++ b/test/measurements.jl @@ -3,10 +3,10 @@ for T in [Float64, Double64, Float128, BigFloat] E = random_povm(Complex{T}, 2, 3) V = dilate_povm(E) - @test [V' * kron(I(2), proj(i, 3)) * V for i = 1:3] ≈ E + @test [V' * kron(I(2), proj(i, 3)) * V for i in 1:3] ≈ E e = sic_povm(Complex{T}, 2) V = dilate_povm(e) - @test [V' * proj(i, 4) * V for i = 1:4] ≈ [ketbra(e[i]) for i = 1:4] + @test [V' * proj(i, 4) * V for i in 1:4] ≈ [ketbra(e[i]) for i in 1:4] end end @testset "SIC POVMs" begin diff --git a/test/nonlocal.jl b/test/nonlocal.jl index 3c072af..a2d0bc9 100644 --- a/test/nonlocal.jl +++ b/test/nonlocal.jl @@ -6,8 +6,8 @@ @test local_bound(chsh(Int64, 3)) == 6 @test local_bound(cglmp(Int64, 4)) == 9 Random.seed!(1337) - @test seesaw(fp2cg(cglmp()),[3,3,2,2],3)[1] ≈ (15+sqrt(33))/24 - @test seesaw(inn22(),[2,2,3,3],2)[1] ≈ 1.25 + @test seesaw(fp2cg(cglmp()), [3, 3, 2, 2], 3)[1] ≈ (15 + sqrt(33)) / 24 + @test seesaw(inn22(), [2, 2, 3, 3], 2)[1] ≈ 1.25 for T in [Float64, Double64, Float128, BigFloat] @test eltype(chsh(T)) <: T @test eltype(cglmp(T)) <: T @@ -19,8 +19,8 @@ end Aax = povm(mub(2)) @test correlation_tensor(state_phiplus(), Aax, 2) ≈ Diagonal([1, 1, -1, 1]) @test correlation_tensor(state_phiplus(), Aax, 2; marg = false) ≈ Diagonal([1, 1, -1]) - scenario = [2,3,4,5] - cg = randn(scenario[3]*(scenario[1]-1) + 1, scenario[4]*(scenario[2]-1) + 1) - @test fp2cg(cg2fp(cg,scenario)) ≈ cg - @test fp2cg(cg2fp(cg,scenario,true),true) ≈ cg + scenario = [2, 3, 4, 5] + cg = randn(scenario[3] * (scenario[1] - 1) + 1, scenario[4] * (scenario[2] - 1) + 1) + @test fp2cg(cg2fp(cg, scenario)) ≈ cg + @test fp2cg(cg2fp(cg, scenario, true), true) ≈ cg end diff --git a/test/norms.jl b/test/norms.jl index 8c68f4a..5af0195 100644 --- a/test/norms.jl +++ b/test/norms.jl @@ -1,7 +1,7 @@ @testset "Norms " begin @testset "Square operators" begin for T in [Float64, Double64, Float128, BigFloat] - X = Complex{T}.([1.34 5im 2.6; 5.2 4+2.1im 0.1; -4.1 0.9-2.2im 1im]) + X = Complex{T}.([1.34 5im 2.6; 5.2 4 + 2.1im 0.1; -4.1 0.9 - 2.2im 1im]) @test isapprox(schatten_norm(X, 1.32), 12.366965628920612) @test isapprox(trace_norm(X), 15.34224630614291) @test isapprox(schatten_norm(X, 2), 10.22133063744638) @@ -15,7 +15,7 @@ @testset "Rectangular operators" begin for T in [Float64, Double64, Float128, BigFloat] - X = Complex{T}.([1.34 5im 2.6; 5.2 4+2.1im 0.1]) + X = Complex{T}.([1.34 5im 2.6; 5.2 4 + 2.1im 0.1]) @test isapprox(schatten_norm(X, 2.43), 8.68672362793173) @test isapprox(trace_norm(X), 11.8442700353432) @test isapprox(schatten_norm(X, 2), 9.00086662494229) diff --git a/test/random.jl b/test/random.jl index d500e4a..9f7d7fa 100644 --- a/test/random.jl +++ b/test/random.jl @@ -63,23 +63,23 @@ for R in [Float64, Double64, Float128, BigFloat] E = random_povm(R, 2, 3) @test test_povm(E) - for i = 1:length(E) + for i in 1:length(E) @test rank(E[i]; rtol = Ket._rtol(R)) == 2 end E = random_povm(R, 2, 3, 1) @test test_povm(E) - for i = 1:length(E) + for i in 1:length(E) @test rank(E[i]; rtol = Ket._rtol(R)) == 1 end T = Complex{R} E = random_povm(T, 2, 3) @test test_povm(E) - for i = 1:length(E) + for i in 1:length(E) @test rank(E[i]; rtol = Ket._rtol(R)) == 2 end E = random_povm(T, 2, 3, 1) @test test_povm(E) - for i = 1:length(E) + for i in 1:length(E) @test rank(E[i]; rtol = Ket._rtol(R)) == 1 end end diff --git a/test/supermaps.jl b/test/supermaps.jl index 7a67e78..0ac009a 100644 --- a/test/supermaps.jl +++ b/test/supermaps.jl @@ -1,10 +1,10 @@ @testset "Supermaps " begin for R in [Float64, Double64] #Float128 and BigFloat take too long din, dout = 2, 3 - K = [randn(R, dout, din) for _ = 1:3] - @test diamond_norm(K) ≈ diamond_norm(choi(K), [din, dout]) atol = 1e-8 rtol = sqrt(Base.rtoldefault(R)) + K = [randn(R, dout, din) for _ in 1:3] + @test diamond_norm(K) ≈ diamond_norm(choi(K), [din, dout]) atol = 1.0e-8 rtol = sqrt(Base.rtoldefault(R)) T = Complex{R} - K = [randn(T, dout, din) for _ = 1:3] - @test diamond_norm(K) ≈ diamond_norm(choi(K), [din, dout]) atol = 1e-8 rtol = sqrt(Base.rtoldefault(R)) + K = [randn(T, dout, din) for _ in 1:3] + @test diamond_norm(K) ≈ diamond_norm(choi(K), [din, dout]) atol = 1.0e-8 rtol = sqrt(Base.rtoldefault(R)) end end