text:
struct NLADefault{T<:AbstractFloat} <: NonLinearAlgorithm end

"""
    NLADefault()
Return the array untouched, default option.
"""
NLADefault() = NLADefault{Float64}()

nla(::NLADefault{T}, x::AbstractArray{T}) where T<: AbstractFloat = x

struct NLAT1{T<:AbstractFloat} <: NonLinearAlgorithm end
"""
    NLAT1()
Apply the \$ \\text{T}_1 \$ transformation algorithm, as defined in [1] and [2]

[1] Chattopadhyay, Ashesh, et al. "Data-driven prediction of a multi-scale Lorenz 96 chaotic system using a hierarchy of deep learning methods: Reservoir computing, ANN, and RNN-LSTM." (2019).

[2] Pathak, Jaideep, et al. "Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach." Physical review letters 120.2 (2018): 024102.
"""
NLAT1() = NLAT1{Float64}()
nla(::NLAT1{T}, x::AbstractArray{T}) where T<: AbstractFloat = _nlat1(x)

struct NLAT2{T<:AbstractFloat} <: NonLinearAlgorithm end
"""
    NLAT2()
Apply the \$ \\text{T}_2 \$ transformation algorithm, as defined in [1]

[1] Chattopadhyay, Ashesh, et al. "Data-driven prediction of a multi-scale Lorenz 96 chaotic system using a hierarchy of deep learning methods: Reservoir computing, ANN, and RNN-LSTM." (2019).
"""
NLAT2() = NLAT2{Float64}()
nla(::NLAT2{T}, x::AbstractArray{T}) where T<: AbstractFloat = _nlat2(x)

struct NLAT3{T<:AbstractFloat} <: NonLinearAlgorithm end
"""
    NLAT3()
Apply the \$ \\text{T}_3 \$ transformation algorithm, as defined in [1]

[1] Chattopadhyay, Ashesh, et al. "Data-driven prediction of a multi-scale Lorenz 96 chaotic system using a hierarchy of deep learning methods: Reservoir computing, ANN, and RNN-LSTM." (2019).
"""
NLAT3() = NLAT3{Float64}()
nla(::NLAT3{T}, x::AbstractArray{T}) where T<: AbstractFloat = _nlat3(x)

function _nlat1(x_old::AbstractArray{T}) where T<: AbstractFloat
    x_new = copy(x_old)
    for i =1:size(x_new, 1)
        if mod(i, 2)!=0
            x_new[i, :] = copy(x_old[i,:].*x_old[i,:])
        end
    end
    return x_new
end

function _nlat2(x_old::AbstractArray{T}) where T<: AbstractFloat
    x_new = copy(x_old)
    for i =2:size(x_new, 1)-1
        if mod(i, 2)!=0
            x_new[i, :] = copy(x_old[i-1,:].*x_old[i-2,:])
        end
    end
    return x_new
end

function _nlat3(x_old::Array{T}) where T<: AbstractFloat
    x_new = copy(x_old)
    for i =2:size(x_new, 1)-1
        if mod(i, 2)!=0
            x_new[i, :] = copy(x_old[i-1,:].*x_old[i+1,:])
        end
    end
    return x_new
end


text:
# position units are in ppm.
mutable struct MetaboliteRegionType{T}
    interval_starts::Vector{Vector{T}} # [spin_group_index][interval_index]
    interval_fins::Vector{Vector{T}}
end

# for a given compound, spin group, interval.
mutable struct RegionInfoType{T}

    # used for diagnostics.
    interval_indices::Vector{Tuple{Int,Int,Int}} # compound, spin group, interval indexing.

    # used for deciding the cost function positions for this region.
    st::T
    fin::T

    # used for updating results.
    spin_group_indices::Vector{Tuple{Int,Int}} # compound, spin group
end

mutable struct ExperimentInfoType{T}

    regions::Vector{RegionInfoType{T}}
    graph::LightGraphs.SimpleGraphs.SimpleGraph{Int}
    connectivity::Vector{Vector{Int}}

    sts_set::Vector{Vector{Vector{T}}}
    fins_set::Vector{Vector{Vector{T}}}
    IDs_set::Vector{Vector{Vector{Tuple{Int,Int,Int}}}}

    sts::Vector{T}
    fins::Vector{T}
    IDs::Vector{Tuple{Int,Int,Int}}
end

ext: jl


text:
using Bennu
using StaticArrays
using CUDA
using Printf
using LinearAlgebra
using WriteVTK
CUDA.allowscalar(false)

# advection velocity
const v = SVector(1, 1)

# Gaussian initial condition
function gaussian(x, t)
    FT = eltype(x)
    xp = mod1.(x .- v .* t, FT(2π))
    xc = SVector{2, FT}(π, π)
    r² = norm(xp .- xc) ^ 2
    exp(-2r²)
end

# sine initial condition
sineproduct(x, t) = prod(sin.(x .- v .* t))

function run(solution, FT, A, N, K; outputvtk=false, vtkdir="output")
    Nq = N + 1

    cell = LobattoCell{FT, A}(Nq, Nq)

    vert1d = range(FT(0), stop=FT(2π), length=K+1)

    function meshwarp(x)
        x₁, x₂ = x
        x̃₁ = x₁ + sin(x₁ / 2) * sin(x₂)
        x̃₂ = x₂ - sin(x₁) * sin(x₂ / 2) / 2
        SVector(x̃₁, x̃₂)
    end

    grid = brickgrid(meshwarp, cell, (vert1d, vert1d); periodic=(true, true))

    g, J = components(metrics(grid))

    g₁₁, g₁₂, g₂₁, g₂₂ = components(g)

    M = mass(cell)
    D₁, D₂ = derivatives(cell)
    MJ = M * J

    facenormal, faceJ = components(facemetrics(grid))
    faceM = facemass(cell)
    faceMJ = faceM * faceJ
    faceix⁻, faceix⁺ = faceindices(grid)

    facedata = (faceix⁻, faceix⁺, faceMJ, facenormal)

    rhs! = function(dq, q)
        # volume term
        dq .-= (v[1] .* g₁₁ .+ v[2] .* g₂₁) .* (D₁ * q) +
               (v[2] .* g₂₂ .+ v[1] .* g₁₂) .* (D₂ * q)
        # surface term
        for (f⁻, f⁺, fMJ, n) in zip(faceviews.(Ref(cell), facedata)...)
            dq[f⁻] .-= fMJ .* dot.(n, Ref(v)) .* (q[f⁺] .- q[f⁻]) ./ 2 ./ MJ[f⁻]
        end
    end

    timeend = FT(2π)

    # crude dt estimate
    cfl = 1 // 10
    dx = Base.step(vert1d)
    dt = cfl * dx / N

    numberofsteps = ceil(Int, timeend / dt)
    dt = timeend / numberofsteps

    RKA = (
           FT(0),
           FT(-567301805773 // 1357537059087),
           FT(-2404267990393 // 2016746695238),
           FT(-3550918686646 // 2091501179385),
           FT(-1275806237668 // 842570457699),
          )
    RKB = (
           FT(1432997174477 // 9575080441755),
           FT(5161836677717 // 13612068292357),
           FT(1720146321549 // 2090206949498),
           FT(3134564353537 // 4481467310338),
           FT(2277821191437 // 14882151754819),
          )

    if outputvtk
        mkpath(vtkdir)
        pvd = paraview_collection(joinpath(vtkdir, "timesteps"))
    end

    do_output = function(step, time, q)
        if outputvtk && step % ceil(Int, timeend / 100 / dt) == 0
            filename = "step$(lpad(step, 6, '0'))"
            vtkfile = vtk_grid(joinpath(vtkdir, filename), grid)
            P = toequallyspaced(cell)
            vtkfile["q"] = vec(Array(P * q))
            vtk_save(vtkfile)
            pvd[time] = vtkfile
        end
    end

    # initialize state
    q = solution.(points(grid), FT(0))

    # storage for rhs
    dq = similar(q)
    dq .= 0

    # initial output
    step = 0
    time = FT(0)
    do_output(step, time, q)

    ### time integration
    for step = 1:numberofsteps
        if time + dt > timeend
            dt = timeend - time
        end
        for stage = 1:length(RKA)
            dq *= RKA[stage]
            rhs!(dq, q)
            q .+= RKB[stage] * dt * dq
        end
        time += dt
        do_output(step, time, q)
    end

    # final output
    do_output(numberofsteps, timeend, q)
    outputvtk && vtk_save(pvd)

    # compute error
    qexact = solution.(points(grid), timeend)
    errf = sqrt(sum(MJ .* (q .- qexact) .^ 2))
end

let
    FT = Float64
    N = 4

    # run on the GPU if possible
    A = CUDA.has_cuda_gpu() ? CuArray : Array

    @info """Configuration:
    precision        = $FT
    polynomial order = $N
    array type       = $A
    """

    # visualize solution of advected gaussian

    K = 16
    vtkdir = "vtk_advection_K$(K)x$(K)"
    @info "Starting Gaussian advection with ($K, $K) elements"
    run(gaussian, FT, A, N, K; outputvtk=true, vtkdir)
    @info "Finished, vtk output written to $vtkdir"

    # run convergence study using a simple sine field

    @info "Starting convergence study"
    numlevels = 5
    err = zeros(FT, numlevels)
    for l in 1:numlevels
        K = 4 * 2 ^ (l - 1)
        err[l] = run(sineproduct, FT, A, N, K)
        @info @sprintf("Level %d, elements = (%2d, %2d), error = %.16e",
                       l, K, K, err[l])
    end
    rates = log2.(err[1:numlevels-1] ./ err[2:numlevels])
    @info "Convergence rates:\n" *
        join(["rate for levels $l → $(l + 1) = $(rates[l])"
              for l in 1:(numlevels - 1)], "\n")
end


text:
function _fix_compound_string(chars::Array{Char,1})

    # init tmp
    tmp_array = Array{Char,1}()

    @show chars

    if (isempty(chars) == true)
        return tmp_array
    end

    # get next char -
    next_char = pop!(chars)
    if (isnumeric(next_char) == true)

        push!(tmp_array, next_char)

        values = _fix_compound_string(chars)
        for value in values
            push!(tmp_array, value)
        end

    elseif (isnumeric(next_char) == false)

        # ok, so pop another char -
        next_next_char = pop!(chars)
        if (isnumeric(next_next_char) == false)

            push!(tmp_array, next_char)
            push!(tmp_array, '1')

            # push the test char back onto the stack -
            push!(chars, next_next_char)
        else
            push!(tmp_array, next_char)
            push!(tmp_array, next_next_char)
        end

        values = _fix_compound_string(chars)
        for value in values
            push!(tmp_array, value)
        end
    end

    # return -
    return tmp_array
end

#test_char_array = "C5H12O8P2" |> collect
#test_char_array = "KOH" |> collect
test_char_array = "H3PO4" |> collect

# add extra 1 if last char is a letter -
if (isnumeric(last(test_char_array)) == false)
    push!(test_char_array, '1')
end

C_tmp = _fix_compound_string(reverse(test_char_array))


text:
"""
Use bisection search to sample efficiently given a vector of objects and the cumulative weights.
"""
function sample_cweighted(rng::AbstractRNG, objects::AbstractVector, cweights::AbstractVector{Float64})
    t = rand(rng)*cweights[end]
    large = length(cweights) # index of cdf value that is bigger than t
    small = 0 # index of cdf value that is smaller than t
    while large > small + 1
        new = div(small + large, 2)
        if t < cweights[new]
            large = new
        else
            small = new
        end
    end
    return objects[large]
end

# XXX Not complete
# """
# Perform low variance resampling from page 110 of Probabilistic Robotics by Thrun, Burgard, and Fox
# """
# function low_var_resample(rng, b::BehaviorParticleBelief, n)
#     for
#     ps = Vector{S}(undef, n)
#     r = rand(rng)*last(cwts)/n
#     c = weight(b,1)
#     i = 1
#     U = r
#     for m in 1:n
#         while U > c
#             i += 1
#             c += weight(b, i)
#         end
#         U += last(cwts)/n
#         ps[m] = particles(b)[i]
#     end
#     return ParticleCollection(ps)
#
# end


text:
module Loggers # Bukdu.Plug

"""
    Plug.Loggers.config

Logging options for System info and error messages.
 - `:controller_pad`
 - `:action_pad`
 - `:path_pad`
 - `:error_stackframes_range`
"""
config = Dict{Symbol, Any}(
    :controller_pad => 20,
    :action_pad => 16,
    :path_pad => 28,
    :error_stackframes_range => :,
)

using ..Deps
using Logging: AbstractLogger

struct DefaultLogger <: AbstractLogger
    stream
    DefaultLogger() = new(IOContext(Core.stdout, :color => have_color()))
end

struct NullLogger <: AbstractLogger
    stream
    NullLogger() = new(IOContext(Core.stdout, :color => have_color()))
end

@generated have_color() = :(2 != Base.JLOptions().color)

current = Dict{Symbol, Union{<:AbstractLogger}}(:logger => NullLogger())

const style_request_action_others  = :red
const style_request_action = Dict{String,Symbol}(
    "GET"     => :normal,
    "POST"    => :yellow,
    "DELETE"  => :magenta,
    "PATCH"   => :light_green,
    "PUT"     => :green,
    "HEAD"    => :light_cyan,
    "OPTIONS" => :cyan,
    "TRACE"   => :light_cyan,
)

const style_response_status_others = :red
const style_response_status = Dict{Int,Symbol}(
    200 => :normal,       # 200 OK
    301 => :light_red,    # 301 Moved Permanently
    302 => :light_red,    # 302 Found
    401 => :magenta,      # 401 Unauthorized
    404 => :light_blue,   # 404 Not Found
    500 => :red,          # 500 Internal Server Error
    503 => :red,          # 503 Service Unavailable
)

function req_method_style(method::String)
    (color=get(style_request_action, method, style_request_action_others),)
end

function resp_status_style(status::Int16)
    (color=get(style_response_status, status, style_response_status_others),)
end

function _regularize_text(str::String, padding::Int)::String
    s = escape_string(str)
    if textwidth(s) < padding
        padded_str = rpad(s, padding)
        if textwidth(padded_str) > padding
        else
            return s
        end
    end
    n = 0
    a = []
    for (idx, x) in enumerate(s)
        n += textwidth(x)
        push!(a, x)
        if n > padding - 1
            break
        end
    end
    newstr = join(a)
    newdiff = textwidth(s) - textwidth(newstr)
    if length(s) == length(a)
        news = newstr
        npad = padding - textwidth(newstr)
    else
        if newdiff > 0 && length(a) >= 2
            newstr = join(a[1:end-2])
            newpad = padding - textwidth(newstr)
            news = string(newstr, fill('.', newpad)...)
        else
            news = newstr
        end
        npad = padding - textwidth(news)
    end
    news
end

function _unescape_request_target(target)
    str = target
    try
        str = Deps.HTTP.URIs.unescapeuri(target)
    catch
    end
    _regularize_text(str, config[:path_pad])
end

function default_print_message(io, args...; kwargs...)
    Base.println(io, args...; kwargs...)
end

function default_print_internal_error(io, name::Symbol, err)
    Base.printstyled(io, name, color=:red)
    r = config[:error_stackframes_range]
    stackframes = err.stackframes[r]
    Base.println(io, " ", err.exception, "\n    ", join(stackframes, "\n    "))
end

function default_info_response(io, conn, route::NamedTuple{(:controller, :action)})
    controller_name = route.controller === nothing ? "" : String(nameof(route.controller))
    action_name = route.action === nothing ? "" : String(nameof(route.action))
    Base.printstyled(io, "INFO:", color=:cyan)
    Base.printstyled(io, ' ')
    Base.printstyled(io, rpad(conn.request.method, 7); req_method_style(conn.request.method)...)
    Base.printstyled(io, ' ')
    controller_color = Sys.iswindows() ? :normal : 248
    if endswith(controller_name, "Controller")
        Base.printstyled(io, controller_name[1:end-10])
        pad_length = config[:controller_pad] - length(controller_name)
        if pad_length > 0
            Base.printstyled(io, "Controller", color=controller_color)
            Base.printstyled(io, repeat(' ', pad_length))
        elseif pad_length < -10
        else
            Base.printstyled(io, "Controller"[1:10+pad_length-1], color=controller_color)
            Base.printstyled(io, ' ')
        end
    else
        Base.printstyled(io, rpad(controller_name, config[:controller_pad]))
    end
    if length(action_name) > config[:action_pad]
        Base.printstyled(io, action_name[1:config[:action_pad]-1], ' ')
    else
        Base.printstyled(io, rpad(action_name, config[:action_pad]))
    end
    Base.printstyled(io, conn.request.response.status; resp_status_style(conn.request.response.status)...)
    Base.printstyled(io, ' ', _unescape_request_target(conn.request.target))
    Base.println(io)
end

function print_message(logger::DefaultLogger, args...; kwargs...)
    io = logger.stream
    default_print_message(io, args...; kwargs...)
end

function print_internal_error(logger::DefaultLogger, name::Symbol, err)
    io = logger.stream
    default_print_internal_error(io, name, err)
end

function info_response(logger::DefaultLogger, conn, route::NamedTuple{(:controller, :action)})
    io = logger.stream
    default_info_response(io, conn, route)
end

function info_response(logger::NullLogger, conn, route::NamedTuple{(:controller, :action)})
end

function print_message(::T, args...; kwargs...) where {T <: AbstractLogger}
    io = IOContext(Core.stdout, :color => have_color())
    default_print_message(io, args...; kwargs...)
end

function print_internal_error(::T, name::Symbol, err) where {T <: AbstractLogger}
    io = IOContext(Core.stdout, :color => have_color())
    default_print_internal_error(io, name, err)
end

function info_response(::T, conn, route::NamedTuple{(:controller, :action)}) where {T <: AbstractLogger}
    io = IOContext(Core.stdout, :color => have_color())
    default_info_response(io, conn, route)
end

function print_message(args...; kwargs...)
    logger = current[:logger]
    print_message(logger, args...; kwargs...)
end

function print_internal_error(name::Symbol, err)
    logger = current[:logger]
    print_internal_error(logger, name, err)
end

function info_response(conn, route::NamedTuple{(:controller, :action)})
    logger = current[:logger]
    info_response(logger, conn, route)
end

end # module Bukdu.Plug.Loggers

"""
    plug(::Type{T}, args...; kwargs...) where {T <: AbstractLogger}
"""
function plug(::Type{T}, args...; kwargs...) where {T <: AbstractLogger}
    Loggers.current[:logger] = T(args...; kwargs...)
end

# module Bukdu.Plug


text:
"""Container for data parsed by PowerModels"""
struct PowerModelsData
    data::Dict{String, Any}
end

"""
Constructs PowerModelsData from a raw file.
Currently Supports MATPOWER and PSSE data files parsed by PowerModels.
"""
function PowerModelsData(file::Union{String, IO}; kwargs...)
    validate = get(kwargs, :pm_data_corrections, true)
    import_all = get(kwargs, :import_all, false)
    pm_dict = parse_file(file; import_all = import_all, validate = validate)
    pm_data = PowerModelsData(pm_dict)
    correct_pm_transformer_status!(pm_data)
    return pm_data
end

"""
Constructs a System from PowerModelsData.

# Arguments
- `pm_data::Union{PowerModelsData, Union{String, IO}}`: PowerModels data object or supported
load flow case (*.m, *.raw)

# Keyword arguments
- `ext::Dict`: Contains user-defined parameters. Should only contain standard types.
- `runchecks::Bool`: Run available checks on input fields and when add_component! is called.
  Throws InvalidValue if an error is found.
- `time_series_in_memory::Bool=false`: Store time series data in memory instead of HDF5.
- `config_path::String`: specify path to validation config file
- `pm_data_corrections::Bool=true` : Run the PowerModels data corrections (aka :validate in PowerModels)
- `import_all:Bool=false` : Import all fields from PTI files

# Examples
```julia
sys = System(
    pm_data, config_path = "ACTIVSg25k_validation.json",
    bus_name_formatter = x->string(x["name"]*"-"*string(x["index"])),
    load_name_formatter = x->strip(join(x["source_id"], "_"))
)
```
"""
function System(pm_data::PowerModelsData; kwargs...)
    runchecks = get(kwargs, :runchecks, true)
    data = pm_data.data
    if length(data["bus"]) < 1
        throw(DataFormatError("There are no buses in this file."))
    end

    @info "Constructing System from Power Models" data["name"] data["source_type"]

    sys = System(data["baseMVA"]; kwargs...)

    bus_number_to_bus = read_bus!(sys, data; kwargs...)
    read_loads!(sys, data, bus_number_to_bus; kwargs...)
    read_loadzones!(sys, data, bus_number_to_bus; kwargs...)
    read_gen!(sys, data, bus_number_to_bus; kwargs...)
    read_branch!(sys, data, bus_number_to_bus; kwargs...)
    read_shunt!(sys, data, bus_number_to_bus; kwargs...)
    read_dcline!(sys, data, bus_number_to_bus; kwargs...)
    read_storage!(sys, data, bus_number_to_bus; kwargs...)
    if runchecks
        check(sys)
    end
    return sys
end

function correct_pm_transformer_status!(pm_data::PowerModelsData)
    for (k, branch) in pm_data.data["branch"]
        if !branch["transformer"] &&
           pm_data.data["bus"][string(branch["f_bus"])]["base_kv"] !=
           pm_data.data["bus"][string(branch["t_bus"])]["base_kv"]
            branch["transformer"] = true
            @warn "Branch $k endpoints have different voltage levels, converting to transformer."
        end
    end
end

"""
Internal component name retreval from pm2ps_dict
"""
function _get_pm_dict_name(device_dict)
    if haskey(device_dict, "name")
        name = device_dict["name"]
    elseif haskey(device_dict, "source_id")
        name = strip(join(string.(device_dict["source_id"]), "-"))
    else
        name = string(device_dict["index"])
    end
    return name
end

"""
Internal branch name retreval from pm2ps_dict
"""
function _get_pm_branch_name(device_dict, bus_f, bus_t)
    index = get(device_dict, "name", device_dict["index"])
    return "$(get_name(bus_f))-$(get_name(bus_t))-i_$index"
end

"""
Creates a PowerSystems.Bus from a PowerSystems bus dictionary
"""
function make_bus(bus_dict::Dict{String, Any})
    bus = Bus(
        bus_dict["number"],
        bus_dict["name"],
        bus_dict["bustype"],
        bus_dict["angle"],
        bus_dict["voltage"],
        bus_dict["voltage_limits"],
        bus_dict["base_voltage"],
        bus_dict["area"],
        bus_dict["zone"],
    )
    return bus
end

function make_bus(bus_name, bus_number, d, bus_types, area)
    bus = make_bus(
        Dict{String, Any}(
            "name" => bus_name,
            "number" => bus_number,
            "bustype" => bus_types[d["bus_type"]],
            "angle" => d["va"],
            "voltage" => d["vm"],
            "voltage_limits" => (min = d["vmin"], max = d["vmax"]),
            "base_voltage" => d["base_kv"],
            "area" => area,
            "zone" => nothing,
        ),
    )
    return bus
end

# "From http://www.pserc.cornell.edu/matpower/MATPOWER-manual.pdf Table B-1"
IS.@scoped_enum(
    MatpowerBusTypes,
    MATPOWER_PQ = 1,
    MATPOWER_PV = 2,
    MATPOWER_REF = 3,
    MATPOWER_ISOLATED = 4,
)

const _BUS_TYPE_MAP = Dict(
    MatpowerBusTypes.MATPOWER_ISOLATED => BusTypes.ISOLATED,
    MatpowerBusTypes.MATPOWER_PQ => BusTypes.PQ,
    MatpowerBusTypes.MATPOWER_PV => BusTypes.PV,
    MatpowerBusTypes.MATPOWER_REF => BusTypes.REF,
)

function Base.convert(::Type{BusTypes}, x::MatpowerBusTypes)
    return _BUS_TYPE_MAP[x]
end

# Disabling this because not all matpower files define areas even when bus definitions
# contain area references.
#function read_area!(sys::System, data; kwargs...)
#    if !haskey(data, "areas")
#        @info "There are no Areas in this file"
#        return
#    end
#
#    for (key, val) in data["areas"]
#        area = Area(string(val["col_1"]))
#        add_component!(sys, area; skip_validation = SKIP_PM_VALIDATION)
#    end
#end

function read_bus!(sys::System, data; kwargs...)
    @info "Reading bus data"
    bus_number_to_bus = Dict{Int, Bus}()

    bus_types = instances(MatpowerBusTypes)
    bus_data = sort!(collect(data["bus"]), by = x -> parse(Int, x[1]))

    if isempty(bus_data)
        @error "No bus data found" # TODO : need for a model without a bus
    end

    _get_name = get(kwargs, :bus_name_formatter, _get_pm_dict_name)
    for (i, (d_key, d)) in enumerate(bus_data)
        # d id the data dict for each bus
        # d_key is bus key
        d["name"] = get(d, "name", string(d["bus_i"]))
        bus_name = strip(_get_name(d))
        bus_number = Int(d["bus_i"])

        area_name = string(d["area"])
        area = get_component(Area, sys, area_name)
        if isnothing(area)
            area = Area(area_name)
            add_component!(sys, area; skip_validation = SKIP_PM_VALIDATION)
        end

        bus = make_bus(bus_name, bus_number, d, bus_types, area)
        has_component(Bus, sys, bus_name) && throw(
            DataFormatError(
                "Found duplicate bus names of $(get_name(bus)), consider formatting names with `bus_name_formatter` kwarg",
            ),
        )

        bus_number_to_bus[bus.number] = bus

        add_component!(sys, bus; skip_validation = SKIP_PM_VALIDATION)
    end

    return bus_number_to_bus
end

function make_load(d, bus, sys_mbase; kwargs...)
    _get_name = get(kwargs, :load_name_formatter, x -> strip(join(x["source_id"])))
    return PowerLoad(;
        name = _get_name(d),
        available = true,
        model = LoadModels.ConstantPower,
        bus = bus,
        active_power = d["pd"],
        reactive_power = d["qd"],
        max_active_power = d["pd"],
        max_reactive_power = d["qd"],
        base_power = sys_mbase,
    )
end

function read_loads!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    if !haskey(data, "load")
        @error "There are no loads in this file"
        return
    end

    sys_mbase = data["baseMVA"]
    for d_key in keys(data["load"])
        d = data["load"][d_key]
        if d["pd"] != 0.0
            bus = bus_number_to_bus[d["load_bus"]]
            load = make_load(d, bus, sys_mbase; kwargs...)
            has_component(PowerLoad, sys, get_name(load)) && throw(
                DataFormatError(
                    "Found duplicate load names of $(get_name(load)), consider formatting names with `load_name_formatter` kwarg",
                ),
            )
            add_component!(sys, load; skip_validation = SKIP_PM_VALIDATION)
        end
    end
end

function make_loadzone(name, active_power, reactive_power; kwargs...)
    _get_name = get(kwargs, :loadzone_name_formatter, _get_pm_dict_name)
    return LoadZone(;
        name = name,
        peak_active_power = sum(active_power),
        peak_reactive_power = sum(reactive_power),
    )
end

function read_loadzones!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    zones = Set{Int}()
    for (i, bus) in data["bus"]
        push!(zones, bus["zone"])
    end

    for zone in zones
        buses = [
            bus_number_to_bus[b["bus_i"]] for b in values(data["bus"]) if b["zone"] == zone
        ]
        bus_names = Set{String}()
        for bus in buses
            push!(bus_names, get_name(bus))
        end

        active_power = Vector{Float64}()
        reactive_power = Vector{Float64}()

        for (key, load) in data["load"]
            load_bus = bus_number_to_bus[load["load_bus"]]
            if get_name(load_bus) in bus_names
                push!(active_power, load["pd"])
                push!(reactive_power, load["qd"])
            end
        end

        load_zone = make_loadzone(string(zone), active_power, reactive_power; kwargs...)
        for bus in buses
            set_load_zone!(bus, load_zone)
        end
        add_component!(sys, load_zone; skip_validation = SKIP_PM_VALIDATION)
    end
end

function make_hydro_gen(gen_name, d, bus, sys_mbase)
    ramp_agc = get(d, "ramp_agc", get(d, "ramp_10", get(d, "ramp_30", abs(d["pmax"]))))
    curtailcost = TwoPartCost(0.0, 0.0)

    base_conversion = sys_mbase / d["mbase"]
    return HydroDispatch( # No way to define storage parameters for gens in PM so can only make HydroDispatch
        name = gen_name,
        available = Bool(d["gen_status"]),
        bus = bus,
        active_power = d["pg"] * base_conversion,
        reactive_power = d["qg"] * base_conversion,
        rating = calculate_rating(d["pmax"], d["qmax"]) * base_conversion,
        prime_mover = parse_enum_mapping(PrimeMovers, d["type"]),
        active_power_limits = (
            min = d["pmin"] * base_conversion,
            max = d["pmax"] * base_conversion,
        ),
        reactive_power_limits = (
            min = d["qmin"] * base_conversion,
            max = d["qmax"] * base_conversion,
        ),
        ramp_limits = (up = ramp_agc, down = ramp_agc),
        time_limits = nothing,
        operation_cost = curtailcost,
        base_power = d["mbase"],
    )
end

function make_renewable_dispatch(gen_name, d, bus, sys_mbase)
    cost = TwoPartCost(0.0, 0.0)
    base_conversion = sys_mbase / d["mbase"]

    rating = calculate_rating(d["pmax"], d["qmax"])
    if rating > d["mbase"]
        @warn "rating is larger than base power for $gen_name, setting to $(d["mbase"])"
        rating = d["mbase"]
    end

    generator = RenewableDispatch(;
        name = gen_name,
        available = Bool(d["gen_status"]),
        bus = bus,
        active_power = d["pg"] * base_conversion,
        reactive_power = d["qg"] * base_conversion,
        rating = rating * base_conversion,
        prime_mover = parse_enum_mapping(PrimeMovers, d["type"]),
        reactive_power_limits = (
            min = d["qmin"] * base_conversion,
            max = d["qmax"] * base_conversion,
        ),
        power_factor = 1.0,
        operation_cost = cost,
        base_power = d["mbase"],
    )

    return generator
end

function make_renewable_fix(gen_name, d, bus, sys_mbase)
    base_conversion = sys_mbase / d["mbase"]
    generator = RenewableFix(;
        name = gen_name,
        available = Bool(d["gen_status"]),
        bus = bus,
        active_power = d["pg"] * base_conversion,
        reactive_power = d["qg"] * base_conversion,
        rating = float(d["pmax"]) * base_conversion,
        prime_mover = parse_enum_mapping(PrimeMovers, d["type"]),
        power_factor = 1.0,
        base_power = d["mbase"],
    )

    return generator
end

function make_generic_battery(storage_name, d, bus)
    storage = GenericBattery(;
        name = storage_name,
        available = Bool(d["status"]),
        bus = bus,
        prime_mover = PrimeMovers.BA,
        initial_energy = d["energy"],
        state_of_charge_limits = (min = 0.0, max = d["energy_rating"]),
        rating = d["thermal_rating"],
        active_power = d["ps"],
        input_active_power_limits = (min = 0.0, max = d["charge_rating"]),
        output_active_power_limits = (min = 0.0, max = d["discharge_rating"]),
        efficiency = (in = d["charge_efficiency"], out = d["discharge_efficiency"]),
        reactive_power = d["qs"],
        reactive_power_limits = (min = d["qmin"], max = d["qmax"]),
        base_power = d["thermal_rating"],
    )
    return storage
end

"""
The polynomial term follows the convention that for an n-degree polynomial, at least n + 1 components are needed.
    c(p) = c_n*p^n+...+c_1p+c_0
    c_o is stored in the  field in of the Econ Struct
"""
function make_thermal_gen(gen_name::AbstractString, d::Dict, bus::Bus, sys_mbase::Number)
    if haskey(d, "model")
        model = GeneratorCostModels(d["model"])
        if model == GeneratorCostModels.PIECEWISE_LINEAR
            cost_component = d["cost"]
            power_p = [i for (ix, i) in enumerate(cost_component) if isodd(ix)]
            cost_p = [i for (ix, i) in enumerate(cost_component) if iseven(ix)]
            cost = [(p, c) for (p, c) in zip(cost_p, power_p)]
            fixed = max(
                0.0,
                cost[1][1] -
                (cost[2][1] - cost[1][1]) / (cost[2][2] - cost[1][2]) * cost[1][2],
            )
            cost = [(c[1] - fixed, c[2]) for c in cost]
        elseif model == GeneratorCostModels.POLYNOMIAL
            if d["ncost"] == 0
                cost = (0.0, 0.0)
                fixed = 0.0
            elseif d["ncost"] == 1
                cost = (0.0, 0.0)
                fixed = d["cost"][1]
            elseif d["ncost"] == 2
                cost = (0.0, d["cost"][1]) ./ sys_mbase
                fixed = d["cost"][2]
            elseif d["ncost"] == 3
                cost = (d["cost"][1] / sys_mbase, d["cost"][2]) ./ sys_mbase
                fixed = d["cost"][3]
            else
                throw(
                    DataFormatError(
                        "invalid value for ncost: $(d["ncost"]). PowerSystems only supports polynomials up to second degree",
                    ),
                )
            end
        end
        startup = d["startup"]
        shutdn = d["shutdown"]
    else
        @warn "Generator cost data not included for Generator: $gen_name"
        tmpcost = ThreePartCost(nothing)
        cost = tmpcost.variable
        fixed = tmpcost.fixed
        startup = tmpcost.start_up
        shutdn = tmpcost.shut_down
    end

    # Ignoring due to  GitHub #148: ramp_agc isn't always present. This value may not be correct.
    ramp_lim = get(d, "ramp_10", get(d, "ramp_30", abs(d["pmax"])))

    operation_cost = ThreePartCost(;
        variable = cost,
        fixed = fixed,
        start_up = startup,
        shut_down = shutdn,
    )

    ext = Dict{String, Any}()
    if haskey(d, "r_source")
        ext["z_source"] = (r = d["r_source"], x = d["x_source"])
    end

    base_conversion = sys_mbase / d["mbase"]
    thermal_gen = ThermalStandard(
        name = gen_name,
        status = Bool(d["gen_status"]),
        available = true,
        bus = bus,
        active_power = d["pg"] * base_conversion,
        reactive_power = d["qg"] * base_conversion,
        rating = calculate_rating(d["pmax"], d["qmax"]) * base_conversion,
        prime_mover = parse_enum_mapping(PrimeMovers, d["type"]),
        fuel = parse_enum_mapping(ThermalFuels, d["fuel"]),
        active_power_limits = (
            min = d["pmin"] * base_conversion,
            max = d["pmax"] * base_conversion,
        ),
        reactive_power_limits = (
            min = d["qmin"] * base_conversion,
            max = d["qmax"] * base_conversion,
        ),
        ramp_limits = (up = ramp_lim, down = ramp_lim),
        time_limits = nothing,
        operation_cost = operation_cost,
        base_power = d["mbase"],
        ext = ext,
    )

    return thermal_gen
end

"""
Transfer generators to ps_dict according to their classification
"""
function read_gen!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading generator data"

    if !haskey(data, "gen")
        @error "There are no Generators in this file"
        return nothing
    end

    generator_mapping = get(kwargs, :generator_mapping, nothing)
    if generator_mapping isa AbstractString || isnothing(generator_mapping)
        generator_mapping = get_generator_mapping(generator_mapping)
    end

    sys_mbase = data["baseMVA"]

    _get_name = get(kwargs, :gen_name_formatter, _get_pm_dict_name)

    for (name, pm_gen) in data["gen"]
        gen_name = _get_name(pm_gen)

        bus = bus_number_to_bus[pm_gen["gen_bus"]]
        pm_gen["fuel"] = get(pm_gen, "fuel", "OTHER")
        pm_gen["type"] = get(pm_gen, "type", "OT")
        @debug "Found generator" _group = IS.LOG_GROUP_PARSING gen_name bus pm_gen["fuel"] pm_gen["type"]

        gen_type = get_generator_type(pm_gen["fuel"], pm_gen["type"], generator_mapping)
        if gen_type == ThermalStandard
            generator = make_thermal_gen(gen_name, pm_gen, bus, sys_mbase)
        elseif gen_type == HydroEnergyReservoir
            generator = make_hydro_gen(gen_name, pm_gen, bus, sys_mbase)
        elseif gen_type == RenewableDispatch
            generator = make_renewable_dispatch(gen_name, pm_gen, bus, sys_mbase)
        elseif gen_type == RenewableFix
            generator = make_renewable_fix(gen_name, pm_gen, bus, sys_mbase)
        elseif gen_type == GenericBattery
            @warn "GenericBattery should be defined as a PowerModels storage... Skipping"
            continue
        else
            @error "Skipping unsupported generator" gen_type
            continue
        end

        has_component(typeof(generator), sys, get_name(generator)) && throw(
            DataFormatError(
                "Found duplicate $(typeof(generator)) names of $(get_name(generator)), consider formatting names with `gen_name_formatter` kwarg",
            ),
        )
        add_component!(sys, generator; skip_validation = SKIP_PM_VALIDATION)
    end
end

function make_branch(name, d, bus_f, bus_t)
    primary_shunt = d["b_fr"]
    alpha = d["shift"]
    branch_type = get_branch_type(d["tap"], alpha, d["transformer"])

    if d["transformer"]
        if branch_type == Line
            throw(DataFormatError("Data is mismatched; this cannot be a line. $d"))
        elseif branch_type == Transformer2W
            value = make_transformer_2w(name, d, bus_f, bus_t)
        elseif branch_type == TapTransformer
            value = make_tap_transformer(name, d, bus_f, bus_t)
        elseif branch_type == PhaseShiftingTransformer
            value = make_phase_shifting_transformer(name, d, bus_f, bus_t, alpha)
        else
            error("Unsupported branch type $branch_type")
        end
    else
        value = make_line(name, d, bus_f, bus_t)
    end

    return value
end

function make_line(name, d, bus_f, bus_t)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)

    return Line(;
        name = name,
        available = d["br_status"] == 1,
        active_power_flow = pf,
        reactive_power_flow = qf,
        arc = Arc(bus_f, bus_t),
        r = d["br_r"],
        x = d["br_x"],
        b = (from = d["b_fr"], to = d["b_to"]),
        rate = d["rate_a"],
        angle_limits = (min = d["angmin"], max = d["angmax"]),
    )
end

function make_transformer_2w(name, d, bus_f, bus_t)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)
    return Transformer2W(;
        name = name,
        available = d["br_status"] == 1,
        active_power_flow = pf,
        reactive_power_flow = qf,
        arc = Arc(bus_f, bus_t),
        r = d["br_r"],
        x = d["br_x"],
        primary_shunt = d["b_fr"],  # TODO: which b ??
        rate = d["rate_a"],
    )
end

function make_tap_transformer(name, d, bus_f, bus_t)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)
    return TapTransformer(;
        name = name,
        available = d["br_status"] == 1,
        active_power_flow = pf,
        reactive_power_flow = qf,
        arc = Arc(bus_f, bus_t),
        r = d["br_r"],
        x = d["br_x"],
        tap = d["tap"],
        primary_shunt = d["b_fr"],  # TODO: which b ??
        rate = d["rate_a"],
    )
end

function make_phase_shifting_transformer(name, d, bus_f, bus_t, alpha)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)
    return PhaseShiftingTransformer(;
        name = name,
        available = d["br_status"] == 1,
        active_power_flow = pf,
        reactive_power_flow = qf,
        arc = Arc(bus_f, bus_t),
        r = d["br_r"],
        x = d["br_x"],
        tap = d["tap"],
        primary_shunt = d["b_fr"],  # TODO: which b ??
        α = alpha,
        rate = d["rate_a"],
    )
end

function read_branch!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading branch data"
    if !haskey(data, "branch")
        @info "There is no Branch data in this file"
        return
    end

    _get_name = get(kwargs, :branch_name_formatter, _get_pm_branch_name)

    for (d_key, d) in data["branch"]
        d["name"] = get(d, "name", d_key)
        bus_f = bus_number_to_bus[d["f_bus"]]
        bus_t = bus_number_to_bus[d["t_bus"]]
        name = _get_name(d, bus_f, bus_t)
        value = make_branch(name, d, bus_f, bus_t)

        add_component!(sys, value; skip_validation = SKIP_PM_VALIDATION)
    end
end

function make_dcline(name, d, bus_f, bus_t)
    return HVDCLine(;
        name = name,
        available = d["br_status"] == 1,
        active_power_flow = get(d, "pf", 0.0),
        arc = Arc(bus_f, bus_t),
        active_power_limits_from = (min = d["pminf"], max = d["pmaxf"]),
        active_power_limits_to = (min = d["pmint"], max = d["pmaxt"]),
        reactive_power_limits_from = (min = d["qminf"], max = d["qmaxf"]),
        reactive_power_limits_to = (min = d["qmint"], max = d["qmaxt"]),
        loss = (l0 = d["loss0"], l1 = d["loss1"]),
    )
end

function read_dcline!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading DC Line data"
    if !haskey(data, "dcline")
        @info "There is no DClines data in this file"
        return
    end

    _get_name = get(kwargs, :branch_name_formatter, _get_pm_branch_name)

    for (d_key, d) in data["dcline"]
        d["name"] = get(d, "name", d_key)
        bus_f = bus_number_to_bus[d["f_bus"]]
        bus_t = bus_number_to_bus[d["t_bus"]]
        name = _get_name(d, bus_f, bus_t)
        dcline = make_dcline(name, d, bus_f, bus_t)
        add_component!(sys, dcline, skip_validation = SKIP_PM_VALIDATION)
    end
end

function make_shunt(name, d, bus)
    return FixedAdmittance(;
        name = name,
        available = Bool(d["status"]),
        bus = bus,
        Y = (-d["gs"] - d["bs"]im),
    )
end

function read_shunt!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading branch data"
    if !haskey(data, "shunt")
        @info "There is no shunt data in this file"
        return
    end

    _get_name = get(kwargs, :shunt_name_formatter, _get_pm_dict_name)

    for (d_key, d) in data["shunt"]
        d["name"] = get(d, "name", d_key)
        name = _get_name(d)
        bus = bus_number_to_bus[d["shunt_bus"]]
        shunt = make_shunt(name, d, bus)

        add_component!(sys, shunt; skip_validation = SKIP_PM_VALIDATION)
    end
end

function read_storage!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading storage data"
    if !haskey(data, "storage")
        @info "There is no storage data in this file"
        return
    end

    _get_name = get(kwargs, :gen_name_formatter, _get_pm_dict_name)

    for (d_key, d) in data["storage"]
        d["name"] = get(d, "name", d_key)
        name = _get_name(d)
        bus = bus_number_to_bus[d["storage_bus"]]
        storage = make_generic_battery(name, d, bus)

        add_component!(sys, storage; skip_validation = SKIP_PM_VALIDATION)
    end
end


text:
using Plots
plotly()

p1 = P2[[9.0, 2.93166], [9.0, 3.66587], [4.19113, 3.27867], [4.20168, 2.97916]]
p2 = P2[[4.07356, 3.05364], [4.19125, 3.27868], [0.0, 2.94463], [0.0, 1.51997], [6.24798, 3.02374], [6.19429, 3.31842]]

x1, y1 = Plots.unzip(p1)
x2, y2 = Plots.unzip(p2)

plot(Shape(x1, y1), fillcolor=:red, label="A")
plot!(Shape(x2, y2), fillcolor=:blue, alpha=0.5, label="B")
xticks!(0.0:0.1:9.0)
yticks!(0.0:0.1:6.0)


text:
function twodBilinear( kernel, ϕ, test, wg )
#= """
  twodBilinear - routine to compute \int{ kernel*ϕ*test }

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    M = twodBilinear(kernel, ϕ, test, wg)
  ```

  Arguments:
  - `kernel`: Kernel function in the integral evaluated at the quadrature points
  - `ϕ`: Matrix of element test functions evaluated at the quadrature points
           (dim: n_gauss, n_dof)
  - `test`: Matrix of test functions evaluated at the quadrature points
            (dim: n_gauss, n_dof)
  - `wg`: Row vector of quadrature weights

  Output argument:
  - `M`: Numerical integral of kernel*phi*test
""" =#

#  M = test'*diagm( wg.*kernel )*ϕ

  wk = wg.*kernel

  M = test'*(LinearAlgebra.Diagonal(wk)*ϕ)

  return M

end


text:
import Base.isless

struct PixelKey{CT, N}
    val::CT
    time_step::Int
    source::CartesianIndex{N}
end
isless(a::PixelKey, b::PixelKey) = (a.val < b.val) || (a.val == b.val && a.time_step < b.time_step)

"""Calculate the euclidean distance between two `CartesianIndex` structs"""
@inline _euclidean(a::CartesianIndex{N}, b::CartesianIndex{N}) where {N} = sqrt(sum(Tuple(a - b) .^ 2))

"""
```
segments                = watershed(img, markers; compactness, mask)
```
Segments the image using watershed transform. Each basin formed by watershed transform corresponds to a segment.
If you are using image local minimas as markers, consider using [`hmin_transform`](@ref) to avoid oversegmentation.

Parameters:
-    img            = input grayscale image
-    markers        = An array (same size as img) with each region's marker assigned a index starting from 1. Zero means not a marker.
                      If two markers have the same index, their regions will be merged into a single region.
                      If you have markers as a boolean array, use `label_components`.
-    compactness    = Use the compact watershed algorithm with the given compactness parameter. Larger values lead to more regularly
                      shaped watershed basins.[^1]
-    mask           = Only segment pixels where the value of `mask` is true, used to restrict segmentation to only areas of interest

[^1]: https://www.tu-chemnitz.de/etit/proaut/publications/cws_pSLIC_ICPR.pdf

# Example

```jldoctest; setup = :(using ImageCore, ImageMorphology, ImageSegmentation)
julia> seeds = falses(100, 100); seeds[50, 25] = true; seeds[50, 75] = true;

julia> dists = distance_transform(feature_transform(seeds)); # calculate distances from seeds

julia> markers = label_components(seeds); # give each seed a unique integer id

julia> results = watershed(dists, markers);

julia> labels_map(results); # labels of segmented image
```
"""
function watershed(img::AbstractArray{T, N},
                   markers::AbstractArray{S,N};
                   mask::AbstractArray{Bool, N}=fill(true, axes(img)),
                   compactness::Real = 0.0) where {T<:ImageCore.NumberLike, S<:Integer, N}

    if axes(img) != axes(markers)
        error("image size doesn't match marker image size")
    elseif axes(img) != axes(mask)
        error("image size doesn't match mask size")
    end

    compact = compactness > 0.0
    segments = copy(markers)
    PK = PixelKey{compact ? floattype(T) : T, N}
    pq = PriorityQueue{CartesianIndex{N}, PK}()
    time_step = 0

    R = CartesianIndices(axes(img))
    Istart, Iend = first(R), last(R)
    for i in R
        if markers[i] > 0
            for j in CartesianIndices(_colon(max(Istart,i-_oneunit(i)), min(i+_oneunit(i),Iend)))
                if segments[j] == 0
                    segments[j] = markers[i]
                    enqueue!(pq, j, PK(img[i], time_step, j))
                    time_step += 1
                end
            end
        end
    end

    while !isempty(pq)
        curr_idx, curr_elem = dequeue_pair!(pq)
        segments_current = segments[curr_idx]

        # If we're using the compact algorithm, we need assign grouping for a given location
        # when it comes off the queue since we could find a better suited watershed later.
        if compact
            if segments_current > 0 && curr_idx != curr_elem.source
                # this is a non-marker location that we've already assigned
                continue
            end
            # group this location with its watershed
            segments[curr_idx] = segments[curr_elem.source]
        end

        img_current = img[curr_idx]
        for j in CartesianIndices(_colon(max(Istart,curr_idx-_oneunit(curr_idx)), min(curr_idx+_oneunit(curr_idx),Iend)))

            # if this location is false in the mask, we skip it
            (!mask[j]) && continue
            # only continue if this is a position that we haven't assigned yet
            if segments[j] == 0
                # if we're doing a simple watershed, we can go ahead and set the final grouping for a new
                # ungrouped position the moment we first encounter it
                if !compact
                    segments[j] = segments_current
                    new_value = img_current
                else
                    # in the compact algorithm case, we don't set the grouping
                    # at push-time and instead calculate a temporary value based
                    # on the weighted sum of the intensity and distance to the
                    # current source marker.
                    new_value = floattype(T)(img_current + compactness * _euclidean(j, curr_elem.source))
                end

                # if this position is in the queue and we're using the compact algorithm, we need to replace
                # its watershed if we find one that it better belongs to
                if compact && j in keys(pq)
                    elem = pq[j]
                    new_elem = PK(new_value, time_step, curr_elem.source)

                    if new_elem < elem
                        pq[j] = new_elem # update the watershed
                        time_step += 1
                    end
                else

                    pq[j] = PK(new_value, time_step, curr_elem.source)
                    time_step += 1
                end
            end
        end
    end

    TM = meantype(T)
    num_segments        = Int(maximum(segments))
    labels              = Array(1:num_segments)
    region_means        = Dict{Int, TM}()
    region_pix_count    = Dict{Int, Int}()

    for i in R
        region_pix_count[segments[i]] = get(region_pix_count, segments[i], 0) + 1
        region_means[segments[i]] = get(region_means, segments[i], zero(TM)) + (img[i] - get(region_means, segments[i], zero(TM)))/region_pix_count[segments[i]]
    end
    return SegmentedImage(segments, labels, region_means, region_pix_count)
end

"""
```
out = hmin_transform(img, h)
```
Suppresses all minima in grayscale image whose depth is less than h.

H-minima transform is defined as the reconstruction by erosion of (img + h) by img. See Morphological image analysis by Soille pg 170-172.
"""
function hmin_transform(img::Array{T, N}, h::Real) where {T<:ImageCore.NumberLike, N}
    out = img.+h
    while true
        temp = max.(img, erode(out))
        if temp == out
            break
        end
        out = temp
    end
    return out
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function testchebx(N::Int)
	x  = chebx(5, N)
	Tx = cos(N*acos(x))
	return abs(Tx)
end

function testchebd(N::Int)
	DR = Float64[chebd(fld(N, 3), j, N) for j in 1:N+1]
	return sum(DR)
end

function testchebdpoly(N::Int)::Float64
    x     = chebgrid(N)
    fx    = x.^5
    edfx  = 5.*x.^4
    chebD = Float64[chebd(i,j,N) for i in 1:N+1, j in 1:N+1]
    ndfx  = chebD*fx
    return maximum(abs.(ndfx - edfx))
end

function testchebw(N::Int)
	w = Float64[chebw(i, N) for i in 1:N+1]
	return sum(w)
end

@test chebx(3, 10) ≈ 0.8090169943749475
@test chebw(9, 12) ≈ 0.2258075258075258
@test vandermonde(10)[2,5] == cheb(4, chebx(2,10))
@test pseudovandermonde(10, chebgrid(14))[3,7] == cheb(6, chebx(3,14))
@test chebgrid(23)[7] == chebx(7,23)
@test chebweights(12)[8] == chebw(8,12)

@test testchebx(rand(10:20)) ≈ 1.0
@test testchebd(rand(10:20)) < 1e-14
@test testchebw(rand(2:20)) ≈ 2.0
@test chebgrid(4, 2, 1) == (chebgrid(4) + 1)/2
@test chebgrid(4, 2, 2) == (chebgrid(4) - 1)/2
@test chebgrid(2, 2, 1)[end] == chebgrid(2, 2, 2)[1]
@test testchebdpoly(10) < 1e-13


text:
module JSOSolverTemplate

# standard lib
using LinearAlgebra

# JSO packages
using NLPModels
using SolverTools

include("uncsolver.jl")

end # module


text:
module GeodesicComplexes
export GeodesicComplex, npoints, nlandmarks, points, landmarks, nonlandmarks

using Distances
using LightGraphs
using NearestNeighbors
using RecipesBase
using SimpleWeightedGraphs
using StaticArrays

using Random
using SparseArrays

struct GeodesicComplex{T, P<:AbstractVector{T}, M<:Metric, K<:NNTree{P, M}} <:
                       AbstractSimpleWeightedGraph{Int, T}
    points        ::Vector{P}
    landmark_idxs ::Vector{Int}
    graph         ::SimpleWeightedGraph{Int, T}
    tree          ::K
    cover         ::Vector{Set{Int}}
    radius        ::T
    metric        ::M
end

function GeodesicComplex(pts::AbstractVector{<:AbstractVector{T}}, r; landmarks = nothing,
                         metric = Euclidean(), tree = KDTree) where {T}
    kdt = tree(pts, metric)
    landmark_idxs, cover = getcover(landmarks, pts, r, kdt)
    I = Int[]
    J = Int[]
    V = T[]
    n = length(landmark_idxs)

    for i in 1:n, j in 1:(i - 1)
        if !isempty(intersect(cover[i], cover[j]))
            d = evaluate(metric, pts[landmark_idxs[i]], pts[landmark_idxs[j]])
            append!(I, (i, j))
            append!(J, (j, i))
            append!(V, (d, d))
        end
    end
    graph = SimpleWeightedGraph(sparse(I, J, V, n, n))

    GeodesicComplex(pts, landmark_idxs, graph, kdt, cover, T(r), metric)
end

function getcover(::Nothing, pts, r, tree)
    covered = falses(length(pts))
    idxs = shuffle(eachindex(pts))
    landmarks = Int[]
    cover = Set{Int}[]

    for i in idxs
        covered[i] && continue
        covered[i] = true
        inball = inrange(tree, pts[i], r)
        covered[inball] .= true
        push!(landmarks, i)
        push!(cover, Set(inball))
    end
    ord = sortperm(landmarks)
    landmarks[ord], cover[ord]
end

function getcover(landmarks, pts, r, tree)
    cover = Set{Int}[]

    for l in landmarks
        inball = inrange(tree, pts[l], r)
        push!(cover, Set(inball))
    end
    ord = sortperm(landmarks)
    landmarks[ord], cover[ord]
end

Base.eltype(::GeodesicComplex{<:Any, P}) where {P} = P

# LightGraphs interface
for f in [:edges, :ne, :nv, :vertices, :is_directed, :weights]
    @eval LightGraphs.$f(gc::GeodesicComplex) = $f(gc.graph)
end
LightGraphs.is_directed(::Type{GeodesicComplex}) = false
LightGraphs.is_directed(::Type{GeodesicComplex{T,P,M,K}}) where {T,P,M,K} = false
LightGraphs.has_edge(gc::GeodesicComplex, u::Integer, v::Integer) = has_edge(gc.graph, u, v)
LightGraphs.has_vertex(gc::GeodesicComplex, v::Integer) = has_edge(gc.graph, v)
LightGraphs.inneighbors(gc::GeodesicComplex, v::Integer) = inneighbors(gc.graph, v)
LightGraphs.outneighbors(gc::GeodesicComplex, v::Integer) = outneighbors(gc.graph, v)

function Base.show(io::IO, gc::GeodesicComplex{T, P}) where {T, P}
    print(io, "GeodesicComplex{$T, $P} with $(npoints(gc)) points, " *
          "$(nlandmarks(gc)) landmarks and $(ne(gc)) edges")
end

npoints(gc::GeodesicComplex) =
    length(gc.points)

nlandmarks(gc::GeodesicComplex) =
    length(gc.landmark_idxs)

points(gc::GeodesicComplex) =
    gc.points
points(gc::GeodesicComplex, idxs) =
    gc.points[idxs]

landmarks(gc::GeodesicComplex, idxs=1:nlandmarks(gc)) =
    gc.points[gc.landmark_idxs[idxs]]

landmark_idxs(gc::GeodesicComplex) =
    gc.landmark_idxs
landmark_idxs(gc::GeodesicComplex, idxs) =
    gc.landmark_idxs[idxs]

nonlandmarks(gc::GeodesicComplex, idxs=1:npoints(gc)-nlandmarks(gc)) =
    gc.points[setdiff(1:npoints(gc), gc.landmark_idxs)[idxs]]

# Plots recipe
function getxyz(pts)
    xs = get.(pts, 1, 0.0)
    ys = get.(pts, 2, 0.0)
    zs = get.(pts, 3, 0.0)
    if all(iszero, zs)
        xs, ys
    else
        xs, ys, zs
    end
end

@recipe function plot(gc::GeodesicComplex;
                      only_landmarks = true, graph = true)
    # edges
    if graph
        @series begin
            label := "edges"
            edgepoints = Vector{Float64}[]
            for e in edges(gc)
                s = landmarks(gc, src(e))
                d = landmarks(gc, dst(e))
                append!(edgepoints, (s, d, fill(NaN, length(s))))
            end
            getxyz(edgepoints)
        end
    end
    # landmarks
    @series begin
        markersize --> 1.0
        seriestype --> :scatter
        label := "landmarks"
        getxyz(landmarks(gc))
    end
    # others
    if !only_landmarks
        @series begin
            markersize --> 0.5
            seriestype --> :scatter
            label := "others"
            getxyz(points(gc, setdiff(1:npoints(gc), landmark_idxs(gc))))
        end
    end
end

end


text:
using Base.Meta: isexpr
using PyCall
using PyCall: pyjlwrap_new
import Conda

const _getproperty = try
    pyimport("sys").executable
    getproperty
catch
    getindex
end
if _getproperty === getindex
    _setproperty!(value, name, x) = setindex!(value, x, name)
else
    const _setproperty! = setproperty!
end

compatattr_imp(x) = x
compatattr_imp(ex::Expr) =
    if ex.head == :. && length(ex.args) == 2 && !isexpr(ex.args[2], :tuple)
        :($_getproperty($(compatattr_imp.(ex.args)...)))
    else
        Expr(ex.head, compatattr_imp.(ex.args)...)
    end

macro compatattr(ex)
    esc(compatattr_imp(ex))
end

julia_exepath() =
    joinpath(VERSION < v"0.7.0-DEV.3073" ? JULIA_HOME : Base.Sys.BINDIR,
             Base.julia_exename())

_start_ipython(name::Symbol) = _start_ipython(getipythonjl(name), PyAny)
# helper function that is used only in tests

function _start_ipython(f::PyObject,
                        returntype = Nothing)
    return pycall(
        f, returntype,
        pyfunctionret(JuliaAPI.eval_str, PyObject, String),
        pyjlwrap_new(JuliaAPI),
    )
end

start_ipython() = _start_ipython(_customized_ipython)

getipythonjl(name) = _getproperty(pyimport("ipython_jl"), name)

const _customized_ipython = PyNULL()

function __init__()
    pushfirst!(PyVector(@compatattr pyimport("sys")."path"), @__DIR__)
    copy!(_customized_ipython, getipythonjl(:customized_ipython))
    afterreplinit(init_repl)
end


text:
# === PRIVATE METHODS THAT ARE NOT EXPORTED ========================================================== #
# ==================================================================================================== #

# === PUBLIC METHODS THAT ARE EXPORTED =============================================================== #
"""
    build_user_model_object(configDictionary::Dict{String,Any})::PSResult
"""
function build_user_model_object(configDictionary::Dict{String,Any})::PSResult

    # initialize -
    # ...

    try

        loginCredentiallDictionary = configDictionary["api"]

        # check for some keys ...
        # username -
        if (haskey(loginCredentiallDictionary,"username") == false)
            throw(ArgumentError("Configuration dictionary error: missing key username"))
        end

        # apikey -
        if (haskey(loginCredentiallDictionary,"apikey") == false)
            throw(ArgumentError("Configuration dictionary error: missing key apikey"))
        end

        # label -
        if (haskey(loginCredentiallDictionary,"label") == false)
            throw(ArgumentError("Configuration dictionary error: missing key label"))
        end

        # ok, we should have all the data to create the user object -
        userName = loginCredentiallDictionary["username"]
        apiKey = loginCredentiallDictionary["apikey"]
        label = loginCredentiallDictionary["label"]

        # build -
        user_object = PSTextMagicAPIUserObject(userName, apiKey, label)

        # return -
        return PSResult(user_object)
    catch error
        return PSResult(error)
    end
end
# ==================================================================================================== #


text:
"""
    Bounded(block, size) <: WrapperBlock

A [`WrapperBlock`](#) for annotating spatial data blocks with
size information for their spatial bounds. As an example,
`Image{2}()` doesn't carry any size information since it
supports variable-size images, but sometimes it can be
useful to have the exact size as information where it can
be known.

Encoding using [`ProjectiveTransforms`](#) returns `Bounded`s since it
crops any input to the same size.

## Examples

```julia
block = Image{2}()  # a 2D-image block with out size informatio
wrapper = Bounded(Image{2}(), (128, 128))  # a 2D-image block with fixed size

@test checkblock(block, rand(10, 10))
@test !checkblock(wrapper, rand(10, 10))  # Expects size `(128, 128)`
```

Wrapping a `Bounded` into another `Bounded` with the same dimensionality
will update the bounds:

```julia
block = Image{2}()
Bounded(Bounded(block, (16, 16)), (8, 8)) == Bounded(block, (8, 8))
```

"""
struct Bounded{N, B<:AbstractBlock} <: WrapperBlock
    block::B
    size::NTuple{N, Int}
end

function Bounded(bounded::Bounded{M}, size::NTuple{N, Int}) where {N, M}
    N == M || error("Cannot rewrap a `Bounded` with different dimensionalities $N and $M")
    Bounded(wrapped(bounded), size)
end

InlineTest.@testset "Bounded [block, wrapper]" begin
    @test_nowarn Bounded(Image{2}(), (16, 16))
    bounded = Bounded(Image{2}(), (16, 16))
    @test Bounded(bounded, (16, 16)) == bounded
end


text:
using BioSequences

#AG
println("running program")

seq = randdnaseq(50) # random seed appears fixed
println(seq)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
# fold(l|r) & mapfold(l|r)
@test foldl(-, 1:5) == -13
@test foldl(-, 10, 1:5) == -5

@test Base.mapfoldl(abs2, -, 2:5) == -46
@test Base.mapfoldl(abs2, -, 10, 2:5) == -44

@test foldr(-, 1:5) == 3
@test foldr(-, 10, 1:5) == -7

@test Base.mapfoldr(abs2, -, 2:5) == -14
@test Base.mapfoldr(abs2, -, 10, 2:5) == -4

# reduce & mapreduce
@test reduce((x,y)->"($x+$y)", [9:11]) == "((9+10)+11)"
@test reduce(max, [8 6 7 5 3 0 9]) == 9
@test reduce(+, 1000, [1:5]) == (1000 + 1 + 2 + 3 + 4 + 5)

@test mapreduce(-, +, [-10 -9 -3]) == ((10 + 9) + 3)
@test mapreduce((x)->x[1:3], (x,y)->"($x+$y)", ["abcd", "efgh", "01234"]) == "((abc+efg)+012)"

# sum

@test sum(Int8[]) === 0
@test sum(Int[]) === 0
@test sum(Float64[]) === 0.0

@test sum(int8(3)) === int8(3)
@test sum(3) === 3
@test sum(3.0) === 3.0

@test sum([int8(3)]) === 3
@test sum([3]) === 3
@test sum([3.0]) === 3.0

z = reshape(1:16, (2,2,2,2))
fz = float(z)
@test sum(z) === 136
@test sum(fz) === 136.0

@test_throws ErrorException sum(sin, Int[])
@test sum(sin, 3) == sin(3.0)
@test sum(sin, [3]) == sin(3.0)
@test sum(sin, z) == sum(sin, fz) == sum(sin(fz))

z = [-4, -3, 2, 5]
fz = float(z)
@test Base.sumabs(Float64[]) === 0.0
@test Base.sumabs([int8(-2)]) === 2
@test Base.sumabs(z) === 14
@test Base.sumabs(fz) === 14.0

@test Base.sumabs2(Float64[]) === 0.0
@test Base.sumabs2([int8(-2)]) === 4
@test Base.sumabs2(z) === 54
@test Base.sumabs2(fz) === 54.0

# check variants of summation for type-stability and other issues (#6069)
sum2(itr) = invoke(sum, (Any,), itr)
plus(x,y) = x + y
sum3(A) = reduce(plus, A)
sum4(itr) = invoke(reduce, (Function, Any), plus, itr)
sum5(A) = reduce(plus, 0, A)
sum6(itr) = invoke(reduce, (Function, Int, Any), plus, 0, itr)
sum7(A) = mapreduce(x->x, plus, A)
sum8(itr) = invoke(mapreduce, (Function, Function, Any), x->x, plus, itr)
sum9(A) = mapreduce(x->x, plus, 0, A)
sum10(itr) = invoke(mapreduce, (Function, Function, Int, Any), x->x,plus,0,itr)
for f in (sum2, sum5, sum6, sum9, sum10)
    @test sum(z) == f(z)
    @test sum(Int[]) == f(Int[]) == 0
    @test sum(Int[7]) == f(Int[7]) == 7
    @test typeof(f(Int8[])) == typeof(f(Int8[1])) == typeof(f(Int8[1 7]))
end
for f in (sum3, sum4, sum7, sum8)
    @test sum(z) == f(z)
    @test_throws ErrorException f(Int[])
    @test sum(Int[7]) == f(Int[7]) == 7
end
@test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7]))

@test sum_kbn([1,1e100,1,-1e100]) == 2

# prod

prod(Int[]) === 0
prod(Int8[]) === 0
prod(Float64[]) === 0.0

prod([3]) === 0
prod([int8(3)]) === 0
prod([3.0]) === 0.0

prod(z) === 120
prod(fz) === 120.0

# check type-stability
prod2(itr) = invoke(prod, (Any,), itr)
@test prod(Int[]) === prod2(Int[]) === 1
@test prod(Int[7]) === prod2(Int[7]) === 7
@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1, 7])) == Int
@test typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7])) == Int

# maximum & minimum & extrema

@test_throws ErrorException maximum(Int[])
@test_throws ErrorException minimum(Int[])

@test maximum(5) == 5
@test minimum(5) == 5
@test extrema(5) == (5, 5)

@test maximum([4, 3, 5, 2]) == 5
@test minimum([4, 3, 5, 2]) == 2
@test extrema([4, 3, 5, 2]) == (2, 5)

@test isnan(maximum([NaN]))
@test isnan(minimum([NaN]))
@test isequal(extrema([NaN]), (NaN, NaN))

@test maximum([NaN, 2., 3.]) == 3.
@test minimum([NaN, 2., 3.]) == 2.
@test extrema([NaN, 2., 3.]) == (2., 3.)

@test maximum([4., 3., NaN, 5., 2.]) == 5.
@test minimum([4., 3., NaN, 5., 2.]) == 2.
@test extrema([4., 3., NaN, 5., 2.]) == (2., 5.)

@test Base.maxabs(Int[]) == 0
@test_throws ErrorException Base.minabs(Int[])

@test Base.maxabs(-2) == 2
@test Base.minabs(-2) == 2
@test Base.maxabs([1, -2, 3, -4]) == 4
@test Base.minabs([-1, 2, -3, 4]) == 1

@test maximum(x->abs2(x), 3:7) == 49
@test minimum(x->abs2(x), 3:7) == 9

# any & all

@test any(Bool[]) == false
@test any([true]) == true
@test any([false, false]) == false
@test any([false, true]) == true
@test any([true, false]) == true
@test any([true, true]) == true
@test any([true, true, true]) == true
@test any([true, false, true]) == true
@test any([false, false, false]) == false

@test all(Bool[]) == true
@test all([true]) == true
@test all([false, false]) == false
@test all([false, true]) == false
@test all([true, false]) == false
@test all([true, true]) == true
@test all([true, true, true]) == true
@test all([true, false, true]) == false
@test all([false, false, false]) == false

@test any(x->x>0, Int[]) == false
@test any(x->x>0, [-3]) == false
@test any(x->x>0, [4]) == true
@test any(x->x>0, [-3, 4, 5]) == true

@test all(x->x>0, Int[]) == true
@test all(x->x>0, [-3]) == false
@test all(x->x>0, [4]) == true
@test all(x->x>0, [-3, 4, 5]) == false

# in

@test in(1, Int[]) == false
@test in(1, Int[1]) == true
@test in(1, Int[2]) == false
@test in(0, 1:3) == false
@test in(1, 1:3) == true
@test in(2, 1:3) == true

# count & countnz

@test count(x->x>0, Int[]) == 0
@test count(x->x>0, -3:5) == 5

@test countnz(Int[]) == 0
@test countnz(Int[0]) == 0
@test countnz(Int[1]) == 1
@test countnz([1, 0, 2, 0, 3, 0, 4]) == 4

## cumsum, cummin, cummax

z = rand(10^6)
let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5])
    @test (es - sum(z)) < es * 1e-13
    cs = cumsum(z)
    @test (es - cs[end]) < es * 1e-13
    @test (es2 - cs[10^5]) < es2 * 1e-13
end

@test isequal(cummin([1, 2, 5, -1, 3, -2]), [1, 1, 1, -1, -1, -2])
@test isequal(cummax([1, 2, 5, -1, 3, -2]), [1, 2, 5, 5, 5, 5])

@test isequal(cummax([1 0; 0 1], 1), [1 0; 1 1])
@test isequal(cummax([1 0; 0 1], 2), [1 1; 0 1])
@test isequal(cummin([1 0; 0 1], 1), [1 0; 0 0])
@test isequal(cummin([1 0; 0 1], 2), [1 0; 0 0])

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# The node order depends on the order of keys in a dict, thus it is not determined and might vary across operating systems
# TODO: this test needs to check the node order and replace the connectivity correspondingly

# @testset "grid refinement" begin
#     #### test grid refinement against Figure 6 from https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620200412

#     # pvd = paraview_collection("test")
#     # test with grid from paper
#     nodes = Node.([(0.0, 0.0), (2., 1.), (-1.5, 2.2), (2.4, 2.1), (2.3, 4.0), (0.3, 4.2), (-1.1, 4.1)])
#     cells = Triangle.([(1,2,3), (2,4,3), (3,4,6), (4,5,6), (3,6,7)])
#     grid = Grid(cells, nodes; topology=Ferrite.ExclusiveTopology(cells))
#     # pvd[0] = vtk_grid("initial_grid", grid)

#     cells_to_split = [5]
#     grid_1 = refine(grid, cells_to_split)
#     @test grid_1.cells == [  Triangle((1, 2, 3)),
#                             Triangle((4, 5, 6)),
#                             Triangle((4, 6, 9)),
#                             Triangle((6, 8, 9)),
#                             Triangle((8, 3, 9)),
#                             Triangle((2, 4, 9)),
#                             Triangle((3, 2, 9)),
#                             Triangle((6, 7, 8)),
#                             Triangle((7, 3, 8))]
#     # pvd[1] = vtk_grid("1st_grid", grid_1)

#     cells_to_split = [8,9]
#     grid_2 = refine(grid_1, cells_to_split)
#     @test grid_2.cells == [ Triangle((1, 2, 3)),
#                             Triangle((4, 5, 6)),
#                             Triangle((4, 6, 9)),
#                             Triangle((6, 8, 9)),
#                             Triangle((8, 3, 9)),
#                             Triangle((2, 4, 9)),
#                             Triangle((3, 2, 9)),
#                             Triangle((7, 8, 10)),
#                             Triangle((8, 6, 10)),
#                             Triangle((3, 8, 11)),
#                             Triangle((8, 7, 11))]
#     # pvd[2] = vtk_grid("2nd_grid", grid_2)

#     cells_to_split = [8,9,10,11]
#     grid_3 = refine(grid_2, cells_to_split)
#     @test grid_3.cells == [ Triangle((2, 4, 9)),
#                             Triangle((2, 9, 17)),
#                             Triangle((9, 16, 17)),
#                             Triangle((16, 3, 17)),
#                             Triangle((6, 10, 13)),
#                             Triangle((10, 8, 13)),
#                             Triangle((4, 5, 18)),
#                             Triangle((5, 6, 18)),
#                             Triangle((8, 10, 12)),
#                             Triangle((10, 7, 12)),
#                             Triangle((1, 2, 17)),
#                             Triangle((3, 1, 17)),
#                             Triangle((6, 13, 15)),
#                             Triangle((13, 8, 15)),
#                             Triangle((8, 9, 15)),
#                             Triangle((7, 11, 12)),
#                             Triangle((11, 8, 12)),
#                             Triangle((8, 11, 14)),
#                             Triangle((11, 3, 14)),
#                             Triangle((8, 14, 16)),
#                             Triangle((14, 3, 16)),
#                             Triangle((9, 8, 16)),
#                             Triangle((6, 15, 18)),
#                             Triangle((15, 9, 18)),
#                             Triangle((9, 4, 18))]
#     # pvd[3] = vtk_grid("3rd_grid", grid_3)

#     # vtk_save(pvd)
# end

@testset "linear to quadratic interpolation" begin

    grid = generate_grid(Triangle, (1,1))

    dh_lin = DofHandler(grid)
    ip_u = Lagrange{2,RefTetrahedron,1}()
    push!(dh_lin, :u, 2, ip_u)
    close!(dh_lin)

    dh_quad = DofHandler(grid)
    ip_u = Lagrange{2,RefTetrahedron,2}()
    push!(dh_quad, :u, 2, ip_u)
    close!(dh_quad)

    a_lin = [1., 1., 2., 2., 2., 2., 4., 4.]
    a_quad = linear_to_quadratic(dh_lin, dh_quad, a_lin)

    @test a_quad == [1., 1., 2., 2., 2., 2., 1.5, 1.5, 2., 2., 1.5, 1.5, 4., 4., 3., 3., 3., 3.]
end


text:
import DSGE: Measurement, get_setting, n_observables, n_shocks_exogenous

"""
```
measurement{T<:AbstractFloat}(m::BondLabor{T}, TTT::Matrix{T}, RRR::Matrix{T},
                              CCC::Vector{T})
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::BondLabor{T}, TTT::Matrix{T},
                     TTT_jump::Matrix{T},
                     RRR::Matrix{T}, CCC::Vector{T}) where {T<:AbstractFloat}
    endo      = m.endogenous_states_unnormalized
    exo       = m.exogenous_shocks
    obs       = m.observables

    _n_model_states = n_model_states(m)
    _n_states       = n_backward_looking_states(m)
    _n_jumps        = n_jumps(m)

    _n_observables = n_observables(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_model_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)

    ν   = m[:ν]
    γ   = m[:γ]
    abar = m[:abar]
    R    = m[:R]

    l   = m[:lstar].value
    μ   = m[:μstar].value
    c   = m[:cstar].value
    η   = m[:ηstar].value
    χss = m[:χstar].value

    xgrid = m.grids[:xgrid].points
    xgrid_total = m.grids[:xgrid_total]
    xwts  = m.grids[:xgrid].weights
    sgrid = m.grids[:sgrid].points
    sgrid_total = m.grids[:sgrid_total]
    swts  = m.grids[:sgrid].weights
    weights_total = m.grids[:weights_total]

    aborrow = abar/R
    chipW = (1+ν)./(ν./(aborrow + χss - xgrid_total) + γ./χss)
    chipR = (1/(R*R))*(ν*abar./(aborrow + χss - xgrid_total))./(ν./(aborrow + χss - xgrid_total) + γ./χss)

    GDP  = 0.0
    GDPZ = 0.0
    GDPR = 0.0

    GDPfn = zeros(1,2*ns*nx +2) # GDP as function of un-normalized MU Z l R
    for ix =1:nx
        for is=1:ns
            i = ix + nx*(is-1)
            GDP += μ[i]*swts[is]*xwts[ix]*η[i]
            GDPZ += sgrid[is]*μ[i]*swts[is]*xwts[ix]*(η[i]/ν - ((γ*η[i])/(ν*c[i]))*chipW[i]*((l[i]^(-1/γ))>χss[i]) )
            GDPR += -(sgrid[is]*μ[i]*γ*η[i]/(ν*c[i]))*swts[is]*xwts[ix]*chipR[i]*((l[i]^(-1/γ))>χss[i])
        end
    end

    # Output in log levels
    GDPfn[1, endo[:μ′_t]] = η'.*sgrid_total'.*weights_total'
    GDPfn[1, endo[:z′_t]] = GDPZ + GDP
    GDPfn[1, endo[:l′_t]] = (1/γ)*μ'.*(sgrid_total'*(γ/ν).*(η./c)').*weights_total'.*( l.^(-1-1/γ) )'.*((l.^(-1/γ)).<=χss)'
    GDPfn[1, endo[:R′_t]] = GDPR

    Qx, Qy, _, _ = compose_normalization_matrices(m)
    gx2  = Qy'*TTT_jump*Qx

    # now we need to create GDP as a function of the normalized states
    ZZ_states = (1/GDP)*GDPfn*[eye(nx*ns+1); gx2]*Qx' # this is for log GDP
                                                      # to use the level of gdp, remove (1/GDP)
    ZZ = Matrix{Float64}(_n_observables, _n_model_states)
    ZZ[1:_n_states] = ZZ_states
    ZZ[_n_states+1:end] = zeros(_n_jumps)

    # Measurement error
    EE[obs[:obs_gdp], obs[:obs_gdp]] = m[:e_y]

    # Variance of innovations
    QQ[exo[:z_sh], exo[:z_sh]] = m[:σ_z]^2

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0)
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
using StatisticalRethinking, StructuralCausalModels
using StatsPlots

ProjDir = @__DIR__

println()
N = 100
df = DataFrame();
df[!, :A] = rand(Normal(), N)
df[!, :M] = [rand(Normal(-a), 1)[1] for a in df[:, :A]]
df[!, :D] = [rand(Normal(a), 1)[1] for a in df[:, :A]]

StatsPlots.cornerplot(Array(df), label=names(df))
savefig("$(ProjDir)/AMD_df_1.png")

println()

for i in 1:3
  include(scm_path("..", "examples", "SR", "fork", "AMD_m$i.jl"))
end

if success(rc)

  r1 = plotcoef([m5_1s, m5_2s, m5_3s], [:bA, :bM], "$(ProjDir)/AMD_fork_1.png",
    "Particles (Normal) estimates")
  display(r1)

end

d = OrderedDict(
  :M => [:A],
  :D => [:A]
);

dag = DAG("fork_correlated", d, df=df);
show(dag)

display(dag.s); println()

shipley_test_dag_3 = shipley_test(dag)
if !isnothing(shipley_test_dag_3)
  display(shipley_test_dag_3)
end
println()

f = [:M]; s = [:D]; sel = vcat(f, s)
c = [:A]

e = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) = $e")

e = d_separation(dag, f, s, c=c)
println("d_separation($(dag.name), $f, $s, c=$c) = $e")
println()
N = 100
df = DataFrame();
df[!, :A] = rand(Normal(), N)
df[!, :M] = [rand(Normal(-a), 1)[1] for a in df[:, :A]]
df[!, :D] = [rand(Normal(a), 1)[1] for a in df[:, :A]] -
  [rand(Normal(m), 1)[1] for m in df[:, :M]]

StatsPlots.cornerplot(Array(df), label=names(df))
savefig("$(ProjDir)/AMD_df_2.png")

println()

for i in 1:3
  include(scm_path("..", "examples", "SR", "fork", "AMD_m$i.jl"))
end

if success(rc)

  r1 = plotcoef([m5_1s, m5_2s, m5_3s], [:bA, :bM], "$(ProjDir)/AMD_fork_2.png",
    "Particles (Normal) estimates")
  display(r1)

end

d = OrderedDict(
  :M => [:A],
  :D => [:A]
);

dag = DAG("fork_correlated", d, df=df);
show(dag)

display(dag.s); println()

shipley_test_dag_3 = shipley_test(dag)
if !isnothing(shipley_test_dag_3)
  display(shipley_test_dag_3)
end
println()

f = :M; s = :D; sel = vcat(f, s)
c = :A

e = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) = $e")

e = d_separation(dag, f, s, c=c)
println("d_separation($(dag.name), $f, $s, c=$c) = $e")

as = adjustment_sets(dag, :A, :D)
as |> display


text:
"""
Original publication:
L. Yin, H. Wang, and W. Fan. Active learning based
support vector data description method for robust
novelty detection. Knowledge-Based Systems, pages
40–52, Aug. 2018.
"""
struct BoundaryNeighborCombinationQs <: HybridQs
    occ::OCClassifier
    nn_dist::Array{Float64,1}
    η::Float64
    p::Float64
    BoundaryNeighborCombinationQs(occ::OCClassifier, x::Array{T, 2}; η=0.7, p=0.15) where T <: Real = new(occ, knn_mean_dist(x; k=1), η, p)
end

function qs_score(qs::BoundaryNeighborCombinationQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    @assert length(qs.nn_dist) == size(x, 2)
    if rand() < qs.p
        return rand(size(x,2))
    end
    prediction = abs.(predict(qs.occ, x))
    d_hyp = (prediction .- minimum(prediction[labels[:U]])) ./ maximum(prediction[labels[:U]])
    d_nn = (qs.nn_dist .- minimum(qs.nn_dist[labels[:U]])) ./ maximum(qs.nn_dist[labels[:U]])
    return qs.η * (-d_hyp) + (1 - qs.η) * (-d_nn)
end


text:
#! /usr/bin/julia

# Rosetta Code, "Sort an array of composite structures"

type Pair{T<:String}
    name::T
    value::T
end

Base.show(a::Pair) = @sprintf "%s => %s" a.name a.value

x = Pair[Pair("gold", "shiny"),
         Pair("neon", "inert"),
         Pair("sulphur", "yellow"),
         Pair("iron", "magnetic"),
         Pair("zebra", "striped"),
         Pair("star", "brilliant"),
         Pair("apple", "tasty"),
         Pair("ruby", "red"),
         Pair("dice", "random"),
         Pair("coffee", "stimulating"),
         Pair("book", "interesting")]

println("x, the original list of pairs: ")
for p in x
    println("    ", show(p))
end

println()
println("x sorted by name:")
x = sort(x, by=a->a.name)
for p in x
    println("    ", show(p))
end

println()
println("x sorted by value:")
x = sort(x, by=a->a.value)
for p in x
    println("    ", show(p))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using jlmie
using Plots

# settings
# structure
nmat = 4  # refractive index of a sphere
radius = 100*1e-9

# environment
nenv = 1.00

# wavelenth
# lbd0 = 618*1e-9  # ED
# lbd0 = 831*1e-9  # MD
lbd0 = 918*1e-9  # 1nd Kerker
# lbd0 = 756*1e-9  # 2nd Kerker
lbdp = lbd0.*1e9

# angular range
theta = range(0,2π,length=73) # per 5 degrees with length 37
phi = 0.0
# phi = pi/2

# other settings
nmax = -1  # -1: namx large enough (determined from x)

Isff = zeros(length(theta))
for i = 1:length(theta)
    Isff[i],_,_ = jlmie_Isff(nmat,radius,lbd0,nenv,theta[i],phi,nmax)
end
Isff = Isff/maximum(Isff)

# show results
plt = plot(theta, Isff,
           proj=:polar,
           legend=false
)
display(plt)

# save
if true
    runningfilename = splitext(splitpath(@__FILE__)[end])[1]
    outfilename = ".\\result\\" * runningfilename
    savefig(plt,outfilename)
    println("Graph has been saved in result directory")
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
using PyCall, ImageMagick, ImageTransformations, FileIO
imgfile = "ILSVRC2012_val_00000001.JPEG"

tf = pyimport("tensorflow")
torch = pyimport("torch")
Image = pyimport("PIL.Image")
transforms = pyimport("torchvision.transforms")
rdims(a) = permutedims(a, ndims(a):-1:1)
interpolation = Image.BILINEAR
t1 = transforms.Compose([
    transforms.Resize(224; interpolation),     # torch resize image
    transforms.ToTensor(),
])
t2 = transforms.Compose([
    transforms.ToTensor(),
    transforms.Resize(224; interpolation),     # torch resize tensor
])

# https://github.com/tensorflow/models/blob/master/official/vision/image_classification/resnet/imagenet_preprocessing.py
t3(image, height_width) = tf.compat.v1.image.resize(
    tf.keras.preprocessing.image.img_to_array(image), # 3/4-D Tensor of shape [batch?, height, width, channels] if default=channels_last
    height_width,
    method=tf.image.ResizeMethod.BILINEAR,
    align_corners=false
)

img = Image.open(imgfile)
img2 = FileIO.load(imgfile)
pilimage = t1(img).numpy()          # (3, 224, 298) = (channels, height, width)
pttensor = t2(img).numpy()
tfresize = t3(img,[224,298]).numpy() ./ 255 |> (x->permutedims(x,(3,1,2)))
jlresize = imresize(img2, (224,298)) |> channelview |> Array{Float32}
@show display(pilimage[1,101:105,101:105])
@show display(pttensor[1,101:105,101:105])
@show display(tfresize[1,101:105,101:105])
@show display(jlresize[1,101:105,101:105])

# resize variants using ImageMagick don't seem to work as well :(
# mobilenet_v2_100_224_pt: 0.7070 (using torchvision.transforms.Resize)
# mobilenet_v2_100_224_pt: 0.6990 (using torchvision.transforms.Resize after ToTensor)
# mobilenet_v2_100_224_pt: 0.7010 (using imresize)
# for i in */*; do convert -resize 256x256^ $i foo.jpg; mv foo.jpg $i; done
# mobilenet_v2_100_224_pt: 0.6930
# for i in */*; do convert -adaptive-resize 256x256^ $i foo.jpg; mv foo.jpg $i; done
# mobilenet_v2_100_224_pt: 0.6920
# for i in */*; do convert -interpolative-resize 256x256^ $i foo.jpg; mv foo.jpg $i; done
# mobilenet_v2_100_224_pt: 0.69
# for i in */*; do convert -interpolate Bilinear -interpolative-resize 256x256^ $i foo.jpg; mv foo.jpg $i; done
# mobilenet_v2_100_224_pt: 0.6940

#= Chasing down source code for transforms.Resize:
/opt/anaconda3/lib/python3.8/site-packages/torchvision/transforms/transforms.py
Resize: interpolation=InterpolationMode.BILINEAR
            img (PIL Image or Tensor): Image to be scaled.
        return F.resize(img, self.size, self.interpolation)
from . import functional as F

/opt/anaconda3/lib/python3.8/site-packages/torchvision/transforms/functional.py
resize:
     if not isinstance(img, torch.Tensor):
        pil_interpolation = pil_modes_mapping[interpolation]
        return F_pil.resize(img, size=size, interpolation=pil_interpolation)
     return F_t.resize(img, size=size, interpolation=interpolation.value)
pil_modes_mapping = {
    InterpolationMode.NEAREST: 0,
    InterpolationMode.BILINEAR: 2,
    InterpolationMode.BICUBIC: 3,
    InterpolationMode.BOX: 4,
    InterpolationMode.HAMMING: 5,
    InterpolationMode.LANCZOS: 1,
}

/opt/anaconda3/lib/python3.8/site-packages/torchvision/transforms/functional_tensor.py
    from torch.nn.functional import grid_sample, conv2d, interpolate, pad as torch_pad
    img, need_cast, need_squeeze, out_dtype = _cast_squeeze_in(img, [torch.float32, torch.float64])
    align_corners = False if interpolation in ["bilinear", "bicubic"] else None
    img = interpolate(img, size=[size_h, size_w], mode=interpolation, align_corners=align_corners)
    if interpolation == "bicubic" and out_dtype == torch.uint8:
        img = img.clamp(min=0, max=255)
    img = _cast_squeeze_out(img, need_cast=need_cast, need_squeeze=need_squeeze, out_dtype=out_dtype)

/opt/anaconda3/lib/python3.8/site-packages/torchvision/transforms/functional_pil.py
def resize(img, size, interpolation=Image.BILINEAR):
return img.resize((ow, oh), interpolation)
img has a resize operation. _is_pil_image(img):

/opt/anaconda3/lib/python3.8/site-packages/PIL/Image.py
resize(self, size, resample=BICUBIC, box=None, reducing_gap=None):
presumably we call this with PIL.Image.BILINEAR
        return self._new(self.im.resize(size, resample, box))
self.im = core.new(mode, size)
    from . import _imaging as core
    # If the _imaging C module is not present, Pillow will not load.
    # Note that other modules should not refer to _imaging directly;
    # import Image and use the Image.core variable instead.
    # Also note that Image.core is not a publicly documented interface,
    # and should be considered private and subject to change.
https://github.com/python-pillow/Pillow/blob/master/src/_imaging.c
        imOut = ImagingResample(imIn, xsize, ysize, filter, box);
filter=bilinear, box=img-size
https://github.com/python-pillow/Pillow/blob/d374015504df2c7f0316154deade1bc1fa5075df/src/libImaging/Resample.c
return ImagingResampleInner(
        imIn, xsize, ysize, filterp, box, ResampleHorizontal, ResampleVertical);
# filter=&BILINEAR, box=img-size
           case IMAGING_TYPE_UINT8:
                ResampleHorizontal = ImagingResampleHorizontal_8bpc;
                ResampleVertical = ImagingResampleVertical_8bpc;
https://github.com/python-pillow/Pillow/blob/d374015504df2c7f0316154deade1bc1fa5075df/src/libImaging/Resample.c#L616
calls ResampleHorizontal, ResampleVertical

=#


text:
function commandOpen(ed::Editor, args::String)
    if args == ""
        editorOpen(ed)
    else
        editorOpen(ed, args)
    end
end

addCommand(:open, commandOpen,
           help="open: open a file")


text:
recipe = Recipe{Assembler}("Production 1", 15)
@RIN "Red Circuit" 5
@RIN "Green Circuit" 5
@ROUT "Production 1" 1

recipe = Recipe{Assembler}("Production 2", 30)
@RIN "Red Circuit" 5
@RIN "Blue Circuit" 5
@RIN "Production 1" 4
@ROUT "Production 2" 1

recipe = Recipe{Assembler}("Production 3", 60)
@RIN "Red Circuit" 5
@RIN "Blue Circuit" 5
@RIN "Production 2" 4
@RIN "Alien Artifact" 1
@ROUT "Production 3" 1

recipe = Recipe{Assembler}("Blue Circuit", 15)
@RIN "Green Circuit" 20
@RIN "Red Circuit" 2
@RIN "Sulphuric Acid" 0.5
@ROUT "Blue Circuit" 1

recipe = Recipe{Assembler}("Red Circuit", 8)
@RIN "Green Circuit" 2
@RIN "Plastic" 2
@RIN "Copper Cable" 4
@ROUT "Red Circuit" 1

recipe = Recipe{Assembler}("Green Circuit", 0.5)
@RIN "Iron Plate" 1
@RIN "Copper Cable" 3
@ROUT "Green Circuit" 1

recipe = Recipe{Assembler}("Copper Cable", 0.5)
@RIN "Copper Plate" 1
@ROUT "Copper Cable" 2

recipe = Recipe{ChemPlant}("Plastic", 1)
@RIN "Coal" 1
@RIN "Petroleum Gas" 3
@ROUT "Plastic" 2

recipe = Recipe{Smelter}("Iron Plate", 3.5)
@RIN "Iron Ore" 1
@ROUT "Iron Plate" 1

recipe = Recipe{Smelter}("Copper Plate", 3.5)
@RIN "Copper Ore" 1
@ROUT "Copper Plate" 1

recipe = Recipe{Drill}("Iron Ore", 1)
@RIN "Rock" 1
@ROUT "Iron Ore" 1

recipe = Recipe{Drill}("Copper Ore", 1)
@RIN "Rock" 1
@ROUT "Copper Ore" 1

recipe = Recipe{Assembler}("Production Science", 14) # 354 pumpjacks per minute
@RIN "Pumpjack" 1
@RIN "Electric Engine" 1
@RIN "Electic Furnace" 1
@ROUT "Production Science" 2

recipe = Recipe{Assembler}("Pumpjack", 20)
@RIN "Steel" 15
@RIN "Iron Gears" 10
@RIN "Green Circuits" 10
@RIN "Pipe" 10
@ROUT "Pumpjack" 1

recipe = Recipe{Assembler}("Pipe", 0.5)
@RIN "Iron Plate" 1
@ROUT "Pipe" 1

recipe = Recipe{Assembler}("Iron Gears", 0.5)
@RIN "Iron Plate" 2
@ROUT "Iron Gears" 1

recipe = Recipe{Assembler}("Electric Furnace", 5)
@RIN "Steel" 15
@RIN "Stone Block" 10
@RIN "Red Circuit" 5
@ROUT "Electric Furnace" 1

recipe = Recipe{Assembler}("Electric Engine", 20)
@RIN "Engine Unit" 1
@RIN "Green Circuit" 2
@RIN "Lubricant" 2
@ROUT "Electric Engine" 1

recipe = Recipe{Assembler}("Electric Unit", 20)
@RIN "Steel" 1
@RIN "Iron Gears" 1
@RIN "Pipe" 1
@ROUT "Electric Unit" 1

recipe = Recipe{Smelter}("Stone Block", 3.5)
@RIN "Stone" 2
@ROUT "Stone Block" 1

recipe = Recipe{ChemPlant}("Lubricant", 1)
@RIN "Heavy Oil" 1
@ROUT "Lubricant" 1

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
"""
Type representing any object drawable on image
"""
abstract type Drawable end

"""
Root type for the polygon filling algorithm in `ImageDraw` packages

Any concrete polygon filling algorithm shall subtype it to support
`draw` and `draw!` APIs.

"""
abstract type AbstractPolyFillAlgorithm end

"""
    p = Point(x,y)
    p = Point(c)

A `Drawable` point on the image
"""
struct Point <: Drawable
    x::Int
    y::Int
end

abstract type Line <: Drawable end
abstract type Circle <: Drawable end
abstract type Rectangle <: Drawable end

"""
    line = LineTwoPoints(p1, p2)

A `Drawable` infinite length line passing through the two points
`p1` and `p2`.
"""
struct LineTwoPoints <: Line
    p1::Point
    p2::Point
end

"""
    line = LineNormal(ρ, θ)

A `Drawable` infinte length line having perpendicular length `ρ` from
origin and angle `θ` between the perpendicular and x-axis

"""
struct LineNormal{T<:Real, U<:Real} <: Line
    ρ::T
    θ::U
end

"""
    circle = CircleThreePoints(p1, p2, p3)

A `Drawable` circle passing through points `p1`, `p2` and `p3`
"""
struct CircleThreePoints <: Circle
    p1::Point
    p2::Point
    p3::Point
end

"""
    circle = CirclePointRadius(center, ρ; thickness = 0, fill = false)

A `Drawable` circle having center `center` and radius `ρ`
with keyword arguments `thickness` and `fill` for drawing hollow circle

# Arguments
- `thickness::Int`: thickness of the circle
- `fill::Bool`: argument to determine whether to fill the circle or not

# Examples
```
using ImageDraw, TestImages, ImageView, ColorVectorSpace
img = testimage("lighthouse")

draw!(img, Ellipse(CirclePointRadius(350,200,100)))
draw!(img, Ellipse(CirclePointRadius(150,200,100; thickness = 70 , fill = false)))
```

"""
struct CirclePointRadius{T<:Real} <: Circle
    center::Point
    ρ::T
    thickness::Int
    fill::Bool
    function CirclePointRadius{T}(center::Point, ρ::T; thickness::Int = 0, fill::Bool = true) where {T<:Real}
        thickness >= ρ && throw(ArgumentError("Thickness $thickness should be smaller than $(ρ)."))
        new{T}(center, ρ, thickness, fill)
    end
end

CirclePointRadius(center::Point, ρ::T; kwargs...) where {T<:Real} = CirclePointRadius{T}(center, ρ; kwargs...)

"""
    ls = LineSegment(p1, p2)

A `Drawable` finite length line between `p1` and `p2`
"""
struct LineSegment <: Drawable
    p1::Point
    p2::Point
end

"""
    path = Path([point])

A `Drawable` sequence of line segments connecting consecutive pairs
of points in `[point]`.
!!! note
    This will create a non-closed path. For a closed path, see `Polygon`
"""
struct Path <: Drawable
    vertices::Vector{Point}
end

"""
    ellipse = Ellipse(center, ρx, ρy; thickness=0, fill=false)

A `Drawable` ellipse with center `center` and parameters `ρx` and `ρy`
with keyword arguments `thickness` and `fill` for drawing hollow ellipse

# Arguments
- `thickness::Int`: thickness of the circle
- `fill::Bool`: argument to determine whether to fill the ellipse or not

# Examples
```
using ImageDraw,TestImages,ImageView, ColorVectorSpace
img = testimage("lighthouse")

draw(img, Ellipse(5, 5, 5, 5), Gray{N0f8}(0.5))
draw(img, Ellipse(CartesianIndex(5,5), 5, 3))
draw(img, Ellipse(CirclePointRadius(Point(6, 6), 5; thickness=4, fill=false)))
```
"""

struct Ellipse{T<:Real, U<:Real} <: Drawable
    center::Point
    ρx::T
    ρy::U
    thickness::Int
    fill::Bool
    function Ellipse{T, U}(center, ρx::T, ρy::U; thickness::Int = 0, fill::Bool = true) where {T<:Real, U<:Real}
        thickness >= min(ρx, ρy) && throw(ArgumentError("Thickness $thickness should be smaller than $(min(ρx, ρy))."))
        new{T, U}(center, ρx, ρy, thickness, fill)
    end
end

Ellipse(center::Point, ρx::T, ρy::U; kwargs...) where {T<:Real, U<:Real} = Ellipse{T,U}(center, ρx, ρy; kwargs...)

"""
    polygon = Polygon([vertex])

A `Drawable` polygon i.e. a closed path created by joining the
consecutive points in `[vertex]` along with the first and last point.
!!! note
    This will create a closed path. For a non-closed path, see `Path`
"""
struct Polygon <: Drawable
    vertices::Vector{Point}
end

"""
    BoundaryFill <: AbstractPolyFillAlgorithm
    BoundaryFill(x::Int, y::Int, fill_value::T, boundary_value::T)

    draw(img, verts, alg::BoundaryFill; closed)
    draw!(img, verts, alg::BoundaryFill; closed)

Applies boundary filling algortithm on image inside the vertices provided in order

# Output

Return the boundary filled image inside the vertices(provided in order).

# Arguments

## `x` && `y`

Boundary fill is a seeded polygon filling algorithm.
So we need to provide the seed point (x,y) inside the image from where the algorithm can start its function.

## `fill_value`

The pixels inside the boundary of vertices of polygon are filled with this colorant parameter.

## `boundary_value`

The pixels to define the boundary using verts is filled with this colorant parameter.
Without this parameter, the whole picture will be filled with the `fill_value`

# `closed`

`closed` keyword specifies whether to connect the polygon edges using the verts provided.
Edges will be filled and connected using `boundary_value`.

```julia
using ImageDraw

img = zeros(RGB, 7, 7)
expected = copy(img)
expected[2:6, 2:6] .= RGB{N0f8}(1)

verts = [CartesianIndex(2, 2), CartesianIndex(2, 6), CartesianIndex(6, 6), CartesianIndex(6, 2), CartesianIndex(2,2)]

fill_method = draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
```
"""
struct BoundaryFill{T<:Colorant} <: AbstractPolyFillAlgorithm
    x::Int
    y::Int
    fill_value::T
    boundary_value::T
    function BoundaryFill(x::Int, y::Int, fill_value::T, boundary_value::T) where {T <: Colorant}
        new{T}(x, y, fill_value, boundary_value)
    end
end

BoundaryFill(x::Int = 1, y::Int = 1; fill_value::Colorant = RGB(1), boundary_value::Colorant = fill_value) = BoundaryFill(x, y, fill_value, boundary_value)
BoundaryFill(p::CartesianIndex{2}; fill_value::Colorant = RGB(1), boundary_value::Colorant = fill_value) = BoundaryFill(p[1], p[2], fill_value, boundary_value)
BoundaryFill(p::Point; fill_value::Colorant = RGB(1), boundary_value::Colorant = fill_value) = BoundaryFill(p.y, p.x, fill_value, boundary_value)

"""
    rectangle = RectanglePoints(p1, p2)
    rectangle = RectanglePoints(x1, y1, x2, y2)

A `Drawable` rectangle i.e. a closed path where parameters `p1` and `p2`
are diagonally opposite vertices of the rectangle.

Parameters `p1` and `p2` can be passed in either `Point` struct format or `CartesianIndex` format as shown in example below.
`(x1, y1, x2, y2)` is interpreted as `(CartesianIndex(x1, y1), CartesianIndex(x2, y2))`.

# Examples

```julia
using TestImages, ImageDraw, ImageCore
img = testimage("lighthouse")

# rectangles drawn on img using different methods of passing parameters
draw!(img, Polygon(RectanglePoints(Point(10, 10), Point(100, 100))), RGB{N0f8}(1))
draw!(img, Polygon(RectanglePoints(CartesianIndex(110, 10), CartesianIndex(200, 200))), RGB{N0f8}(1))
draw!(img, Polygon(RectanglePoints(220, 10, 300, 300)), RGB{N0f8}(1))
```

"""
struct RectanglePoints <: Rectangle
    p1::Point
    p2::Point
end

"""
    rp = RegularPolygon(center, side_count, side_length, θ)

A `Drawable` regular polygon.

#Arguments
* `center::Point` : the center of the polygon
* `side_count::Int` : number of sides of the polygon
* `side_length::Real` : length of each side
* `θ::Real` : orientation of the polygon w.r.t x-axis (in radians)

"""
struct RegularPolygon{T<:Real, U<:Real} <: Drawable
    center::Point
    side_count::Int
    side_length::T
    θ::U
end

"""
    cross = Cross(c, range::UnitRange{Int})
A `Drawable` cross passing through the point `c` with arms ranging across `range`.
"""
struct Cross <: Drawable
    c::Point
    range::UnitRange{Int}
end

"""
    img = draw!(img, drawable, color)
    img = draw!(img, drawable)

Draws `drawable` on `img` using color `color` which
defaults to `oneunit(eltype(img))`
"""
draw!(img::AbstractArray{T,2}, object::Drawable) where {T<:Colorant} = draw!(img, object, oneunit(T))

"""
    img = draw!(img, [drawable], [color])
    img = draw!(img, [drawable] ,color)
    img = draw!(img, [drawable])

Draws all objects in `[drawable]` in the given order on `img` using
corresponding colors from `[color]` which defaults to `oneunit(eltype(img))`
If only a single color `color` is specified then all objects will be
colored with that color.
"""
function draw!(img::AbstractArray{T,2}, objects::AbstractVector{U}, colors::AbstractVector{V}) where {T<:Colorant, U<:Drawable, V<:Colorant}
    colors = copy(colors)
    while length(colors) < length(objects)
        push!(colors, oneunit(T))
    end
    foreach((object, color) -> draw!(img, object, color), objects, colors)
    img
end

draw!(img::AbstractArray{T,2}, objects::AbstractVector{U}, color::T = oneunit(T)) where {T<:Colorant, U<:Drawable} =
    draw!(img, objects, [color for i in 1:length(objects)])

"""
    img_new = draw(img, drawable, color)
    img_new = draw(img, [drawable], [color])

Draws the `drawable` object on a copy of image `img` using color
`color`. Can also draw multiple `Drawable` objects when passed
as a `AbstractVector{Drawable}` with corresponding colors in `[color]`
"""
draw(img::AbstractArray{T,2}, args...) where {T<:Colorant} = draw!(copy(img), args...)

Point(τ::Tuple{Int, Int}) = Point(τ...)
Point(p::CartesianIndex) = Point(p[2], p[1])

function draw!(img::AbstractArray{T,2}, point::Point, color::T) where T<:Colorant
    drawifinbounds!(img, point, color)
end

"""

    img_new = drawifinbounds!(img, y, x, color)
    img_new = drawifinbounds!(img, Point, color)
    img_new = drawifinbounds!(img, CartesianIndex, color)

Draws a single point after checkbounds() for coordinate in the image.
Color Defaults to oneunit(T)

"""
drawifinbounds!

drawifinbounds!(img::AbstractArray{T,2}, p::Point, color::T = oneunit(T)) where {T<:Colorant} = drawifinbounds!(img, p.y, p.x, color)
drawifinbounds!(img::AbstractArray{T,2}, p::CartesianIndex{2}, color::T = oneunit(T)) where {T<:Colorant} = drawifinbounds!(img, Point(p), color)

function drawifinbounds!(img::AbstractArray{T,2}, y::Int, x::Int, color::T) where {T<:Colorant}
    if checkbounds(Bool, img, y, x) img[y, x] = color end
    img
end