using LinearAlgebra
using CVChannel
using Convex
using Test
using DelimitedFiles

"""
This code first verifies the non-multiplicativity of the Werner-Holevo over the PPT
cone when tensored with the identity. This is different than over the separable cone.
We verify this using an a linear program (LP) that we scan over.
"""
#This script looks at multiplicativity of WH with identity over PPT cone
println("\nHere we initialize the identity channel and the solver. (One moment please...)")
identChan(X) = X
eaCVDual(choi(identChan,2,2),2,2)

println("\nNow we look at tensoring the Werner-Holevo channel with the identity.")
println("\nWhat makes this interesting is it is multiplcative over the separable")
println("cone, but over the PPT it is non-multiplicative. This exemplifies a separation.\n\n")

function nonMultWHID()
    println("\nNow we look at the Werner-Holevo channels tensored with the identity")
    println("\nfor dimension [2,8] for both channels for values of λ in [0.4,0.5,...,1]")

    table_vals = zeros(343,7)
    λ_vals = 0:0.1:0.6
    d1_range = 2:8
    d2_range = 2:8
    λ_ctr = 1
    table_ctr = 1
    for λ in λ_vals
        println("\nNow evaluating for Werner-Holevo channel with parameter = ", round(λ, digits = 2))
        #Note the LP is parameterized backwards with respect to how our SDP is
        for d1 in d1_range
            println("Now evaluating for d = ", d1, "...")
            #Here we get the communication value of the Werner-Holevo channel itself
            A,B,g,a = generalWHLPConstraints(1,3,(1-λ)*ones(1))
            cv_WH, v1 =wernerHolevoCVPPT(1,3,A,B,g,a)

            for d2 in d2_range
                cv = WHIDLP(d1,d2,λ)
                table_vals[table_ctr,:] = [
                    λ,            #parameter
                    d1,             #WH dimension
                    d2,             #identity channel dimension
                    cv_WH,          #comm val of WH
                    cv_WH*d2,       #cv(𝒩_WH ⊗ id)
                    cv,             #cv^{PPT}(𝒩_WH ⊗ id)
                    cv - cv_WH*d2   #cv^{PPT} - cv
                ];
                table_ctr = table_ctr + 1;
            end
        end
        λ_ctr = λ_ctr + 1
    end

    header = ["param" "WH Dim" "id Dim" "cv(WH)" "cv" "cv_ppt" "diff"]
    data_to_save = vcat(header,table_vals)
    info_vec = Vector{Union{Nothing,String}}(nothing, size(data_to_save)[1])
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by werner-holevo-channel-with-identity.jl"
    info_vec[3] = "WH dim range: " * string(d1_range)
    info_vec[4] = "id dim range: " * string(d2_range)
    info_vec[5] = "WH parameter range: " * string(λ_vals)
    data_to_save = hcat(data_to_save,info_vec)

    print("\nPlease name the file you'd like to write the results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')
    if all(table_vals -> table_vals < 1e-4 , table_vals[2:end,7])
        print("\nSomething went awry. Please check your results.")
        return false
    else
        #Here we construct the highlights of the data
        println("Here are the highlights.")
        println("There is no multiplicativity for qubit to qubit, but qutrit to qutrit there is.")
        qubit_to_qutrit = vcat(header,table_vals[6:9,:])
        show(stdout, "text/plain", qubit_to_qutrit)
        println("For p >= 0.5, there is multiplicativity. This is when Werner-Holevo becomes PPT.")
        p_geq_half = vcat(header,
                        [
                            table_vals[49,:]';table_vals[98,:]';table_vals[147,:]';
                            table_vals[196,:]';table_vals[245,:]';table_vals[294,:]';
                            table_vals[343,:]'
                        ])
        show(stdout, "text/plain", p_geq_half)
        println("\nFor the full data, please see the saved results.")
        return true
    end
    end

@testset "non-multiplicativity" begin
    @test nonMultWHID()
end
