using XLSX

nM = XLSX.readxlsx("NetMatrix.xlsx")["Sheet1"] |> Matrix
fM = XLSX.readxlsx("FlowMatrix.xlsx")["Sheet1"] |> Matrix
N, nFlows, OFQM, qLM = OptiFlows(nM, fM)
nl = flowDir(fM)

if N > 0
    println("Optimization completed.")
    for n in 1:nFlows
        XLSX.writetable("OptiFlows.xlsx", (nl[1, n], "⟶", nl[2, n]), OFQM[:, :, n], sheetname = string(nl[1, n], "⟶", nl[2, n]))
    end
else 
    println("Optimization impossible. Need to increase the bandwidth of the channels.")
end

XLSX.writetable("LoadFactors.xlsx", qLM)

function flowDir(nM)
    qN = size(nM, 1)
    nfl = count(x -> x > 0, nM)
    nl = zeros(Int, 2, nfl)
    c = 1
    for a in 1:qN
        for b in 1:qN
            if nM[a, b] > 0
                nl[1, c] = a
                nl[2, c] = b
                c += 1
            end
        end
    end
    return nl
end

using Optim
function OptiFlows(nM, fM)
    qn = length(nM)
    nEdges = count(!iszero, nM)
    nFlows = count(!iszero, fM)
    source, target, Mu = opti1(nM, qn, nEdges)
    iFlows = opti2(qn, nEdges, source, target)
    fsource, ftarget, flows = opti3(fM, qn, nFlows)
    flowSourceTarget = opti4(qn, nFlows, fsource, ftarget, flows)
    stLimit = opti5(nEdges, nFlows, source, target, fsource, ftarget)
    QFlows = opti6(qn, iFlows, nEdges, nFlows, stLimit, flowSourceTarget)
    Aeq, Beq, Lambda0, A, B = opti7(qn, nFlows, nEdges, QFlows, flowSourceTarget, Mu)
    Lambda, N = opti8(Aeq, Beq, Lambda0, A, B, Mu, nFlows)
    OFQM, qLM = opti9(nFlows, nEdges, Lambda, qn, source, target, Mu)


function opti1(nM, qn, nEdges)
    target = zeros(Int, nEdges)
    source = zeros(Int, nEdges)
    Mu = zeros(Int, nEdges)
    u = 1
    for p in 1:qn
        for q in 1:qn
            if nM[p, q] > 0
                source[u] = p
                target[u] = q
                Mu[u] = nM[p, q]
                u += 1
            end
        end
    end
    return source, target, Mu
end

function opti2(qn, nEdges, source, target)
    iFlows = zeros(Int, qn, nEdges)
    for l in 1:nEdges
        iFlows[source[l], l] = -1
        iFlows[target[l], l] = 1
    end
    return iFlows
end

function opti3(fM, qn, nFlows)
    ftarget = zeros(Int, nFlows)
    fsource = zeros(Int, nFlows)
    flows = zeros(Int, nFlows)
    u = 1
    for p in 1:qn
        for q in 1:qn
            if fM[p, q] > 0
                fsource[u] = p
                ftarget[u] = q
                flows[u] = fM[p, q]
                u += 1
            end
        end
    end
    return fsource, ftarget, flows
end

function opti4(qn, nFlows, fsource, ftarget, flows)
    fST = zeros(Int, qn, nFlows)
    for l in 1:nFlows
        fST[fsource[l], l] = -flows[l]
        fST[ftarget[l], l] = flows[l]
    end
    flowSourceTarget = zeros(Int, qn + 2, nFlows)
    flowSourceTarget[1:qn, :] .= fST
    return flowSourceTarget
end

function opti5(nEdges, nFlows, source, target, fsource, ftarget)
    stLimit = zeros(Int, 2, nEdges, nFlows)
    for v in 1:nFlows
        for w in 1:nEdges
            if source[w] == ftarget[v]
                stLimit[1, w, v] = 1
            elseif target[w] == fsource[v]
                stLimit[2, w, v] = 1
            end
        end
    end
    return stLimit
end

function opti6(qn, iFlows, nEdges, nFlows, stLimit, flowSourceTarget)
    QFlows = zeros(Int, qn + 2, nEdges, nFlows)
    for o in 1:nFlows
        QFlows[1:qn, :, o] .= iFlows
        QFlows[qn + 1:qn + 2, :, o] .= stLimit[:, :, o]
    end
    for k in 1:nFlows
        for i in 1:qn
            if flowSourceTarget[i, k] < 0
                for j in 1:nEdges
                    if QFlows[i, j, k] == 1
                        QFlows[i, j, k] = 0
                    end
                end
            elseif flowSourceTarget[i, k] > 0
                for j in 1:nEdges
                    if QFlows[i, j, k] == -1
                        QFlows[i, j, k] = 0
                    end
                end
            end
        end
    end
    return QFlows
end

function opti7(qn, nFlows, nEdges, QFlows, flowSourceTarget, Mu)
    Aeq = zeros(Int, (qn + 2) * nFlows, nEdges * nFlows)
    Beq = zeros(Int, (qn + 2) * nFlows)
    for j in 0:nFlows - 1
        Aeq[(qn + 2) * j + 1:(qn + 2) * (j + 1), nEdges * j + 1:nEdges * (j + 1)] .= QFlows[:, :, j + 1]
        Beq[(qn + 2) * j + 1:(qn + 2) * (j + 1)] .= flowSourceTarget[:, j + 1]
    end
    Lambda0 = zeros(Int, nFlows * nEdges)
    A = zeros(Int, nEdges, nEdges * nFlows)
    for z in 1:nEdges
        for y in 0:nFlows - 1
            A[z, z + y * nEdges] = 1
        end
    end
    B = zeros(Int, nEdges, 1)
    B[:, 1] .= Mu
    return Aeq, Beq, Lambda0, A, B
end

function opti8(Aeq, Beq, Lambda0, A, B, Mu, nFlows)
    options = Optim.Options(show_trace=false, iterations=2^18, tol=2^-32)
    result = optimize(ObjF, Lambda0, A, B, Aeq, Beq, options)
    Lambda = result.minimizer
    N = result.minimum

    function ObjF(LambdaFlows)
        fS = zeros(length(Mu))
        AvgQ = 0.0
        for i in 1:length(Mu)
            for j in 0:nFlows-1
                fS[i] += LambdaFlows[length(Mu) * j + i]
            end
        end
        for k in 1:length(Mu)
            AvgQ += (fS[k] / Mu[k]) / (1 - fS[k] / Mu[k])
        end
        return AvgQ
    end
end

function opti9(nFlows, nEdges, Lambda, qn, source, target, Mu)
    oFlows = zeros(nFlows, nEdges)
    for n in 0:nFlows-1
        oFlows[n + 1, :] .= Lambda[n * nEdges + 1:nEdges * (n + 1)]
    end
    for i in 1:nFlows
        for j in 1:nEdges
            if oFlows[i, j] < 10^-10
                oFlows[i, j] = 0
            end
        end
    end
    OFQM = zeros(qn, qn, nFlows)
    for k in 1:nFlows
        for x in 1:nEdges
            OFQM[source[x], target[x], k] = oFlows[k, x]
        end
    end
    qL = sum(oFlows, dims=1) ./ Mu
    qLM = zeros(qn, qn)
    for v in 1:nEdges
        qLM[source[v], target[v]] = qL[v]
    end
    return OFQM, qLM
end

    return N, nFlows, OFQM, qLM
end

using ExcelReaders
using Graphs
using Plots
using LightGraphs


nM = readdlm("NetMatrix.xlsx", ",")
fM = readdlm("FlowMatrix.xlsx", ",")
lM = readdlm("LoadFactors.xlsx", ",")
nfl, nl, fG = flowGraph(fM)
nG = netGraph(nM)
sG = simpGraph(nM)
fnums = reshape(nl, (2, nfl))
oM = zeros(size(nM, 1), size(nM, 2), nfl)

for n in 1:nfl
    oM[:, :, n] = readdlm("OptiFlows.xlsx", n, ",")
end

plot(sG, node_color=[245/256 64/256 33/256], edge_font_size=10, edge_color=[127/256 255/256 212/256], linewidth=1.5, edge_alpha=1, node_font_size=12, node_label_color=[1 1 1])
title!("Simplified network graph")
plot!(size=(800, 600), legend=false)
savefig("OptiFlows.pdf", append=true)

plot(fG, node_label=nl, edge_label=weights(fG), node_color=[245/256 64/256 33/256], arrow_size=13, edge_font_size=12, edge_color=[127/256 255/256 212/256], linewidth=1.5, edge_alpha=1, node_font_size=12, node_label_color=[1 1 1], edge_label_color=[0 1 0])
title!("Flow directions")
plot!(size=(800, 600), legend=false)
savefig("OptiFlows.pdf", append=true)

plot(nG, layout=:force, edge_label=weights(nG), node_color=[245/256 64/256 33/256], edge_font_size=10, edge_color=[127/256 255/256 212/256], arrow_size=10, node_font_size=12, node_label_color=[1 1 1], linewidth=1.5, marker_size=4, edge_label_color=[0 1 0], edge_alpha=1)
title!("Network graph")
plot!(size=(800, 600), legend=false)
savefig("OptiFlows.pdf", append=true)

plot(lG, layout=:force, edge_label=weights(lG), node_color=[245/256 64/256 33/256], edge_font_size=10, edge_color=[127/256 255/256 212/256], arrow_size=10, node_font_size=12, node_label_color=[1 1 1], linewidth=1.5, marker_size=4, edge_label_color=[0 1 0], edge_alpha=1)
title!("Load factors of communication lines")
plot!(size=(800, 600), legend=false)
savefig("OptiFlows.pdf", append=true)

for m in nfl:-1:1
    optiname = "Optimal flows from node $(fnums[1, m]) to node $(fnums[2, m])"
    oG = netGraph(oM[:, :, m])
    plot(oG, layout=:force, edge_label=weights(oG), node_color=[245/256 64/256 33/256], edge_font_size=10, edge_color=[127/256 255/256 212/256], arrow_size=10, node_font_size=12, node_label_color=[1 1 1], linewidth=1.5, marker_size=4, edge_label_color=[1 1 0], edge_alpha=1)
    title!(optiname)
    plot!(size=(800, 600), legend=false)
    savefig("OptiFlows.pdf", append=true)
end

function netGraph(nM)
    qN = size(nM, 1)
    t = zeros(Int, nnz(nM))
    s = t
    bW = t
    u = 1
    for p in 1:qN
        for q in 1:qN
            if nM[p, q] > 0
                s[u] = p
                t[u] = q
                bW[u] = nM[p, q]
                u += 1
            end
        end
    end
    return SimpleDiGraph(s, t, bW)
end

function simpGraph(nM)
    qN = size(nM, 1)
    rM = nM ./ nM

    for a in 1:qN
        for b in 1:qN
            if isnan(rM[a, b])
                rM[a, b] = 0
            end
        end
    end

    ruM = triu(rM)
    trM = rM - ruM
    rrM = zeros(qN, qN)
    for p in 1:qN
        for q in 1:qN
            rrM[q, p] = trM[p, q]
        end
    end
    grM = rrM + ruM

    t = zeros(Int, nnz(grM))
    s = t
    bW = t
    u = 1
    for p in 1:qN
        for q in 1:qN
            if grM[p, q] > 0
                s[u] = p
                t[u] = q
                bW[u] = grM[p, q]
                u += 1
            end
        end
    end
    return SimpleGraph(s, t)
end


function flowGraph(nM)
    qN = size(nM, 1)
    nfl = count(!iszero, nM)
    nl = zeros(Int, 2 * nfl)
    flG = zeros(Int, 2, nfl)
    flG[1, :] .= 1:2:2*nfl-1
    flG[2, :] .= 2:2:2*nfl
    flI = zeros(Int, nfl)
    c = 1
    for a in 1:qN
        for b in 1:qN
            if nM[a, b] > 0
                nl[2*c - 1] = a
                nl[2*c] = b
                flI[c] = nM[a, b]
                c += 1
            end
        end
    end
    fG = DiGraph(flG[1, :], flG[2, :], flI)

    return nfl, nl, fG
end