PlotPFC.jl

#* Functions that allow to compare different power flow solutions
#*------------------------------------------------------------------------------

#* NOTE: Utility functions used are defined in PlotUtils.jl

#*------------------------------------------------------------------------------

#=
Scatters branch loads contained in y_ndd versus branch loads in x_ndd. Possible plot settings can be seen in PlotUtils.jl.
=#
function scatter_branchloads(
        x_ndd::Dict{String,<:Any},
        y_ndd::Dict{String,<:Any},
        settings::Dict{String,<:Any}; # plot settings
        figpath = "" # where to save the figure
    )

    ### Set plot settings
    settings = _recursive_merge(
        _default_settings(:scatter_branchloads), settings
    )
    
    ### Get branch loadings from both network data dictionaries
    x_bl_dict = Dict(
        i => br[settings["x_pf_type"]] for (i, br) in x_ndd["branch"] 
        if br["br_status"] == 1
    )
    y_bl_dict = Dict(
        i => br[settings["y_pf_type"]] for (i, br) in y_ndd["branch"]
        if br["br_status"] == 1
    )

    L = length(y_bl_dict)
    x_bl, y_bl = zeros(L), zeros(L) # arrays for branch loads
    todelete = Int64[]
    ### Write branch loads with correct ordering into arrays    
    for (n, key) in enumerate(keys(y_bl_dict))
        ### Check whether branch is also active in x_ndd
        if key ∉ keys(x_bl_dict) # branch inactive in x_ndd
            push!(todelete, n) # n-th entries will be deleted later
            println("Branch $key active in y_ndd and inactive in x_ndd!")
        else # branch also active in y_ndd
            x_bl[n] = x_bl_dict[key]
            y_bl[n] = y_bl_dict[key]
        end
    end
    ### Remove branches that are inactive in one of the ndd's
    deleteat!(x_bl, todelete)
    deleteat!(y_bl, todelete)

    ### Plot straight lines with slope 1 as reference
    x = [i for i in 0:0.01:1]
    plt.plot(x, x, color="g", alpha=1.) # plot line with slope 1

    ### Scatter branch loads versus each other
    bl_diff = y_bl - x_bl # deviation from ratio 1 for coloring

    ### Adjust the range of the colorbar
    if settings["cmap_range"] == "relative" # range according to data
        mcolors = pyimport("matplotlib.colors")
        offset = mcolors.TwoSlopeNorm(
            vcenter=0, vmin=minimum(bl_diff), vmax=maximum(bl_diff)
        ) # match diff to values between 0 and 1
        vmin, vmax = minimum(offset(bl_diff)), maximum(offset(bl_diff))
        coloring = offset(bl_diff)
        coloring[coloring .== offset(0.0)] .= -10 # for unique coloring of 0.0
    elseif settings["cmap_range"] == "full" # range from -1. to 1.
        vmin, vmax = -1., 1.
        offset = norm=plt.Normalize(vmin, vmax)
        coloring = bl_diff
        coloring[coloring .== 0.0] .= -10 # for unique coloring of 0.0
    else
        throw(ArgumentError("Unknown cmap_range $(settings["cmap_range"])."))
    end

    ### Colormap
    pycopy = pyimport("copy")
    cmap = pycopy.copy(plt.get_cmap(settings["cmap"]))
    cmap.set_under(color="tab:green") # color perfect matches green
    
    ### Scatter plot
    sc = plt.scatter(x_bl, y_bl,       
        s = settings["size"],
        c = coloring,
        cmap = cmap,
        vmin = vmin,
        vmax = vmax,
        edgecolors = settings["ec"],
        linewidths = settings["lw"],
        alpha = settings["alpha"]
    )

    ### Plot colorbar
    sm = plt.cm.ScalarMappable(cmap=cmap, offset)
    cbar = plt.colorbar(sm)
    cbar.ax.plot([-1, 1], [0.0, 0.0], "tab:green") # perfect match
    cbar.ax.set_ylabel(
        L"Difference $y-x$", rotation=-90, va="bottom"
    )

    ### Plot labels
    plt.xlabel(settings["xlabel"])
    plt.ylabel(settings["ylabel"])

    ### Plot horizontal and vertical lines to indicate overloads
    plt.axhline(1.0, linestyle="--", color="k", alpha=.4) # horizontal line
    plt.axvline(1.0, linestyle="--", color="k", alpha=.4) # vertical line

    plt.xlim(left=0.)
    plt.ylim(bottom=0.)

    plt.savefig(figpath, bbox_inches="tight") # save figure
    plt.close("all") # close figure
    return nothing
end

#*------------------------------------------------------------------------------

#=
Plots a power grid similar to plot_pg in PlotPG.jl and colors the edges according to the difference in branch loads in ndd and ref_ndd. Possible plot settings can be seen in PlotUtils.jl.
=#
function plot_pg_pf_diff(
        ndd::Dict{String,<:Any}, # NDD with (new) power flow
        ref_ndd::Dict{String,<:Any}, # NDD with reference power flow
        settings = Dict{String,Any}(); # plot settings
        figpath = "" # where to save the figure
    )
    
    ### Setup figure
    figure::Figure, ax::PyObject = plt.subplots()::Tuple{Figure,PyObject}
    w::Float64, h::Float64 = plt.figaspect(2/3)::Vector{Float64}
    figure.set_size_inches(1.5w, 1.5h)
    ax.set_aspect("equal")

    ### Set plot settings and plot power grid
    settings = _recursive_merge(_default_settings(:plot_pg_pf_diff), settings)
    _plot_pg_pf_diff!(ax, ndd, ref_ndd, settings, figpath)

    return nothing
end

function _plot_pg_pf_diff!(
        ax::PyObject, # axes to draw power grid onto
        ndd::Dict{String,<:Any}, # NDD with (new) power flow
        ref_ndd::Dict{String,<:Any}, # NDD with reference power flow
        settings = Dict{String,Any}(), # plot settings
        figpath = "" # where to save the figure
    )
    
    ### Draw power grid graph
    nx = pyimport("networkx")
    G::PyObject = nx.Graph() # empty graph

    ### Plot buses into graph
    G, bus_markers, bus_labels = _draw_buses!(G, ndd, settings)

    ### Plot branches into graph
    G, br_markers, br_labels, br_cbar = _draw_br_pf_diff!(
        G, ndd, ref_ndd, settings
    )

    ### Check for a predefined area to show
    if haskey(settings, "area")
        area = settings["area"]
        plt.xlim(area[1], area[2])
        plt.ylim(area[3], area[4])
    end

    ### Axes settings
    ax.tick_params(
        left = settings["draw_ticks"][1],
        bottom = settings["draw_ticks"][2], 
        labelleft = settings["draw_ticks"][3], 
        labelbottom = settings["draw_ticks"][4]
    )
    plt.xlabel(settings["xlabel"])
    plt.ylabel(settings["ylabel"], rotation=90)

    ### Draw legend, if wanted
    if settings["draw_legend"] == true
        all_markers = vcat(bus_markers, br_markers)
        all_labels = vcat(bus_labels, br_labels)
        plt.legend(all_markers, all_labels)
    end

    plt.savefig(figpath, bbox_inches="tight") # save figure
    plt.close("all") # close figure
    return ax, G
end

function _draw_br_pf_diff!(
        G::PyObject, # power grid graph
        ndd::Dict{String,<:Any}, # NDD with (new) power flow
        ref_ndd::Dict{String,<:Any}, # NDD with reference power flow
        settings::Dict{String,<:Any} # dictionary containing plot settings
    )
        
    pos = Dict(
        b["index"] => (b["bus_lon"], b["bus_lat"]) 
        for b in collect(values(ndd["bus"]))
    ) # geographic bus locations
    branches = collect(values(ndd["branch"])) # branch dictionaries

    ### Get edges contained in the NDD and the loading difference
    edges = Array{Tuple{Int64,Int64},1}() # array for edges
    bl_diff = Array{Float64,1}() # array for branch load differences
    for (i, br) in ndd["branch"]
        ### Check whether branch is active in both ndd's
        if br["br_status"] == 1 && ref_ndd["branch"][i]["br_status"] == 1
            push!(edges, (br["f_bus"], br["t_bus"]))
            pf_type = settings["Branches"]["pf_type"]
            ref_pf_type = settings["Branches"]["ref_pf_type"]
            push!(bl_diff, br[pf_type] - ref_ndd["branch"][i][ref_pf_type])
        end
    end

    ### Draw edges and color them according to the difference in loading
    ### Adjust the range of the colorbar
    if settings["Branches"]["cmap_range"] == "relative" # same range as data
        mcolors = pyimport("matplotlib.colors")
        offset = mcolors.TwoSlopeNorm(
            vcenter=0, vmin=minimum(bl_diff), vmax=maximum(bl_diff)
        ) # match diff to values between 0 and 1
        vmin, vmax = minimum(offset(bl_diff)), maximum(offset(bl_diff))
        coloring = offset(bl_diff)
        coloring[coloring .== offset(0.0)] .= -10 # for unique coloring of 0.0
    elseif settings["Branches"]["cmap_range"] == "full" # range from -1. to 1.
        vmin, vmax = -1., 1.
        offset = norm=plt.Normalize(vmin, vmax)
        coloring = bl_diff
        coloring[coloring .== 0.0] .= -10 # for unique coloring of 0.0
    else
        cmap_range = settings["Branches"]["cmap-range"]
        throw(ArgumentError("Unknown cmap_range $cmap_range."))
    end

    pycopy = pyimport("copy")
    cmap = pycopy.copy(plt.get_cmap(settings["Branches"]["cmap"]))
    cmap.set_under(color="tab:green") # color perfect matches green

    nx = pyimport("networkx")
    drawnedges = nx.draw_networkx_edges(
        G, pos, 
        width = settings["Branches"]["br_lw"],
        edgelist = edges,
        edge_color = coloring,
        edge_vmin = vmin,
        edge_vmax = vmax,
        edge_cmap = cmap
    )

    ### Add colorbar
    sm = plt.cm.ScalarMappable(cmap=cmap, offset)
    cbar = plt.colorbar(sm)
    cbar.ax.plot([-1, 1], [0.0, 0.0], "tab:green") # perfect match
    cbar.ax.set_ylabel(
        settings["Branches"]["cbar_label"], rotation=-90, va="bottom"
    )

    return G, [], [], cbar
end