export norm_diamond, diamond_distance

const MOI = Convex.MOI

"""
$(SIGNATURES)
- `Î¦`: DynamicalMatrix

Return [diamond norm](https://arxiv.org/pdf/1207.5726.pdf) of dynamical matrix `Î¦`.
"""
norm_diamond(Î¦::AbstractQuantumOperation, method::Symbol=:primal, eps=1e-7) = norm_diamond(Î¦, Val(method), eps)

function norm_diamond(Î¦::DynamicalMatrix, ::Val{:primal}, eps)
    J = Î¦.matrix
    # TODO: compare d1, d2 with idim, odim
    d1 = Î¦.idim
    d2 = Î¦.odim
    X = ComplexVariable(d1*d2, d1*d2)
    t = 0.5*inner_product(X, J) + 0.5*inner_product(X', J')

    Ï�â‚€ = ComplexVariable(d1, d1)
    Ï�â‚� = ComplexVariable(d1, d1)

    constraints = [Ï�â‚€ in :SDP, Ï�â‚� in :SDP]
    constraints += tr(Ï�â‚€) == 1
    constraints += tr(Ï�â‚�) == 1
    constraints += [ð�•€(d2) âŠ— Ï�â‚€ X; X' ð�•€(d2) âŠ— Ï�â‚�] in :SDP

    problem = maximize(t, constraints)
    solve!(
        problem,
        MOI.OptimizerWithAttributes(SCS.Optimizer, "eps_abs" => eps);
        silent_solver = true
    )
    problem.optval
end


function norm_diamond(Î¦::DynamicalMatrix, ::Val{:dual}, eps)
    J = Î¦.matrix
    # TODO: compare d1, d2 with idim, odim
    d1 = Î¦.idim
    d2 = Î¦.odim
    Yâ‚€ = ComplexVariable(d1*d2, d1*d2)
    Yâ‚� = ComplexVariable(d1*d2, d1*d2)

	t = 0.5*sigmamax(partialtrace(Yâ‚€, 1, [d2,d1])) +
		0.5*sigmamax(partialtrace(Yâ‚�, 1, [d2,d1]))
	Z = [Yâ‚€ -J; -J' Yâ‚� ]

    constraints = [Yâ‚€ in :SDP, Yâ‚� in :SDP]
    constraints += Z+Z' in :SDP

    problem = minimize(t, constraints)
    solve!(
        problem,
        MOI.OptimizerWithAttributes(SCS.Optimizer, "eps_abs" => eps);
        silent_solver = true
    )
    problem.optval
end



# function norm_diamond(Î¦::AbstractQuantumOperation{T}, args...) where T<:AbstractMatrix{<:Number}
#     norm_diamond(DynamicalMatrix{T}(Ï•), args...)
# end

"""
$(SIGNATURES)
- `Î¦1`: DynamicalMatrix
- `Î¦2`: DynamicalMatrix

Return [diamond distance](https://arxiv.org/pdf/1004.4110.pdf) between dynamical matrices `Î¦1` and `Î¦2`.
"""
function diamond_distance(Î¦1::DynamicalMatrix{T}, Î¦2::DynamicalMatrix{T}, args...) where T<:AbstractMatrix{<:Number}
    J1 = Î¦1.matrix
    J2 = Î¦2.matrix
    # TODO: Test dimnesions
    Î¦ = DynamicalMatrix{T}(J1-J2, Î¦1.idim, Î¦1.odim)
    norm_diamond(Î¦, args...)
end

function diamond_distance(Î¦1::AbstractQuantumOperation{T}, Î¦2::AbstractQuantumOperation{T}, args...) where T<:AbstractMatrix{<:Number}
    diamond_distance(convert(DynamicalMatrix{T},Î¦1), convert(DynamicalMatrix{T},Î¦2), args...)
end
