text:
# Specify model parameters and/or initial values to optimize
function get_search_index()::Tuple{Array{Int64,1},Array{Int64,1}}
    # parameters
    search_idx_params::Vector{Int} = [
        C.VmaxPY
        C.KmPY
        C.kdeg
        C.kf47
        C.Vmaxr47
        C.Kmf47
        C.Kmr47
        C.kf48
        C.Kmf48
        C.Kmr48
        #C.PTEN
        C.kf49
        C.kr49
        C.Kmf49
        C.Kmr49
        C.Kmr49b
        C.kr49b
        C.kf51
        C.Vmaxr51
        C.Kmf51
        C.Kmr51
        C.Kmrb51
        C.kf52
        C.Vmaxr52
        C.Kmf52
        C.Kmr52
        C.kf54
        C.Vmaxr54
        C.Kmf54
        C.Kmr54
        C.kf55
        C.Vmaxr55
        C.Kmf55
        C.Kmr55
        C.kf38
        C.kf39
        C.kf50
        C.a98
        C.b98
        C.koff46
        C.EGF_off
        C.HRGoff_3
        C.HRGoff_4
        C.koff4
        C.koff5
        C.koff6
        C.koff7
        C.koff8
        C.koff9
        C.koff61
        C.koff62
        C.koff16
        C.koff17
        C.koff18
        C.koff19
        C.koff20
        C.koff21
        C.koff22
        C.koff23
        C.koff24
        C.koff25
        C.koff26
        C.koff27
        C.koff28
        C.koff29
        C.koff30
        C.koff31
        C.koff32
        C.koff33
        C.koff34
        C.koff35
        C.koff36
        C.koff37
        C.koff65
        C.koff66
        C.koff67
        C.koff68
        C.koff69
        C.koff70
        C.koff71
        C.koff72
        C.koff40
        C.koff41
        C.koff42
        C.koff43
        C.koff44
        C.koff45
        C.koff57
        C.koff58
        C.koff59
        C.koff60
        C.kPTP10
        C.kPTP11
        C.kPTP12
        C.kPTP13
        C.kPTP14
        C.kPTP15
        C.kPTP63
        C.kPTP64
        C.koff73
        C.koff74
        C.koff75
        C.koff76
        C.koff77
        C.koff78
        C.koff79
        C.koff80
        C.kPTP38
        C.kPTP39
        C.koff88
        C.kPTP50
        C.kf81
        C.Vmaxr81
        C.Kmf81
        C.Kmr81
        C.kf82
        C.Vmaxr82
        C.Kmf82
        C.Kmr82
        C.kf83
        C.Vmaxr83
        C.Kmf83
        C.Kmr83
        C.kf84
        C.Vmaxr84
        C.Kmf84
        C.Kmr84
        C.kf85
        C.Vmaxr85
        C.Kmf85
        C.Kmr85
        C.kcon49
        C.kon1
        C.kon86
        C.kon2
        C.kon3
        C.kon87
        C.kon4
        C.kon5
        C.kon6
        C.kon7
        C.kon8
        C.kon9
        C.kon61
        C.kon62
        C.kf10
        C.kf11
        C.kf12
        C.kf13
        C.kf14
        C.kf15
        C.kf63
        C.kf64
        C.kon16
        C.kon17
        C.kon18
        C.kon73
        C.kon19
        C.kon20
        C.kon21
        C.kon74
        C.kon22
        C.kon23
        C.kon24
        C.kon25
        C.kon75
        C.kon26
        C.kon27
        C.kon28
        C.kon29
        C.kon76
        C.kon30
        C.kon31
        C.kon32
        C.kon33
        C.kon77
        C.kon34
        C.kon35
        C.kon36
        C.kon37
        C.kon78
        C.kon65
        C.kon66
        C.kon67
        C.kon68
        C.kon79
        C.kon69
        C.kon70
        C.kon71
        C.kon72
        C.kon80
        C.kon40
        C.kon41
        C.kon42
        C.kon43
        C.kon44
        C.kon45
        C.kon88
        C.kon46
        C.kon57
        C.kon58
        C.kon59
        C.kon60
        #
        C.V1
        C.Km1
        C.V5
        C.Km5
        C.V10
        C.Km10
        #C.n10
        C.p11
        C.p12
        C.p13
        C.V14
        C.Km14
        C.V15
        C.Km15
        C.KimDUSP
        C.KexDUSP
        C.V20
        C.Km20
        C.V21
        C.Km21
        C.V24
        C.Km24
        C.V25
        C.Km25
        C.KimRSK
        C.KexRSK
        C.V27
        C.Km27
        C.V28
        C.Km28
        C.V29
        C.Km29
        C.V30
        C.Km30
        C.V31
        C.Km31
        #C.n31
        C.p32
        C.p33
        C.p34
        C.V35
        C.Km35
        C.V36
        C.Km36
        C.V37
        C.Km37
        C.KimFOS
        C.KexFOS
        C.V42
        C.Km42
        C.V43
        C.Km43
        C.V44
        C.Km44
        C.p47
        C.m47
        C.p48
        C.p49
        C.m49
        C.p50
        C.p51
        C.m51
        C.V57
        C.Km57
        #C.n57
        C.p58
        C.p59
        C.p60
        C.p61
        C.KimF
        C.KexF
        C.p63
        C.KF31
        #C.nF31
        #
        C.scale_E1
        C.scale_E2
        C.scale_E3
        C.scale_E4
        C.scale_G
        #C.scale_S
        C.a_SHC1
        C.a_SHC2
        C.a_SHC3
        C.a_SHC4
        #C.scale_I
        C.a_PIK3CA
        C.a_PIK3CB
        C.a_PIK3CD
        C.a_PIK3CG
        #ver.5
        C.scale_PTEN
        #C.scale_R
        C.a_RASA1
        C.a_RASA2
        C.a_RASA3
        #C.scale_O
        C.a_SOS1
        C.a_SOS2
        C.scale_A
        #C.scale_Akt
        C.a_AKT1
        C.a_AKT2
        #C.scale_RsD
        C.a_HRAS
        C.a_KRAS
        C.a_NRAS
        #C.scale_Raf
        C.a_ARAF
        C.a_BRAF
        C.a_RAF1
        #C.scale_MEK
        C.a_MAP2K1
        C.a_MAP2K2
        C.scale_T
        C.scale_CREBn
        #C.scale_ERKc
        C.a_MAPK1
        C.a_MAPK3
        C.scale_Elk1n
        #C.scale_RSKc
        C.a_RPS6KA1
        C.a_RPS6KA2
        C.a_RPS6KA3
        C.a_RPS6KA6
    ]

    # initial values
    search_idx_initials::Vector{Int} = [
        V.P2,
    ]

    return search_idx_params, search_idx_initials
end

function get_search_region()::Matrix{Float64}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()
    search_param::Vector{Float64} = init_search_param(search_idx, p, u0)

    search_rgn::Matrix{Float64} = zeros(2, length(p)+length(u0))

    # Default: 0.1 ~ 10x
    for (i,j) in enumerate(search_idx[1])
        search_rgn[1,j] = search_param[i]*0.1  # lower bound
        search_rgn[2,j] = search_param[i]*10.0  # upper bound
    end

    # Default: 0.5 ~ 2x
    for (i,j) in enumerate(search_idx[2])
        search_rgn[1,j+length(p)] = search_param[i+length(search_idx[1])]*0.5  # lower bound
        search_rgn[2,j+length(p)] = search_param[i+length(search_idx[1])]*2.0  # upper bound
    end

    # search_rgn[:, C.param_name] = [lower_bound, upper_bound]
    # search_rgn[:, V.var_name+length(p)] = [lower_bound, upper_bound]

    #search_rgn[:,C.PTEN] = [1.0,1000.0]
    search_rgn[:,V.P2+length(p)] = [1.0,1000.0]

    search_rgn[:,C.scale_E1] = [0.1, 100.0]
    search_rgn[:,C.scale_E2] = [0.1, 100.0]
    search_rgn[:,C.scale_E3] = [0.1, 100.0]
    search_rgn[:,C.scale_E4] = [0.1, 100.0]
    search_rgn[:,C.scale_G] = [0.1, 100.0]
    #search_rgn[:,C.scale_S] = [0.1, 100.0]
    search_rgn[:,C.a_SHC1] = [0.1, 100.0]
    search_rgn[:,C.a_SHC2] = [0.1, 100.0]
    search_rgn[:,C.a_SHC3] = [0.1, 100.0]
    search_rgn[:,C.a_SHC4] = [0.1, 100.0]
    #search_rgn[:,C.scale_I] = [0.1, 100.0]
    search_rgn[:,C.a_PIK3CA] = [0.1, 100.0]
    search_rgn[:,C.a_PIK3CB] = [0.1, 100.0]
    search_rgn[:,C.a_PIK3CD] = [0.1, 100.0]
    search_rgn[:,C.a_PIK3CG] = [0.1, 100.0]
    # ver.5
    search_rgn[:,C.scale_PTEN] = [0.1, 100.0]
    #search_rgn[:,C.scale_R] = [0.1, 100.0]
    search_rgn[:,C.a_RASA1] = [0.1, 100.0]
    search_rgn[:,C.a_RASA2] = [0.1, 100.0]
    search_rgn[:,C.a_RASA3] = [0.1, 100.0]
    #search_rgn[:,C.scale_O] = [0.1, 100.0]
    search_rgn[:,C.a_SOS1] = [0.1, 100.0]
    search_rgn[:,C.a_SOS2] = [0.1, 100.0]
    search_rgn[:,C.scale_A] = [0.1, 100.0]
    #search_rgn[:,C.scale_Akt] = [0.1, 100.0]
    search_rgn[:,C.a_AKT1] = [0.1, 100.0]
    search_rgn[:,C.a_AKT2] = [0.1, 100.0]
    #search_rgn[:,C.scale_RsD] = [0.1, 100.0]
    search_rgn[:,C.a_HRAS] = [0.1, 100.0]
    search_rgn[:,C.a_KRAS] = [0.1, 100.0]
    search_rgn[:,C.a_NRAS] = [0.1, 100.0]
    #search_rgn[:,C.scale_Raf] = [0.1, 100.0]
    search_rgn[:,C.a_ARAF] = [0.1, 100.0]
    search_rgn[:,C.a_BRAF] = [0.1, 100.0]
    search_rgn[:,C.a_RAF1] = [0.1, 100.0]
    #search_rgn[:,C.scale_MEK] = [0.1, 100.0]
    search_rgn[:,C.a_MAP2K1] = [0.1, 100.0]
    search_rgn[:,C.a_MAP2K2] = [0.1, 100.0]
    search_rgn[:,C.scale_T] = [0.1, 100.0]
    search_rgn[:,C.scale_CREBn] = [0.1, 100.0]
    #search_rgn[:,C.scale_ERKc] = [0.1, 100.0]
    search_rgn[:,C.a_MAPK1] = [0.1, 100.0]
    search_rgn[:,C.a_MAPK3] = [0.1, 100.0]
    search_rgn[:,C.scale_Elk1n] = [0.1, 100.0]
    #search_rgn[:,C.scale_RSKc] = [0.1, 100.0]
    search_rgn[:,C.a_RPS6KA1] = [0.1, 100.0]
    search_rgn[:,C.a_RPS6KA2] = [0.1, 100.0]
    search_rgn[:,C.a_RPS6KA3] = [0.1, 100.0]
    search_rgn[:,C.a_RPS6KA6] = [0.1, 100.0]

    search_rgn = conv_lin2log!(search_rgn, search_idx)

    return search_rgn
end

function update_param(indiv::Vector{Float64})::Tuple{Array{Float64,1},Array{Float64,1}}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()

    for (i,j) in enumerate(search_idx[1])
        @inbounds p[j] = indiv[i]
    end
    for (i,j) in enumerate(search_idx[2])
        @inbounds u0[j] = indiv[i+length(search_idx[1])]
    end

    # scaling
    (p, u0) = tpm2ival!(p, u0, "MCF7")
    #
    # constraints --------------------------------------------------------------
    p[C.V6] = p[C.V5]
    p[C.Km6] = p[C.Km5]
    p[C.KimpDUSP] = p[C.KimDUSP]
    p[C.KexpDUSP] = p[C.KexDUSP]
    p[C.KimpcFOS] = p[C.KimFOS]
    p[C.KexpcFOS] = p[C.KexFOS]
    p[C.p52] = p[C.p47]
    p[C.m52] = p[C.m47]
    p[C.p53] = p[C.p48]
    p[C.p54] = p[C.p49]
    p[C.m54] = p[C.m49]
    p[C.p55] = p[C.p50]
    p[C.p56] = p[C.p51]
    p[C.m56] = p[C.m51]
    # --------------------------------------------------------------------------

    return p, u0
end

function decode_gene2val(indiv_gene::Vector{Float64})::Vector{Float64}
    search_rgn::Matrix{Float64} = get_search_region()
    indiv::Vector{Float64} = zeros(length(indiv_gene))

    for i in eachindex(indiv_gene)
        indiv[i] = 10^(
            indiv_gene[i] * (
                search_rgn[2,i] - search_rgn[1,i]
            ) + search_rgn[1,i]
        )
    end

    return round.(indiv,sigdigits=7)
end

function encode_val2gene(indiv::Vector{Float64})
    search_rgn::Matrix{Float64} = get_search_region()
    indiv_gene::Vector{Float64} = zeros(length(indiv))

    for i in eachindex(indiv)
        indiv_gene[i] = (
            log10(indiv[i]) - search_rgn[1,i]
        ) / (
            search_rgn[2,i] - search_rgn[1,i]
        )
    end

    return indiv_gene
end

function encode_bestIndivVal2randGene(j::Int64, best_indiv::Vector{Float64},
                                        p0_bounds::Vector{Float64})::Float64
    search_rgn::Matrix{Float64} = get_search_region()
    rand_gene::Float64 = (
        log10(
            best_indiv[j]*10^(
                rand() * log10(p0_bounds[2]/p0_bounds[1]) + log10(p0_bounds[1])
            )
        ) - search_rgn[1,j]
    ) / (
        search_rgn[2,j] - search_rgn[1,j]
    )

    return rand_gene
end

function init_search_param(search_idx::Tuple{Array{Int64,1},Array{Int64,1}},
                            p::Vector{Float64}, u0::Vector{Float64})::Vector{Float64}
    search_param = zeros(
        length(search_idx[1]) + length(search_idx[2])
    )
    for (i,j) in enumerate(search_idx[1])
        @inbounds search_param[i] = p[j]
    end
    for (i,j) in enumerate(search_idx[2])
        @inbounds search_param[i+length(search_idx[1])] = u0[j]
    end

    if any(x -> x == 0.0, search_param)
        message::String = "search_param must not contain zero."
        for (_, idx) in enumerate(search_idx[1])
            if p[idx] == 0.0
                error(
                    @sprintf(
                        "`C.%s` in search_idx_params: ", C.parameters[idx]
                    ) * message
                )
            end
        end
        for (_, idx) in enumerate(search_idx[2])
            if u0[idx] == 0.0
                error(
                    @sprintf(
                        "`V.%s` in search_idx_initials: ", V.species[idx]
                    ) * message
                )
            end
        end
    end

    return search_param
end

function conv_lin2log!(search_rgn::Matrix{Float64},
                        search_idx::Tuple{Array{Int64,1},Array{Int64,1}}
                        )::Matrix{Float64}
    for i=1:size(search_rgn,2)
        if minimum(search_rgn[:,i]) < 0.0
            message = "search_rgn[lb,ub] must be positive.\n"
            if i <= C.n_parameters
                error(
                    @sprintf(
                        "`C.%s` ", C.parameters[i]
                    ) * message
                )
            else
                error(
                    @sprintf(
                        "`V.%s` ", V.species[i-C.n_parameters]
                    ) * message
                )
            end
        elseif minimum(search_rgn[:,i]) == 0.0 && maximum(search_rgn[:,i]) != 0.0
            message = "lower_bound must be larger than 0.\n"
            if i <= C.n_parameters
                error(
                    @sprintf(
                        "`C.%s` ", C.parameters[i]
                    ) * message
                )
            else
                error(
                    @sprintf(
                        "`V.%s` ", V.species[i-C.n_parameters]
                    ) * message
                )
            end
        elseif search_rgn[2,i] - search_rgn[1,i] < 0.0
            message = "lower_bound < upper_bound\n"
            if i <= C.n_parameters
                error(
                    @sprintf(
                        "`C.%s` ", C.parameters[i]
                    ) * message
                )
            else
                error(
                    @sprintf(
                        "`V.%s` ", V.species[i-C.n_parameters]
                    ) * message
                )
            end
        end
    end

    nonzero_idx::Vector{Int} = []
    for i=1:size(search_rgn,2)
        if search_rgn[:,i] != [0.0,0.0]
            push!(nonzero_idx,i)
        end
    end
    difference::Vector{Int} = collect(
        symdiff(
            Set(nonzero_idx),
            Set(
                append!(
                    search_idx[1], C.n_parameters .+ search_idx[2]
                )
            )
        )
    )
    if length(difference) > 0
        for (i,j) in enumerate(difference)
            if j <= C.n_parameters
                print(
                    @sprintf(
                        "`C.%s`\n", C.parameters[Int(j)]
                    )
                )
            else
                print(
                    @sprintf(
                        "`V.%s`\n", V.species[Int(j)-C.n_parameters]
                    )
                )
            end
        end
        error(
            "Set these search_params in both search_idx and search_rgn."
        )
    end

    search_rgn = search_rgn[:,nonzero_idx]

    return log10.(search_rgn)
end


text:
module Clothoid

using Gtk
import Graphics
using LinearAlgebra

# GUI parameters
curvature_scaling = -2000
curvature_sparsity = 10
tangent_length = 80
show_tangents = false
show_endpoints = false
point_size = 5
rectangle_scale = 1.4

# GUI variables
closed_curve = false
subsampling = 25
iterations = 100
alpha = -1.0
precise = false
precision = 1.0e-5

# Global variables
points = []
curve = []
indices = []
curvature_comb = []

distance(p, q) = norm(p - q)

function perp2d(p)
    u = [-p[2], p[1]]
    unorm = norm(u)
    unorm > 0 ? u / unorm : u
end

"Like `linspace`, but the last value is not included."
linspace_nolast(start, stop, n) = collect(range(start, stop=stop, length=n+1))[1:n]

"The return value is the discrete curvature at the central point."
function curvature(prev, p, next)
    denom = distance(prev, p) * distance(p, next) * distance(prev, next)
    2 * det(hcat(p - prev, next - p)) / denom
end

"""
Subsamples the given data by `subsampling`, cycling through if the second argument is `true`.
"""
function subsample(original, closedp)
    if length(original) < 2
        return []
    end

    # Compute the average distance & sampling density
    distances = map(norm, original[2:end] - original[1:end-1])
    if closedp
        push!(distances, norm(original[1] - original[end]))
    end
    avg_distance = sum(distances) / length(distances)
    density = subsampling / avg_distance

    n = length(original)
    result = []
    global indices = [1]
    for i in 2:n
        resolution = Int(round(distances[i-1] * density))
        append!(result, linspace_nolast(original[i-1], original[i], resolution))
        push!(indices, length(result) + 1)
    end
    if closedp
        resolution = Int(round(distances[end] * density))
        append!(result, linspace_nolast(original[end], original[1], resolution))
        push!(indices, length(result) + 1)
    else
        push!(result, original[end])
    end

    result
end

function target_curvature(indices, curvatures, closedp)
    result = []
    curvatures = map(x -> sign(x)*abs(x)^(-alpha), curvatures)
    for i in 2:length(curvatures)
        resolution = indices[i] - indices[i-1]
        append!(result, linspace_nolast(curvatures[i-1], curvatures[i], resolution))
    end
    if closedp
        resolution = indices[end] - indices[end-1]
        append!(result, linspace_nolast(curvatures[end], curvatures[1], resolution))
    else
        push!(result, curvatures[end])
    end
    map(x -> sign(x)*abs(x)^(-1/alpha), result)
end

"Updates one point, given its neighbors and the target curvature."
function update(prev, p, next, target)
    # Rearrange the equation
    target *= distance(prev, p) * distance(p, next) * distance(prev, next) / 2.0
    target += prev[1] * next[2] - prev[2] * next[1]
    beta_x, beta_y = next[2] - prev[2], prev[1] - next[1]
    # Now we need: x beta_x + y beta_y = target
    start = (next + prev) / 2
    dir = perp2d(next - prev)
    # [x,y] is of the form start + alpha * dir
    target -= start[1] * beta_x + start[2] * beta_y
    alpha = target / (dir[1] * beta_x + dir[2] * beta_y)
    start + alpha * dir
end

"""
Discrete clothoid curve generation, as in R. Schneider, L. Kobbelt,
_Discrete Fairing of Curves and Surfaces Based on Linear Curvature Distribution_,
In Curve and Surface Design: Saint Malo, pp. 371-380, University Press, 2000.

Uses the indirect approach, with 0 curvature at the ends for open curves.
"""
function generate_curve()
    # Generate initial curve points
    global curve = subsample(points, closed_curve)

    n = length(points)
    if n < 3
        return
    end

    # Approximate clothoid curve
    curve_curv = []
    iter = 0
    while precise || iter < iterations
        # Generate curvature values at the seed points
        point_curv = []
        if closed_curve
            push!(point_curv, curvature(curve[end], curve[1], curve[2]))
        else
            push!(point_curv, 0.0)
        end
        for i in 2:length(points)-1
            j = indices[i]
            push!(point_curv, curvature(curve[j-1], curve[j], curve[j+1]))
        end
        if closed_curve
            j = indices[n]
            push!(point_curv, curvature(curve[j-1], curve[j], curve[j+1]))
        else
            push!(point_curv, 0.0)
        end

        # Generate target curvature values on the curve
        curve_curv = target_curvature(indices, point_curv, closed_curve)

        # Update curve points
        tmp = similar(curve)
        tmp[1] = curve[1]
        for i in 2:n
            tmp[indices[i]] = curve[indices[i]]
            for j in indices[i-1]+1:indices[i]-1
                tmp[j] = update(curve[j-1], curve[j], curve[j+1], curve_curv[j])
            end
        end
        if closed_curve
            for j in indices[end-1]+1:indices[end]-2
                tmp[j] = update(curve[j-1], curve[j], curve[j+1], curve_curv[j])
            end
            tmp[end] = update(curve[end-1], curve[end], curve[1], curve_curv[end])
        end

        # Find maximal distance and exit when small
        change = maximum(norm, curve - tmp)
        curve = tmp
        if change < precision
            break
        end
        iter += 1
    end

    # Generate curvature comb
    normals = similar(curve)
    for i in eachindex(curve)
        if i == 1
            normals[1] = perp2d(curve[2] - curve[1])
        elseif i == length(curve)
            normals[end] = perp2d(curve[end] - curve[end-1])
        else
            normals[i] = perp2d(curve[i+1] - curve[i-1])
        end
    end
    if isempty(curve_curv)
        curve_curv = zeros(length(curve))
    end
    global curvature_comb = curve .+ normals .* curve_curv * curvature_scaling
end

# Graphics

function draw_polygon(ctx, poly, closep = false)
    if isempty(poly)
        return
    end
    Graphics.new_path(ctx)
    Graphics.move_to(ctx, poly[1][1], poly[1][2])
    for p in poly[2:end]
        Graphics.line_to(ctx, p[1], p[2])
    end
    if closep && length(poly) > 2
        Graphics.line_to(ctx, poly[1][1], poly[1][2])
    end
    Graphics.stroke(ctx)
end

draw_callback = @guarded (canvas) -> begin
    ctx = Graphics.getgc(canvas)

    # White background
    Graphics.rectangle(ctx, 0, 0, Graphics.width(canvas), Graphics.height(canvas))
    Graphics.set_source_rgb(ctx, 1, 1, 1)
    Graphics.fill(ctx)

    # Input polygon
    # Graphics.set_source_rgb(ctx, 1, 0, 0)
    # Graphics.set_line_width(ctx, 1.0)
    # draw_polygon(ctx, points, closed_curve)

    # Generated curve
    Graphics.set_source_rgb(ctx, 0.8, 0.3, 0)
    Graphics.set_line_width(ctx, 2.0)
    draw_polygon(ctx, curve, closed_curve)

    # Approximate curvature comb
    if length(points) > 2
        Graphics.set_source_rgb(ctx, 0, 0, 1)
        Graphics.set_line_width(ctx, 1.0)
        draw_polygon(ctx, curvature_comb, closed_curve)
        for i in 1:length(curvature_comb)
            i % curvature_sparsity != 0 && continue
            Graphics.new_path(ctx)
            Graphics.move_to(ctx, curve[i][1], curve[i][2])
            Graphics.line_to(ctx, curvature_comb[i][1], curvature_comb[i][2])
            Graphics.stroke(ctx)
        end
    end

    if subsampling <= 10
        Graphics.set_source_rgb(ctx, 0.4, 0.4, 0.4)
        for p in curve
            Graphics.arc(ctx, p[1], p[2], 3, 0, 2pi)
            Graphics.fill(ctx)
        end
    end

    # Approximated tangents
    if show_tangents
        Graphics.set_source_rgb(ctx, 1, 0, 1)
        Graphics.set_line_width(ctx, 1.0)
        for i in 2:length(points)-1
            j = indices[i]
            dir = curve[j-1] - curve[j+1]
            len = norm(dir)
            if len > 0
                dir /= len
            end
            p1 = points[i] - dir * tangent_length / 2
            p2 = points[i] + dir * tangent_length / 2
            Graphics.new_path(ctx)
            Graphics.move_to(ctx, p1[1], p1[2])
            Graphics.line_to(ctx, p2[1], p2[2])
            Graphics.stroke(ctx)
        end
    end

    # Input points
    for p in points[1:end]
        Graphics.set_source_rgb(ctx, 0, 1, 0)
        Graphics.arc(ctx, p[1], p[2], point_size, 0, 2pi)
        Graphics.fill(ctx)
        Graphics.set_source_rgb(ctx, 0, 0, 0)
        rect = [p + [-point_size, -point_size] * rectangle_scale,
                p + [-point_size,  point_size] * rectangle_scale,
                p + [ point_size,  point_size] * rectangle_scale,
                p + [ point_size, -point_size] * rectangle_scale]
        draw_polygon(ctx, rect, true)
    end
    if show_endpoints && !isempty(points)
        Graphics.set_source_rgb(ctx, 0, 0.8, 0.8)
        for p in (points[1], points[end])
            Graphics.arc(ctx, p[1], p[2], point_size, 0, 2pi)
            Graphics.fill(ctx)
        end
    end
end

# GUI

mousedown_handler = @guarded (canvas, event) -> begin
    p = [event.x, event.y]
    global clicked = findfirst(points) do q
        distance(p, q) < 10
    end
    if clicked === nothing
        push!(points, p)
        clicked = length(points)
        generate_curve()
        draw(canvas)
    end
end

mousemove_handler = @guarded (canvas, event) -> begin
    global clicked
    points[clicked] = [event.x, event.y]
    generate_curve()
    draw(canvas)
end

function setup_gui()
    win = GtkWindow("Discrete Clothoid")
    vbox = GtkBox(:v)

    # Canvas widget
    canvas = GtkCanvas(500, 500)
    canvas.mouse.button1press = mousedown_handler
    canvas.mouse.button1motion = mousemove_handler
    draw(draw_callback, canvas)
    push!(win, vbox)
    push!(vbox, canvas)

    # Reset button
    reset = GtkButton("Start Over")
    signal_connect(reset, "clicked") do _
        global points = []
        global curve = []
        draw(canvas)
    end
    hbox = GtkBox(:h)
    set_gtk_property!(hbox, :spacing, 10)
    push!(vbox, hbox)
    push!(hbox, reset)

    # Closed Checkbox
    closedp = GtkCheckButton("Closed curve")
    set_gtk_property!(closedp, :active, closed_curve)
    signal_connect(closedp, "toggled") do cb
        global closed_curve = get_gtk_property(cb, :active, Bool)
        generate_curve()
        draw(canvas)
    end
    push!(hbox, closedp)

    # Subsampling Spinbutton
    sss = GtkSpinButton(5:5:100)
    set_gtk_property!(sss, :value, subsampling)
    signal_connect(sss, "value-changed") do sb
        global subsampling = get_gtk_property(sb, :value, Int)
        generate_curve()
        draw(canvas)
    end
    push!(hbox, GtkLabel("Subsampling:"))
    push!(hbox, sss)

    # Iterations Spinbutton
    its = GtkSpinButton(0:20:1000)
    set_gtk_property!(its, :value, iterations)
    signal_connect(its, "value-changed") do sb
        global iterations = get_gtk_property(sb, :value, Int)
        generate_curve()
        draw(canvas)
    end
    push!(hbox, GtkLabel("# of iterations:"))
    push!(hbox, its)

    # Alpha Choices
    hbox = GtkBox(:h)
    push!(vbox, hbox)
    push!(hbox, GtkLabel("Alpha: "))
    choices = [-2.5, -1.0, -0.8, -0.5, -0.2, -0.1]
    radios = [GtkRadioButton(string(choice)) for choice in choices]
    for r in radios
        set_gtk_property!(r, :group, radios[1])
        signal_connect(r, "toggled") do _
            global alpha = choices[findfirst(rb -> rb === r, radios)]
            generate_curve()
            draw(canvas)
        end
        push!(hbox, r)
    end
    set_gtk_property!(radios[2], :active, true)

    # Precise Checkbox
    precisep = GtkCheckButton("Precision:")
    set_gtk_property!(precisep, :active, precise)
    signal_connect(precisep, "toggled") do cb
        global precise = get_gtk_property(cb, :active, Bool)
        generate_curve()
        draw(canvas)
    end
    push!(hbox, precisep)

    # Precision Spinbutton
    prs = GtkSpinButton(1:1:10)
    set_gtk_property!(prs, :value, Int(round(-log10(precision))))
    signal_connect(prs, "value-changed") do sb
        global precision = exp10(-get_gtk_property(sb, :value, Int))
        generate_curve()
        draw(canvas)
    end
    push!(hbox, prs)

    generate_curve()
    showall(win)
end

run() = begin setup_gui(); nothing end

end # module

ext: jl

max_stars_repo_name: salvipeter/discrete-clothoids

max_issues_repo_name: salvipeter/discrete-clothoids

max_forks_repo_name: salvipeter/discrete-clothoids


text:
#=
Portion of the source code below, the _lanczos4_opencv method in particular,
are copyrights of the below and are licensed under the terms below with the following
disclaimers. The source code has been significantly modified and ported to Julia.

By downloading, copying, installing or using the software you agree to this license.
If you do not agree to this license, do not download, install,
copy or use the software.

                           License Agreement
                For Open Source Computer Vision Library

Copyright (C) 2000-2008, 2017, Intel Corporation, all rights reserved.
Copyright (C) 2009, Willow Garage Inc., all rights reserved.
Copyright (C) 2014-2015, Itseez Inc., all rights reserved.
Copyright (C) 2020-2021 Miles Lucas, Tim Holy, and Mark Kittisopikul, all rights reserved
Third party copyrights are property of their respective owners.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

  * Redistribution's of source code must retain the above copyright notice,
    this list of conditions and the following disclaimer.

  * Redistribution's in binary form must reproduce the above copyright notice,
    this list of conditions and the following disclaimer in the documentation
    and/or other materials provided with the distribution.
  * The name of the copyright holders may not be used to endorse or promote products
    derived from this software without specific prior written permission.

This software is provided by the copyright holders and contributors "as is" and
any express or implied warranties, including, but not limited to, the implied
warranties of merchantability and fitness for a particular purpose are disclaimed.
In no event shall the Intel Corporation or contributors be liable for any direct,
indirect, incidental, special, exemplary, or consequential damages
including, but not limited to, procurement of substitute goods or services;
loss of use, data, or profits; or business interruption) however caused
and on any theory of liability, whether in contract, strict liability,
or tort (including negligence or otherwise) arising in any way out of
the use of this software, even if advised of the possibility of such damage.
=#
using StaticArrays

export Lanczos4OpenCV

"""
    Lanczos4OpenCV()

Alternative implementation of Lanczos resampling using algorithm `lanczos4` function of OpenCV:
https://github.com/opencv/opencv/blob/de15636724967faf62c2d1bce26f4335e4b359e5/modules/imgproc/src/resize.cpp#L917-L946
"""
struct Lanczos4OpenCV <: AbstractLanczos end

degree(::Lanczos4OpenCV) = 4

value_weights(::Lanczos4OpenCV, δx) = _lanczos4_opencv(δx)

# s45 = sqrt(2)/2
const s45 = 0.70710678118654752440084436210485

# l4_2d_cs is a lookup table that could be generated by
# x = (0:7)*45*5
# l4_2d_cs = [cosd.(x) sind.(x)]'
const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45]

function _lanczos4_opencv(δx)
    p_4 = π / 4
    y0 = -(δx + 3) * p_4
    s0, c0 = sincos(y0)
    cs = ntuple(8) do i
        y = (δx + 4 - i) * p_4
        # Improve precision of Lanczos OpenCV4 #451, avoid NaN
        if iszero(y)
            y = eps(oneunit(y))/8
        end
        # Numerator is the sin subtraction identity
        # It is equivalent to the following
        # f(δx,i) = sin( π/4*( 5*(i-1)-δx-3 ) )
        (l4_2d_cs[i, 1] * s0 + l4_2d_cs[i, 2] * c0) / y^2
    end
    sum_cs = sum(cs)
    normed_cs = ntuple(i -> cs[i] / sum_cs, Val(8))
    return normed_cs
end


text:
using Test
using LinearAlgebra
include("module-3-exercise-3.jl")

@testset "$(get_module_3_excersise_3_author_name()): problem 2 tests" begin
    @test get_module_3_excersise_3_author_name() != "YOUR NAME HERE"
    @test problem_3([1 2 3;2 1 4;3 4 1]) == true
    @test problem_3([1 2 59;2 1 4;3 4 1]) == false
    @test problem_3([1 2 3 4; 2 3 4 5;3 4 5 6; 4 5 6 7]) == true
    @test problem_3([1 2 3 34; 2 3 4 5;3 4 5 6; 4 5 6 7]) == false
end


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

"""
    share_pce!(d::Dict)
`pce`: Regional shares of final consumption
"""
function share_pce!(d::Dict)
    println("  Calculating regional shares of final consumption")
    d[:pce] /= transform_over(d[:pce], :r)
    verify_over(d[:pce],:r) !== true && @error("PCE shares don't sum to 1.")
    return d[:pce]
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
import Base

"""
    define_from(::Type{Cascade{Data}}, df)
    define_from(::Type{Plot{Data}}, cascade; kwargs...)
    define_from(::Type{T}, cascade::Cascade{Data}; kwargs...) where T<:Axis

# Keyword Arguments
- `ylabel::String`, y-axis label
-
"""
function define_from(::Type{Cascade{Data}}, df::DataFrames.DataFrame;
    permutation = missing,
    correlate = true,
    permute = true,
    kwargs...,
)
    rows = size(df,1)
    maxstep = min(rows-2,length(COLORCYCLE))

    # Check for numbers in here lager than the colorcycle length?
    permutation = coalesce(permutation, collect(1:maxstep))

    idx = [1; sort(permutation).+1; rows]
    data = define_from(Vector{Data}, df[idx,:]; kwargs...)

    cascade = Cascade(
        start = first(data),
        stop = last(data),
        steps = data[2:end-1],
        correlation = random_rotation(length(permutation);
            nrandom=correlate,
            maxdim=maxstep,
            permutation=permutation,
            kwargs...,
        ),
        permutation = permutation,
        iscorrelated = false,
        ispermuted = false,
    )

    correlate!==false && correlate!(cascade)
    permute && (cascade = permute!(cascade))

    return cascade
end

function define_from(::Type{Data}, sdf::DataFrames.SubDataFrame;
    label,
    sublabel,
    value=VALUE_COL,
    order=[],
    kwargs...,
)
    return Data( ;
        label = sdf[1,label],
        sublabel = sdf[1,sublabel],
        value = sdf[:,value],
        order = isempty(order) ? collect(1:size(sdf,1)) : order,
    )
end

function define_from(::Type{Vector{Data}}, df::DataFrames.DataFrame; kwargs...)
    gdf = fuzzify(df; kwargs...)
    data = [define_from(Data, sdf; kwargs...) for sdf in gdf]

    # Sort "start" samples by magnitude. Ensure consistent ordering in all subsequent steps.
    iiorder = sortperm(get_value(data[1]))
    [set_order!(data[ii], iiorder) for ii in 1:length(data)]

    return data
end

function define_from(::Type{Plot{T}}, cascade::Cascade{Data}; kwargs...) where T<:Geometry
    cascade_geometry = set_geometry(cascade, T; kwargs...)
    legend = _define_legend(cascade_geometry; kwargs...)

    return Plot( ;
        cascade = cascade_geometry,
        xaxis = define_from(XAxis, cascade; kwargs...),
        yaxis = define_from(YAxis, cascade; kwargs...),
        title = _define_title(cascade; kwargs...),
        path = _define_path(cascade, T; kwargs...),
        # Add other annotations to the legend.
        legend = _push!(legend, cascade, Horizontal, Statistics.mean; kwargs...),
    )
end

function define_from(::Type{XAxis}, cascade::Cascade{Data}; kwargs...)
    label_type = Vector{Label{Vector{String}}}
    pad = SEP/2

    return XAxis( ;
        ticklabels = _define_from(label_type, cascade, XAxis, :label;
            yshift = pad,
            kwargs...,
        ),
        ticksublabels = _define_from(label_type, cascade, XAxis, :sublabel;
            yshift = pad+FONTSIZE,
            scale = 0.8,
            kwargs...,
        ),
        frame = _define_from(Line, (Luxor.Point(0,HEIGHT), Luxor.Point(WIDTH+2*SEP,HEIGHT)))
    )
end

function define_from(::Type{YAxis}, cascade::Cascade{Data};
    ylabel,
    x = 0,
    y = HEIGHT,
    kwargs...,
)
    label_type = Label{String}

    return YAxis( ;
        ticks = _define_ticks(cascade, YAxis; kwargs...),
        label = _define_from(label_type, ylabel, Luxor.Point(-(LEFT_BORDER-SEP/2), HEIGHT/2);
            angle = -pi/2,
            valign = :top,
        ),
        ticklabels = _define_from(Vector{label_type}, cascade, YAxis;
            halign = :right,
            xshift = -SEP,
            scale = 0.9,
            kwargs...,
        ),
        frame = _define_from(Arrow, (Luxor.Point(0,HEIGHT), Luxor.Point(0,0))),
    )
end

function define_from(::Type{T}, x; kwargs...) where T <: Luxor.Colorant
    return scale_hsv(parse(T, _define_colorant(x)); kwargs...)
end

"""
    _define_from(Shape, Color, pos, sgn; kwargs...)
    _define_from(::Type{T}, sign, pos; kwargs...) where T<:Vector{Blending}
    _define_from(::Type{T}, sign; kwargs...) where T<:Union{Vector{Coloring},Coloring}

# Keyword Arguments:
- `style`

    _define_from(::Type{Vector{L}}, ::Type{A}, cascade; kwargs...) where {L<:Label, A<:Axis}
    _define_from(::Type{Vector{L}}, ::Type{XAxis}, cascade, fun; kwargs...) where L<:Label
These methods define LABELS for the input `Axis` subtype (XAxis or YAxis).

# Arguments
- `cascade::Cascade`
- `field::Symbol`, if defining ticks for the XAxis, include the Data field name that will
serve as the x-axis tick label or sublabel. Allowed options: `:label`, `:sublabel`.

# Returns
- `ticklabels::Vector{L}` of tick LABELS
- `ticksublabels::Vector{L}` (returned if `A=XAxis` and `fun` unspecified), of tick SUBLABELS

    _define_from(Handle, cascade, fun, args...)

    _define_from(Annotation)
"""
function _define_from(Shape, Color, pos::Vector{T1}, sgn::Vector{T2};
    kwargs...,
) where {T1<:AbstractArray, T2<:AbstractArray}
    return _define_from.(Shape, Color, pos, sgn; kwargs...)
end

function _define_from(Shape, Color, pos::Vector{T}, sgn;
    style,
    kwargs...,
) where T <: Luxor.Point
    return Shape(pos, _define_from(Color, sgn, pos; kwargs...), style)
end

function _define_from(Shape, Color, pos::Vector{T}, sgn;
    style,
    kwargs...,
) where T <: Tuple
    return Shape.(pos, _define_from(Vector{Color}, sgn, pos; kwargs...), style)
end

function _define_from(::Type{Coloring}, sign, args...;
    alpha = 0.8,
    kwargs...,
)
    hue = define_from(Luxor.RGB, sign; kwargs...)
    alpha = _define_alpha(alpha; kwargs...)
    return Coloring(hue, alpha)
end

function _define_from(::Type{Vector{Coloring}}, args...; kwargs...)
    return _define_from.(Coloring, args...; kwargs...)
end

function _define_from(::Type{Vector{Blending}}, sgn, pos; kwargs...)
    hue = _define_gradient.(sgn)

    # Calculate direction.
    xmid = mid(pos; dims=1)

    ii = Dict(k => sgn.==k for k in [-1,1])
    ymax = Dict(k => maximum(pos[v]; dims=2) for (k,v) in ii)
    ymin = Dict(k => minimum(pos[v]; dims=2) for (k,v) in ii)

    begin y1 = fill(0.0, size(sgn)); y2 = fill(0.0, size(sgn)) end
    begin y1[ii[-1]] .= ymin[-1]; y2[ii[-1]] .= ymax[-1] end
    begin y1[ii[+1]] .= ymax[+1]; y2[ii[+1]] .= ymin[+1] end

    direction = tuple.(Luxor.Point.(xmid,y1), Luxor.Point.(xmid,y2))
    return Blending.(direction, hue)
end

## Define Label; lists of Labels

function _define_from(::Type{Vector{T}}, cascade::Cascade{Data}, args...; kwargs...) where T <: Label
    txt = _define_text(cascade, args...; kwargs...)
    pos = _define_position(cascade, args...; kwargs...)
    return _define_from(Vector{T}, txt, pos; kwargs...)
end

function _define_from(::Type{T}, text, pos::Luxor.Point;
    scale = 1,
    halign = :center,
    valign = :middle,
    leading = 1,
    angle = 0,
    kwargs...,
) where T <: Label
    return Label( ;
        text = text,
        scale = scale,
        position = pos,
        halign = halign,
        valign = valign,
        angle = angle,
        leading = leading,
    )
end

function _define_from(::Type{T}, text::String, pos::Luxor.Point, width::Float64;
    scale = 1,
    kwargs...,
) where T <: Label
    return _define_from(T, wrap_to(text, width; scale), pos; scale=scale, kwargs...)
end

function _define_from(::Type{Vector{Label{T}}}, txt::AbstractArray, pos::AbstractArray;
    kwargs...,
) where T <: Vector{String}
    N = length(pos)
    wid = width(N; space=2)
    return [_define_from(Label, txt[ii], pos[ii], wid; kwargs...) for ii in 1:N]
end

function _define_from(::Type{Vector{Label{T}}}, text::AbstractArray, pos::AbstractArray;
    kwargs...,
) where T <: Any
    return [_define_from(Label{T}, text[ii], pos[ii]; kwargs...) for ii in 1:length(pos)]
end

function _define_from(::Type{T}, x; kwargs...) where T<:Label
    # If no point is given, initialize label position to (0,0)
    return _define_from(Label, x, Luxor.Point(0,0); kwargs...)
end

## Define Arrow, Line (for plot axes)

function _define_from(::Type{T}, x::Tuple{Luxor.Point,Luxor.Point};
    hue = "black",
    alpha = 1.0,
    style = :stroke,
    kwargs...,
) where T <: Union{Arrow,Line}
    return T(x, _define_from(Coloring, hue; alpha=alpha), style)
end

## Define Handle

function _define_from(::Type{Handle}, cascade::Cascade{T}, str::String;
    scale = 0.8,
    kwargs...,
) where T <: Geometry
    shape = get_shape(cascade)
    label = _define_from(Label, uppercase(str); halign=:left, scale=scale)
    return set_position!(Handle(shape, label); scale=scale, kwargs...)
end

function _define_from(::Type{Handle}, cascade, fun::Function, args...; kwargs...)
    return _define_from(Handle, cascade, string(fun); kwargs...)
end

function _define_from(::Type{Handle}, cascade, str::String, hue; kwargs...)
    h = _define_from(Handle, cascade, str; kwargs...)
    return set_hue!(h, hue)
end

function _define_from(::Type{Handle}, cascade; colorcycle, kwargs...)
    #= !!!!
    Should probably set a lower limit on alpha for large numbers of samples.
    This gets really light and hard to see.
    =#
    return if colorcycle
        [_define_from(Handle, cascade, "SAMPLE", "black"; idx=1, kwargs...)]
    else
        [_define_from(Handle, cascade, str, hue; idx=idx, kwargs...)
            for (idx, str, hue) in zip(1:2, ["GAIN","LOSS"], [HEX_GAIN,HEX_LOSS])]
    end
end

## Define Annotation

function _define_from(::Type{Annotation}, cascade::Cascade{Data}, Geometry::DataType,
    args...;
    scale = 0.8,
    kwargs...,
)
    return Annotation(
        label = _define_label(cascade, Geometry, args...; scale=scale, kwargs...),
        cascade = set_geometry(cascade, Geometry, args...; alpha=1.0, style=:stroke, kwargs...),
    )
end

## Define values for DataTypes PROPERTIES

"""
    _define_alpha(N::Int; kwargs...)
This function calculates an alpha for `N` samples to scale transparency for overlays:

```math
\\alpha = \\min\\left\\lbrace\\dfrac{w}{f(N)},\\, 0\\right\\rbrace
```

# Keyword Arguments
- `factor::Real=0.25`, scaling weighting ``w``
- `fun::Function=log`, scaling function ``\\ln``
"""
function _define_alpha(N::Int; factor::Real=0.25, fun::Function=log, kwargs...)
    return min(factor/fun(N) ,1.0)
end

_define_alpha(x; kwargs...) = x

"""
    _define_colorant(idx::Int)
This method returns the color defined at index `idx` of `COLORCYCLE`

    _define_colorant(sgn::Float64)
This method returns
- **Red**, given a negative value or
- **Blue**, given a positive value.

    _define_colorant(lst::AbstractVector)
This method returns an average of the colors defined in the `lst`
"""
function _define_colorant(lst::AbstractVector)
    rgb = _define_colorant.(lst)
    rgb_average = [Statistics.mean(getproperty.(rgb, f)) for f in fieldnames(Luxor.RGB)]
    return Luxor.Colors.RGB(rgb_average...)
end

_define_colorant(idx::Int) = COLORCYCLE[idx]
_define_colorant(sgn::Float64) = sign(sgn)>0 ? HEX_GAIN : HEX_LOSS
_define_colorant(x) = parse(Luxor.Colorant, x)

"""
    _define_gradient(x; kwargs...)
This method defines the starting and stopping colors for a color gradient, such that the
starting color is defined using [`Waterfall._define_colorant`](@ref) and the stopping color
is 50% "lighter" than this color.

# Returns
- `x::Tuple{Luxor.RGB,Luxor.RGB}`, color gradient from the starting to stopping pos of
    each waterfall in the cascade.
"""
function _define_gradient(x; kwargs...)
    lightness = 0.5
    h1 = define_from(Luxor.RGB, x)
    h2 = define_from(Luxor.RGB, x; lightness=lightness)
    return tuple(h1, h2)
end

function _define_gradient(x1, x2; kwargs...)
    return tuple(define_from(Luxor.RGB, x1; kwargs...), define_from(Luxor.RGB, x2; kwargs...))
end

"""
    _define_label(cascade, Geometry::DataType, args...; kwargs...)
"""
function _define_label(cascade::Cascade{Data}, Geometry::DataType, args...; kwargs...)
    txt = _define_text(cascade, args...; kwargs...)
    pos = _define_position(cascade, Geometry; kwargs...)
    return _define_from(Vector{Label{String}}, txt, pos; valign=:bottom, kwargs...)
end

"""
    _define_legend(cascade; kwargs)
"""
function _define_legend(cascade::Cascade{T}; kwargs...) where T <: Geometry
    legend = _define_from(Handle, cascade; kwargs...) .=> missing
    return legend |> Vector{Pair{Handle,Any}}
end

"""
    _push!(legend, cascade, Geometry, args...; kwargs...)

# Arguments
- `legend`
- `cascade::Cascade{Data}`
- `fun::Function`
"""
function _push!(
    legend::Vector{Pair{Handle,T}},
    cascade::Cascade{Data},
    Geometry::DataType,
    args...;
    kwargs...,
) where T<:Any
    a =  _define_from(Annotation, cascade, Geometry, args...; kwargs...)
    h = _define_from(Handle, a.cascade, args...; idx=length(legend)+1, kwargs...)
    push!(legend, h => a)
    return legend
end

"""
    _define_path(cascade::Cascade, Geometry; kwargs...)
This method defines the path to which to write the file:
    Geometry/Geometry_n<nsample>_colorcycle<0/1>_corr<0/1>_<permutation>.png
"""
function _define_path(x::Cascade, Geometry;
    colorcycle,
    ext = ".png",
    figdir = "fig",
    maxsample = 50,
    kwargs...,
)
    separator = "_"
    Geometry = lowercase(string(Geometry))

    # Define the path name and create the directory if it doesn't already exist.
    path = joinpath(WATERFALL_DIR, figdir, Geometry)

    if !isdir(path)
        @info("Creating directory: $path")
        mkpath(path)
    end

    str = joinpath(path, join([
        Geometry,
        "n" * _define_path(length(x); maxchar=get_order(maxsample)+1), # number of samples
        _define_path(get_label.(x.steps); kwargs...),                  # labels, in order
        "corr"       * _define_path(x.iscorrelated),                   # is it correlated?
        "colorcycle" * _define_path(colorcycle),                       # colorscheme?
        ], separator,
    ) * ext)

    println("Writing plot: $str")
    return str
end

function _define_path(lst::Vector{String}; maxchar=missing, kwargs...)
    maxchar = coalesce(maxchar, maximum(length.(lst)))
    lst = _define_path.(tryinteger(lst); maxchar=maxchar)
    return join(lst, maxchar>1 ? "-" : "")
end

_define_path(x::Int; maxchar=1) = string(Printf.@sprintf("%010.0f", x)[end-(maxchar-1):end])
_define_path(x::Bool; kwargs...) = string(convert(Int, x))
_define_path(x::String; kwargs...) = x

"""
    _define_position(cascade, ::Type{T}; kwargs...) where T <: Axis
This method defines the TICK POSITIONS for the input Axis subtype (XAxis or YAxis)
`XAxis` tick labels are centered below the associated cascade waterfall.

    _define_position(cascade, Geometry::DataType; kwargs...)
This method defines a point above the top of each cascade waterfall, so that an annotation
label can be bottom/center-aligned at these points.

# Arguments
- `cascade`
- ``

# Keyword Arguments
- `x=0`, x-coordinate of y-axis.
- `y=HEIGHT`, y-coordinate of x-axis.
- `xshift`
- `yshift`
"""
function _define_position(cascade, ::Type{YAxis};
    x = 0,
    xshift = 0,
    kwargs...,
)
    y = collect_ticks(cascade; kwargs...)
    y = scale_y(y; vlim(cascade; kwargs...)...)
    return Luxor.Point.(x+xshift, y)
end

function _define_position(cascade, ::Type{XAxis}, args...;
    y = HEIGHT,
    yshift = 0,
    kwargs...,
)
    x = scale_x( ; nrow=length(collect_data(cascade)))
    return Luxor.Point.(x, y+yshift)
end

function _define_position(cascade, Geometry; kwargs...)
    yshift = -SEP/2
    x = scale_x( ; nrow=length(collect_data(cascade)))
    pos = vectorize.(scale_for(cascade, Geometry; kwargs...))
    y = minimum.(pos; dims=2)
    return Luxor.Point.(x, y.+yshift)
end

"""
    _define_text(x::Float64; kwargs...)
This method formats `x` as a string.

# Keyword Arguments
- `digits=2`, fractional digits to include
- `sgn=true`, indicates whether to show a leading `+`

    _define_text(cascade, fun::Function, args...; kwargs...)
This method calculates a value
"""
function _define_text(x::Float64; sgn=true, digits=2)
    x = round(x; digits=digits)
    sgn = sgn ? sgn : Base.sign(x)<0

    len = length(string(abs(convert(Int, round(x))))) + 1 + sgn
    str = sgn ? @Printf.sprintf("%+2.10f", x) : @Printf.sprintf("%2.10f", x)

    return string(str[1:len+digits])
end

function _define_text(cascade, fun::Function, args...; kwargs...)
    v = collect_value(cascade)
    v = calculate(v, fun, args...)

    digits = abs(minimum(get_order.(v))) + 1

    str = _define_text.(v; digits=digits)
    [str[ii] = str[ii][2:end] for ii in [1,size(str,1)]]

    return str
end

function _define_text(cascade, ::Type{YAxis}; scale=0.9, kwargs...)
    v = collect_ticks(cascade; kwargs...)
    digits = abs(minimum(get_order.(v)))
    return _define_text.(v; digits=digits, sgn=false)
end

function _define_text(cascade, ::Type{XAxis}, field::Symbol; kwargs...)
    return getfield.(collect_data(cascade), field)
end

_define_text(x) = x

"""
"""
function _define_ticks(cascade, ::Type{T};
    x = 0,
    y = HEIGHT,
    len = SEP,
    kwargs...,
) where T <: Axis

    p1 = _define_position(cascade, T; x=x-0.5*len, y=y-0.5*len, kwargs...)
    p2 = _define_position(cascade, T; x=x+0.5*len, y=y+0.5*len, kwargs...)

    return _define_from.(Line, tuple.(p1,p2))
end

"""
    _define_title(cascade::Cascade{Data}; kwargs...)
This method formats a title to describe, where applicable:
1. The number of samples,
2. Range of correlation coefficients (if a correlation was applied), and
3. The distribution function (normal or uniform) used to create samples
    (if more than one sample).
"""
function _define_title(cascade::Cascade{Data}; nsample,
    titlescale = 1.1,
    titleleading = 1.2,
    kwargs...,
)
    # Select which rows of the title to show:
    # (1) Always show the number of samples,
    # (2) Show correlation range if one is applied,
    # (3) Show sample distribution if multiple samples are given.
    idx = [true, cascade.iscorrelated, nsample>1]

    str = [
        "$nsample SAMPLE" * (nsample>1 ? "S" : ""),
        _title_correlation( ; kwargs...),
        _title_distribution( ; kwargs...)
    ][idx]

    return Label( ;
        text = str,
        scale = titlescale,
        position = Luxor.Point(WIDTH/2, -TOP_BORDER),
        halign = :center,
        valign = :middle,
        angle = 0,
        leading = titleleading,
    )
end

"Format title string for normal distribution function"
_title_normal(x) = "f(x; $(x[1]) < s < $(x[2]))"

"Format title string for uniform distribution function"
function _title_uniform(x; abs::Bool=false)
    abs = abs ? "|" : ""
    return "f($(abs)x$(abs); a>$(x[1]), b<$(x[2]))"
end

"Format title string for SAMPLE distribution"
function _title_distribution( ; distribution, fuzziness, kwargs...)
    str = uppercase("$distribution sample distribution: ")

    return if distribution==:normal; str * _title_normal(fuzziness)
    elseif distribution==:uniform;   str * _title_uniform(fuzziness)
    else; ""
    end
end

"Format title string for correlation coefficient"
function _title_correlation( ; interactivity, kwargs...)
    return "CORRELATION COEFFICIENT: " * _title_uniform(interactivity; abs=true)
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using ACE.SphericalHarmonics: SHBasis, index_y
using StaticArrays
using JuLIP: AbstractCalculator, Atoms, JVec
using JuLIP.Potentials: SitePotential, SZList, ZList
using NeighbourLists: neigs

import JuLIP, JuLIP.MLIPs

export PIPotential

# from pibasis:
# StandardEvaluator, DAGEvaluator, graphevaluator, standardevaluator

"""
`struct PIPotential` : specifies a PIPotential, which is basically defined
through a PIBasis and its coefficients
"""
mutable struct PIPotential{T, NZ, TPI, TEV} <: SitePotential
   pibasis::TPI
   coeffs::NTuple{NZ, Vector{T}}
   dags::NTuple{NZ, CorrEvalGraph{T, Int}}
   evaluator::TEV
end

cutoff(V::PIPotential) = cutoff(V.pibasis)

==(V1::PIPotential, V2::PIPotential) =
      (V1.pibasis == V2.pibasis) && (V1.coeffs == V2.coeffs)

# TODO: this doesn't feel right ... should be real(T)?
fltype(::PIPotential{T}) where {T} = real(T)

zlist(V::PIPotential) = zlist(V.pibasis)

graphevaluator(V::PIPotential) =
   PIPotential(V.pibasis, V.coeffs, V.dags, DAGEvaluator())

standardevaluator(V::PIPotential) =
   PIPotential(V.pibasis, V.coeffs, V.dags, StandardEvaluator())

maxorder(V::PIPotential) = maxorder(V.pibasis)

# ------------------------------------------------------------
#   Initialisation code
# ------------------------------------------------------------

combine(basis::PIBasis, coeffs::AbstractVector) =
      PIPotential(basis, identity.(collect(coeffs)))

# assemble from basis with global coeff vector
function PIPotential(basis::PIBasis, coeffs::Vector{<: Number})
   coeffs_t = ntuple(iz0 -> coeffs[basis.inner[iz0].AAindices], numz(basis))
   return PIPotential(basis, coeffs_t)
end

# assemble from basis with coeff vectors separated into individual species
function PIPotential(basis::PIBasis, coeffs_t::Tuple)
   dags = ntuple(iz0 -> _getdagfrombasis(basis.inner[iz0], coeffs_t[iz0]), numz(basis))
   return PIPotential(basis, coeffs_t, dags, DAGEvaluator())
end

function _getdagfrombasis(inner, c)
   nodes = inner.dag.nodes
   vals = inner.dag.vals
   dag = DAG.CorrEvalGraph(nodes, zeros(eltype(c), length(nodes)),
                           inner.dag.num1, inner.dag.numstore)
   for (i, idx) in enumerate(vals)
      if idx != 0  # idx = 0 means this is an auxiliary basis with 0 coefficient!
         dag.vals[i] = c[idx]
      end
   end
   return dag
end

# ------------------------------------------------------------
#   experimental sparsification code

function deletezeros(V::PIPotential)
   idx0 = 0
   innernew = []
   coeffs = []
   for iz0 = 1:length(V.coeffs)
      # zeros
      Iz = findall(V.coeffs[iz0] .== 0)
      # non-zeros
      Inz = findall(V.coeffs[iz0] .!= 0)
      # inverse mapping from old to new (local) indices
      invInz = Dict{Int, Int}()
      for (i, val) in enumerate(Inz)
         invInz[val] = i
      end
      # extract the inner basis to start rebuilding it
      inner = V.pibasis.inner[iz0]     # old
      iAA2iA = inner.iAA2iA[Inz, :]    # new
      # construct the new b2iAA mapping
      b2iAA = Dict{ACE.PIBasisFcn, Int}()   # new
      for (b, iAA) in inner.b2iAA
         if haskey(invInz, iAA)
            b2iAA[b] = invInz[iAA]
         end
      end
      # create the new inner basis
      push!(innernew,
            ACE.InnerPIBasis( inner.orders[Inz],
                              iAA2iA,
                              b2iAA,
                              inner.b2iA,
                              (idx0+1):(idx0+length(Inz)),
                              inner.z0,
                              ACE.DAG.CorrEvalGraph{Int, Int}() ) )
      # and the new coefficients
      push!(coeffs, V.coeffs[iz0][Inz])
   end
   pibasis = PIBasis( V.pibasis.basis1p, V.pibasis.zlist,
                      tuple( innernew... ), V.pibasis.evaluator )
   return standardevaluator(PIPotential( pibasis, tuple(coeffs...) ))
end

# ------------------------------------------------------------
#   FIO code
# ------------------------------------------------------------

write_dict(V::PIPotential) = Dict(
      "__id__" => "ACE_PIPotential",
     "pibasis" => write_dict(V.pibasis),
      "coeffs" => [ write_dict.(V.coeffs)... ] )

#    if ntests > 0
#       tests = ACE.Random.rand_nhd(Nat, J::ScalarBasis, species = :X)
#       Pr = V.pibasis.basis1p.J
#       for ntest = 1:5
#
#          r = ACE.rand_radial(V.pibasis.basis1p.J)
#          Pr = evaluate(B3.J, r)
#          push!(rtests, Dict("r" => r, "Pr" => Pr))
#       end
#
#    end
#    return D
# end

read_dict(::Val{:SHIPs_PIPotential}, D::Dict; tests = true) =
   read_dict(Val{:ACE_PIPotential}(), D)

read_dict(::Val{:ACE_PIPotential}, D::Dict; tests = true) =
   PIPotential( read_dict(D["pibasis"]),
                tuple( read_dict.( D["coeffs"] )... ) )

# function test_compat(V::PIPotential, rtests, tests)
#
# end

# ------------------------------------------------------------
#   Dispatching the Evaluation codes
# ------------------------------------------------------------

alloc_temp(V::PIPotential, maxN::Integer) =  alloc_temp(V.evaluator, V, maxN)
alloc_temp_d(V::PIPotential, maxN::Integer) =  alloc_temp_d(V.evaluator, V, maxN)
evaluate!(tmp, V::PIPotential, args...) = evaluate!(tmp, V, V.evaluator, args...)
evaluate_d!(dEs, tmp, V::PIPotential, args...) = evaluate_d!(dEs, tmp, V, V.evaluator, args...)

# ------------------------------------------------------------
#   Standard Evaluation code
# ------------------------------------------------------------

# TODO: generalise the R, Z, allocation
alloc_temp(::StandardEvaluator, V::PIPotential{T}, maxN::Integer) where {T} =
   (
      R = zeros(JVec{real(T)}, maxN),
      Z = zeros(AtomicNumber, maxN),
      tmp_pibasis = alloc_temp(V.pibasis, maxN),
   )

# compute one site energy
function evaluate!(tmp, V::PIPotential, ::StandardEvaluator,
                   Rs::AbstractVector{JVec{T}},
                   Zs::AbstractVector{<:AtomicNumber},
                   z0::AtomicNumber) where {T}
   iz0 = z2i(V, z0)
   A = evaluate!(tmp.tmp_pibasis.A, tmp.tmp_pibasis.tmp_basis1p,
                 V.pibasis.basis1p, Rs, Zs, z0)
   inner = V.pibasis.inner[iz0]
   c = V.coeffs[iz0]
   Es = zero(T)
   @inbounds for iAA = 1:length(inner)
      Esi = one(Complex{T}) * c[iAA]    # TODO: OW - NASTY!!!
      for α = 1:inner.orders[iAA]
         Esi *= A[inner.iAA2iA[iAA, α]]
      end
      Es += real(Esi)
   end
   return Es
end

# TODO: generalise the R, Z, allocation
alloc_temp_d(::StandardEvaluator, V::PIPotential{T}, N::Integer) where {T} =
      (
      dAco = zeros(fltype(V.pibasis),
                   maximum(length(V.pibasis.basis1p, iz) for iz=1:numz(V))),
       tmpd_pibasis = alloc_temp_d(V.pibasis, N),
       dV = zeros(JVec{real(T)}, N),
        R = zeros(JVec{real(T)}, N),
        Z = zeros(AtomicNumber, N),
        dAAt = zero(MVector{maxorder(V), fltype(V.pibasis.basis1p)})
      )

# compute one site energy
function evaluate_d!(dEs, tmpd, V::PIPotential, ::StandardEvaluator,
                     Rs::AbstractVector{<: JVec{T}},
                     Zs::AbstractVector{AtomicNumber},
                     z0::AtomicNumber
                     ) where {T}
   iz0 = z2i(V, z0)
   basis1p = V.pibasis.basis1p
   tmpd_1p = tmpd.tmpd_pibasis.tmpd_basis1p
   Araw = tmpd.tmpd_pibasis.A

   # stage 1: precompute all the A values
   A = evaluate!(Araw, tmpd_1p, basis1p, Rs, Zs, z0)

   # stage 2: compute the coefficients for the ∇A_{klm} = ∇ϕ_{klm}
   dAco = tmpd.dAco
   c = V.coeffs[iz0]
   inner = V.pibasis.inner[iz0]
   dAAt = tmpd.dAAt
   fill!(dAco, 0)
   for iAA = 1:length(inner)
      c_ = c[iAA]
      ord = inner.orders[iAA]
      Afwd = one(eltype(A))
      dAAt[1] = 1
      for α = 1:ord-1
         Afwd *= A[inner.iAA2iA[iAA, α]]
         dAAt[α+1] = Afwd
      end
      Abwd = one(eltype(A))
      for α = ord:-1:2
         Abwd *= A[inner.iAA2iA[iAA, α]]
         dAAt[α-1] *= Abwd
      end
      for α = 1:ord
         iAα = inner.iAA2iA[iAA, α]
         dAco[iAα] += c_ * dAAt[α]
      end

      # for α = 1:inner.orders[iAA]
      #    CxA_α = c[iAA]
      #    for β = 1:inner.orders[iAA]
      #       if β != α
      #          CxA_α *= A[inner.iAA2iA[iAA, β]]
      #       end
      #    end
      #    iAα = inner.iAA2iA[iAA, α]
      #    dAco[iAα] += CxA_α
      # end
   end

   # stage 3: get the gradients
   fill!(dEs, zero(JVec{T}))
   dAraw = tmpd.tmpd_pibasis.dA
   for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
      evaluate_d!(Araw, dAraw, tmpd_1p, basis1p, R, Z, z0)
      iz = z2i(basis1p, Z)
      zinds = basis1p.Aindices[iz, iz0]
      for iA = 1:length(basis1p, iz, iz0)
         dEs[iR] += real(dAco[zinds[iA]] * dAraw[zinds[iA]])
      end
   end

   return dEs
end

# ------- Evaluation code

_maxstore(V::PIPotential) = maximum( dag.numstore for dag in V.dags )

alloc_temp(::DAGEvaluator, V::PIPotential{T}, maxN::Integer) where {T} =
   (
   R = zeros(JVec{real(T)}, maxN),
   Z = zeros(AtomicNumber, maxN),
   tmp_basis1p = alloc_temp(V.pibasis.basis1p),
   AA = zeros(fltype(V.pibasis.basis1p), _maxstore(V)),
   A = alloc_B(V.pibasis.basis1p)
   )

function evaluate!(tmp, V::PIPotential, ::DAGEvaluator, Rs, Zs, z0)
   AAdag = tmp.AA
   A = tmp.A
   iz0 = z2i(V, z0)
   dag = V.dags[iz0]
   nodes = dag.nodes
   vals = dag.vals
   @assert length(A) >= dag.num1
   @assert length(AAdag) >= dag.numstore

   evaluate!(A, tmp.tmp_basis1p, V.pibasis.basis1p, Rs, Zs, z0)

   Es = zero(fltype(V))
   @inbounds for i = 1:dag.num1
      AAdag[i] = a = A[i]
      Es = muladd(vals[i], real(a), Es)
   end

   @inbounds for i = (dag.num1+1):dag.numstore
      n1, n2 = nodes[i]
      AAdag[i] = a = AAdag[n1] * AAdag[n2]
      Es = muladd(vals[i], real(a), Es)
   end

   @inbounds for i = (dag.numstore+1):length(dag)
      n1, n2 = nodes[i]
      a = AAdag[n1] * AAdag[n2]
      Es = muladd(vals[i], real(a), Es)
   end

   return Es
end

alloc_temp_d(::DAGEvaluator, V::PIPotential{T}, maxN::Integer) where {T} =
   (
   R = zeros(JVec{real(T)}, maxN),
   Z = zeros(AtomicNumber, maxN),
    dV = zeros(JVec{real(T)}, maxN),
   tmpd_basis1p = alloc_temp_d(V.pibasis.basis1p),
   AA = zeros(fltype(V.pibasis.basis1p), _maxstore(V)),
   B = zeros(fltype(V.pibasis.basis1p), _maxstore(V)),
   A = alloc_B(V.pibasis.basis1p),
   dA = alloc_dB(V.pibasis.basis1p)
    )

function evaluate_d!(dEs, tmpd, V::PIPotential{T}, ::DAGEvaluator, Rs, Zs, z0) where {T}
   iz0 = z2i(V, z0)
   AA = tmpd.AA
   B = tmpd.B
   tmpd_basis1p = tmpd.tmpd_basis1p
   basis1p = V.pibasis.basis1p
   dag = V.dags[iz0]
   nodes = dag.nodes
   coeffs = dag.vals

   # we start from the representation
   #     V = sum B[i] AA[i]
   # i.e. this vector represents the constributions c[i] ∂AA[i]
   copy!(B, coeffs)

   # FORWARD PASS
   # ------------
   # evaluate the 1-particle basis
   # this puts the first `dag.num1` 1-b (trivial) correlations into the
   # storage array, and from these we can build the rest
   evaluate!(AA, tmpd_basis1p, basis1p, Rs, Zs, z0)

   # Stage 2 of evaluate!
   # go through the dag and store the intermediate results we need
   @inbounds @fastmath for i = (dag.num1+1):dag.numstore
      n1, n2 = nodes[i]
      AA[i] = muladd(AA[n1], AA[n2], AA[i])
   end

   # BACKWARD PASS
   # --------------
   # fill the B array -> coefficients of the derivatives
   #  AA_i = AA_{n1} * AA_{n2}
   #  ∂AA_i = AA_{n1} * ∂AA_{n2} + AA_{n1} * AA_{n2}
   #  c_{n1} * ∂AA_{n1} <- (c_{n1} + c_i AA_{n2}) ∂AA_{n1}
   @inbounds @fastmath for i = length(dag):-1:(dag.numstore+1)
      c = coeffs[i]
      n1, n2 = nodes[i]
      B[n1] = muladd(c, AA[n2], B[n1])
      B[n2] = muladd(c, AA[n1], B[n2])
   end
   # in stage 2 c = C[i] is replaced with b = B[i]
   @inbounds @fastmath for i = dag.numstore:-1:(dag.num1+1)
      n1, n2 = nodes[i]
      b = B[i]
      B[n1] = muladd(b, AA[n2], B[n1])
      B[n2] = muladd(b, AA[n1], B[n2])
   end

   # stage 3: get the gradients
   fill!(dEs, zero(JVec{T}))
   A = tmpd.A
   dA = tmpd.dA
   for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
      evaluate_d!(A, dA, tmpd_basis1p, basis1p, R, Z, z0)
      iz = z2i(basis1p, Z)
      inds = basis1p.Aindices[iz, iz0]
      for iA = 1:length(basis1p, iz, iz0)
         dEs[iR] += real(B[inds[iA]] * dA[inds[iA]])
      end
   end

   return dEs
end


text:
module SimplePosets

using SimpleGraphs, Primes

import Base.show, Base.isequal, Base.hash
import Base.inv, Base.intersect #, Base.zeta
import Base.adjoint, Base.*, Base.+, Base./ #, Base.\
import Base.==

import SimpleGraphs.add!, SimpleGraphs.has, SimpleGraphs.delete!
import SimpleGraphs.relabel

export SimplePoset, IntPoset, check, hash, element_type
export elements, relations, incomparables
export card, show, add!, has, delete!
export above, below, interval
export maximals, minimals
export relabel
export zeta_matrix, zeta, mobius_matrix, mobius

export Chain, Antichain, Divisors, BooleanLattice, StandardExample
export RandomPoset, PartitionLattice

export ComparabilityGraph, CoverDigraph

export inv, intersect, stack, height

"""
`SimplePoset()` creates a new partially ordered set (poset) in which
the elements can be of `Any` type.

Use `SimplePoset(T)` or `SimplePlot{T}()` to create a new poset in
which the elements are of type `T`.
"""
mutable struct SimplePoset{T}
    D::SimpleDigraph{T}
    function SimplePoset{T}() where T
        D = SimpleDigraph{T}()
        forbid_loops!(D)
        new(D)
    end
end

# Create a new poset whose elements have a specific type (default Any)
SimplePoset(T::DataType=Any) = SimplePoset{T}()

# Validation check. This should not be necessary to ever use if the
# poset was properly built.
function check(P::SimplePoset)

    # cycle detection
    PP = deepcopy(P)
    while true
        bottoms = minimals(PP)
        if length(bottoms)==0
            break
        end
        for b in bottoms
            delete!(PP,b)
        end
    end
    if card(PP)>0
        warn("Cycles detected")
        return false
    end

    # transitive closure check
    Z = zeta_matrix(P)
    if countnz(Z) != countnz(Z*Z)
        warn("Not transitively closed")
        return false
    end
    return true
end

"""
`element_type(P)` returns the type of elements in this `SimplePoset`.
"""
element_type(P::SimplePoset{T}) where T = T

# Check if two posets are the same
isequal(P::SimplePoset, Q::SimplePoset) = isequal(P.D,Q.D)
==(P::SimplePoset, Q::SimplePoset) = isequal(P,Q)

# hash function for this class based on P.D
hash(P::SimplePoset, h::UInt64 = UInt64(0)) = hash(P.D,h)

# return a list of the elements in P
"""
`elements(P)` returns a list of the elements of `P`.
"""
elements(P::SimplePoset) = vlist(P.D)

# return a list of all the < relations in P
"""
`relations(P)` returns a list of ordered pairs `(u,v)` where `u` and
`v` are elements of the poset that satisfy `u<v` in `P`.
"""
relations(P::SimplePoset) = elist(P.D)

# list all pairs of elements that are incomparable to each other

"""
`incomparables(P)` returns a list of ordered pairs `(u,v)` such that
`u` and `v` are incomparable in `P`. Note that if `(u,v)` appears in
the list, we do not also include `(v,u)`.
"""
function incomparables(P::SimplePoset{T}) where T
    els = elements(P)
    n   = length(els)

    pairs = Tuple{T,T}[]
    for j=1:n-1
        for k=j+1:n
            push!(pairs, (els[j],els[k]) )
        end
    end

    filter( p -> !has(P,p[1],p[2]) && !has(P,p[2],p[1]) , pairs)
end

# return the cardinality of this poset
"""
`card(P)` returns the cardinality (number of elements) of `P`.
"""
card(P::SimplePoset) = NV(P.D)

# How we print posets to the terminal
function show(io::IO, P::SimplePoset)
    print(io, "SimplePoset{$(element_type(P))} ($(card(P)) elements)")
end

display(P::SimplePoset) = show(P)

# Add an element to the groundset of this poset
"""
`add!(P,x)` adds the element `x` to the poset with no relations to any
other elements.

`add!(P,x,y)` adds elements `x` and `y` to the poset (if they are not
already present) and, more importantly, adds the relation `x<y` as
well. This fails if adding this relation would violate transitivity.
"""
function add!(P::SimplePoset{T}, x) where T
    return add!(P.D, x)
end

# Add x<y as a relation in this poset
function add!(P::SimplePoset{T}, x, y) where T
    # start with some basic checks
    if !has(P,x)
        add!(P,x)
    end
    if !has(P,y)
        add!(P,y)
    end
    if x==y || has(P,y,x) || has(P,x,y)
        return false
    end

    U = above(P,y)
    push!(U,y)

    D = below(P,x)
    push!(D,x)

    for u in U
        for d in D
            add!(P.D,d,u)
        end
    end
    return true
end

# Delete an element from P
"""
`delete!(P,x)` deletes the element `x` from the poset `P`.

`delete!(P,x,y)` deletes the relation `x<y` from `P` (assuming that
such a relation exists in `P`). It then deletes other relations as
needed to ensure what remains is still a poset (i.e., still
transitive).
"""
delete!(P::SimplePoset, x) = delete!(P.D,x)

# Delete a relation from P (see Doc folder in github for explanation)
function delete!(P::SimplePoset, x, y)
    if !has(P,x) || !has(P,y) || x==y || !has(P,x,y)
        return false
    end

    delete!(P.D, x, y)

    for z in P.D.V
        if z==x || z==y
            continue
        end
        if has(P,x,z) && has(P,z,y)
            delete!(P.D,x,z)
            delete!(P.D,z,y)
        end
    end
    return true
end

# Check if a particular element is in the ground set

"""
`has(P,x)` checks if `x` is an element of the poset `P`.

`has(P,x,y)` checks if `x<y` is a relation in the poset `P`.
"""
has(P::SimplePoset{T}, x)  where T = has(P.D, x)

# Check if x<y holds in this poset
has(P::SimplePoset{T}, x, y) where T = has(P.D, x, y)

# return a list of all elements > xe

"""
`above(P,x)` returns a list of all elements `y` in the poset `P` for
which `x<y`.
"""
function above(P::SimplePoset, x)
    if !has(P,x)
        error("This poset does not contain ", x)
    end
    return collect(P.D.N[x])
end

"""
`below(P,x)` returns a list of all elements `y` in the poset `P` for
which `y<x`.
"""
function below(P::SimplePoset, x)
    if !has(P,x)
        error("This poset does not contain ", x)
    end
    return collect(P.D.NN[x])
end

# return a list of all elements z with x < z < y

"""
`interval(P,x,y)` returns a list of all elements `z` in the poset `P`
for which `x<z<y`.
"""
function interval(P::SimplePoset, x, y)
    A = Set(above(P,x))
    B = Set(below(P,y))
    return collect(intersect(A,B))
end

# Construct an antichain with n elements 1,2,...,n
"""
`AntiChain(n)` creates a new poset with `Int` elements `1:n` with no
relations between those elements.

`AntiChain(list)` creates a new poset with elements from `list` with
no relations between those elements.
"""
function Antichain(n::Int)
    if n < 0
        error("Number of elements must be nonnegative")
    end
    P = SimplePoset(Int)
    for e = 1:n
        add!(P,e)
    end
    return P
end

"""
`IntPoset(n)` creates a new `SimplePoset{Int}()` prepopulated with
elements `1:n` but no relations.
"""
IntPoset(n::Int) = Antichain(n)

# Construction an antichain from a list of elements
function Antichain(els::Array{T,1}) where T
    P = SimplePoset(T)
    for e in els
        add!(P,e)
    end
    return P
end

# Construct a chain 1<2<3<...<n
"""
`Chain(n)` creates a new poset whose elements are the `Int` values
`1:n` with the relations `1<2<3<...<n`.

`Chain(list)` creates a new poset whose elements are the values in
`list` with `list[1]<list[2]<...`.
"""
function Chain(n::Int)
    P = Antichain(n)
    for k=1:n
        add!(P,k)
    end
    if n > 1
        for k=1:n-1
            add!(P,k,k+1)
        end
    end
    return P
end

# Construct a chain given a list of elements
function Chain(els::Array{T,1}) where T
    P = Antichain(els)
    n = length(els)
    for k=1:n-1
        add!(P,els[k],els[k+1])
    end
    return P
end

# requires n>1, but we don't check. gives first prime factor. this is
# not exposed.
function first_prime_factor(n::Int)
    if isprime(n)
        return n
    end

    for k=2:n
        if n%k == 0
            return k
        end
    end
end

# creates the set of divisors of a positive integer. should we expose?
function divisors(n::Int)
    if n<1
        error("divisors only works on positive integers")
    end

    if n==1
        return BitSet(1)
    end

    p = first_prime_factor(n)
    if n==p
        return BitSet([1,p])
    end

    A = divisors(div(n,p))
    Alist = collect(A)

    Blist = [ p*x for x in Alist ]

    B = BitSet(Blist)

    return union(A,B)
end

# Create the poset of the divisors of a positive integer
"""
`Divisors(n)` creates a new poset whose elements are the positive
divisor of `n` (including `1` and `n` itself) in which we have the
relations `(u,v)` precisely when `u` is a factor of `v`.
"""
function Divisors(n::Int)
    if n<1
        error("Argument must be a positive integer")
    end

    A = divisors(n)
    P = SimplePoset(Int)

    for a in A
        add!(P,a)
    end

    for a in A
        for b in A
            if a!=b && b%a == 0
                add!(P,a,b)
            end
        end
    end
    return P
end

# Create the Boolean lattice poset. Elements are n-long binary
# strings.
"""
`BooleanLattice(n)` creates the Boolean lattice whose elements are `n`-long
character strings of 0s and 1s. Ordering is coordinatewise.
"""
function BooleanLattice(n::Int)
    if n<1
        error("Argument must be a positive integer")
    end

    P = SimplePoset(String)

    NN = (1<<n) - 1
    for e = 0:NN
        # add!(P,bin(e,n))
        add!(P,string(e,base=2,pad=n))
    end

    for e = 0:NN
        for f=0:NN
            if e!=f && e|f == f
                # add!(P.D,bin(e,n), bin(f,n))
                add!(P.D,string(e,base=2,pad=n),string(f,base=2,pad=n))
            end
        end
    end

    return P
end

# Helper function for RandomPoset
function vec_less(x::Array{Float64,1}, y::Array{Float64,1})
    n = length(x)
    return all([ x[k] <= y[k] for k=1:n ])
end

# Create a random d-dimensional poset with n elements

"""
`RandomPoset(n,d)` creates a random `d`-dimensional poset with
elements `1:n`.
"""
function RandomPoset(n::Int, d::Int)
    if n<1 || d<1
        error("Require n and d positive in RandomPoset(n,d)")
    end
    vectors = [ rand(d) for k=1:n ]

    P = SimplePoset(Int)
    for k=1:n
        add!(P,k)
    end

    for i=1:n
        for j=1:n
            if i!=j && vec_less(vectors[i],vectors[j])
                add!(P,i,j)
            end
        end
    end
    return P
end

# Create standard example poset. Lower level named by negatives and
# upper level by positives.

"""
`StandardExample(n)` creates a new poset with `2n` elements in two
levels. Each element on the lower level is below exactly `n-1`
elements from the upper level.

We name the lower elements `-1,-2,...,-n` and the upper elements
`1,2,...,n`. We have `-i` below all the positive elements *except*
`+i`.

```
julia> P = StandardExample(4)
SimplePoset{Int64} (8 elements)

julia> elements(P)
8-element Array{Int64,1}:
 -4
 -3
 -2
 -1
  1
  2
  3
  4

julia> relations(P)
12-element Array{Tuple{Int64,Int64},1}:
 (-4,1)
 (-4,2)
 (-4,3)
 (-3,1)
 (-3,2)
 (-3,4)
 (-2,1)
 (-2,3)
 (-2,4)
 (-1,2)
 (-1,3)
 (-1,4)
```
"""
function StandardExample(n::Int)
    if n<1
        error("Argument must be a positive integer")
    end

    P = SimplePoset(Int)
    for e=1:n
        add!(P,e)
        add!(P,-e)
        for f=1:n
            if e!=f
                add!(P.D,-e,f)
            end
        end
    end
    return P
end

# maximal and minimal elements

"""
`maximals(P)` returns a list of maximal elements of `P`.
"""
maximals(P::SimplePoset) = filter(x->out_deg(P.D,x)==0, elements(P))

"""
`minimals(P)` returns a list of minimal elements of `P`.
"""
minimals(P::SimplePoset) = filter(x->in_deg(P.D,x)==0, elements(P))

# The inverse of a poset is a new poset with the order reversed
"""
`inv(P)` creates a new poset with the same elements as `P` in which
all of `P`'s relations have been reversed. That is `(u,v)` is a
relation of `P` iff `(v,u)` is a relation of `inv(P)`.

Use `P'` as a shortcut for `inv(P)`.
"""
function inv(P::SimplePoset{T}) where T
    Q = SimplePoset(T)
    for e in P.D.V
        add!(Q,e)
    end
    for r in relations(P)
        x,y = r[1],r[2]
        add!(Q.D,y,x)
    end
    return Q
end

# We can use P' to mean inv(P) also
adjoint(P::SimplePoset) = inv(P)

# Create the intersection of two posets (must be of same element
# type). Ideally, the two posets have the same set of elements, but
# this is not necessary; if they don't we just intersect the element
# sets first.

"""
`intersect(P,Q)` constructs a new poset that is the intersection of
the two given posets (which must contain elements of the same
datatype).
"""
function intersect(P::SimplePoset{T}, Q::SimplePoset{T}) where T
    R = SimplePoset(T)
    elist = filter(x -> has(P,x), elements(Q))
    for e in elist
        add!(R,e)
    end

    rlist = filter( r -> has(P,r[1],r[2]), relations(Q))
    for r in rlist
        add!(R.D, r[1],r[2])
    end
    return R
end

# Produce the cartesian product of two posets
"""
`P*Q` is the Cartesian product of the two posets. The elements need
not be of the same type.
"""
function (*)(P::SimplePoset{S}, Q::SimplePoset{T}) where {S,T}
    PQ = SimplePoset{Tuple{S,T}}()

    for a in P.D.V
        for b in Q.D.V
            add!(PQ,(a,b))
        end
    end

    elist = elements(PQ)
    for alpha in elist
        for beta in elist
            if alpha != beta
                if has(P,alpha[1],beta[1]) || alpha[1]==beta[1]
                    if has(Q,alpha[2],beta[2]) || alpha[2]==beta[2]
                        add!(PQ.D, alpha, beta)
                    end
                end
            end
        end
    end
    return PQ
end

# Disjoint union of posets

# This seems to fix some warnings. Can't say I understand why.
function +()
    nothing
end

"""
`P+Q` is the disjoint union of the two (or more) posets. The poset
elements must all be of the same type.
"""
function (+)(x::SimplePoset{T}...) where T
    PP = SimplePoset{Tuple{T,Int}}()

    for i=1:length(x)
        P = x[i]
        for e in P.D.V
            add!(PP, (e,i))
        end

        for r in relations(P)
            a = r[1]
            b = r[2]
            add!(PP.D, (a,i), (b,i))
        end
    end

    return PP
end

# Stack a bunch of posets one atop the next. The first one in the
# argument list is at the bottom.

"""
`stack(P...)` stacks its poset arguments each on top of each
other. The first poset in the list is on the bottom.
"""
function stack(x::SimplePoset{T}...) where T
    np = length(x)
    PP = +(x...)

    for i=1:np-1
        P = x[i]
        for j=i+1:np
            Q = x[j]
            for a in P.D.V
                for b in Q.D.V
                    add!(PP.D, (a,i),(b,j))
                end
            end
        end
    end
    return PP
end

# Binary operator version of stack P/Q puts P on top while P\Q puts Q
# on top.

"""
`P/Q` stacks poset `P` on top of poset `Q`.
"""
(/)(P::SimplePoset{T}, Q::SimplePoset{T}) where T = stack(Q,P)
#
# """
# `P\Q` stacks poset `P` under poset `Q`.
# """
# (\){T}(P::SimplePoset{T}, Q::SimplePoset{T}) = stack(P,Q)
#
# # Zeta function as a matrix

"""
`zeta_matrix(P)` creates the zeta matrix of a poset `P`. The
rows/columns of this matrix are indexed by the elements of `P`. We
have a `1` in the `i,j`-position exactly when the `i`th element of `P`
is `<=` the `j`th element of `P`.

The ordering of the rows/columns does not necessarily lead to an upper
triangular matrix.
"""
function zeta_matrix(P::SimplePoset)
 elist = elements(P)
    n = length(elist)
    Z = zeros(Int, n, n)

    for i=1:n
        for j=1:n
            if i==j || has(P,elist[i],elist[j])
                Z[i,j] = 1
            end
        end
    end
    return Z
end

# Mobius function as a matrix

"""
`mobius_matrix(P)` returns the inverse of `zeta_matrix(P)`.
"""
mobius_matrix(P::SimplePoset) = round.(Int,inv(zeta_matrix(P)))

# Zeta function as a dictionary

"""
`zeta(P)` returns the zeta function of `P` as a `Dict`. If `d=zeta(P)`
then `d[u,v]==1` when `u` and `v` are elements of the poset for which
`u==v` or `(u,v)` is a relation of `P`.
"""
function zeta(P::SimplePoset{T}) where T
    z = Dict{Tuple{T,T},Int}()
    els = elements(P)
    for a in els
        for b in els
            if a==b || has(P,a,b)
                z[a,b] = 1
            else
                z[a,b] = 0
            end
        end
    end
    return z
end

# Mobius function of this poset.

"""
`mobius(P)` returns the Mobius function of the poset `P` as a `Dict`.
See `zeta`.
"""
function mobius(P::SimplePoset{T}) where T
    mu = Dict{Tuple{T,T},Int}()
    els = elements(P)
    M = mobius_matrix(P)
    n = length(els)

    for i=1:n
        a = els[i]
        for j=1:n
            b = els[j]
            mu[a,b] = M[i,j]
        end
    end
    return mu
end

# The comparability graph of a poset

"""
`ComparabilityGraph(P)` returns the comparability graph of the poset
`P`. The result is a `SimpleGraph` whose vertices are the elements of
`P` in which there is an edge `(u,v)` provided `(u,v)` or `(v,u)` is a
relation of `P`.
"""
ComparabilityGraph(P::SimplePoset) = simplify(P.D)

# The CoverDigraph of a poset P is a directed graph that has the same
# vertices as P, in which (x,y) is an edge iff x<y and there is no z with x<z<y
"""
`CoverDigraph(P)` creates a `SimpleDirectedGraph` `G` whose vertices
are the elements of `P` and in which we have the edge `(u,v)`
provided: (a) `u<v` is a relation of `P` and (b) there is no `w` in
`P` with `u<w<v`.
"""
function CoverDigraph(P::SimplePoset{T}) where T
    CD = SimpleDigraph{T}()
    for v in P.D.V
        add!(CD,v)
    end

    for r in relations(P)
        x = r[1]
        y = r[2]
        add_flag::Bool = true
        for z in P.D.V
            if has(P,x,z) && has(P,z,y)
                add_flag = false
                break
            end
        end
        if add_flag
            add!(CD,x,y)
        end
    end
    return CD
end

# Relabel the vertics of a poset based on a dictionary mapping old
# element names to new

"""
`relabel(P)` creates a new poset `P` in which the elements are
relabeled using the values `1:n` (where `n=card(P)`).

`relabel(P,d)` relabels elements according to the dictionary `d`.
"""
function relabel(P::SimplePoset{S}, label::Dict{S,T}) where {S,T}
    Q = SimplePoset{T}()
    Q.D = relabel(P.D, label)
    return Q
end

# Relabel the elements with the integers 1:n
function relabel(P::SimplePoset{S}) where S
    verts = vlist(P.D)
    n = length(verts)
    label = Dict{S,Int}()
    sizehint!(label,n)

    for idx = 1:n
        label[verts[idx]] = idx
    end

    return relabel(P,label)
end

# Compute the height by stripping off minimals repeatedly

"""
`height(P)` gives the size of a largest chain in the poset `P`.
"""
function height(P::SimplePoset)
    PP = deepcopy(P)
    result = 0

    while card(PP) > 0
        result += 1
        M = minimals(PP)
        for x in M
            delete!(PP,x)
        end
    end

    return result
end

using SimplePartitions

"""
`PartitionLattice(n)` creates the poset whose elements are the
partitions of the set `{1,2,...,n}` ordered by refinement.
"""
function PartitionLattice(n::Int)
    X = all_partitions(n)
    P = SimplePoset{Partition{Int}}()
    for x in X
        for y in X
            if x <= y
                add!(P,x,y)
            end
        end
    end
    return P
end

include("linear_extensions.jl")
include("average_height.jl")

end # end of module SimplePosets


text:
export PointForecast, ScenariosForecast, QuantilesForecast

export PointForecastMetrics, ScenariosForecastMetrics

export forecast_metrics

abstract type Forecast end
abstract type ProbabilisticForecast <: Forecast end

"""
    PointForecast

Define the results of point forecasts.
"""
mutable struct PointForecast{T<:Real} <: Forecast
    name::String
    timestamps::Vector{DateTime}
    forecast::Vector{T}

    function PointForecast(
        name::String,
        timestamps::Vector{DateTime},
        forecast::Vector{T},
    ) where {T<:Real}

        if length(timestamps) != length(forecast)
            throw(DimensionMismatch("timestamps and forecast do not have the same length."))
        end

        if !allunique(timestamps)
            throw(ErrorException("timestamps must be unique."))
        end

        if !issorted(timestamps)
            throw(ErrorException("timestamps must be sorted."))
        end

        return new{T}(name, timestamps, forecast)
    end
end

# Additional constructors
function PointForecast(
    name::String,
    timestamps::Vector{Date},
    forecast::Vector{T},
) where {T<:Real}

    return PointForecast(name, DateTime.(timestamps), forecast)
end

## Evaluation Metrics for Point Forecast
struct PointForecastMetrics{T<:Real}
    errors::Vector{T}
    absolute_percentage_errors::Vector{T}
end

function forecast_metrics(
    point_forecast::PointForecast{T},
    real_ts::TimeSeries{T},
) where {T}

    if length(point_forecast.forecast) != length(real_ts.vals)
        throw(DimensionMismatch("real values and forecast do not have the same length."))
    end

    err = error(real_ts.vals, point_forecast.forecast)
    if observations_close_to_zero(real_ts)
        @warn(
            "The real observations have values too close to zero. This makes the " *
            "absolute percentage error impractical, NaNs will be returned."
        )
        abs_percentage_err = [NaN]
    else
        abs_percentage_err =
            absolute_percentage_error(real_ts.vals, point_forecast.forecast)
    end
    return PointForecastMetrics{T}(err, abs_percentage_err)
end

error(real::Vector{T}, forecast::Vector{T}) where {T} = real .- forecast
absolute_percentage_error(real::Vector{T}, forecast::Vector{T}) where {T} =
    abs.(error(real, forecast) ./ real)

"""
    ScenariosForecast

Define the probabilistic forecast results calculated from scenarios.
"""
mutable struct ScenariosForecast{T<:Real} <: ProbabilisticForecast
    name::String
    timestamps::Vector{DateTime}
    scenarios::Matrix{T}
    quantiles_probabilities::Vector{T}
    quantiles::Matrix{T}

    function ScenariosForecast(
        name::String,
        timestamps::Vector{DateTime},
        scenarios::Matrix{T},
        quantiles_probabilities::Vector{T},
        quantiles::Matrix{T},
    ) where {T<:Real}

        if length(timestamps) != size(scenarios, 1)
            throw(DimensionMismatch("timestamps and scenarios do not have the same length."))
        end

        if length(timestamps) != size(quantiles, 1)
            throw(DimensionMismatch("timestamps and quantiles do not have the same length."))
        end

        if length(quantiles_probabilities) != size(quantiles, 2)
            throw(DimensionMismatch("quantiles_probabilities and quantiles do not have the same length."))
        end

        if !allunique(timestamps)
            throw(ErrorException("timestamps must be unique."))
        end

        if !issorted(timestamps)
            throw(ErrorException("timestamps must be sorted."))
        end

        return new{T}(name, timestamps, scenarios, quantiles_probabilities, quantiles)
    end
end

# Additional constructors
function ScenariosForecast(
    name::String,
    timestamps::Vector{Date},
    scenarios::Matrix{T},
    quantiles_probabilities::Vector{T},
    quantiles::Matrix{T},
) where {T<:Real}

    return ScenariosForecast(
        name,
        DateTime.(timestamps),
        scenarios,
        quantiles_probabilities,
        quantiles,
    )
end

## Evaluation Metrics for Point Forecast
struct ScenariosForecastMetrics
    # A Vector of dictionaaries for each lead time
    # where each key is a quantile and each value
    # is wheather the observation is smaller than
    # the quantile.
    probabilistic_calibration::Vector{Dict{Float64,Bool}}
    # A Vector of dictionaaries for each lead time
    # where each key is the confidence interval
    # and each value is the width of the interval in the
    # unit of the dependent variable.
    interval_width::Vector{Dict{Float64,Float64}}
    # A Vector that contains the crps score for each
    # lead time
    crps::Vector{Float64}
end

function forecast_metrics(
    scenarios_forecast::ScenariosForecast,
    real_ts::TimeSeries{T},
) where {T}

    if size(scenarios_forecast.scenarios, 1) != length(real_ts.vals)
        throw(DimensionMismatch("real values and forecast do not have the same length."))
    end

    probabilistic_calibration =
        evaluate_probabilistic_calibration(scenarios_forecast.scenarios, real_ts.vals)
    interval_width = evaluate_interval_width(scenarios_forecast.scenarios)
    crps = evaluate_crps(scenarios_forecast.scenarios, real_ts.vals)

    return ScenariosForecastMetrics(probabilistic_calibration, interval_width, crps)
end

function get_quantiles(quantile_probs::Vector{T}, scenarios::Matrix{T}) where {T}
    quantiles = mapslices(x -> quantile(x, quantile_probs), scenarios; dims = 2)
    return quantiles
end

## probabilistic_calibration functions
hitted_below_quantile(val::T, quantile::T) where {T} = val <= quantile

function evaluate_probabilistic_calibration(scenarios::Matrix{T}, vals::Vector{T}) where {T}
    quantile_probs = collect(0.025:0.05:0.975)
    quantiles = get_quantiles(quantile_probs, scenarios)
    probabilistic_calibration = Vector{Dict}(undef, size(scenarios, 1))
    for k = 1:size(scenarios, 1)
        probabilistic_calibration[k] = Dict{Float64,Bool}()
        for (i, q) in enumerate(quantile_probs)
            probabilistic_calibration[k][q] =
                hitted_below_quantile(vals[k], quantiles[k, i])
        end
    end
    return probabilistic_calibration
end

lower_quantile(interval_prob::T) where {T} = (1 - interval_prob) / 2
upper_quantile(interval_prob::T) where {T} = (1 + interval_prob) / 2

width_of_interval(upper_quantile::T, lower_quantile::T) where {T} =
    upper_quantile - lower_quantile

function evaluate_interval_width(scenarios::Matrix{T}) where {T}
    intervals_probs = collect(0.05:0.05:0.95)
    interval_width = Vector{Dict}(undef, size(scenarios, 1))
    upper_quantiles = get_quantiles(upper_quantile.(intervals_probs), scenarios)
    lower_quantiles = get_quantiles(lower_quantile.(intervals_probs), scenarios)
    for k = 1:size(scenarios, 1)
        interval_width[k] = Dict{Float64,Float64}()
        for (i, q) in enumerate(intervals_probs)
            interval_width[k][q] =
                width_of_interval(upper_quantiles[k, i], lower_quantiles[k, i])
        end
    end
    return interval_width
end

# crps functions
discrete_crps_indicator_function(val::T, z::T) where {T} = val < z

function crps(scenarios::Vector{T}, val::T) where {T}
    sorted_scenarios = sort(scenarios)
    m = length(scenarios)
    crps_score = zero(T)

    for i = 1:m
        crps_score +=
            (sorted_scenarios[i] - val) *
            (m * discrete_crps_indicator_function(val, sorted_scenarios[i]) - i + 0.5)
    end

    return (2 / m^2) * crps_score
end

function evaluate_crps(scenarios::Matrix{T}, vals::Vector{T}) where {T}
    crps_scores = Vector{Float64}(undef, length(vals))

    for k = 1:length(vals)
        crps_scores[k] = crps(scenarios[k, :], vals[k])
    end

    return crps_scores
end

function mean_of_metrics(forecast::Vector{ScenariosForecastMetrics})
    # TODO assert all forecast metrics have the same length.
    mean_of_probabilistic_calibration = mean_probabilistic_calibration(forecast)
    mean_of_interval_width = mean_interval_width(forecast)
    mean_of_crps = mean_crps(forecast)
    return mean_of_probabilistic_calibration, mean_of_interval_width, mean_of_crps
end

function mean_probabilistic_calibration(forecast::Vector{ScenariosForecastMetrics})
    m_probabilistic_calibrations =
        Vector{Dict{Float64,Float64}}(undef, length(forecast[1].probabilistic_calibration))
    for t = 1:length(m_probabilistic_calibrations)
        m_probabilistic_calibrations[t] = Dict{Float64,Float64}()
        for forec in forecast
            for (k, v) in forec.probabilistic_calibration[t]
                if !haskey(m_probabilistic_calibrations[t], k)
                    m_probabilistic_calibrations[t][k] = v / length(forecast)
                else
                    m_probabilistic_calibrations[t][k] += v / length(forecast)
                end
            end
        end
    end
    return m_probabilistic_calibrations
end

function mean_interval_width(forecast::Vector{ScenariosForecastMetrics})
    m_interval_width =
        Vector{Dict{Float64,Float64}}(undef, length(forecast[1].interval_width))
    for t = 1:length(m_interval_width)
        m_interval_width[t] = Dict{Float64,Float64}()
        for forec in forecast
            for (k, v) in forec.interval_width[t]
                if !haskey(m_interval_width[t], k)
                    m_interval_width[t][k] = v / length(forecast)
                else
                    m_interval_width[t][k] += v / length(forecast)
                end
            end
        end
    end
    return m_interval_width
end

function mean_crps(forecast::Vector{ScenariosForecastMetrics})
    m_crps = Matrix{Float64}(undef, length(forecast[1].crps), length(forecast))
    for (i, forec) in enumerate(forecast)
        m_crps[:, i] = forec.crps
    end
    return vec(mean(m_crps, dims = 2))
end

"""
    QuantilesForecast

Define the probabilistic forecast results calculated from distributions.
"""
mutable struct QuantilesForecast{T<:Real} <: ProbabilisticForecast
    name::String
    timestamps::Vector{DateTime}
    quantiles_probabilities::Vector{T}
    quantiles::Matrix{T}

    function QuantilesForecast(
        name::String,
        timestamps::Vector{DateTime},
        quantiles_probabilities::Vector{T},
        quantiles::Matrix{T},
    ) where {T<:Real}

        if length(timestamps) != size(quantiles, 1)
            throw(DimensionMismatch("timestamps and quantiles do not have the same length."))
        end

        if length(quantiles_probabilities) != size(quantiles, 2)
            throw(DimensionMismatch("quantiles_probabilities and quantiles do not have the same length."))
        end

        if !allunique(timestamps)
            throw(ErrorException("timestamps must be unique."))
        end

        if !issorted(timestamps)
            throw(ErrorException("timestamps must be sorted."))
        end

        return new{T}(name, timestamps, quantiles_probabilities, quantiles)
    end
end

# Additional constructors
function QuantilesForecast(
    name::String,
    timestamps::Vector{Date},
    quantiles_probabilities::Vector{T},
    quantiles::Matrix{T},
) where {T<:Real}

    return QuantilesForecast(
        name,
        DateTime.(timestamps),
        quantiles_probabilities,
        quantiles,
    )
end


text:
"""
Authors: Yaroslav Salii, 2020+
         Kara Ignatenko, 2021

This is Oboe-Mangle, or maybe Arshyn, with a view to generate instances for testing
computational methods for networked epidemic models with data from
    (a) FluTE/US 2010 Census Tracts (population, coordinates); github.com/dlchao/FluTE
    (b) US 2019+ Domestic Flights, available from BTS

Reading the FluTE data and transforming it to my input format.
Might separate my input format definitions later.

Oboe.jl v.0.1: "From scripts to proper code" edition
        v.0.2: switch to `module`
        v.0.3: got as far as reading FluTE tracts
        v.0.4: added basic by-county and by-state IV aggregation
        v.0.5: join (intersect) BTS with OpenFlights
        v.0.6: added designated AP assignment
'20-12-18   v.0.7: wrote pairs-to-matrix xform by hand (vs. `unstack`, which was unpredictable)
'20-12-30   v.0.8: added a beta node-to-node daily air passenger computation
'21-01-03   v.0.8.1: a half-baked Main(), look in bit bucket. Tested on 3K by-county!
'21-01-22   v.0.9: air travel and commute matrices are in, tested on 2K NW (by-tract)
'21-02-02   v.0.9.1: hotfix, removing deprecated `by`, `names!`, and the like (half-done)
'21-02-02   v.0.9.2: removed all the deprecated stuff, and it all works again
'21-02-14   v.0.9.3: add node partitioning, streamline mkFlightMx2() and mkPsgMx2() with aux
'21-02-14   v.0.9.4: add partition-to-partition flight matrix
'21-02-18   v.0.9.5: add partition-to-partition commute matrix
'21-03-31   v.0.9.6: min,median,max for nonzero elements in talkDnsy
            v.0.A: add by-AP agg with recomputation bypass in mkPsgMx()
'21-07-14   v.0.A.1: add checkIfFilesExist and censorFluteTractByFIPS
'21-07-19   v.1.0: fixed issue where unmatched tract IDs led to a crash
'21-08-10   v.1.1: significantly improved performance of mkPsgMx
'21-08-17   v.2.0: separated functionality into submodules;
                   API changes to reduce coupling between submodules
'21-08-24   v.2.0.1: renamed io.jl to data-io.jl, no external API change
"""
module Oboe
using FromFile

const callsign = "This is Oboe v.2.0.1" * begin
    #display the name of the current branch after the version if it's not master
    #useful for debugging/notebooks/etc
    try
        branch = readchomp(`git branch --show-current`)
        branch == "master" ? "" : "-" * branch
    catch err
        ""
    end
end

@from "data-io.jl"     using DataIO
@from "airports.jl"    using Airports
@from "tracts.jl"      using Tracts
@from "aggregation.jl" using Aggregation
@from "travel.jl"      using Travel

end


text:
# Resampling schemes for particle filters
# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

# default resampling scheme
function resample(w::AbstractVector{<:Real}, num_particles::Integer=length(w))
    return resample_systematic(w, num_particles)
end

# More stable, faster version of rand(Categorical)
function randcat(p::AbstractVector{T}) where T<:Real
    r, s = rand(T), 1
    for j in eachindex(p)
        r -= p[j]
        if r <= zero(T)
            s = j
            break
        end
    end
    return s
end

function resample_multinomial(w::AbstractVector{<:Real}, num_particles::Integer)
    return rand(Distributions.sampler(Categorical(w)), num_particles)
end

function resample_residual(w::AbstractVector{<:Real}, num_particles::Integer)

    M = length(w)

    # "Repetition counts" (plus the random part, later on):
    Ns = floor.(length(w) .* w)

    # The "remainder" or "residual" count:
    R = Int(sum(Ns))

    # The number of particles which will be drawn stocastically:
    M_rdn = num_particles - R

    # The modified weights:
    Ws = (M .* w - floor.(M .* w)) / M_rdn

    # Draw the deterministic part:
    indx1, i = Array{Int}(undef, R), 1
    for j in 1:M
        for k in 1:Ns[j]
            indx1[i] = j
            i += 1
        end
    end

    # And now draw the stocastic (Multinomial) part:
    return append!(indx1, rand(Distributions.sampler(Categorical(w)), M_rdn))
end

function resample_stratified(w::AbstractVector{<:Real}, num_particles::Integer)

    Q, N = cumsum(w), num_particles

    T = Array{Float64}(undef, N + 1)
    for i=1:N,
        T[i] = rand() / N + (i - 1) / N
    end
    T[N+1] = 1

    indx, i, j = Array{Int}(undef, N), 1, 1
    while i <= N
        if T[i] < Q[j]
            indx[i] = j
            i += 1
        else
            j += 1
        end
    end
    return indx
end

function resample_systematic(w::AbstractVector{<:Real}, num_particles::Integer)

    Q, N = cumsum(w), num_particles

    T = collect(range(0, stop = maximum(Q)-1/N, length = N)) .+ rand()/N
    push!(T, 1)

    indx, i, j = Array{Int}(undef, N), 1, 1
    while i <= N
        if T[i] < Q[j]
            indx[i] = j
            i += 1
        else
            j += 1
        end
    end
    return indx
end