text:
"""
Plot the results for a 2d solution.
"""
function plot_2d_verification_results(results_path; num_dfa_states=1, min_threshold=0.95, extents_dict=Dict(), extents_labels_colors=Dict(), prob_plots=false, lazy_flag=false, plot_tag="", state_outlines_flag=false, trajectories=nothing, dfa_init_state=1)
    extents, _ = deserialize_region_data("$results_path/regions.bson", lazy_flag=lazy_flag)
    num_regions = extents.count - 1
    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    mat_files = filter(x -> occursin(r"globally-safe-.*.mat", x), readdir(results_path))

    for mat_file in mat_files
        res = matread("$results_path/$mat_file")
        base_str = mat_file[1:end-4]
        maxPrs = res["indVmax"][dfa_init_state:num_dfa_states:end]
        minPrs = res["indVmin"][dfa_init_state:num_dfa_states:end]

        if prob_plots
            # Plot the maximumprobabilities
            plt_max = plot(aspect_ratio=1,
                        size=(300,300), dpi=300,
                        xlims=[minx, maxx], ylims=[miny, maxy],
                        xtickfont=font(10),
                        ytickfont=font(10),
                        titlefont=font(10),
                        xticks = [minx, 0, maxx],
                        yticks = [miny, 0, maxy],
                        grid=false)
            [plot_cell(extents[i], maxPrs[i], state_outlines_flag=state_outlines_flag) for i in 1:num_regions]
            plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
            savefig(plt_max, "$results_path/$base_str$plot_tag-max-heatmap.png")

            # Plot the minimum probabilities
            plt_min = plot(aspect_ratio=1,
                        size=(300,300), dpi=300,
                        xlims=[minx, maxx], ylims=[miny, maxy],
                        xtickfont=font(10),
                        ytickfont=font(10),
                        titlefont=font(10),
                        xticks = [minx, 0, maxx],
                        yticks = [miny, 0, maxy],
                        grid=false,
                        backgroundcolor=128)
            [plot_cell(extents[i], minPrs[i], state_outlines_flag=state_outlines_flag) for i in 1:num_regions]
            plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
            savefig(plt_min, "$results_path/$base_str$plot_tag-min-heatmap.png")
        end

        plt_verification = plot(aspect_ratio=1,
                                size=(300,300), dpi=300,
                                xlims=[minx, maxx], ylims=[miny, maxy],
                                xtickfont=font(10),
                                ytickfont=font(10),
                                titlefont=font(10),
                                xticks = [minx, 0, maxx],
                                yticks = [miny, 0, maxy],
                                grid=false,
                                backgroundcolor=128)

        # Plot the cells
        [plot_cell_verify(extents[i], minPrs[i], maxPrs[i], min_threshold, state_outlines_flag=state_outlines_flag) for i in 1:num_regions]
        plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")

        # Plot extents with labels
        for key in keys(extents_dict)
            labelled_extents = extents_dict[key]
            for extent in labelled_extents
                plot_labelled_extent_outline(extent, key, color=extents_labels_colors[key])
            end
        end

        savefig(plt_verification, "$results_path/$base_str$plot_tag-verification.png")

        if !isnothing(trajectories)
            for trajectory in trajectories
                [plot!([trajectory[i][1], trajectory[i+1][1]], [trajectory[i][2], trajectory[i+1][2]], color=:black, label="", linewith=2) for i=1:length(trajectory)-1]
                scatter!([trajectory[end][1]], [trajectory[end][2]], color=:black, markershape=:star5, label="")
            end
            savefig(plt_verification, "$results_path/$base_str$plot_tag-simulations.png")
        end
    end
end

"""
Plot a labelled extent
"""
function plot_labelled_extent_outline(extent, label; color=:white)
    xvals = extent["x1"]
    yvals = extent["x2"]
    plot!(extent_dict_to_shape(extent), fillalpha=0, linecolor=:black, linewidth=1, label="")
    annotate!((minimum(xvals)+maximum(xvals))/2, (minimum(yvals)+maximum(yvals))/2, text(label, color, :center, 10,))
end

"""
Plot a cell with alpha associated with probability value
"""
function plot_cell(extent, prob_value; state_outlines_flag=false)
    x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
    y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
    shape = Plots.Shape(x, y)
    linealpha = state_outlines_flag ? 1.0 : 0.0
    plot!(shape, color=:black, fillalpha=1-prob_value, linewidth=1, linecolor=:black, linealpha=linealpha, foreground_color_border=:white, foreground_color_axis=:white, label="")
end

"""
Determine the cell color based on the verification probability threshold.
"""
function plot_cell_verify(extent, min_prob_value, max_prob_value, threshold; state_outlines_flag=false)
    x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
    y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
    shape = Plots.Shape(x, y)

    linealpha = state_outlines_flag ? 1.0 : 0.0

    if min_prob_value >= threshold
        plot!(shape, color=:purple, fillalpha=0.10, linewidth=1, linecolor=:black, linealpha=linealpha, foreground_color_border=:white, foreground_color_axis=:white, label="")
    elseif max_prob_value < threshold
        plot!(shape, color=:purple, fillalpha=0.99, linewidth=1, linecolor=:black, linealpha=linealpha, foreground_color_border=:white, foreground_color_axis=:white, label="")
    else
        plot!(shape, color=:purple, fillalpha=0.50, linewidth=1, linecolor=:black, linealpha=linealpha, foreground_color_border=:white, foreground_color_axis=:white, label="")
    end
end

function plot_gamma_value(results_path, min_ResMat, max_ResMat; num_dfa_states=1)
    function plot_cell(extent, prob_value)
        x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
        y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
        shape = Plots.Shape(x, y)
        plot!(shape, color=:black, fillalpha=1-prob_value, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
    end

    extents, post_extents = deserialize_region_data("$results_path/regions.bson")
    num_regions = extents.count - 1
    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    minPrs = min_ResMat[1:num_dfa_states:end, 3]
    maxPrs = max_ResMat[1:num_dfa_states:end, 4]

    base_str = "gamma-values"
    plt = plot(#title=L"$k Step Minimum Safety",
                       aspect_ratio=1,
                       size=(300,300), dpi=300,
                       xlims=[minx, maxx], ylims=[miny, maxy],
                       #xlabel="X Units", ylabel="Y Units",
                       xtickfont=font(10),
                       ytickfont=font(10),
                       titlefont=font(10),
                       grid=false,
                       backgroundcolor=128)
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
    [plot_cell(extents[i], abs(minPrs[i]-maxPrs[i])) for i in 1:num_regions]
    savefig(plt, "$results_path/$base_str.png")
end

function plot_gp_fields(results_path, dyn_fn)
    files = readdir("$results_path/gps")
    gp_file = filter(x->occursin(".bin", x), files)[1]
    f = open("$results_path/gps/$gp_file")
    gp_save_info = deserialize(f)
    close(f)
    extents, post_extents = deserialize_region_data("$results_path/regions.bson")

    if length(gp_save_info) == 1
        plot_gp_1d(gp_set, region_data, results_path, dyn_fn)
        return
    elseif length(gp_save_info) > 2
        plot_gp_field_slice(results_path, dyn_fn, 0.0; nl_flag=false)
        return
    end

    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    x_range = minx:0.25:maxx
    y_range = miny:0.25:maxy
    meshgrid(x, y) = (repeat(x, outer=length(y)), repeat(y, inner=length(x)))
    x, y = meshgrid(x_range, y_range)
    gp_x1 = gp_save_info[:gp_set]["x1"]
    gp_x2 = gp_save_info[:gp_set]["x2"]
    delta_x_gp = [x[i] - predict_y(gp_x1, hcat([x[i], y[i]]))[1][1] for i in 1:length(x)]
    delta_y_gp = [y[i] - predict_y(gp_x2, hcat([x[i], y[i]]))[1][1] for i in 1:length(x)]
    u_gp = delta_x_gp*0.5
    v_gp = delta_y_gp*0.5
    plt1 = plot(size=(300,300), dpi=600, aspect_ratio=1, grid=false, foreground_color_border=:white, foreground_color_axis=:white,
                xlims = [minx, maxx], ylims = [miny, maxy],
                xticks = [minx, 0, maxx], yticks = [miny, 0, maxy],
                xtickfont=font(10), ytickfont=font(10), titlefont=font(10))
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="", linealpha=0.1)
    quiver!(x, y, quiver=(-u_gp, -v_gp), color=:black, linewidth=1.5)
    plt2 = plot(size=(300,300), dpi=600, aspect_ratio=1, grid=false, foreground_color_border=:white, foreground_color_axis=:white,
                xlims = [minx, maxx], ylims = [miny, maxy],
                xticks = [minx, 0, maxx], yticks = [miny, 0, maxy],
                xtickfont=font(10), ytickfont=font(10), titlefont=font(10))
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="", linealpha=0.1)

    u_true = 0.5*[x[i] - dyn_fn([x[i], y[i]])[1] for i in 1:length(x)]
    v_true = 0.5*[y[i] - dyn_fn([x[i], y[i]])[2] for i in 1:length(x)]
    quiver!(x, y, quiver=(-u_true, -v_true), color=:black, linewidth=1.5)
    plt = plot(plt1, plt2, layout=(1,2))
    savefig(plt1, "$results_path/gp-vector-fields.png")
    savefig(plt2, "$results_path/true-vector-fields.png")

    plt = plot(gp_x1, var=true, title="Variance of GP1", fill=true)
    savefig(plt, "$results_path/gp1-var.pdf")
    plt = plot(gp_x2, var=true, title="Variance of GP2", fill=true)
    savefig(plt, "$results_path/gp2-var.pdf")

    # Plot where the data was collected and a convex hull
    data_shape = VPolygon(convex_hull([gp_x1.x[:, i] for i=1:gp_x1.nobs]))
    plt = plot(data_shape, fillalpha=0.5)
    [scatter!([gp_x1.x[1, i]], [gp_x1.x[2,i]], label="", color=:black) for i=1:gp_x1.nobs]
    savefig(plt, "$results_path/data-hull.pdf")
end

function plot_gp_1d(gp_set, region_data, results_path, dyn_fn)
    extents = region_data[:extents]
    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    x = collect(minx:0.25:maxx)
    gp_x1 = gp_set["x1"]
    delta_x_gp = [x[i] - predict_y(gp_x1,[x[i]])[1][1] for i in 1:length(x)]

    plt1 = plot(gp_x1)
    savefig(plt1, "$results_path/gp-vector-fields.png")
    plt = plot(gp_x1, var=true, title="Variance of GP1")
    savefig(plt, "$results_path/gp1-var.pdf")
end

function plot_gp_field_slice(results_path, dyn_fn, x3; nl_flag=false)
    f = open("$results_path/gaussian_processes_set.bin")
    gp_set = deserialize(f)
    close(f)
    #A_mats = collect(keys(gp_set))
    #A = A_mats[1]
    #A = [0.9 0.4; -0.4 0.5] # Rotation
    #gp_dict = gp_set[A_mats[1]]

    extents, post_extents = deserialize_region_data("$results_path/regions.bson")
    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    x_range = minx:0.25:maxx
    y_range = miny:0.25:maxy
    meshgrid(x, y) = (repeat(x, outer=length(y)), repeat(y, inner=length(x)))
    y, x = meshgrid(x_range, y_range)
    gp_x1 = gp_set["x1"]
    gp_x2 = gp_set["x2"]
    delta_x_gp = [x[i] - predict_y(gp_x1, hcat([x[i], y[i], x3]))[1][1] for i in 1:length(x)]
    delta_y_gp = [y[i] - predict_y(gp_x2, hcat([x[i], y[i], x3]))[1][1] for i in 1:length(x)]
    u_gp = delta_x_gp*0.5
    v_gp = delta_y_gp*0.5
    plt1 = plot(size=(300,300), dpi=600, aspect_ratio=1, grid=false, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], xtickfont=font(10), ytickfont=font(10), titlefont=font(10))
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
    quiver!(x, y, quiver=(-u_gp, -v_gp), color=:black, linewidth=1.5)
    plt2 = plot(size=(300,300), dpi=600, aspect_ratio=1, grid=false, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], xtickfont=font(10), ytickfont=font(10), titlefont=font(10))
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="", linealpha=0.1)

    u_true = 0.5*[x[i] - dyn_fn([x[i], y[i], x3])[1] for i in 1:length(x)]
    v_true = 0.5*[y[i] - dyn_fn([x[i], y[i], x3])[2] for i in 1:length(x)]
    quiver!(x, y, quiver=(-u_true, -v_true), color=:black, linewidth=1.5)
    plt = plot(plt1, plt2, layout=(1,2))
    savefig(plt1, "$results_path/gp-vector-fields-$x3.png")
    savefig(plt2, "$results_path/true-vector-fields-$x3.png")

    # plt = plot(gp_x1, var=true, title="Variance of GP1", fill=true)
    # savefig(plt, "$results_path/gp1-var.pdf")
    # plt = plot(gp_x2, var=true, title="Variance of GP2", fill=true)
    # savefig(plt, "$results_path/gp2-var.pdf")

    # Plot where the data was collected and a convex hull
    # data_shape = VPolygon(convex_hull([gp_x1.x[:, i] for i=1:gp_x1.nobs]))
    # plt = plot(data_shape, fillalpha=0.5)
    # [scatter!([gp_x1.x[1, i]], [gp_x1.x[2,i]], label="", color=:black) for i=1:gp_x1.nobs]
    # savefig(plt, "$results_path/data-hull.pdf")
end

function plot_synthesis_results(results_path, res_mat, imdp, dfa, pimdp; plot_gp=false, min_threshold=0.95, trajectories=nothing, filename=nothing, old_mins=nothing)
    extents, post_extents = deserialize_region_data("$results_path/regions.bson")
    num_regions = extents.count - 1
    X = extents[extents.count]

    if X.count == 2
        plot_2d_verification_results(results_path, num_dfa_states=length(dfa.states), trajectories=trajectories)
        return
    end

    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    # mat_files = filter(x -> occursin(r"globally-safe-.*.mat", x), readdir(results_path))

    # if plot_gp
    #     f = open("$results_path/gaussian_processes_set.bin")
    #     gp_set = deserialize(f)
    #     gp_x1 = gp_set["x1"]
    #     close(f)
    # end

    # for mat_file in mat_files

    # res = matread("$results_path/$mat_file")
    # k = occursin("--1", mat_file) ? -1 : parse(Int, split(mat_file, "-")[3])
    # base_str = mat_file[1:end-4]
    # Plot the maximum results
    # Probabilities for each cell from 0 to 1 corresponding to regions 1 to N
    # TODO: Need to somehow generalize this skipping step
    Qsize = length(dfa.states)
    maxPrs = res_mat[1:Qsize:end, 4]
    minPrs = res_mat[1:Qsize:end, 3]
    min_minPr = minimum(minPrs)
    min_maxPr = minimum(maxPrs)

    # Get the number of x3 extents
    x3_extents = []

    for region_id in keys(extents)
        if !in(extents[region_id]["x3"], x3_extents)
            push!(x3_extents, extents[region_id]["x3"])
        end
    end

    function plot_cell(extent, prob_value, norm_val)
        x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
        y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
        shape = Plots.Shape(x, y)
        plot!(shape, color=:black, fillalpha=1-prob_value, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
    end

    for (j, x3_extent) in enumerate(x3_extents)
        indeces = []
        for i=1:num_regions
            if extents[i]["x3"] == x3_extent
                push!(indeces, i)
            end
        end
        base_str = "synthesis-res-x3$j"

        plt_max = plot(#title="$k Step Maximum Safety",
                        aspect_ratio=1,
                        size=(300,300), dpi=300,
                        xlims=[minx, maxx], ylims=[miny, maxy],
                        #xlabel="X Units", ylabel="Y Units",
                        xtickfont=font(10),
                        ytickfont=font(10),
                        titlefont=font(10),
                        grid=false)
        plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
        [plot_cell(extents[i], maxPrs[i], min_maxPr) for i in indeces]
        savefig(plt_max, "$results_path/$base_str-max-heatmap.png")
        # Plot the minimum results
        plt_min = plot(#title=L"$k Step Minimum Safety",
                        aspect_ratio=1,
                        size=(300,300), dpi=300,
                        xlims=[minx, maxx], ylims=[miny, maxy],
                        #xlabel="X Units", ylabel="Y Units",
                        xtickfont=font(10),
                        ytickfont=font(10),
                        titlefont=font(10),
                        grid=false,
                        backgroundcolor=128)
        plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
        [plot_cell(extents[i], minPrs[i], min_minPr) for i in indeces]
        savefig(plt_min, "$results_path/$base_str-min-heatmap.png")

        plt_min = plot(#title=L"$k Step Minimum Safety",
                        aspect_ratio=1,
                        size=(300,300), dpi=300,
                        xlims=[minx, maxx], ylims=[miny, maxy],
                        #xlabel="X Units", ylabel="Y Units",
                        xtickfont=font(10),
                        ytickfont=font(10),
                        titlefont=font(10),
                        grid=false,
                        backgroundcolor=128)

        function plot_cell_verify(extent, min_prob_value, max_prob_value, threshold)
            x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
            y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
            shape = Plots.Shape(x, y)

            if min_prob_value >= threshold
                plot!(shape, color=:purple, fillalpha=0.10, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
            elseif max_prob_value < threshold
                plot!(shape, color=:purple, fillalpha=0.99, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
            else
                plot!(shape, color=:purple, fillalpha=0.50, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
            end
        end
        # Plot the cells
        plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
        [plot_cell_verify(extents[i], minPrs[i], maxPrs[i], min_threshold) for i in indeces]
        savefig(plt_min, "$results_path/$base_str-verification.png")

        if !isnothing(trajectories)
            for trajectory in trajectories
                [plot!([trajectory[i][1], trajectory[i+1][1]], [trajectory[i][2], trajectory[i+1][2]], color=:black, label="", linewith=2) for i=1:length(trajectory)-1]
                scatter!([trajectory[end][1]], [trajectory[end][2]], color=:black, markershape=:star5)
            end
            savefig(plt_min, "$results_path/$filename")
        end
    end
    # end
end

function extent_dict_to_shape(extent)
    xmin = minimum(extent["x1"])
    xmax = maximum(extent["x1"])
    ymin = minimum(extent["x2"])
    ymax= maximum(extent["x2"])

    x = [xmin, xmin, xmax, xmax]
    y = [ymin, ymax, ymax, ymin]

    return Plots.Shape(x, y)
end

"""
Plot the results from file using the legacy format.
"""
function plot_online_results(extents, res_mat, plots_basename, num_dfa_states; min_threshold=0.95, trajectories=nothing, trajectory_dict=nothing, extents_dict=nothing, online_data_dict=nothing)
    num_regions = extents.count - 1
    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    # if plot_gp
    #     files = readdir("$results_path/gps")
    #     gp_file = filter(x->occursin(".bin", x), files)[1]
    #     f = open("$results_path/gps/$gp_file")
    #     gp_set = deserialize(f)
    #     gp_x1 = gp_set["x1"]
    #     close(f)
    # end

    # Plot the maximum results
    maxPrs = res_mat[1:num_dfa_states:end,4]
    minPrs = res_mat[1:num_dfa_states:end,3]
    min_minPr = minimum(minPrs)
    min_maxPr = minimum(maxPrs)

    function plot_cell(extent, prob_value, norm_val)
        x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
        y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
        shape = Plots.Shape(x, y)
        plot!(shape, color=:black, fillalpha=1-prob_value, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
    end

    plt_max = plot(#title="$k Step Maximum Safety",
                    aspect_ratio=1,
                    size=(300,300), dpi=300,
                    xlims=[minx, maxx], ylims=[miny, maxy],
                    #xlabel="X Units", ylabel="Y Units",
                    xtickfont=font(10),
                    ytickfont=font(10),
                    titlefont=font(10),
                    grid=false)
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
    [plot_cell(extents[i], maxPrs[i], min_maxPr) for i in 1:num_regions]
    savefig(plt_max, "$plots_basename-max-heatmap.png")
    # Plot the minimum results
    plt_min = plot(#title=L"$k Step Minimum Safety",
                    aspect_ratio=1,
                    size=(300,300), dpi=300,
                    xlims=[minx, maxx], ylims=[miny, maxy],
                    #xlabel="X Units", ylabel="Y Units",
                    xtickfont=font(10),
                    ytickfont=font(10),
                    titlefont=font(10),
                    grid=false,
                    backgroundcolor=128)
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
    [plot_cell(extents[i], minPrs[i], min_minPr) for i in 1:num_regions]

    if !isnothing(extents_dict)
        for key in keys(extents_dict)
            for extent in extents_dict[key]
                shape = extent_dict_to_shape(extent)
                if key == "!a∧!b∧c"
                    plot!(shape, color=:black, label="", fillalpha=1.0)
                else
                    plot!(shape, color=:black, label="", fillalpha=0.0)
                end
            end
        end
    end

    savefig(plt_min, "$plots_basename-min-heatmap.png")

    plt_min = plot(#title=L"$k Step Minimum Safety",
                    aspect_ratio=1,
                    size=(300,300), dpi=300,
                    xlims=[minx, maxx], ylims=[miny, maxy],
                    #xlabel="X Units", ylabel="Y Units",
                    xtickfont=font(10),
                    ytickfont=font(10),
                    titlefont=font(10),
                    grid=false,
                    backgroundcolor=128)

    function plot_cell_verify(extent, min_prob_value, max_prob_value, threshold)
        x = [extent["x1"][1], extent["x1"][1], extent["x1"][2], extent["x1"][2]]
        y = [extent["x2"][1], extent["x2"][2], extent["x2"][2], extent["x2"][1]]
        shape = Plots.Shape(x, y)

        # if min_prob_value >= threshold
        plot!(shape, color=:black, fillalpha=min_prob_value/2, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
        # elseif max_prob_value < threshold
            # plot!(shape, color=:purple, fillalpha=0.99, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
        # else
            # plot!(shape, color=:purple, fillalpha=0.50, linealpha=0.05, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], label="")
        # end
    end
    # Plot the cells
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, linewidth=2, label="")
    [plot_cell_verify(extents[i], minPrs[i], maxPrs[i], min_threshold) for i in 1:num_regions]

    # savefig(plt_min, "$plots_basename-verification.png")

    if !isnothing(trajectories)
        for trajectory in trajectories
            if length(trajectory) > 1
                [plot!([trajectory[i][1], trajectory[i+1][1]], [trajectory[i][2], trajectory[i+1][2]], color=:black, label="", linewith=2) for i=1:length(trajectory)-1]
                scatter!([trajectory[end][1]], [trajectory[end][2]], color=:black, markershape=:star5, label="")
            end
        end

        savefig(plt_min, "$plots_basename-online-sim.png")
    end

    if !isnothing(extents_dict)
        for key in keys(extents_dict)
            for extent in extents_dict[key]
                shape = extent_dict_to_shape(extent)
                if key == "!a∧!b∧c"
                    plot!(shape, color=:black, label="", fillalpha=1.0)
                else
                    plot!(shape, color=:black, label="", fillalpha=0.0)
                end
            end
        end
    end

    if !isnothing(trajectory_dict)
        for (jj, key) in enumerate(keys(trajectory_dict))
            trajectories = trajectory_dict[key]
            for trajectory in trajectories
                if length(trajectory) > 1
                    [plot!([trajectory[i][1], trajectory[i+1][1]], [trajectory[i][2], trajectory[i+1][2]],  color=jj, label="", linewith=2) for i=1:length(trajectory)-1]
                    scatter!([trajectory[end][1]], [trajectory[end][2]], color=jj, markershape=:star5, label="")
                end
            end
            scatter!([trajectories[1][1][1]], [trajectories[1][1][2]], markershape=:circle, color=jj, label=key)
        end
        savefig(plt_min, "$plots_basename-online-sim-dict.png")
    end

    if !isnothing(online_data_dict)
        # Plot where the data was collected and a convex hull
        x_data = []
        for i=1:online_data_dict.count
            for dp in online_data_dict[i]["x_online"]
                push!(x_data, dp)
            end
        end
        # data_shape = VPolygon(convex_hull(x_data))
        # plot!(data_shape, fillalpha=0.25)
        [scatter!([x_data[i][1]], [x_data[i][2]], label="", color=:black) for i=1:length(x_data)]
        savefig(plt_min, "$plots_basename-online-data-hull.png")
    end
end

function plot_extent_outlines(extents, extents_dict, plots_basename)
    num_regions = extents.count - 1
    X = extents[extents.count]
    minx = minimum(X["x1"])
    maxx = maximum(X["x1"])
    miny = minimum(X["x2"])
    maxy = maximum(X["x2"])

    plt_min = plot(#title=L"$k Step Minimum Safety",
                    aspect_ratio=1,
                    size=(300,300), dpi=300,
                    xlims=[minx, maxx], ylims=[miny, maxy],
                    #xlabel="X Units", ylabel="Y Units",
                    xtickfont=font(10),
                    ytickfont=font(10),
                    titlefont=font(10),
                    grid=false,
                    backgroundcolor=128)
    plot!(Plots.Shape([minx, minx, maxx, maxx], [miny, maxy, maxy, miny]), fillalpha=0, linecolor=:black, foreground_color_border=:white, foreground_color_axis=:white, xticks = [minx, 0, maxx], yticks = [miny, 0, maxy], linewidth=2, label="")

    for key in keys(extents_dict)
        for extent in extents_dict[key]
            shape = extent_dict_to_shape(extent)
            if key == "!a∧!b∧c"
                plot!(shape, color=:black, label="", fillalpha=1.0)
            else
                plot!(shape, color=:black, label="", fillalpha=0.0)
            end
        end
    end
    savefig(plt_min, "$plots_basename-extent-dict.png")

end


text:
using Splines
using MathGL
function visualize3d{T}(song::AbstractArray{T, 1})
#Draw the sonogram
	sonog = sonogram(song, transformation=STFT.powSpec);
	sonogramMax = maximum(sonog)
	ops = plotOpStack()
	push!(ops, gr->Zoom(gr, -.3, 1.3, -.3, 1.3))
	push!(ops, gr->SetRanges(gr, 0, size(sonog, 2), 0, size(sonog, 1), minimum(sonog), maximum(sonog)))
	push!(ops, gr->SetRange(gr, 'c', minimum(sonog), maximum(sonog)))
	push!(ops, gr->Surf(gr, sonog), "Surface")
#	ops = Surf(sonog)

	insert!(ops, 2, gr->Alpha(gr, true))
#	insert!(ops, 2, gr->SetAlphaDef(gr, .09), "Set transparency value")
	insert!(ops, 2, gr->SetDefScheme(gr, "Wrcbk")) #Set the color scheme

#Draw just a color map
	push!(ops, gr->Dens(gr, sonog), "Color map")

#Draw the chirps.
	function drawChirpsFlat(gr::mglGraph)
		for i in 1:size(phr, 1)
			Plot(gr, [Float64(i) for i in phrb[i][1]:phrb[i][2]], phr[i], [sonogramMax for i in phrb[i][1]:phrb[i][2]], "b=")
		end
	end

	push!(ops, drawChirpsFlat, "Draw chirps")

#Draw the pitch as a function of time, with amplituide information
	sonogramTrace = centroidfrequencies(song, STFT.gaussian(windowlength), windowstep, lowharmonic, highharmonic, STFT.powSpec)
	env = envelope_along_centroid(song, STFT.gaussian(windowlength), windowstep, lowharmonic, highharmonic, STFT.powSpec, 4)
	push!(ops, gr->Area(gr, [i for i in 1:size(sonogramTrace, 1)], sonogramTrace, env, "y", "legend 'Pitch and volume vs time'"), "Pitch and volume vs. time")
	push!(ops, gr->Plot(gr, [i for i in 1:size(sonogramTrace, 1)], sonogramTrace, [sonogramMax for i in 1:size(sonogramTrace, 1)], "u:", "legend 'Pitch vs. time'"), "Pitch vs time")
	push!(ops, gr->Plot(gr, [i for i in 1:size(sonogramTrace, 1)], [0 for i in 1:size(sonogramTrace, 1)], env, "n", "legend Volume"), "Volume vs time")

#Get the phrases in the song
	phr = phrases(song, 10);
	phrb = phraseBoundaries(song, 10);
	if size(phrb, 1) == 0
		return ops
	end
	songStart = phrb[1][1]
	songEnd = phrb[end][2]
	push!(ops, gr->Mark(gr, [Float64(mean(i)) for i in phrb], [Float64(mean(i)) for i in phr], [sonogramMax for i in phrb], [Float64(size(i, 1))/100 for i in phr], "ro", "legend 'Chirp mean frequencies'"), "Mark chirp means")

boundPlaneHeight = maximum(sonog)/20 #The height of the planes drawn at the chirp boundaries

#Draw planes at the boundaries of the chirps
	function drawBoundsPlanes(gr::mglGraph)
		Alpha(gr, true)
		SetTranspType(gr, 1)
		for i in phrb
			FaceX(gr, mglPoint(i[1], 0, 0), 54, boundPlaneHeight, "H")
			FaceX(gr, mglPoint(i[2], 0, 0), 54, boundPlaneHeight, "H")
		end
	end

	push!(ops, drawBoundsPlanes, "Draw planes at chirp boundaries")

#Draw lines at the boundaries of the chirps
	lineBoundStyle = "r="
	function drawBoundsLines(gr::mglGraph)
		Alpha(gr, true)
		SetTranspType(gr, 1)
		for i in phrb
			Line(gr, mglPoint(i[1], 0, 0), mglPoint(i[1], 54, 0), lineBoundStyle)
			Line(gr, mglPoint(i[2], 0, 0), mglPoint(i[2], 54, 0), lineBoundStyle)
		end
	end

	push!(ops, drawBoundsLines, "Draw lines at chirp boundaries")
	push!(ops, gr->AddLegend(gr, "Chirp Boundaries", lineBoundStyle), "Set linestyle for lines at chirp boundaries")

#Draw smoothed pitch as a function of time
	chirpTimes = [Float64(mean(i)) for i in phrb]
	tone = [Float64(mean(i)) for i in phr] #Pitch as a function of time, only defined in the chirps.
	if length(chirpTimes) > 1 #Handle the case where chirpTimes is empty
		toneSpline = Spline(chirpTimes, tone)
		toneStart = mean(phrb[1])
		toneEnd = mean(phrb[end])
		toneSmooth = [Float64(toneSpline(Float64(i))) for i in toneStart:10:toneEnd]
		push!(ops, gr->Plot(gr, [ i for i in toneStart:(toneEnd-toneStart)/(size(toneSmooth, 1)-1):toneEnd], toneSmooth, [sonogramMax for i in toneSmooth], "g2", "legend 'Smoothed pitch vs time'"), "Draw smooth tone")
	else
		warn("Empty chirpTimes; will not plot smoothed pitch")
	end

#Legend
	push!(ops, gr->Legend(gr, 3, "Wbk#"), "Legend")

#Axes, colorbar, and labels
	push!(ops, gr->Label(gr, 'x', "Time (seconds)"))
	push!(ops, gr->Label(gr, 'y', "Frequency (Hz)"))
	push!(ops, gr->Label(gr, 'z', "Amplitude"))

	push!(ops, gr->SetRanges(gr, 0, size(sonog, 2)*windowstep/samplerate, (samplerate/windowlength)*lowharmonic, (samplerate/windowlength)*highharmonic), )
	push!(ops, gr->MathGL.Box(gr))
	push!(ops, gr->Axis(gr))
	push!(ops, gr->Colorbar(gr));

	return ops
end

export visualize3d

ext: jl


text:
"""
```
krylov_reduction(Γ0, Γ1, Ψ, Π, C; F = x -> x, n_Z = 0)
valuef_reduction(grids, params, n_g_red, Γ0, Γ1, Ψ, Π, C, n_v)
```

Perform Krylov reduction on the state space system defined by Γ0 and Γ1
and/or value function reduction using a spline approximation across some
subset of the state space.

Returns the new matrices, the change of basis matrix and its inverse,
and the new size of the state space.

"""
function krylov_reduction(m::AbstractCTModel, Γ0::SparseMatrixCSC{Float64,Int64},
                          Γ1::SparseMatrixCSC{Float64,Int64},
                          Ψ::SparseMatrixCSC{Float64,Int64}, Π::SparseMatrixCSC{Float64,Int64},
                          C::SparseMatrixCSC{Float64,Int64})

    # Get Krylov subspace function
    if !in(:F, keys(m.settings))
        m <= Setting(:F, identity, "Function applied during Krylov reduction")
    end
    F = get_setting(m, :F)

    # Grab Dimensions
    n_state_vars_unreduce = Int64(get_setting(m, :n_state_vars_unreduce))
    n_state_vars          = Int64(get_setting(m, :n_state_vars))
    n_jump_vars           = Int64(get_setting(m, :n_jump_vars))
    krylov_dim            = Int64(get_setting(m, :krylov_dim))
    n_total               = n_jump_vars + n_state_vars
    n_vars                = n_states(m)
    n_state_vars         -= n_state_vars_unreduce

    # Slice Γ1 into different parts. See Why Inequality matters paper.
    Γ1_full = Array{Float64}(Γ1)
    B_pv    = -Γ1_full[n_total .+ 1:n_vars, n_total + 1:n_vars] \ Γ1_full[n_total .+ 1:n_vars, 1:n_jump_vars]
    B_pg    = -Γ1_full[n_total .+ 1:n_vars, n_total .+ 1:n_vars] \
           Γ1_full[n_total .+ 1:n_vars, n_jump_vars .+ 1:n_jump_vars + n_state_vars]
    B_pZ    = -Γ1_full[n_total .+ 1:n_vars, n_total .+ 1:n_vars] \
        Γ1_full[n_total+1:n_vars, n_jump_vars+n_state_vars .+ 1:n_jump_vars+n_state_vars+n_state_vars_unreduce]
    B_gg    = Γ1_full[n_jump_vars .+ 1:n_jump_vars + n_state_vars, n_jump_vars .+ 1:n_jump_vars + n_state_vars]
    B_gv    = Γ1_full[n_jump_vars .+ 1:n_jump_vars + n_state_vars, 1:n_jump_vars]
    B_gp    = Γ1_full[n_jump_vars .+ 1:n_jump_vars .+ n_state_vars, n_total .+ 1:n_vars]
    # Drop redundant equations
    obs        = B_pg
    ~, d0, V_g = svd(obs)
    aux        = d0/d0[1]
    n_Bpg      = Int64(sum(aux .> 10*eps()))
    V_g        = V_g[:, 1:n_Bpg] .* aux[1:n_Bpg]'
    # Compute Krylov subspace
    A(x::Matrix{Float64}) = (F == identity ? 0. : F(B_gv' * x)) .+ B_gg' * x + B_pg' * (B_gp' * x)
    V_g, ~, ~             = deflated_block_arnoldi(A, V_g, krylov_dim)
    n_state_vars_red      = size(V_g, 2) # number of state variables after reduction
    # Build state space reduction transform
    reduced_basis = spzeros(Float64, n_jump_vars+n_state_vars_red,n_vars)
    reduced_basis[1:n_jump_vars,1:n_jump_vars] = SparseMatrixCSC{Float64}(I, n_jump_vars, n_jump_vars)
    reduced_basis[n_jump_vars .+ 1:n_jump_vars+n_state_vars_red,n_jump_vars .+ 1:n_jump_vars+n_state_vars] = V_g'
    reduced_basis[n_jump_vars + n_state_vars_red .+ 1:n_jump_vars + n_state_vars_red + n_state_vars_unreduce,
                  n_jump_vars + n_state_vars .+ 1:n_jump_vars + n_state_vars + n_state_vars_unreduce] =
                      Matrix{Float64}(I, n_state_vars_unreduce,
                                      n_state_vars_unreduce)

    # Build inverse transform
    inv_reduced_basis                              = spzeros(Float64, n_vars, n_jump_vars + n_state_vars_red)
    inv_reduced_basis[1:n_jump_vars,1:n_jump_vars] = SparseMatrixCSC(I, n_jump_vars, n_jump_vars)
    inv_reduced_basis[n_jump_vars+1:n_jump_vars+n_state_vars,n_jump_vars+1:n_state_vars_red+n_jump_vars] = V_g
    inv_reduced_basis[n_total+1:n_vars,1:n_jump_vars] = B_pv
    inv_reduced_basis[n_total+1:n_vars,n_jump_vars+1:n_jump_vars .+ n_state_vars_red] = B_pg * V_g
    inv_reduced_basis[n_total+1:n_vars,n_jump_vars+n_state_vars_red .+ 1:n_jump_vars+n_state_vars_red+n_state_vars_unreduce] = B_pZ
    inv_reduced_basis[n_jump_vars+n_state_vars .+ 1:n_total,n_jump_vars+n_state_vars_red .+ 1:n_jump_vars+n_state_vars_red+n_state_vars_unreduce] = SparseMatrixCSC{Float64}(I, n_state_vars_unreduce, n_state_vars_unreduce)
    m <= Setting(:n_state_vars_red, n_state_vars_red + n_state_vars_unreduce,
                 "Number of state variables after reduction")

    # change basis
    Γ0_kry, Γ1_kry, Ψ_kry, Π_kry, C_kry =
        change_basis(reduced_basis, inv_reduced_basis, Γ0, Γ1, Ψ, Π, C; ignore_Γ0 = true)

    return Γ0_kry, Γ1_kry, Ψ_kry, Π_kry, C_kry, reduced_basis, inv_reduced_basis
end

function valuef_reduction(m::AbstractCTModel,
                          Γ0::SparseMatrixCSC{Float64,Int64}, Γ1::SparseMatrixCSC{Float64,Int64},
                          Ψ::SparseMatrixCSC{Float64,Int64},
                          Π::SparseMatrixCSC{Float64,Int64}, C::SparseMatrixCSC{Float64,Int64})

    n_jump_vars      = Int64(get_setting(m, :n_jump_vars))
    n_state_vars_red = Int64(get_setting(m, :n_state_vars_red))
    spline_grid      = get_setting(m, :spline_grid)
    knots_dict       = get_setting(m, :knots_dict)
    n_prior          = Int64(get_setting(m, :n_prior))
    n_post           = Int64(get_setting(m, :n_post))

    # Function calls to create basis reduction
    knots_dim       = sum(map(i -> length(knots_dict[i]) + 1, 1:length(keys(knots_dict))))
    spline_grid_dim = prod(size(spline_grid))
    from_spline     = spzeros(Float64, spline_grid_dim, knots_dim)
    to_spline       = spzeros(Float64, knots_dim, spline_grid_dim)

    # create spline basis
    from_spline[:,:], to_spline[:,:] = spline_basis(spline_grid, knots_dict, 2)
    # extend to other dimensions along which we did not use a spline approximation.
    # In this case, we computed a spline basis for wealth but not for income
    from_spline, to_spline = extend_to_nd(from_spline, to_spline, n_prior, n_post)
    # extra jump variables besides value function
    extra_jv = Int64(n_jump_vars - length(m.endogenous_states[:value_function]))

    if extra_jv > 0
        # Add additional dimensions for any jump variables we did not reduce b/c
        # spline reduction only for value function
        dim1_from, dim2_from = size(from_spline)
        from_spline = sparse([from_spline spzeros(dim1_from, extra_jv)]) # add extra columns
        from_spline = sparse([from_spline; spzeros(extra_jv, dim2_from + extra_jv)]) # add extra zeros
        to_spline = sparse([to_spline spzeros(dim2_from, extra_jv)])     # same as above
        to_spline = sparse([to_spline; spzeros(1, dim1_from + extra_jv)])
        for i = 0:(extra_jv - 1)
            from_spline[end - i, end - i] = 1 # add 1 along diagonal of added dimensions
            to_spline[end - i, end - i] = 1
        end
    end
    n_splined = Int64(size(from_spline, 2))

    # Create projection matrix that projects value function onto spline basis
    from_spline, to_spline = projection_for_subset(from_spline, to_spline, 0, n_state_vars_red)

    # Reduce the decision vector
    Γ0_spl, Γ1_spl, Ψ_spl, ~, C_spl = change_basis(to_spline, from_spline, Γ0, Γ1, Ψ, Π, C; ignore_Γ0 = true)

    Π_spl = to_spline * Π * from_spline[1:n_jump_vars, 1:n_splined]
    m  <= Setting(:n_splined, n_splined, "Dimension of jump variables after spline basis reduction")

    return Γ0_spl, Γ1_spl, Ψ_spl, Π_spl, C_spl, to_spline, from_spline
 end

# Just think through how matrices can be used to do a change of basis
function change_basis(basis, inv_basis, Γ0::SparseMatrixCSC{Float64,Int64},
                      Γ1::SparseMatrixCSC{Float64,Int64},
                      Ψ::SparseMatrixCSC{Float64,Int64}, Π::SparseMatrixCSC{Float64,Int64},
                      C::SparseMatrixCSC{Float64,Int64};
                      ignore_Γ0::Bool = false)

    g1 = basis * Γ1 * inv_basis

    if ignore_Γ0
        g0 = SparseMatrixCSC{Float64}(I, size(g1))
    else
        g0 = basis * Γ0 * inv_basis
    end
    c    = basis * C
    psi  = basis * Ψ
    Pi   = basis * Π

    return g0, g1, psi, Pi, c
end

# Cleans Γ0 matrix to be identity by solving out static conditions
function solve_static_conditions(Γ0::SparseMatrixCSC{Float64,Int64},
                                 Γ1::SparseMatrixCSC{Float64,Int64},
                                 Ψ::SparseMatrixCSC{Float64,Int64},
                                 Π::SparseMatrixCSC{Float64,Int64}, C::Array{Float64})

    redundant     = maximum(abs.([Γ0, Ψ]), 2) .== 0  # Find rows of both Γ0 & Ψ with only zeros
    inv_state_red = null(Γ1(redundant,:)) # Compute orthonormal basis of null space of Γ1
    state_red     = inv_state_red'

    g0  = state_red * Γ0 * inv_state_red # adjust Γ0 to new basis
    g1  = state_red * Γ1 * inv_state_red # Zero out redundant rows/columns in Γ1
    g1  = g0 \ g1                        # Force Γ1 and other matrices
    Psi = g0 \ (state_red * Ψ)           # to be the ones corresponding to when
    Pi  = g0 \ (state_red * Π)           # Γ0 is the identity
    c   = g0 \ (state_red * C)

    return g0, g1, c, Pi, Psi, state_red, inv_state_red
end

# When we have sparse matrices
#=
function solve_static_conditions(Γ0::SparseMatrixCSC{Float64, Int64}, Γ1::SparseMatrixCSC{Float64, Int64},
                  Ψ::SparseMatrixCSC{Float64, Int64},
                  Π::SparseMatrixCSC{Float64, Int64}, C::SparseMatrixCSC{Float64, Int64})

    redundant     = maximum(abs.([Γ0, Ψ]), 2) .== 0
    inv_state_red = null(full(Γ1(redundant,:)))
    state_red     = inv_state_red'

    g0  = state_red * full(Γ0) * inv_state_red
    g1  = state_red * full(Γ1) * inv_state_red
    g1  = g0 \ g1
    Psi = g0 \ (state_red * full(Ψ))
    Pi  = g0 \ (state_red * full(Π))
    c   = g0 \ (state_red * full(C))

    return g0, g1, c, Pi, Psi, state_red, inv_state_red
end
=#


text:
# %%-----------Koch Curve---------------------

function pointskoch(points, maxk, α = sqrt(3)/2)
  Q = [0 -1; 1 0]
  for k = 1:maxk
    n = length(points)
    new_points = Vector{Float64}[]
    for i = 1:n-1
      p1, p2 = points[i], points[i+1]
      v = (p2 - p1) / 3
      q1 = p1 + v
      q2 = p1 + 1.5v + α * Q * v
      q3 = q1 + v
      append!(new_points, [p1, q1, q2, q3])
    end
    push!(new_points, points[end])
    points = new_points
  end
  return points
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using LinearAlgebra: LinearAlgebra, Symmetric, SymTridiagonal, I

using NQCBase: NQCBase, PeriodicCell, InfiniteCell

struct RingPolymerParameters{T<:AbstractFloat}
    n_beads::Int
    ω_n::T
    ω_n²::T
    springs::Symmetric{T,Matrix{T}}
    normal_mode_springs::Vector{T}
    U::Matrix{T}
    quantum_atoms::Vector{Int}
    tmp::Vector{T}
    function RingPolymerParameters{T}(
        n_beads::Integer, temperature::Real,
        quantum_atoms::Vector{<:Integer}) where {T<:AbstractFloat}

        ω_n = n_beads * temperature
        new(n_beads, ω_n, ω_n^2,
            get_spring_matrix(n_beads, ω_n),
            get_normal_mode_springs(n_beads, ω_n),
            get_normal_mode_transformation(n_beads),
            quantum_atoms,
            zeros(T, n_beads))
    end
end

Base.broadcastable(beads::RingPolymerParameters) = Ref(beads)

"""Constructor for choosing specific elements to be quantum."""
function RingPolymerParameters{T}(n_beads::Integer, temperature::Real, atom_types::AbstractVector{Symbol}, quantum_nuclei::Vector{Symbol}) where {T}
    quantum_atoms = findall(in(quantum_nuclei), atom_types)
    RingPolymerParameters{T}(n_beads, temperature, quantum_atoms)
end

"""Constructor for the case where all nuclei are quantum."""
function RingPolymerParameters{T}(n_beads::Integer, temperature::Real, n_atoms::Integer) where {T}
    RingPolymerParameters{T}(n_beads, temperature, collect(1:n_atoms))
end

"""
    get_L(n_beads, mass, ω_n)

Return the Circulant symmetric matrix for the ring polymer springs.
"""
function get_spring_matrix(n_beads::Integer, ω_n::Real)::Symmetric
    if n_beads == 1
        spring = zeros(1, 1)
    elseif n_beads == 2
        spring = [2 -2; -2 2]
    else
        spring = SymTridiagonal(fill(2, n_beads), fill(-1, n_beads-1))
        spring = convert(Matrix, spring)
        spring[end,1] = spring[1, end] = -1
    end
    Symmetric(spring .*  ω_n^2 / 2)
end

get_normal_mode_springs(n_beads::Integer, ω_n::Real) = get_matsubara_frequencies(n_beads, ω_n) .^2 / 2
get_matsubara_frequencies(n::Integer, ω_n::Real) = 2ω_n*sin.((0:n-1)*π/n)

# """
#     get_normal_mode_transformation(n::Int)::Matrix

# Get the transformation matrix that converts to normal mode coordinates.
# """
# function get_normal_mode_transformation(n::Int)::Matrix
#     a = ([exp(2π*im/n * j * k) for j=0:n-1, k=0:n-1])
#     (real(a) + imag(a)) / sqrt(n)
# end

"""
Creates normal mode transformation for `n` beads.
"""
function get_normal_mode_transformation(n::Integer)::Matrix
    # Real and imaginary parts of the discrete Fourier transform matrix.
    U = sqrt(2/n) .* hcat([cos(2π * j * k / n) for j=0:n-1, k=0:n÷2],
                        [sin(2π * j * k / n) for j=0:n-1, k=n÷2+1:n-1])

    # Normalisation
    U[:, 1] ./= sqrt(2)
    iseven(n) && (U[:, n÷2+1] ./= sqrt(2))

    U
end

function transform_to_normal_modes!(p::RingPolymerParameters, R::AbstractArray{T,3}) where {T}
    @views for i in p.quantum_atoms
        for j in axes(R, 1)
            LinearAlgebra.mul!(p.tmp, p.U', R[j,i,:])
            R[j,i,:] .= p.tmp
        end
    end
end

function transform_from_normal_modes!(p::RingPolymerParameters, R::AbstractArray{T,3}) where {T}
    @views for i in p.quantum_atoms
        for j in axes(R, 1)
            LinearAlgebra.mul!(p.tmp, p.U, R[j,i,:])
            R[j,i,:] .= p.tmp
        end
    end
end

function transform!(p::RingPolymerParameters, A::RingPolymerArray)
    if A.normal
        transform_from_normal_modes!(p, A)
    else
        transform_to_normal_modes!(p, A)
    end
    A.normal = !A.normal
end

nbeads(beads::RingPolymerParameters) = beads.n_beads
Base.length(beads::RingPolymerParameters) = nbeads(beads)
Base.range(beads::RingPolymerParameters) = range(1; length=length(beads))

"""
    cayley_propagator(beads::RingPolymerParameters{T}, dt::Real; half::Bool=true) where {T}

J. Chem. Phys. 151, 124103 (2019); doi: 10.1063/1.5120282
"""
function cayley_propagator(beads::RingPolymerParameters{T}, dt::Real; half::Bool=true) where {T}

    cay(dtA::Matrix)::Matrix = LinearAlgebra.inv(I - dtA/2) * (I + dtA/2)

    ω_k = get_matsubara_frequencies(length(beads), beads.ω_n)
    prop = [Array{T}(undef, 2, 2) for i=1:length(beads)]
    for (i, ω) in enumerate(ω_k)
        A = [0 1; -ω^2 0]
        prop[i] .= half ? real.(sqrt(cay(dt.*A))) : cay(dt.*A)
    end
    prop
end

"""
    get_spring_energy(beads::RingPolymerParameters, masses, R)

Calculate the ring polymer spring potential.
"""
function get_spring_energy(beads::RingPolymerParameters, masses, R)
    E = zero(eltype(R))
    for bead=1:nbeads(beads)-1
        for i in beads.quantum_atoms # Only for quantum nuclei
            for j in axes(R, 1)
                E += masses[i] * (R[j, i, bead] - R[j, i, bead+1])^2
            end
        end
    end
    for i in beads.quantum_atoms
        for j in axes(R, 1)
            E += masses[i] * (R[j, i, end] - R[j, i, 1])^2
        end
    end
    E * beads.ω_n^2/2
end

function NQCBase.apply_cell_boundaries!(cell::PeriodicCell, R::AbstractArray{T,3}, beads::RingPolymerParameters) where {T}
    transform_to_normal_modes!(beads, R)
    R[:,beads.quantum_atoms,1] ./= sqrt(length(beads))
    @views NQCBase.apply_cell_boundaries!(cell, R[:,:,1])
    R[:,beads.quantum_atoms,1] .*= sqrt(length(beads))
    transform_from_normal_modes!(beads, R)
end
NQCBase.apply_cell_boundaries!(::InfiniteCell, ::AbstractArray{T,3}, ::RingPolymerParameters) where {T} = nothing


text:
# example 1: basic string hashing

include("../src/Murmur3.jl")
using Murmur3

value = "Hello world"

hash = Murmur3.x86.hash32(value)
println("x86  32-bit: \"$value\" => $(hex(hash))")

hash = Murmur3.x86.hash128(value)
println("x86 128-bit: \"$value\" => $(hex(hash))")

hash = Murmur3.x64.hash128(value)
println("x64 128-bit: \"$value\" => $(hex(hash))")


text:
using DataFrames
using PrettyTables
import ModelParameters: Param
import MetaGraphs

abstract type NetworkNode end

@def network_node begin
    name::String
    area::Float64
end

Base.@kwdef mutable struct GenericNode <: NetworkNode
    @network_node
end

"""Create an arbitrary node."""
function create_node(node::Type{<:NetworkNode}, name::String, area::Float64)
    return node{Param, Float64}(; name=name, area=area)
end

function GenericNode(name::String, spec::Dict)
    mod_spec = copy(spec)
    delete!(mod_spec, "node_type")
    delete!(mod_spec, "name")
    delete!(mod_spec, "area")
    delete!(mod_spec, "func")
    n = GenericNode(; name=name, area=spec["area"], mod_spec...)

    return n
end

Base.@kwdef mutable struct GenericDirectNode <: NetworkNode
    @network_node

    outflow = []
end

function GenericDirectNode(name::String, spec::Dict)
    mod_spec = copy(spec)
    area = spec["area"]
    delete!(mod_spec, "node_type")
    delete!(mod_spec, "name")
    delete!(mod_spec, "area")
    delete!(mod_spec, "func")
    n = GenericDirectNode(; name=name, area=area, mod_spec...)

    return n
end

"""
    param_info(node::NetworkNode)

Generic parameter information extractor.

Extracts parameter names, values, and bounds
"""
function param_info(node::NetworkNode; kwargs...)::Tuple
    tmp = Model(node)
    values = collect(tmp[:val])
    bounds = collect(tmp[:bounds])
    param_names = collect(tmp[:fieldname])

    return param_names, values, bounds
end

"""Retrieve network node_id for a given gauge (by name)."""
function get_node_id(mg::MetaDiGraph, node_name::String)::Int64
    v = collect(MetaGraphs.filter_vertices(mg, :name, node_name))

    @assert length(v) > 0 || error("Node with name '$(node_name)' not found")
    @assert length(v) == 1 || error("Found more than 1 node with name '$(node_name)'")
    return v[1]
end

"""
    get_node(sn::StreamfallNetwork, node_name::String)

Retrieve node_id and node property for a specified gauge.

# Arguments
- sn : Streamfall Network
- node_name : name of node of interest

# Returns
Tuple, node position id and node object
"""
get_node(sn::StreamfallNetwork, node_name::String) = get_node(sn.mg, node_name)
function get_node(mg, node_name::String)::Tuple
    v_id = get_node_id(mg, node_name)
    return v_id, MetaGraphs.get_prop(mg, v_id, :node)
end

get_node(sn::StreamfallNetwork, v_id::Int)::NetworkNode = MetaGraphs.get_prop(sn.mg, v_id, :node)
get_node_name(sn::StreamfallNetwork, v_id::Int) = MetaGraphs.get_prop(sn.mg, v_id, :name)

get_node(mg, v_id::Int) = MetaGraphs.get_prop(mg, v_id, :node)
get_node_name(mg, v_id::Int) = MetaGraphs.get_prop(mg, v_id, :name)

"""
    subcatchment_data(node::NetworkNode, data::DataFrame)

Extract all data for a given node from a dataframe.
"""
function subcatchment_data(node::NetworkNode, data::DataFrame, partial::String=nothing)::DataFrame
    name = node.name
    cols = filter(x -> occursin(name, string(x)), names(data))

    if !isnothing(partial)
        cols = filter(x -> occursin(name, x), cols)
    end

    return data[:, cols]
end

function extract_node_spec(node::NetworkNode)
    area = node.area
    param_names, x0, _ = param_info(node)
    params = Dict(zip(param_names, x0))

    if hasproperty(node, :storage)
        params[:initial_storage] = node.storage[1]
    end

    node_type = typeof(node)
    spec = Dict(
        "node_type" => string(nameof(node_type)),
        "area" => area,
        "parameters" => params
    )

    return spec
end

Base.show(io::IO, ::MIME"text/plain", n::NetworkNode) = show(io, n)
function Base.show(io::IO, n::NetworkNode)
    ntype_name = nameof(typeof(n))
    println(io, "Name: $(n.name) [$(ntype_name)]")
    println(io, "Area: $(n.area)")

    param_names, x0, bounds = param_info(n; with_level=false)
    lb, ub = zip(bounds...)
    details = hcat(param_names, x0, [lb...], [ub...])

    pretty_table(io, details, ["Parameter", "Value", "Lower Bound", "Upper Bound"])
    print("\n")
end


text:
module GenomePermutations

import BED
import Distributions
import GenomicFeatures
import Statistics
import HypothesisTests

include("Overlaps.jl")
include("Dist.jl")

export anyoverlapping, alloverlapping, isin, countoverlapping
export dist
export generatedistribution
export randomiseregions, randomisegenome
export simpleP, SimplePTest, pvalue
export permtest, overlappermtest, PermTestResult

"""
	vec_getcollection(collection::GenomicFeatures.IntervalCollection{T}, sequence::String)

(inefficiently) retun a interval collection with only features in the specified sequence

See Also [`GenomePermutations.iter_getcollection`](@ref)
"""
function vec_getcollection(collection, sequence::String)
	GenomicFeatures.IntervalCollection([GenomicFeatures.Interval(i.seqname,
	i.first, i.last, i.strand, i.metadata ) for i in collection.trees[sequence]])
end

"""
	iter_getcollection(collection::GenomicFeatures.IntervalCollection{T}, sequence::String)

(inefficiently) retun a interval collection with only features in the specified sequence.

see also [`GenomePermutations.vec_getcollection`](@ref)
"""
function iter_getcollection(
	collection::GenomicFeatures.IntervalCollection{T},
	sequence::String
) where {T}
	tmp = GenomicFeatures.IntervalCollection{T}()
	for i in collection.trees[sequence]
		push!(tmp, i)
	end
	return(tmp)
end

"""
	generatedistribution(collection)

Generate a Distributions.MixtureModel from the collection

The model assumes each position of the genome has equal probability of being sampled.
Practically, it is modelled as a mixture model of individual discrete uniform distributions,
one for each interval in the collection.
The component distributions are distribuited according to a categorical distribution, with
probabilities proportional to the length of each interval.

# Note

At the moment it only supports collections with single sequence.
Generating a distribution from an overlapping collection is undocumented
(I assume it would add up the probability but this is untested)
"""
function generatedistribution(collection)

	d = Vector{Distributions.DiscreteUniform}(undef, 0)
	l = Vector{Int}(undef,0)
	for i in collection
		push!(d, Distributions.DiscreteUniform(i.first, i.last))
		push!(l, i.last-i.first)
	end
	l = l/sum(l)
	return Distributions.MixtureModel(d, Distributions.Categorical(l))
end

"""
	_randominterval(interval::GenomicFeatures.Interval{T},
		distribution::Distributions.Sampleable, regions::GenomicFeatures.IntervalCollection{S};
		collection::GenomicFeatures.IntervalCollection{T} = GenomicFeatures.IntervalCollection{T}(),
		allow_overlap::Bool = true, max_tries::Int = 1000) where {T, S}}

This is the private version of randominterval, it skips sequence assertions in the name of performance
Does it actually improve performance? Debeatable, in theory it prevents 100s of IFs. In practice,
I haven't tested it. Indeed, branch prediction on asymmetrical branches should mean
that the ifs are not really that damaging to performance.
"""
function _randominterval(
	interval::GenomicFeatures.Interval{T},
	distribution::Distributions.Sampleable,
	regions::GenomicFeatures.IntervalCollection{S};
	collection::GenomicFeatures.IntervalCollection{T} = GenomicFeatures.IntervalCollection{T}(),
	allow_overlap::Bool = true,
	max_tries::Int = 1000,
	onfail = :throw
) where {T,S}

	il = interval.last - interval.first
	new = GenomicFeatures.Interval{T}

	# Can you turn this into branchless code for performance?
	for i in 1:max_tries
		newfirst = rand(distribution,1)
		new = GenomicFeatures.Interval(interval.seqname, newfirst[1], newfirst[1]+il,
				interval.strand, interval.metadata)
		if allow_overlap
				isin(new, regions) && return new
		else
				isin(new,regions) && !(anyoverlapping(new, collection)) && return new
		end
	end
	# allow passsing of the orignal interval if a random interval cannot be generated.
	if onfail == :orig
		if !(allow_overlap)
			@warn "Returning the original interval can currently cause overlapping segements - this is a known, currently untested issue!"
		end
		# if we allow overlapping, we can just return the original interval
		return interval
	elseif onfail == :skip
		return nothing
	else
		error("""
			Error: can't randomise regions with $max_tries tries,
			Please consider increasing max_iter,
			Please also check the intervals come from the same sequnce used to generate
			the distribution from the regions
		""")
	end
end

"""
	randominterval(interval, distribution, regions, collection; allow_overlap, max_tries, onfail)

Randomises the position of an interval according to a given distribution.

# Arguments

- `interval::GenomicFeatures.IntervalCollection`: The interval to randomise
- `distribution::Distributions.Sampleable`: distribution of possible interval start location.
Should be a discrete distribution, although right now this is not enforced at the signature
level
- `regions::`: the regions the new interval should be fully contained in (checked by`isin`)
- `collecton::GenomicFeatures.IntervalCollection =  GenomicFeatures.IntervalCollection{T}()`: Collection to avoid overlapping
to if allow_overlap is set to false, defaults to an empty collection of metadata .
- `allow_overlap::Bool = true`: whether to allow overlaps to collection.
- `max_tries::Int = 1000`: the maximum number of attempts to draw an interval that fits in
the region from the distribution
- `onfail::Symbol = :throw`: what to do if the interval fails to fit in the region.
Will return the orignal interval if :orig, returns nothing if :skip or throw an error if :throw.

# Note

To prevent intervals in a loop from overlapping each other, use
`allow_overlap = false` and add the result of each iteration to the collection
paseed to the function.

See Also (_randominterval)[@ref]
"""
function randominterval(
	interval::GenomicFeatures.Interval{T},
	distribution::Distributions.Sampleable,
	regions::GenomicFeatures.IntervalCollection{S};
	collection::GenomicFeatures.IntervalCollection{T} = GenomicFeatures.IntervalCollection{T}(),
	allow_overlap::Bool = true,
	max_tries::Int = 1000,
	onfail = :throw
) where {S,T}

	# check interval key is in regions
	interval.seqname in keys(regions.trees) || throw(KeyError("$interval.seqname not found in $regions"))
	# check distribution is valid, given the regions
	if !(distribution == generatedistribution(iter_getcollection(regions, interval.seqname)))
		throw(KeyError("$distribution does not match the regions for $interval.seqname"))
	end

	return _randominterval(
		interval, distribution, regions;
		collection = collection,
		allow_overlap = allow_overlap,
		max_tries = max_tries,
		onfail = onfail
	)
end

"""
	randomiseregions(collection, distribution, regions; allow_overlap, max_tries, onfail)

Generates randomised regions from an interval collection according to a start postion distribution.
It guarantees that the new intervals are contained in the given regions.

# Arguments

- `collection::GenomicFeatures.IntervalCollection{T}`: collection to randomise,
must have only 1 sequence.
- `distribution::Distributions.Sampleable`: distribution of possible interval
start locations, must match collection.
- `regions::GenomicFeatures.IntervalCollection{S}`: regions to randomise the intervals in.
- `allow_overlap::Bool = false.
- `max_tries::Int = 1000`: maximum number of attempts to draw an interval
that fits in the region from the distribution.
- `onfail::Symbol = :throw`: what to do if an interval fails to fit in the region.
Will return the orignal interval if :orig, skip the interval if :skip, otherwise throws an error.

# Notes

this functions assumes that the interval collection contains intervals belonging
to a single sequence that is the same sequence used for the distribution.
The distribution should be derived from the regions.
If the distribution and the regions don't match it will most likely fail safely
New intervals overlapping with other intervals in the nascient collection are disallowed by default.
To allow overlaps use `allow_overlap = true
"""
function randomiseregions(
	collection::GenomicFeatures.IntervalCollection{T},
	distribution::Distributions.Sampleable,
	regions::GenomicFeatures.IntervalCollection{S};
	allow_overlap = false, max_tries::Int = 1000, onfail = :throw
) where {T, S}

	new_collection = GenomicFeatures.IntervalCollection{T}()
	for i in collection
		# we no longer handle errors here. If onfail is set to :skip, we skip
		# before adding nothing to the collection
		new_interval =  GenomePermutations._randominterval(i, distribution, regions;
				collection = new_collection, allow_overlap, max_tries, onfail = onfail)
		# Should be a quite asymmetrical branch so not too much performance drag
		isnothing(new_interval) || push!(new_collection, new_interval)
	end
	return new_collection
end

"""
	randomisegenome(col, regions; allow_overlap, max_tries)

Randomises a collection, sequence by sequence along the specified regions, returning the new collection

# rguments

- `col::GenomicFeatrues.IntervalCollection{T}`: collection to randomise
- `regions::GenomicFeatures.IntervalCollection{S}`: regions to randomise the collection in.
- `allow_overlap::Bool = false`: whether to allow overlapping intervals in the result.
- `max_tries::Int = 1000`: maximum number of attempts to draw an interval
that fits in the regions.
- `onfail::Symbol = :throw`: what to do if an interval fails to fit in the regions.
Will return the orignal interval if :orig, skip the interval if :skip, otherwise throws an error.

Note: it will skip any sequences in the collection not found in the target regions.
If no sequences are found returns an empty collection.
"""
function randomisegenome(
	col::GenomicFeatures.IntervalCollection{T},
	regions::GenomicFeatures.IntervalCollection{S};
	allow_overlap = false, max_tries::Int = 1000,
	onfail=:throw
) where {S,T}

	# iterate per chromosome, generating all distributions
	ran = GenomicFeatures.IntervalCollection{T}()

	D = Dict{String, Distributions.MixtureModel}()
	for seq in keys(regions.trees)
		push!(D, seq => generatedistribution(regions.trees[seq]))
	end

	a = Dict{String,
	GenomicFeatures.IntervalCollection{T}}()
	for seq in keys(col.trees)
			push!(a, seq => GenomePermutations.iter_getcollection(col, seq))
	end

	r = Dict{String,
	GenomicFeatures.IntervalCollection{S}}()
	for seq in keys(regions.trees)
		push!(r, seq => GenomePermutations.iter_getcollection(regions, seq))
	end

	for seq in keys(col.trees)

		if !(seq in keys(regions.trees))
			@warn "sequence $seq is not found in the target regions, skipping it..."
			continue
		else
			# This is ugly but it works
			for interval in randomiseregions(a[seq],D[seq],r[seq]; allow_overlap, max_tries, onfail)
				push!(ran, interval)
			end
		end
	end
	return ran::GenomicFeatures.IntervalCollection{T}
end

"""
	simpleP(obs, ran, iterations; alternative :Auto)

Calculates the p value for the alternative hypothesis that the observed value
is either more or less than what we expect from the Vecor of random values.

If no altrantive is provided, will be automatically selected
on the basis of the data.

See also: [regioneR](https://doi.org/doi:10.18129/B9.bioc.regioneR), [`SimplePTest`](@ref)
Please see regioneR for full details
"""
function simpleP(obs, ran, iterations; alternative = :Auto)

	if alternative == :Less
		p = (length(filter(x-> x <= obs, ran)) +
		1)/(length(ran) + 1)
	elseif alternative == :More
		p = ((length(filter(x-> x >= obs, ran)) +
		1)/(length(ran) + 1))
	else
		# the median should be the same as mean for large numbers but is more robust for small numbers
		obs <= Statistics.median(ran) && return(simpleP(obs, ran, iterations; alternative = :Less))
		return(simpleP(obs, ran, iterations, alternative = :More))
	end
	return(SimplePTest(iterations, p, alternative))
end

"""
	SimplePTest(iterations::Int, p_val::, alternative::)

A structure to save RegioneR - style simple P tests. Implements `pvalue`
"""
struct SimplePTest

	iterations::Int
	p_val::Float64
	alternative::Symbol
	min_p_val::Float64
	function SimplePTest(iterations, p_val, alternative)
		new(iterations, p_val, alternative, (1/(iterations+1)))
	end
end

"""
	PermTestResult{S <: Union{Real, Vector{<:Real}}, N <: Union{Real, Vector{<:Real}}, T <: Any}

A structure to store results of any permutation test
implements pvalue.
"""
struct PermTestResult{S <: Union{Real, Vector{<:Real}}, N <: Union{Real, Vector{<:Real}}, T <: Any}

	iterations::Int
	obs::S
	ran::Vector{N}
	test::T
	tested_regions::Union{String, Nothing}
	randomised_regions::Union{String, Nothing}
	eval_function::Function
	random_strategy::String
end

"""
	HypothesisTests.pvalue(test::GenomePermutations.PermTestResult)

Extends HypothesisTests.pvalue to get the P value of a GenomePermutations.PermTestResult.
"""
function HypothesisTests.pvalue(test::GenomePermutations.PermTestResult)
	return HypothesisTests.pvalue(test.test)
end

"""
	HypothesisTests.pvalue(test::GenomePermutations.SimplePTest)

Extends HypothesisTests.pvalue to get the P value of a GenomePermutations.SimplePTest.
"""
function HypothesisTests.pvalue(test::GenomePermutations.SimplePTest)
	return(test.p_val)
end

"""
	overlappermtest(col1, col2, regions, iterations; bed_names, allow_overlap = false,
	max_tries = 1000, onfail = :throws)

Legacy function to run an overlap permutation test across 2 regions.
will be replaced by `permtest` using countoverlapping as `f` argument.
Will need a `g` function to summarise overlaps across regions. (likely `sum`).

# Arguments

- `col1::GenomicFeatures.IntervalCollection{T}`: the collection to randomise.
- `col2::GenomicFeatures.IntervalCollection{S}`: the constant collection to test.
- `regions::GenomicFeatures.IntervalCollection{S}`: the regions col1 is randomised in.
- `iterations::Int`: the number of iterations to run.
- `bed_names::Union{String, Nothing}`: the names of the bed files to write to.
- `allow_overlap::Bool = false`: whether to allow overlaps in the randomised regions.
- `max_tries::Int`: the maximum number of tries to generate a randomised collection.
- `onfail::Symbol = :throw`: what to do if a random interval fails to fit in the region.
Will return the orignal interval if :orig, skip the interval if :skip, otherwise throws an error.
"""
function overlappermtest(col1, col2, regions, iterations;
	bed_names = (nothing,nothing),
	allow_overlap = false,
	max_tries::Int = 1000,
	onfail = :throw
)

	# count base overlaps
	obs = countoverlapping(col1,col2)
	if !(isin(col2,regions))
		@warn "Intervals to randomise are not a subset of target regions for the randomisation"
	end
	if !(keys(col1.trees) ⊆ keys(regions.trees) && keys(col2.trees) ⊆ keys(regions.trees))
		throw(ErrorException("Sequences in either collection does not match selected regions"))
	end
	# iterate per chromosome, generating all distributions
	D = Dict{String, Distributions.MixtureModel}()
	for seq in keys(regions.trees)
		push!(D, seq => generatedistribution(regions.trees[seq]))
	end

	# iterate per chromosome get all regions

	r = Dict{String, GenomicFeatures.IntervalCollection{}}()
	for seq in keys(regions.trees)
		push!(r, seq => GenomePermutations.iter_getcollection(regions, seq))
	end

	a = Dict{String, GenomicFeatures.IntervalCollection{}}()
	for seq in keys(col1.trees)
		push!(a, seq => GenomePermutations.iter_getcollection(col1, seq))
	end

	b = Dict{String, GenomicFeatures.IntervalCollection{}}()
	for seq in keys(col2.trees)
		push!(b, seq => GenomePermutations.iter_getcollection(col2, seq))
	end

	# over the number of iterations
	ran = Array{Int}(undef, iterations)

	# Because we made an array of right size, no need for bounds checks
	Threads.@threads for i in 1:iterations
		tmprand = Array{Int}(undef,0)
		for seq in keys(col2.trees)
			# check if the sequence is in collection 1 too
			seq in keys(col1.trees) || continue

			# run chromosome by chromosome.
			rand_a = randomiseregions(a[seq],D[seq],r[seq]; allow_overlap,
									  max_tries, onfail)
			push!(tmprand, countoverlapping(rand_a, b[seq]))
		end
		ran[i] = sum(tmprand)
	end
	# get p value for overlap test

	test = simpleP(obs, ran, iterations)

	return PermTestResult(iterations, obs, ran, test, bed_names[1], bed_names[2],
						  countoverlapping, "guess and check")
end

"""
	_randomgenome(regions, n, C, L, names; allow_overlap = false)

Generates a new random genome with n segments from a named sequence distribution

# Arguments

- `regions`: the regions the genome to generate belongs to.
- `n::Int`: the number of segments to generate.
- `C::Distributions.Categorical`: a distibution with each cateogy representing
one of the named sequences. Lenght and order must match names.
- `L::Distributions.DiscreteUniform": the segment lenght distribution.
- `names::Vector{Strimg}`: the names of the sequences. must match the categories in C.
- `allow_overlap::Bool = false`: whether to allow overlapping segments in the result.

This is currently private as I have not tested this and I don't like the categorical distribution
squence names situation.
"""
function _randomgenome(
	regions,
	n::Int, C::Distributions.Categorical,
	L::Distributions.DiscreteDistribution,
	names::Vector{String};
	allow_overlap = false
)

	#generate start distrubutions - Soon to be spun off into a new function
	S = Dict{String, Distributions.MixtureModel}()
	for seq in keys(regions.trees)
		push!(S, seq => generatedistribution(regions.trees[seq]))
	end

	# Check sequence distribution
	Distributions.ncategories(C) == length(names) || throw(KeyError("Number of Sequence in distribution $C does not match names $names"))
	length(S) == length(names) ||  throw(KeyError("Number of sequences in $regions does not match names $names"))

	col = GenomicFeatures.IntervalCollection{Nothing}()
	for i in 1:n
		int = GenomicFeatures.Interval(names[rand(C)],1, rand(L)+1)
		push!(col, _randominterval(int, S[int.seqname], regions; collection = col, allow_overlap = allow_overlap, onfail = :throw))
	end
	return(col)
end

"""
	permtest(col1, col2, regions, iterations, f; allow_overlap = false, max_tries::Int = 1000, onfail = :throw)

Runs a flexible permuation test, using a custsom evaluation function `f`.

# Arguments

- `col1::GenomicFeatures.IntervalCollection{}`: the collection to randomise.
- `col2::GenomicFeatures.IntervalCollection{}`: the constant collection to test.
- `regions::GenomicFeatures.IntervalCollection{}`: the regions col1 is randomised in.
- `iterations::Int`: the number of iterations to run.
- `f::Function = GenomePermutations.dist`: the evaluation function.
Defaults to GenomePermutations.dist.
- `bed_names::Union{String, Nothing}`: the names of the bed files to write to.
- `allow_overlap::Bool = false`: whether to allow overlaps in the randomised regions.
- `max_tries::Int = 1000`: the maximum number of tries to generate a randomised collection.
- `onfail::Symbol = :throw`: what to do if a random interval fails to fit in the region.
Will return the orignal interval if :orig, skips it if :skip, otherwise throws an error.
"""
function permtest(col1, col2, regions, iterations, f::Function = GenomePermutations.dist;
	bed_names = (nothing,nothing), allow_overlap = false, max_tries::Int = 1000, onfail = :throw)

	# count base overlaps
	obs = f(col1, col2)

	# should spin this off into a function, rather than copy paste from the other function
	if !(isin(col2, regions))
		@warn "Intervals to randomise are not a subset of target regions for the randomisation"
	end
	if !(keys(col1.trees) ⊆ keys(regions.trees) && keys(col2.trees) ⊆ keys(regions.trees))
		throw(ErrorException("Sequences in either collection does not match selected regions"))
	end

	# iterate per chromosome, generating all distributions
	D = Dict{String, Distributions.MixtureModel}()
	for seq in keys(regions.trees)
		push!(D, seq => generatedistribution(regions.trees[seq]))
	end

	# iterate per chromosome get all regions
	r = Dict{String, GenomicFeatures.IntervalCollection{}}()
	for seq in keys(regions.trees)
		push!(r, seq => GenomePermutations.iter_getcollection(regions, seq))
	end

	a = Dict{String, GenomicFeatures.IntervalCollection{}}()
	for seq in keys(col1.trees)
		push!(a, seq => GenomePermutations.iter_getcollection(col1, seq))
	end

	b = Dict{String, GenomicFeatures.IntervalCollection{}}()
	for seq in keys(col2.trees)
		push!(b, seq => GenomePermutations.iter_getcollection(col2, seq))
	end

	# over the number of iterations

	ran = Array{typeof(obs)}(undef, iterations)
	Threads.@threads for i in 1:iterations

		tmprand = Array{typeof(obs)}(undef,0)
		for seq in keys(col2.trees)
			# check if the sequence is in collection 1 too
			seq in keys(col1.trees) || continue

			# run chromosome by chromosome.
			rand_a = randomiseregions(a[seq],D[seq],r[seq]; allow_overlap,
									  max_tries, onfail)
			push!(tmprand, f(rand_a, b[seq]))
		end
		ran[i] = reduce(vcat, tmprand)
	end

	# get p value for overlap test, if we only have 1 observaton per iteration
	if length(obs) == 1
		test = simpleP(obs, ran, iterations)
	else
		ran = map(Statistics.mean, ran)
		obs = Statistics.mean(reduce(vcat, obs))
		#obs = map(Statistics.mean, map(Statistics.mean, obs))
		test = simpleP(obs, ran, iterations)
	end

	return PermTestResult(iterations, obs, ran, test, bed_names[1], bed_names[2], f,
						  "guess and check")

end

end


text:
using Base.Test
using MathProgBase
using MathProgBase.SolverInterface

function linprogsolvertest(solver::AbstractMathProgSolver, eps = Base.rtoldefault(Float64))

    m = LinearQuadraticModel(solver)

    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0
    loadproblem!(m, [1. 1.], [0.,0.], [Inf, Inf], [-1.,0.], [-Inf],[1.], :Min)
    @test numvar(m) == 2
    @test numconstr(m) == 1
    @test getvarLB(m) == [0.,0.]
    @test getvarUB(m)[1] >= 1e20
    @test getvarUB(m)[2] >= 1e20
    @test getconstrLB(m)[1] < -1e30
    @test getconstrUB(m)[1] == 1.0
    @test getobj(m) == [-1.,0.]
    @test getsense(m) == :Min

    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) -1 eps
    @test_approx_eq_eps getsolution(m) [1.0, 0.0] eps
    @test_approx_eq_eps getconstrsolution(m) [1.0] eps
    @test_approx_eq_eps getconstrduals(m) [-1.0] eps
    @test_approx_eq_eps getreducedcosts(m) [0.0, 1.0] eps

    setsense!(m, :Max)
    # max x
    setobj!(m, [1.0,0.0])
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 1 eps
    @test_approx_eq_eps getsolution(m) [1.0, 0.0] eps
    @test_approx_eq_eps getconstrsolution(m) [1.0] eps
    @test_approx_eq_eps getconstrduals(m) [1.0] eps
    @test_approx_eq_eps getreducedcosts(m) [0.0, -1.0] eps

    # add new variable to get:
    # max x + 2z
    # s.t. x + y + z <= 1
    # x,y,z >= 0
    addvar!(m, [1], [1.0], 0, Inf, 2.0)

    @test numvar(m) == 3
    @test numconstr(m) == 1

    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 2 eps
    @test_approx_eq_eps getsolution(m) [0.0, 0.0, 1.0] eps
    @test_approx_eq_eps getconstrsolution(m) [1.0] eps
    @test_approx_eq_eps getconstrduals(m) [2.0] eps
    @test_approx_eq_eps getreducedcosts(m) [-1.0, -2.0, 0.0] eps

    setvarLB!(m, [-1.0,0.0,0.0])
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 3 eps

    # fix z to zero
    setvarLB!(m, [0.0,0.0,0.0])
    setvarUB!(m, [Inf,Inf,0.0])
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 1 eps

    setconstrUB!(m, [2.0])
    setconstrLB!(m, [2.0])
    optimize!(m)
    @test_approx_eq_eps getobjval(m) 2 eps

    setobj!(m, [1.0, 2.0, 0.0])
    optimize!(m)
    @test_approx_eq_eps getobjval(m) 4 eps
    @test_approx_eq_eps getsolution(m) [0.0, 2.0, 0.0] eps

    # we now have:
    # max x+2y
    # s.t. x + y + z == 2
    # x,y >= 0, z = 0

    # add constraint x - y >= 0
    addconstr!(m, [1,2], [1.0,-1.0], 0.0, Inf)
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 3 eps
    @test_approx_eq_eps getsolution(m) [1.0, 1.0, 0.0] eps
    @test_approx_eq_eps getconstrduals(m) [1.5, -0.5] eps
    @test_approx_eq_eps getreducedcosts(m) [0.0, 0.0, -1.5] eps

    # test addvar! interface

    m = LinearQuadraticModel(solver)

    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0
    addvar!(m, 0, Inf, -1)
    addvar!(m, 0, Inf, 0)
    addconstr!(m, [1, 2], [1.0, 1.0], -Inf, 1.0)
    @test numvar(m) == 2
    @test numconstr(m) == 1
    @test getvarLB(m) == [0.,0.]
    @test getvarUB(m)[1] >= 1e20
    @test getvarUB(m)[2] >= 1e20
    @test getconstrLB(m)[1] <= -1e20
    @test getconstrUB(m)[1] == 1.0
    @test getobj(m) == [-1.,0.]
    @test getsense(m) == :Min

    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) -1 eps
    @test_approx_eq_eps getsolution(m) [1.0, 0.0] eps
    @test_approx_eq_eps getconstrsolution(m) [1.0] eps
    @test_approx_eq_eps getconstrduals(m) [-1.0] eps
    @test_approx_eq_eps getreducedcosts(m) [0.0, 1.0] eps

    ####################################
    # test setconstrLB! and setconstrUB!
    # Test that:
    #   Modifying lower bound works
    #   Modifying upper bound works
    #   Setting upper and lower bound to same value works
    #

    m = LinearQuadraticModel(solver)
    # Min  x - y
    # s.t. 0.0 <= x
    #             y <= 0.0
    # x,y unbounded

    loadproblem!(m, [ 1.0 0.0 ; 0.0 1.0 ], [-Inf, -Inf], [Inf,Inf], [1.0, -1.0], [0.0,-Inf], [Inf,0.0], :Min)

    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 0.0 eps
    @test_approx_eq_eps getsolution(m) [ 0.0, 0.0 ] eps

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= 0.0
    # x,y unbounded
    setconstrLB!(m,[100.0,-Inf])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 100.0 eps
    @test_approx_eq_eps getsolution(m) [ 100.0, 0.0 ] eps

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= -100.0
    # x,y unbounded
    setconstrUB!(m,[Inf,-100.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq_eps getobjval(m) 200.0 eps
    @test_approx_eq_eps getsolution(m) [ 100.0, -100.0 ] eps

    # Test issue #40 from Gurobi.jl
    # min  x
    # s.t. x >= 0
    #      x >= 3
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [1.0 1.0]', [-Inf], [Inf], [1.0], [0.0, 3.0], [Inf, Inf], :Min)
    optimize!(m)
    for i = 1:length(getconstrLB(m))
        @test getconstrLB(m)[i] <= getconstrsolution(m)[i] + eps
        @test getconstrsolution(m)[i] <= getconstrUB(m)[i] + eps
    end

    # min  x
    # s.t. x <= 0
    #      x <= 3
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [1.0 1.0]', [-Inf], [Inf], [1.0], [-Inf, -Inf], [0.0, 3.0], :Max)
    optimize!(m)
    for i = 1:length(getconstrLB(m))
        @test getconstrLB(m)[i] <= getconstrsolution(m)[i] + eps
        @test getconstrsolution(m)[i] <= getconstrUB(m)[i] + eps
    end
end

function linprogsolvertestextra(solver::AbstractMathProgSolver)
    ####################################
    # test setconstrLB! and setconstrUB!
    # Test that:
    #   Modifying lower bound works
    #   Modifying upper bound works
    #   Setting upper and lower bound to same value works
    #

    m = LinearQuadraticModel(solver)
    # Min  x - y
    # s.t. 0.0 <= x <= 0.0
    #      0.0 <= y <= 0.0
    # x,y unbounded
    loadproblem!(m, [ 1.0 0.0 ; 0.0 1.0 ], [-Inf, -Inf], [Inf,Inf], [1.0, -1.0], [0.0,0.0], [0.0,0.0], :Min)

    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 0.0
    @test_approx_eq getsolution(m) [ 0.0, 0.0 ]

    # Min  x - y
    # s.t. 0.0 <= x <= 100.0
    #      0.0 <= y <= 100.0
    # x,y unbounded
    setconstrUB!(m,[100.0,100.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) -100.0
    @test_approx_eq getsolution(m) [ 0.0, 100.0 ]

    # Min  x - y
    # s.t. -100.0 <= x <= 100.0
    #      -100.0 <= y <= 100.0
    # x,y unbounded
    setconstrLB!(m,[-100.0,-100.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) -200.0
    @test_approx_eq getsolution(m) [ -100.0, 100.0 ]

    # Min  x - y
    # s.t. -100.0 <= x <= 100.0
    #      -100.0 <= y <= 100.0
    # x,y unbounded
    setconstrLB!(m,[10.0,10.0])
    setconstrUB!(m,[10.0,10.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 0.0
    @test_approx_eq getsolution(m) [ 10.0, 10.0 ]

    # Min  x - y
    # s.t. 0.0  <= x <= Inf
    #      -Inf <= y <= 0.0
    # x,y unbounded
    setconstrLB!(m,[0.0,-Inf])
    setconstrUB!(m,[Inf,0.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 0.0
    @test_approx_eq getsolution(m) [ 0.0, 0.0 ]

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= 0.0
    # x,y unbounded
    setconstrLB!(m,[-Inf,-Inf])
    setconstrUB!(m,[Inf,0.0])
    optimize!(m)
    @test status(m) == :Unbounded
    @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8

    # Min  x - y
    # s.t. 0.0 <= x <= Inf
    #      -Inf <= y <= Inf
    # x,y unbounded
    setconstrLB!(m,[0.0,-Inf])
    setconstrUB!(m,[Inf,Inf])
    optimize!(m)
    @test status(m) == :Unbounded
    @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8

    ####################################
    # test setvarLB! and setvarUB!
    # Test that:
    #   Modifying lower bound works
    #   Modifying upper bound works
    #   Setting upper and lower bound to same value works
    #

    m = LinearQuadraticModel(solver)
    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    # 0.0 <= x <= 0.0
    # 0.0 <= y <= 0.0
    loadproblem!(m, [ 1.0 0.0 ; 0.0 1.0 ], [0.0,0.0], [0.0,0.0], [1.0, -1.0], [-Inf, -Inf], [Inf,Inf], :Min)

    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 0.0
    @test_approx_eq getsolution(m) [ 0.0, 0.0 ]

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    # 0.0 <= x <= 100.0
    # 0.0 <= y <= 100.0

    # x,y unbounded
    setvarUB!(m,[100.0,100.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) -100.0
    @test_approx_eq getsolution(m) [ 0.0, 100.0 ]

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    # -100.0 <= x <= 100.0
    # -100.0 <= y <= 100.0
    setvarLB!(m,[-100.0,-100.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) -200.0
    @test_approx_eq getsolution(m) [ -100.0, 100.0 ]

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    # -100.0 <= x <= 100.0
    # -100.0 <= y <= 100.0
    setvarLB!(m,[10.0,10.0])
    setvarUB!(m,[10.0,10.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 0.0
    @test_approx_eq getsolution(m) [ 10.0, 10.0 ]

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    # 0.0  <= x <= Inf
    # -Inf <= y <= 0.0
    # x,y unbounded
    setvarLB!(m,[0.0,-Inf])
    setvarUB!(m,[Inf,0.0])
    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 0.0
    @test_approx_eq getsolution(m) [ 0.0, 0.0 ]

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    # -Inf <= x <= Inf
    # -Inf <= y <= 0.0
    setvarLB!(m,[-Inf,-Inf])
    setvarUB!(m,[Inf,0.0])
    optimize!(m)
    @test status(m) == :Unbounded
    @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8

    # Min  x - y
    # s.t. -Inf <= x <= Inf
    #      -Inf <= y <= Inf
    #  0.0 <= x <= Inf
    # -Inf <= y <= Inf
    # x,y unbounded
    setvarLB!(m,[0.0,-Inf])
    setvarUB!(m,[Inf,Inf])
    optimize!(m)
    @test status(m) == :Unbounded
    @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8
end


text:
abstract type GraphicsArea end

function nrows end
function ncols end
function printrow end

function Base.print(io::IO, c::GraphicsArea)
    for row in 1:nrows(c)
        printrow(io, c, row)
        row < nrows(c) && print(io, '\n')
    end
end

function Base.show(io::IO, c::GraphicsArea)
    b = border_solid
    border_length = ncols(c)
    print_border_top(io, "", border_length, :solid, :light_black)
    print(io, '\n')
    for row in 1:nrows(c)
        print_color(:light_black, io, b[:l])
        printrow(io, c, row)
        print_color(:light_black, io, b[:r])
        print(io, '\n')
    end
    print_border_bottom(io, "", border_length, :solid, :light_black)
end