text:
"""
Custom JSON serializer for log entries.
Handles Module types for now, more can be added later.
"""
struct LogEntrySerialization <: CommonSerialization end

show_json(io::StructuralContext, ::LogEntrySerialization, m::Module) = show_json(io, LogEntrySerialization(), string(m))
show_json(io::StructuralContext, ::LogEntrySerialization, ptr::Ptr) = show_json(io, LogEntrySerialization(), string(ptr))
show_json(io::StructuralContext, ::LogEntrySerialization, sv::Core.SimpleVector) = show_json(io, LogEntrySerialization(), [sv...])
show_json(io::StructuralContext, ::LogEntrySerialization, typ::DataType) = show_json(io, LogEntrySerialization(), string(typ))

function show_json(io::StructuralContext, ::LogEntrySerialization, level::Logging.LogLevel)
    levelstr = (level == Logging.Debug) ? "Debug" :
               (level == Logging.Info)  ? "Info" :
               (level == Logging.Warn)  ? "Warn" :
               (level == Logging.Error) ? "Error" :
               "LogLevel($(level.level))"
    show_json(io, LogEntrySerialization(), levelstr)
end

function show_json(io::StructuralContext, ::LogEntrySerialization, exception::Tuple{Exception,Any})
    iob = IOBuffer()
    Base.show_exception_stack(iob, [exception])
    show_json(io, LogEntrySerialization(), String(take!(iob)))
end

as_text(str::String) = str
function as_text(obj)
    iob = IOBuffer()
    lim = LimitIO(iob, 4*1024)  # fixing it as of now to a large enough size for most use cases
    try
        show(lim, "text/plain", obj)
    catch ex
        if isa(ex, LimitIOException)
            # ignore and take what was printed
            print(iob, "...")
        else
            rethrow()
        end
    end

    String(take!(iob))
end

"""
IndexedLogEntry represents a log entry as a dictionary and its
indexable attributes in a form that is useful to a logging sink.

The index part contains metadata that are to be indexed. Event metadata
consists of attributes like level, module, filepath, line, job id,
process id, user id, etc. It also includes application specific
keywords that the originating piece of code wishes to index.

Keywords that should be considered as metadata are indicated via the
`indexable` constructor parameter.

What metadata can be indexed depends on the type of sink and whether
it has support to index certain types of attributes. Attributes that
the sink can not index are made part of the message itself for storage.

The message part can contain the following keys unless they are empty:
- `metadata`: event metadata that could not be indexed
- `message`: the log message string
- `keywords`: any keywords provided

Constructor parameters:
- `log`: Named tuple containing args to the handle_message method, e.g.: (level, message, _module, group, id, file, line, kwargs)
- `indexable`: list of names from `log` and `log.kwargs` that should be included in the index
"""
struct IndexedLogEntry
    index::Dict{Symbol,Any}
    message::Dict{Symbol,Any}
end

function IndexedLogEntry(log, indexable::Vector{Symbol}=[:level, :module, :filepath, :line])
    index = Dict{Symbol, Any}()
    metadata = Dict{Symbol,Any}()
    keywords = Dict(log.kwargs)

    log_prop_names = propertynames(log)
    for name in log_prop_names
        (name === :kwargs) && continue # skip the kwargs, deal with that separately
        (name === :message) && continue # skip message, we are dealing with that separately
        ((name in indexable) ? index : metadata)[name] = getproperty(log, name)
    end
    for name in keys(keywords)
        (name in log_prop_names) && continue # avoid clobbering reserved names
        if name in indexable
            index[name] = keywords[name]
            delete!(keywords, name)
        end
    end

    message = Dict{Symbol,Any}()
    messagestr = as_text(log.message)
    isempty(messagestr)     || (message[:message]  = messagestr)
    isempty(metadata)       || (message[:metadata] = metadata)
    isempty(keywords)       || (message[:keywords] = keywords)

    IndexedLogEntry(index, message)
end

message_string(entry::IndexedLogEntry, size_limit::Int, newline::Bool=false) = message_string(entry.message, size_limit, newline)
function message_string(message::Dict{Symbol,Any}, size_limit::Int, newline::Bool=false)
    iob = IOBuffer()
    lim = LimitIO(iob, size_limit)
    try
        JSON.show_json(lim, LogEntrySerialization(), message)
        newline && write(lim, '\n')
    catch ex
        if isa(ex, LimitIOException)
            if haskey(message, :keywords)
                # strip off keywords (if any) and retry
                delete!(message, :keywords)
                return message_string(message, size_limit, newline)
            elseif haskey(message, :message)
                # strip off the message and retry (we retain only the metadata in the worst case)
                delete!(message, :message)
                return message_string(message, size_limit, newline)
            end
        end
        rethrow(ex)
    end
    String(take!(iob))
end


text:
export use_pyplot

global plot_system = :none

################################################################################
#                         Selection of plotting system
################################################################################

"""
### function use_pyplot(pyplot_gui)

Select the PyPlot as the plotting system.

##### Args

* pyplot_gui: PyPlot GUI (see PyPlot).

"""

function use_pyplot(pyplot_gui)
    # List obtained from PyPlot.jl package.
    gui2matplotlib = Dict(:wx=>"WXAgg",:gtk=>"GTKAgg",:gtk3=>"GTK3Agg",
                          :qt_pyqt4=>"Qt4Agg",
                          :qt_pyqt5=>"Qt5Agg",
                          :qt_pyside=>"Qt4Agg",
                          :qt4=>"Qt4Agg",
                          :qt5=>"Qt5Agg",
                          :qt=>"Qt4Agg",:tk=>"TkAgg")

    # Import matplotlib.
    const matplotlib = pyimport("matplotlib")

    # Select the GUI and set to interactive mode.
    pygui_start(pyplot_gui)
    matplotlib[:use](gui2matplotlib[pyplot_gui])
    matplotlib[:interactive](true)

    # Import pyplot.
    pltm = pyimport("matplotlib.pyplot") # raw Python module
    global plt = pywrap(pltm)

    # Set the plotting system as PyPlot.
    global plot_system = :pyplot

    nothing
end

################################################################################
#                                    Plots
################################################################################

# Definitions

const blw   = 2.0   # Border line width.
const plw   = 2.5   # Plot line width.
const glw   = 1.3   # Grid line width.
const sglw  = 1.0   # Subgrid line width.
const tifsz = 14    # Title font size.
const mksz  = 80    # Marker size.
const arrsz = 20    # Arrow height.
const alw   = 1.5   # Axes line width.

################################################################################
#                                   Includes
################################################################################

include("plot_bode.jl")
include("plot_lsim.jl")
include("plot_margin.jl")
include("plot_nyquist.jl")
include("plot_rlocus.jl")

ext: jl


text:
function getfxpwidths(x::Real)
    xabs = abs(x)
    xint = floor(xabs)
    xfrac = xabs - xint
    intwidth = iszero(xint) ? 0 : Int(ceil(log2(xint)))
    fracwidth = iszero(xfrac) ? 0 : Int(ceil(-log2(xfrac)))

    return (integral = intwidth, fractional = fracwidth)
end
function getfxpwidths(x::AbstractArray{<:Real})
    widths = getfxpwidths.(x)

    return (integral = maximum(first.(widths)), fractional = maximum(last.(widths)))
end
function getfxpwidths(x::NTuple{N, <:Real}) where N
    widths = map(getfxpwidths, x)

    return (integral = maximum(first.(widths)), fractional = maximum(last.(widths)))
end

function getfixedpoint(x::Real, width)
    xabs = abs(x)
    intrep = Int(floor(xabs * 2^width.fractional))
    bitwidth = width.integral + width.fractional + 1
    binstr = "$(x < 0 ? "-" : "")$bitwidth'b$(string(intrep; pad = bitwidth, base = 2))"

    return binstr
end
getfixedpoint(x::AbstractArray{<:Real}, width) = getfixedpoint.(x, Ref(width))
getfixedpoint(x::NTuple{<:Any, <:Real}, width) = map(x) do xi
    getfixedpoint(xi, width)
end

function constantreplacement!(m::Module)
    maxintwidth = m.bitwidth.integral
    maxfracwidth = m.bitwidth.fractional

    # find require widths
    for v in vertices(m.dfg)
        inputs = getinputs(m.dfg, v)
        for input in inputs
            if (isconstant(input) || isparameter(input)) &&
               (jltypeof(input) <: Union{Real, AbstractArray{<:Real}, NTuple{<:Any, <:Real}})
                width = getfxpwidths(value(input))
                maxintwidth = max(maxintwidth, width.integral)
                maxfracwidth = max(maxfracwidth, width.fractional)
            end
        end
    end

    # replace all constants
    m.bitwidth = (integral = maxintwidth, fractional = maxfracwidth)
    # constreplacements = []
    for v in vertices(m.dfg)
        inputs = getinputs(m.dfg, v)
        for (i, input) in enumerate(inputs)
            if jltypeof(input) <: Union{Real, AbstractArray{<:Real}, NTuple{<:Any, <:Real}}
                if isconstant(input)
                    conststr = getfixedpoint(value(input), m.bitwidth)
                    inputs[i] = setwidth(setname(input, conststr), sum(m.bitwidth) + 1)
                elseif isparameter(input)
                    conststr = getfixedpoint(value(input), m.bitwidth)
                    inputs[i] = setwidth(input, sum(m.bitwidth) + 1)
                    m.parameters[name(input)] = conststr
                end
            end
        end
        set_prop!(m.dfg, v, :inputs, inputs)
    end

    return m
end


text:
module FASTQDemultiplexer

using BioSequences.FASTQ
using YAML

include("utils.jl")
include("protocol.jl")
include("barcodes.jl")
include("record.jl")
include("input.jl")
include("output/output.jl")
include("demultiplexer.jl")

end # module

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text: @nanosoldier `runbenchmarks("array" && !("simd" || "linalg"), vs = ":release-0.4")`

max_stars_repo_name: Seo-Yeon-jl/Magnum

max_issues_repo_name: Seo-Yeon-jl/Magnum

max_forks_repo_name: Seo-Yeon-jl/Magnum


text:
using Distributions, ProgressMeter

"""
    bet(dealer, self_strat, opp_strat, self_prob, opp_prob, self_return,
        opp_return, pot, big_blind)

Processes the bets for the players, given thresholds at which each player will
not put more money into the game, and the calculated odds of each player
winning the showdown. Also returns whether the whole get_return function should
return

"""
function bet(dealer::String, self_strat::Float64, opp_strat::Float64,
  self_prob::Float64, opp_prob::Float64, self_return::Float64,
  opp_return::Float64, pot::Float64, big_blind::Float64)
    to_return = false
    if dealer == "opp"
        if opp_prob > opp_strat
            if self_prob < self_strat
                opp_return += pot
                to_return = true
            else
                opp_return -= 2 * big_blind
                self_return -= 2 * big_blind
                pot += 4 * big_blind
            end
        else
            if self_prob > self_strat
                self_return += pot
                to_return = true
            end
        end
    else
        if self_prob > self_strat
            if opp_prob < opp_strat
                self_return += pot
                to_return = true
            else
                opp_return -= 2 * big_blind
                self_return -= 2 * big_blind
                pot += 4 * big_blind
            end
        else
            if opp_prob > opp_strat
                opp_return += pot
                to_return = true
            end
        end
    end
    return to_return, self_return, opp_return, pot
end

"""
    get_return(self_strat, opp_strat[, preflop_odds_table])

Returns the self return after 1 simulated game, given the singular thresholds
at which either player will not put more money into the game. Optional
optimization is available by specifying a pre-calculated preflop odds table

"""
function get_return(self_strat::Float64, opp_strat::Float64;
  preflop_odds_table::Array{Float64, 2} = zeros(52, 52))
    using_preflop_odds = preflop_odds_table == zeros(52, 52)
    big_blind, small_blind = 1.0, 0.5
    self_return, opp_return, pot = 0.0, 0.0, 0.0

    deck = fresh_deck
    table, self_hand, opp_hand = Hand([]), Hand([]), Hand([])
    dealer = rand(["self", "opp"])
    if dealer == "opp" # Opponent deals
        self_return -= big_blind
        opp_return -= small_blind
        self_hand, deck = deal(deck, 2)
        opp_hand, deck = deal(deck, 2)
    else # Self deals
        self_return -= small_blind
        opp_return -= big_blind
        opp_hand, deck = deal(deck, 2)
        self_hand, deck = deal(deck, 2)
    end
    pot += big_blind + small_blind
    if using_preflop_odds
        self_prob, opp_prob = pre_flop_odds_table[
            findfirst(isequal(self_hand[1]), fresh_deck_cards),
            findfirst(isequal(self_hand[2]), fresh_deck_cards)],
            pre_flop_odds_table[
            findfirst(isequal(opp_hand[1]), fresh_deck_cards),
            findfirst(isequal(opp_hand[2]), fresh_deck_cards)]
    else
        self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
            calculate_odds(opp_hand, table, N = 100)
    end

    # First round of betting
    to_return, self_return, opp_return, pot = bet(dealer, self_strat, opp_strat,
        self_prob, opp_prob, self_return, opp_return, pot, big_blind)
    to_return && return self_return

    # Flop
    dealt, deck = deal(deck, 3)
    table += dealt
    self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
        calculate_odds(opp_hand, table, N = 100)

    # Second round of betting
    to_return, self_return, opp_return, pot = bet(dealer, self_strat, opp_strat,
        self_prob, opp_prob, self_return, opp_return, pot, big_blind)
    to_return && return self_return

    # Turn
    dealt, deck = deal(deck, 1)
    table += dealt
    self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
        calculate_odds(opp_hand, table, N = 100)

    # Third round of betting
    to_return, self_return, opp_return, pot = bet(dealer, self_strat, opp_strat,
        self_prob, opp_prob, self_return, opp_return, pot, big_blind)
    to_return && return self_return

    # River
    dealt, deck = deal(deck, 1)
    table += dealt
    self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
        calculate_odds(opp_hand, table, N = 100)

    # Fourth and final round of betting
    to_return, self_return, opp_return, pot = bet(dealer, self_strat, opp_strat,
        self_prob, opp_prob, self_return, opp_return, pot, big_blind)
    to_return && return self_return

    # Showdown
    if (self_hand + table) < (opp_hand + table)
        opp_return += pot
    elseif (opp_hand + table) < (self_hand + table)
        self_return += pot
    else
        self_return += pot / 2
        opp_return += pot / 2
    end
    return self_return
end

"""
    get_return(self_strat, opp_strat[, preflop_odds_table])

Returns the self return after 1 simulated game, given the multiple thresholds
(one for each round of betting; ie, 4) at which either player will not put more
money into the game. Optional optimization is available by specifying a
pre-calculated preflop odds table

"""
function get_return(self_strat::Array{Float64}, opp_strat::Array{Float64};
  preflop_odds_table::Array{Float64, 2} = zeros(52, 52))
    using_preflop_odds = preflop_odds_table == zeros(52, 52)
    big_blind, small_blind = 1.0, 0.5
    self_return, opp_return, pot = 0.0, 0.0, 0.0

    deck = fresh_deck
    table, self_hand, opp_hand = Hand([]), Hand([]), Hand([])
    dealer = rand(["self", "opp"])
    if dealer == "opp" # Opponent deals
        self_return -= big_blind
        opp_return -= small_blind
        self_hand, deck = deal(deck, 2)
        opp_hand, deck = deal(deck, 2)
    else # Self deals
        self_return -= small_blind
        opp_return -= big_blind
        opp_hand, deck = deal(deck, 2)
        self_hand, deck = deal(deck, 2)
    end
    pot += big_blind + small_blind
    if using_preflop_odds
        self_prob, opp_prob = pre_flop_odds_table[
            findfirst(isequal(self_hand[1]), fresh_deck_cards),
            findfirst(isequal(self_hand[2]), fresh_deck_cards)],
            pre_flop_odds_table[
            findfirst(isequal(opp_hand[1]), fresh_deck_cards),
            findfirst(isequal(opp_hand[2]), fresh_deck_cards)]
    else
        self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
            calculate_odds(opp_hand, table, N = 100)
    end

    # First round of betting
    @inbounds to_return, self_return, opp_return, pot = bet(dealer,
        self_strat[1], opp_strat[1], self_prob, opp_prob, self_return,
        opp_return, pot, big_blind)
    to_return && return self_return

    # Flop
    dealt, deck = deal(deck, 3)
    table += dealt
    self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
        calculate_odds(opp_hand, table, N = 100)

    # Second round of betting
    @inbounds to_return, self_return, opp_return, pot = bet(dealer,
        self_strat[2], opp_strat[2], self_prob, opp_prob, self_return,
        opp_return, pot, big_blind)
    to_return && return self_return

    # Turn
    dealt, deck = deal(deck, 1)
    table += dealt
    self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
        calculate_odds(opp_hand, table, N = 100)

    # Third round of betting
    @inbounds to_return, self_return, opp_return, pot = bet(dealer,
        self_strat[3], opp_strat[3], self_prob, opp_prob, self_return,
        opp_return, pot, big_blind)
    to_return && return self_return

    # River
    dealt, deck = deal(deck, 1)
    table += dealt
    self_prob, opp_prob = calculate_odds(self_hand, table, N = 100),
        calculate_odds(opp_hand, table, N = 100)

    # Fourth and final round of betting
    @inbounds to_return, self_return, opp_return, pot = bet(dealer,
        self_strat[4], opp_strat[4], self_prob, opp_prob, self_return,
        opp_return, pot, big_blind)
    to_return && return self_return

    # Showdown
    if (self_hand + table) < (opp_hand + table)
        opp_return += pot
    elseif (opp_hand + table) < (self_hand + table)
        self_return += pot
    else
        self_return += pot / 2
        opp_return += pot / 2
    end
    return self_return
end

function get_returns(self_strat::Float64, opp_strat::Float64;
  preflop_odds_table::Array{Float64, 2} = zeros(52, 52), N = 5_000)
    return mean(map(x -> get_return(self_strat, opp_strat,
        preflop_odds_table = preflop_odds_table), 1:N))
end

"""
    get_returns(self_strat, opp_strat[, preflop_odds_table, N])

Wrapping of get_return, but the mean of N (default 5,000) times

"""
function get_returns(self_strat::Array{Float64}, opp_strat::Array{Float64};
  preflop_odds_table::Array{Float64, 2} = zeros(52, 52), N = 5_000)
    return mean(map(x -> get_return(self_strat, opp_strat,
        preflop_odds_table = preflop_odds_table), 1:N))
end

"""
    get_returns(self_strat, opp_strat[, preflop_odds_table, N])

Wrapping of get_return, but the mean of N (default 100) times

"""
function get_returns(self_strat::Array{Float64};
  preflop_odds_table::Array{Float64, 2} = zeros(52, 52), N = 100)
    return mean(map(x -> get_return(self_strat,
        [rand(), rand(), rand(), rand()],
        preflop_odds_table = preflop_odds_table), 1:N))
end

"""
    get_returns_table([preflop_odds_table, N])

Calculation of a table of average returns across a gamut of singular threshold
strategies, with the optional optimization of pre-calculating the preflop odds
table. Default N is 5,000 simulations per pair of strategies (total 50,000,000
games)

"""
function get_returns_table(;
  preflop_odds_table::Array{Float64, 2} = zeros(52, 52), N = 5_000)
    returns_table = zeros(100, 100)
    @showprogress for i = 1:100, j = 1:100
        if i != j
            returns_table[i, j] = get_returns(i / 100, j / 100,
                preflop_odds_table = preflop_odds_table, N = N)
        end
    end
    return returns_table
end


text:
include("help.jl")
  function chess!(r::Robot)
    num_hor = moves!(r,West)
    num_vert = moves!(r,Sud)
   #Левый нижний угол
   side = Ost
   if ((num_hor % 2) == 0) #если слево расстояние четна
       if ((num_vert % 2) == 0) # если вертикаль четна
           putmarks2!(r,side)
       else # если вертикаль не четна
           putmarks1!(r,side)
       end
   else
       putmarks1!(r,side)
   end
   #алгоритм - выполнен
   movements!(r,West)
   movements!(r,Sud)
   movements!(r,Ost,num_hor)
   movements!(r,Nord,num_vert)
end
function putmarks1!(r::Robot,side::HorizonSide)
   movements!(r,inverse(side))
   move!(r,side)
   while isborder(r,side) == false
       putmarker!(r)
       for _ in 1:2
           if (isborder(r,side) == false)
               move!(r,side)
           else
               if (isborder(r,Nord) == fasle)
                   move!(r,Nord)
                   side = inverse(side)
                   move!(r,side)
               end
           end
       end
   end
   putmarker!(r)
   if (isborder(r,Nord) == false)
       move!(r,Nord)
       putmarks2!(r,Ost)
   end
end
function putmarks2!(r::Robot,side::HorizonSide)
   movements!(r,inverse(side))
   while isborder(r,side) == false
       putmarker!(r)
       for _ in 1:2
           if (isborder(r, side) == false)
               move!(r,side)
           else
               if (isborder(r,Nord) == false)
                   move!(r,Nord)
                   side = inverse(side)
               end
           end
       end
   end
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
"""
    struct DBSCAN{T<:AbstractFloat, KD<:KDTree}
        df::Matrix{T}
        ϵ::T
        min_pts::Int
        labels::Vector{Int}
        tree::KD
        clusters::Vector{Vector{Int}}
    end

Struct that contains all the relevant information about the model. The `data`, `ϵ`, `min_pts`, `labels`, `KDTree` and `clusters`.
The labels and the clusters are defined during the algorithm inside the routine inside `fit!()` function, which iterate over every observation p in the dataset.

Internal/External Constructors

* Int - DBSCAN(df::Matrix{T}, ϵ::T, min_pts::Int) where {T<:AbstractFloat}
* Ext - DBSCAN(table, ϵ::Real, min_pts::Int)
* Ext - DBSCAN(df::Matrix{T}, ϵ::Real, min_pts::Int) where {T<:AbstractFloat} = DBSCAN(df, T(ϵ), min_pts)
* Ext - DBSCAN(df::Matrix{T}, ϵ::Real, min_pts::Int) where {T} = DBSCAN(Matrix{Float64}(df), ϵ, min_pts)
"""
struct DBSCAN{T<:AbstractFloat, KD<:KDTree}
    df::Matrix{T}
    ϵ::T
    min_pts::Int
    labels::Vector{Int}
    tree::KD
    clusters::Vector{Vector{Int}}

    function DBSCAN(df::Matrix{T}, ϵ::T, min_pts::Int) where {T<:AbstractFloat}
        labels = fill(-1, size(df,1))
        tree = KDTree(df', leafsize=20)
        clusters = Vector{Vector{Int}}()
        KD = typeof(tree)

        new{T, KD}(df, ϵ, min_pts, labels, tree, clusters)
    end
end

# External Constructors
DBSCAN(df::Matrix{T}, ϵ::Real, min_pts::Int) where {T<:AbstractFloat} = DBSCAN(df, T(ϵ), min_pts)
DBSCAN(df::Matrix{T}, ϵ::Real, min_pts::Int) where {T} = DBSCAN(Matrix{Float64}(df), ϵ, min_pts)

function DBSCAN(table, ϵ::Real, min_pts::Int)
    Tables.istable(table) ? (df = Tables.matrix(table)) : throw(ArgumentError("The df argument passed does not implement the Tables.jl interface."))
    return DBSCAN(df, ϵ, min_pts)
end

"""
    dbscan(df, ϵ::Real, min_pts::Int)

Classify data observations in clusters and noises by using a density concept obtained with the parameters input (`ϵ`, `min_pts`).

The number of clusters are obtained during the execution of the model, therefore, initially the user don't know how much clusters it will obtain.
The algorithm use the `KDTree` structure from `NearestNeighbors.jl` to calculate the `RangeQuery` operation more efficiently.

For more detailed explanation of the algorithm, check the [`Algorithm's Overview of DBSCAN`](https://github.com/AugustoCL/ClusterAnalysis.jl/blob/main/algo_overview/dbscan_overview.md)
"""
function dbscan(df, ϵ::Real, min_pts::Int)
    model = DBSCAN(df, ϵ, min_pts)
    fit!(model)
    return model
end

function fit!(model::DBSCAN)
    nlines = size(model.df, 1)                  # Number of observations/points in dataset
    visited = falses(nlines)                    # Array to control if point was visited
    clusters_selected = falses(nlines)          # Array to identify if point was labeled in for-loop
    FIFO = Int[]                                # FIFO data structure
    nb_list = Int[]                             # Array with the neighbors's point

    # for each point p in data
    for p ∈ 1:nlines
        visited[p] && continue                  # Check if point p was visited
        push!(FIFO, p)                          # Add p in FIFO (First In First Out) to execute RangeQuery
        fill!(clusters_selected, false)

        # start RangeQuery all points in FIFO until finished
        # during loop the FIFO got updated with more points
        while !isempty(FIFO)
            newp = popfirst!(FIFO)              # Remove and return the first element in FIFO
            visited[newp] && continue           # Check if point `newp` was visited
            visited[newp] = true                # assign p as visited

            nbs = inrange(model.tree,           # get all the neighbors (index) from point `newp`
                          model.df[newp, :],    # it's the RangeQuery from NearestNeighbors.jl
                          model.ϵ)

            # later, try remove nb_list
            append!(nb_list, nbs)               # save neighbor's indexes in nb_list
            clusters_selected[nb_list] .= true

            if length(nb_list) < model.min_pts  # check if point is core
                empty!(nb_list)
                continue
            end

            for i in nb_list                    # update the FIFO wih nbs
                visited[i] && continue
                push!(FIFO, i)
            end

            empty!(nb_list)                     # clean list to evaluate new point
        end

        # save the indices from each cluster
        ind_cluster = findall(clusters_selected)
        model.min_pts ≤ sum(clusters_selected) && push!(model.clusters, ind_cluster)
    end

    # write labels to struct and disregards duplicate indices
    nc = 1
    for cl in model.clusters
        for i in cl
            model.labels[i] = nc
        end
        nc += 1
    end
end


text:
mutable struct HDF5DiskArray{T, N, C, D, R} <: AbstractDiskArray{T, N}
    ds::HDF5Dataset
    cs::C
    lo::NTuple{N, Int}
    hi::NTuple{N, Int}
    cache::Array{T, D}
end

Base.size(x::HDF5DiskArray{T, N}) where {T, N} = size(x.ds)::NTuple{N, Int}

haschunks(x::HDF5DiskArray{<:Any, <:Any, Nothing}) = Chunked()
haschunks(x::HDF5DiskArray) = Unchunked()

eachchunk(x::HDF5DiskArray{<:Any, <:Any, <:GridChunks}) = x.cs

readblock!(x::HDF5DiskArray, aout, r::OrdinalRange...) = aout .= x.ds[r...]

function readblock!(A::HDF5DiskArray, A_ret, r::AbstractVector...)
    r′ = map(i -> i isa OrdinalRange ? i : minimum(i):maximum(i), r)
    r′′ = map(i -> i isa OrdinalRange ? (:) : i .- i[1] .+ 1, r)
    A_temp = similar(A_ret, length.(r′))
    readblock!(A, A_temp, r′...)
    A_ret .= view(A_temp, r′′...)
    nothing
end

writeblock!(x::HDF5DiskArray, v, r::OrdinalRange...) = x.ds[r...] = v

function writeblock!(A::HDF5DiskArray, A_ret, r::AbstractVector...)
    r′ = map(i -> i isa OrdinalRange ? i : minimum(i):maximum(i), r)
    r′′ = map(i -> i isa OrdinalRange ? (:) : i .- i[1] .+ 1, r)
    A_temp = similar(A_ret, length.(r′))
    A_temp[r′′...] = A_ret
    writeblock!(A, A_temp, r′...)
    nothing
end

const _cache_size = Ref(10 * 1024^2)

function set_cache_size(cache_size)
    _cache_size[] = cache_size
end

get_cache_size() = _cache_size[]

get_cache_size(ds::HDF5Dataset) = _cache_size[] ÷ sizeof(eltype(ds))

function HDF5DiskArray(ds::HDF5Dataset)
    cs = try
        disable_dag()
        GridChunks(ds, get_chunk(ds))
        enable_dag()
    catch
        nothing
    end
    T, N, C = eltype(ds), ndims(ds), typeof(cs)
    strides = cumprod(collect(size(ds)))
    D = findlast(strides .< get_cache_size(ds))
    D = min(something(D, 1) + 1, N)
    if D == 1
        R = min(get_cache_size(ds), size(ds, 1))
    else
        R = min(ceil(Int, get_cache_size(ds) / strides[D - 1]), size(ds, D))
    end
    lo, hi = ntuple(zero, N), ntuple(zero, N)
    cache = zeros(T, size(ds)[1:(D - 1)]..., 0)
    HDF5DiskArray{T, N, C, D, R}(ds, cs, lo, hi, cache)
end

@generated function _getindex(x::HDF5DiskArray{T, N, C, D, R}, r::Integer...) where {T, N, C, D, R}
    colons = fill(:(:), D - 1)
    rl = [:(r[$d]) for d in 1:(D - 1)]
    rr = [:(r[$d]) for d in (D + 1):N]
    cond = :(r[$D] < x.lo[$D] || r[$D] > x.hi[$D])
    for d in (D + 1):N
        cond = :($cond || r[$d] != x.lo[$d])
    end
    ex = quote
        @inbounds if $cond
            x.lo, x.hi = r, min.(r .+ $R, size(x))
            x.cache = x[$(colons...), x.lo[$D]:x.hi[$D], $(rr...)]
        end
        @inbounds v = x.cache[$(rl...), r[$D] - x.lo[$D] + 1]
        return v
    end
    return ex
end

Base.getindex(x::HDF5DiskArray, r::CartesianIndex) = _getindex(x, Tuple(r)...)

Base.getindex(x::HDF5DiskArray, r::Integer...) = _getindex(x, r...)

Base.getindex(x::HDF5DiskArray, i::Integer) =  getindex(x, CartesianIndices(x)[i])

Base.getindex(x::HDF5DiskArray{T, 1}, i::Integer) where T = _getindex(x, i)

Base._reshape(x::HDF5DiskArray, dims::NTuple{N, Int}) where N = Base.__reshape((x, IndexStyle(x)), dims)

Base.Array(x::HDF5DiskArray) = read(x.ds)


text:
using ApproxFun, SpectralMeasures

a = [3/4,-1/4,1/2]; b = [1,3/4]
Q,L = ql(a,b,2.1,0.5);

Q,L = ql(J-5*I)


text:
function get_name(obj::TCGScale)
    return jcall(obj, "getName", JString, ())
end

function offset_from_tai(obj::TCGScale, arg0::AbsoluteDate)
    return jcall(obj, "offsetFromTAI", jdouble, (AbsoluteDate,), arg0)
end

function offset_from_tai(obj::TCGScale, arg0::FieldAbsoluteDate)
    return jcall(obj, "offsetFromTAI", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function to_string(obj::TCGScale)
    return jcall(obj, "toString", JString, ())
end


text:
# This file holds the definition of the functions pertaining to the Retirement
#   type.

export  setRetirementCareerLength!,
        setRetirementAge!,
        setRetirementSchedule!,
        setRetirementIsEither!

"""
```
setRetirementCareerLength!(
    retirement::Retirement,
    careerLength::Real )
```
This function sets the maximum career length of the retirement scheme `retirement` to `careerLength`, where a career length of 0 indicates there is no limit.

If a career length < 0 is entered, the function gives a warning and makes no changes.

The function returns `true` if the maximum career length has been successfully set, and `false` if the entered career length is < 0.
"""
function setRetirementCareerLength!( retirement::Retirement,
    careerLength::Real )::Bool

    if careerLength < 0.0
        @warn "Max career length must be ⩾ 0, not making any changes."
        return false
    end  # if careerLength < 0.0

    retirement.maxCareerLength = careerLength
    return true

end  # setRetirementCareerLength!( retirement, careerLength )

"""
```
setRetirementAge!(
    retirement::Retirement,
    age::Real )
```
This function sets the retirement age of the retirement scheme `retirement` to `age`, where an age of 0 indicates there is no limit.

If a retirement age < 0 is entered, the function gives a warning and makes no changes.

The function returns `true` if the retirement age has been successfully set, and `false` if the entered age is < 0.
"""
function setRetirementAge!( retirement::Retirement, age::Real )::Bool

    if age < 0.0
        @warn "Retirement age must be ⩾ 0, not making any changes."
        return false
    end  # if age < 0.0

    retirement.retirementAge = age
    return true

end  # setRetirementAge!( retirement, age )

"""
```
setRetirementSchedule!(
    retirement::Retirement,
    freq::Real,
    offset::Real = 0.0 )
```
This function sets the schedule of the retirement scheme `retirement` one one check every `freq` time units, with an offset of `offset` w.r.t. the start of the simulation.

If the entered period is ⩽ 0, the function issues a warning and makes no changes.

The function returns `true` if the schedule has been successfully set, and `false` if the entered schedule period is ⩽ 0.
"""
function setRetirementSchedule!( retirement::Retirement, freq::Real, offset::Real = 0.0 )::Bool

    if freq <= 0.0
        @warn "Frequency of retirement schedule must be > 0, not making any changes."
        return false
    end  # if freq <= 0.0

    retirement.freq = freq
    retirement.offset = ( offset % freq ) + ( offset < 0.0 ? freq : 0.0 )
    return true

end  # setRetirementSchedule!( retirement, freq, offset )

"""
```
setRetirementIsEither!(
    retirement::Retirement,
    isEither::Bool )
```
This function sets the `isEither` flag of the retirement scheme `retirement` to `isEither`.

This function returns `true`, indicating that the flag has been successfully set.
"""
function setRetirementIsEither!( retirement::Retirement, isEither::Bool )::Bool

    retirement.isEither = isEither
    return true

end  # setRetirementIsEither!( retirement, isEither )

function Base.show( io::IO, retirement::Retirement )::Nothing

    print( io, "Retirement scheme" )

    if ( retirement.retirementAge == 0.0 ) &&
        ( retirement.maxCareerLength == 0.0 )
        print( io, " undefined" )
        return
    end  # if ( retirement.retirementAge == 0.0 ) && ...

    print( io, "\n  Retirement check occurs with a period of ",
        retirement.freq, " (offset ", retirement.offset, ")" )

    print( io, "\n  Retirement occurs " )

    if retirement.retirementAge == 0.0
        print( io, "after a career of length ", retirement.maxCareerLength )
    elseif retirement.maxCareerLength == 0.0
        print( io, "at age ", retirement.retirementAge )
    else
        print( io, "at age ", retirement.retirementAge,
            retirement.isEither ? " or" : " and", " after a career of length ",
            retirement.maxCareerLength )
    end  # if retirement.retirementAge == 0.0

    return

end  # show( io, retirement )


text:
struct Simplex{D,T<:Number,U<:Complex,V}
  vertices::Vector{Vertex{T,U,V}}
  permabs::Vector{Int}
  function Simplex(vertices::Vector{Vertex{T,U,V}},
      permabs::Vector{Int}=zeros(Int64, length(vertices))
      ) where {T,U<:Complex,V}
    D = length(vertices) - 1
    output = new{D,T,U,V}(vertices, permabs)
    sort!(output)
    return output
  end
end

function Simplex(f::T, ic::U, initial_step::Number
    ) where {T, U<:AbstractVector{<:Number}}
  return Simplex(f, ic, fill!(similar(ic), initial_step))
end

function Simplex(f::T, ic::AbstractVector{U}, initial_steps::AbstractVector{V}
    ) where {T, U<:Number, V<:Number}
  return Simplex(f, vertexpositions(ic, initial_steps))
end
function Simplex(f::T, positions::U
    ) where {T, W<:Number, V<:AbstractVector{W}, U<:AbstractVector{V}}
  if length(unique(length.(positions))) != 1
    throw(ArgumentError("All entries in positions $positions must be the same
                        length"))
  end
  vertices = [Vertex(p, f(p)) for p in positions]
  return Simplex(vertices)
end

import Base: length, iterate, iterate, getindex
import Base: eachindex, sort!, hash, extrema
Base.length(s::Simplex{D}) where {D} = D + 1
Base.iterate(s::Simplex) = iterate(s.vertices)
Base.iterate(s::Simplex, counter) = iterate(s.vertices, counter)
Base.getindex(s::Simplex, index) = s.vertices[index]

sortby(s::Simplex) = v->angle(value(v))

function sortby(s::Simplex{2})
  return function output(v)
    c = centre(s)
    pv = position(v)
    return @inbounds atan(pv[2] - c[2], pv[1] - c[1])
  end
end

function Base.sort!(s::Simplex)
  sort!(s.vertices, by=sortby(s), alg=InsertionSort)
  sortperm!(s.permabs, s.vertices, by=x->abs(value(x)), alg=InsertionSort)
  return nothing
end

issortedbyangle(s::Simplex) = issorted(s, by=sortby(s))

function Base.extrema(s::Simplex)
  return [extrema(position(v)[i] for v in s) for i in 1:dimensionality(s)]
end

dimensionality(s::Simplex{D}) where {D} = D

selectabs(s, index) = @inbounds s.vertices[s.permabs[index]]

bestvertex(s::Simplex) = selectabs(s, 1)
worstvertex(s::Simplex) = selectabs(s, length(s))
secondworstvertex(s::Simplex) = selectabs(s, length(s) - 1)

function centroidposition(s::Simplex, ignoredvertex=worstvertex(s))
  g(v) = isequal(v, ignoredvertex) ? zero(position(v)) : position(v)
  return mapreduce(g, +, s) / (length(s) - 1)
end

centre(s::Simplex) = mapreduce(position, +, s) / length(s)

function hypervolume(s::Simplex)
  m = hcat((vcat(position(v), 1) for v in s)...)
  d = dimensionality(s)
  return abs(det(m)) / factorial(d)
end

function swap!(s::Simplex, this::Vertex, forthat::Vertex)
  @assert this ∈ s.vertices
  s.vertices[findfirst(x -> isequal(x, this), s.vertices)] = forthat
  sort!(s)
  return nothing
end

swapworst!(s::Simplex, forthis::Vertex) = swap!(s, worstvertex(s), forthis)

function closestomiddlevertex(s::Simplex)
  mid = mapreduce(position, +, s) ./ length(s)
  _, index = findmin(map(v->sum(abs, position(v) - mid), s))
  return s[index]
end

struct AssessConvergenceGraph
  toprocess::BitSet
  processed::BitSet
  connectedto::BitSet
end
AssessConvergenceGraph() = AssessConvergenceGraph(BitSet(), BitSet(), BitSet())

function assessconvergence(simplex, config, asg=AssessConvergenceGraph())

  if abs(value(bestvertex(simplex))) <= config[:stopval]
    return :STOPVAL_REACHED
  end

  toprocess = asg.toprocess # toprocess is empty to exit the while loop
  processed = empty!(asg.processed)
  connectedto = asg.connectedto # connectedto is emptied at top of while loop
  push!(toprocess, 1) # now only contains 1
  @inbounds while !isempty(toprocess)
    vi = pop!(toprocess)
    v = simplex.vertices[vi]
    pv = position(v)
    empty!(connectedto)
    for (qi, q) ∈ enumerate(simplex)
      thisxtol = true
      pq = position(q)
      @inbounds for i ∈ 1:dimensionality(simplex)
        thisxtol &= isapprox(pv[i], pq[i],
          rtol=config[:xtol_rel][i], atol=config[:xtol_abs][i])
        thisxtol || break
      end
      thisxtol && push!(connectedto, qi)
      thisxtol && for i in connectedto if i ∉ processed push!(toprocess, i) end end
    end
    push!(processed, vi)
  end
  allxtol = all(i ∈ processed for i ∈ 1:length(simplex))
  allxtol && return :XTOL_REACHED

  allftol = true
  @inbounds for (vi, v) ∈ enumerate(simplex)
    for qi ∈ vi+1:length(simplex)
      q = simplex.vertices[qi]
      allftol &= all(isapprox(value(v), value(q),
                              rtol=config[:ftol_rel], atol=config[:ftol_abs]))
      position(v) == position(q) && return :XTOL_DEGENERATE_SIMPLEX
    end
  end
  allftol && return :FTOL_REACHED

  return :CONTINUE
end

function _πtoπ(ϕ::T) where {T}
  ϕ < -T(π) && return _πtoπ(ϕ + 2π)
  ϕ >= T(π) && return _πtoπ(ϕ - 2π)
  return ϕ
end

function windingangle(s::Simplex{D,T,U}) where {D,T,U}
  θ = zero(real(U))
  @inbounds for i in 1:length(s)
    θ += _πtoπ(angle(value(s.vertices[mod1(i+1, length(s))])) -
               angle(value(s.vertices[i])))
  end
  return θ
end

function windingnumber(s::Simplex)
  radians = windingangle(s)
  return isfinite(radians) ? Int64(round(radians / 2π)) : Int64(0)
end


text:
"""
    test direct policy optimization
"""
prob_dpo = prob_dpo.prob
# objective
z0 = rand(prob_dpo.num_var)
tmp_o(z) = eval_objective(prob_dpo,z)
∇j = zeros(prob_dpo.num_var)
eval_objective_gradient!(∇j, z0, prob_dpo)
@assert norm(ForwardDiff.gradient(tmp_o, z0) - ∇j) < 1.0e-10

# sample dynamics
con_dynamics = sample_dynamics_constraints(prob_dpo.prob,
    prob_dpo.N, prob_dpo.N + prob_dpo.D)
c0 = zeros(con_dynamics.n)
_c(c,z) = constraints!(c, z, con_dynamics,
    prob_dpo.prob, prob_dpo.idx, prob_dpo.N, prob_dpo.D,
    prob_dpo.dist, prob_dpo.sample)
∇c_fd = zeros(con_dynamics.n, prob_dpo.num_var)
FiniteDiff.finite_difference_jacobian!(∇c_fd, _c, z0)
# ∇c_fd = ForwardDiff.jacobian(_c, c0, z0)
spar = constraints_sparsity(con_dynamics, prob_dpo.prob, prob_dpo.idx,
    prob_dpo.N, prob_dpo.D)
∇c_vec = zeros(length(spar))
∇c = zeros(con_dynamics.n, prob_dpo.num_var)
constraints_jacobian!(∇c_vec, z0, con_dynamics,
    prob_dpo.prob, prob_dpo.idx, prob_dpo.N, prob_dpo.D,
    prob_dpo.dist, prob_dpo.sample)
for (i,k) in enumerate(spar)
    ∇c[k[1],k[2]] = ∇c_vec[i]
end
@assert norm(vec(∇c) - vec(∇c_fd)) < 1.0e-5
@assert sum(∇c) - sum(∇c_fd) < 1.0e-5

# policy
prob_dpo.con[2]
model.m * (T - 1) * (N + 1)
con_policy = policy_constraints(prob_dpo.prob, prob_dpo.N)
c0 = zeros(con_policy.n)
_c(c,z) = constraints!(c, z, con_policy, prob_dpo.policy,
    prob_dpo.prob, prob_dpo.idx, prob_dpo.N)

_c(c0,z0)
# ∇c_fd = zeros(con_policy.n, prob_dpo.num_var)
# FiniteDiff.finite_difference_jacobian!(∇c_fd, _c, z0)
∇c_fd = ForwardDiff.jacobian(_c, c0, z0)
spar = constraints_sparsity(con_policy, prob_dpo.policy, prob_dpo.prob,
    prob_dpo.idx, prob_dpo.N)
∇c_vec = zeros(length(spar))
∇c = zeros(con_policy.n, prob_dpo.num_var)
constraints_jacobian!(∇c_vec, z0, con_policy, prob_dpo.policy,
    prob_dpo.prob, prob_dpo.idx, prob_dpo.N)

for (i,k) in enumerate(spar)
    ∇c[k[1],k[2]] = ∇c_vec[i]
end
@assert norm(vec(∇c) - vec(∇c_fd)) < 1.0e-10
@assert sum(∇c) - sum(∇c_fd) < 1.0e-10

# constraints
c0 = zeros(prob_dpo.num_con)
_c(c,z) = eval_constraint!(c, z, prob_dpo)

∇c_fd = zeros(prob_dpo.num_con, prob_dpo.num_var)
FiniteDiff.finite_difference_jacobian!(∇c_fd, _c, z0)
spar = sparsity_jacobian(prob_dpo)
∇c_vec = zeros(length(spar))
∇c = zero(∇c_fd)
eval_constraint_jacobian!(∇c_vec, z0, prob_dpo)
for (i,k) in enumerate(spar)
    ∇c[k[1],k[2]] = ∇c_vec[i]
end
@assert norm(vec(∇c) - vec(∇c_fd)) < 1.0e-5
@assert sum(∇c) - sum(∇c_fd) < 1.0e-5

# sample dynamics
t = 10
xt = view(z0_dpo, prob_dpo.idx.xt[t])
ut = view(z0_dpo, prob_dpo.idx.ut[t])
μ = view(z0_dpo, prob_dpo.idx.mean[prob_dpo.prob.mean.idx.x[t]])
ν = view(z0_dpo, prob_dpo.idx.mean[prob_dpo.prob.mean.idx.u[t]])

propagate_dynamics(model, rand(model.n), rand(model.m), rand(model.d), h, t)
propagate_dynamics_jacobian(model, rand(model.n), rand(model.m), rand(model.d), h, t)

sample_dynamics(model, xt, ut, μ, ν, prob_dpo.dist.w, h, t,
	prob_dpo.sample.β)

a1, a2, a3, a4 = sample_dynamics_jacobian(model, xt, ut, μ, ν, prob_dpo.dist.w, h, t,
	prob_dpo.sample.β)

sdx(y) = sample_dynamics(model, y, ut, μ, ν, prob_dpo.dist.w, h, t,
	prob_dpo.sample.β)[1]
sdx(xt)

sdu(y) = sample_dynamics(model, xt, y, μ, ν, prob_dpo.dist.w, h, t,
	prob_dpo.sample.β)[1]
sdu(ut)

sdμ(y) = sample_dynamics(model, xt, ut, y, ν, prob_dpo.dist.w, h, t,
	prob_dpo.sample.β)[1]
sdμ(μ)

sdν(y) = sample_dynamics(model, xt, ut, μ, y, prob_dpo.dist.w, h, t,
	prob_dpo.sample.β)[1]
sdν(ν)

@assert norm(FiniteDiff.finite_difference_jacobian(sdx, xt) - a1) < 1.0e-5
@assert norm(FiniteDiff.finite_difference_jacobian(sdu, ut) - a2) < 1.0e-5
@assert norm(FiniteDiff.finite_difference_jacobian(sdμ, μ) - a3) < 1.0e-5
@assert norm(FiniteDiff.finite_difference_jacobian(sdν, ν) - a4) < 1.0e-5

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
"""
	similarity_order(wv::WordVectors, alpha=-0.65)

Post-processes the word embeddings `wv` so that the embeddings capture more
information than directly apparent through a linear transformation that
adjusts the similarity order of the model. The function returns a new
`WordVectors` object containing the processed embeddings.

# Arguments
  * `wv::WordVectors` the word embeddings
  # `alpha::AbstractFloat` the `α` parameter of the algorithm (default -0.65)

# References:
  * [Artetxe et al. \"Uncovering divergent linguistic information in
     word embeddings with lessons for intrinsic and extrinsic evaluation\",
     2018](https://arxiv.org/pdf/1809.02094.pdf)
"""
function similarity_order(wv::WordVectors{S,T,H},
                          alpha::T=T(-0.65)
                         ) where {S<:AbstractString, T<:Real, H<:Integer}
    n = length(wv.vocab)
    X = wv.vectors
    M = X*X'
    L, Q = eigen(M)
    Wₐ= Q .* (L.^alpha)
    return WordVectors{S,T,H}(wv.vocab, Wₐ*X, wv.vocab_hash)
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
"""
`estimate_entrainment_area(arnold_tongue, I_max)`

Estimate entrainment area of the Arnold tongue.

**Arguments**
- `arnold_tongue`: A DataFrame representing an Arnold tongue. The DataFrame has
    two columns names `amplitudes` and `periods`. Periods indicate a position
    at x-axis and amplitudes indices a maximal height at y-axis at which was the
    system entrained.
- `I_max`: Maximal input period that was considered for the Arnold tongue.
"""
function estimate_entrainment_area(arnold_tongue, I_max)
    amplitudes = arnold_tongue[:, "amplitudes"]
    entrainment_area = 1 - mean(amplitudes ./ I_max)
    return entrainment_area
end

"""
`plot_arnold_tongue(arnold_tongue, I_max; ax=gca())`

Plot the Arnold tongue.

**Arguments**
- `arnold_tongue`: A DataFrame representing an Arnold tongue. The DataFrame has
    two columns names `amplitudes` and `periods`. Periods indicate a position
    at x-axis and amplitudes indices a maximal height at y-axis at which was the
    system entrained.
- `I_max`: Maximal input period that was considered for the Arnold tongue.

**Keyword Arguments**
- `ax`: PyPlot axes.
"""
function plot_arnold_tongue(arnold_tongue, I_max; ax=gca())
    periods = arnold_tongue[:, "periods"]
    amplitudes = arnold_tongue[:, "amplitudes"]
    n = length(amplitudes)
    arnold_tongue_matrix = fill(0.0, n, n)

    for i_period = 1:n
        for i_amplitude = 1:n
            amplitude = amplitudes[i_period]
            entrained = (I_max * (i_amplitude/n) ) < amplitude
            arnold_tongue_matrix[i_amplitude, i_period] = entrained
        end
    end
    amplitude_ticks = range(0, I_max, length=n)
    ax.pcolor(periods, amplitude_ticks, arnold_tongue_matrix, rasterized=true)
    ax.set_xticks(range(minimum(periods), maximum(periods), length=3))
    ax.set_yticks(range(minimum(amplitude_ticks), maximum(amplitude_ticks),
        length=4))
    ax.set_xlabel("Input period")
    ax.set_ylabel("Input amplitude")
    ax.set_title("Arnold tongue")

    return ax
end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text: @test UncertainValueDataset(rand(30)) isa UncertainValueDataset


text:
# DP formulation for additively separable constraints
# State is Array{Real}: right-hand sides of constraints

include("mathprog_eval.jl")

"""^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Transition function for additively separable constraints

Input:
evals: An array of decomposed terms of a constraint
doms: Array of range representing domain of variables
ords: Function for ordering of the variables

Output: Transition function
"""
function create_constraint_transition_function(evals::SeparableFunctionEvaluation, doms::Array{Range}, ords::Ordering)
	@assert(length(doms) == length(ords)) # number of variables are consistent
	n = length(doms)	# number of variables
	min_ahead = Array{Real}(undef, n)		# For each variable, stores the minimum value that can be obtained starting from the next layer to the terminal
	max_ahead = Array{Real}(undef, n)		# For each variable, stores the maximum value that can be obtained starting from the next layer to the terminal

	#initialize values for the last layer
	min_ahead[get_var(ords,n)] = 0
	max_ahead[get_var(ords,n)] = 0
	for i in n-1:-1:1
		var = get_var(ords,i)		# gets the variable of current layer i
		var_next = get_var(ords,i+1)		# gets the variable in layer i+1
		if isassigned(evals.separate_functions, var_next)		#checks if the function is defined for the variable
			min_v, max_v = min_max_calculator(evals.separate_functions[var_next], doms[var_next])		#computes the min and max of the function over varibale domain
			min_ahead[var] = min_v + min_ahead[var_next]
			max_ahead[var] = max_v + max_ahead[var_next]
		else
			min_ahead[var] = min_ahead[var_next]	#if the function is not defined for the next layer, it does not affect the min max values
			max_ahead[var] = max_ahead[var_next]
		end
	end
	return function (state::Real, variable::Int, value::Int)
		@assert(value in doms[variable])		#ensures that the value belongs to the variable domain

		if isassigned(evals.separate_functions, variable)
			new_state = state + evals.separate_functions[variable](value)		#computes the new state value
		else
			new_state = state 		#if the function is not defined for the variable, the state value remains unchanged
		end

		#To determine feasibility, it is assumed that the constraint form is f(x) <= 0
		if new_state <= -max_ahead[variable]			#new state is feasible no matter what variable assignments are chosen from now on
			new_state = -max_ahead[variable]
		elseif new_state > -min_ahead[variable]			#new state is infeasible no matter what variable assignments are chosen
			new_state = nothing
		end

		return new_state
	end
end

"""^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Initial state for additively separable constraints

Input:
evals: An array of decomposed terms of a constraint

"""
function create_constraint_initial_state(evals::SeparableFunctionEvaluation)::Real
    return eval.constant  	# constant value of the constraint as moved to lhs of <= constraints
end

"""^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Domain function for additively separable constraints

Input:
m: math prog model

Output:
domain function that takes a variable index and returns the domain as range
"""
function create_constraint_domain_function(m::Model)
	nvars = MathProgBase.numvar(m)
	ranges = Array{Range}(undef, nvars)
	for v in 1:nvars
		v_obj = Variable(m, v)
		ranges[v] = ceil(Int, getlowerbound(v_obj)):floor(Int, getupperbound(v_obj))
	end
	return ranges
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra, Random
using ArraysOfArrays, Distributions, PDMats, StatsBase

struct _TestDensityStruct <: AbstractDensity
end

@testset "abstract_density" begin

    BAT.nparams(td::_TestDensityStruct) = Int(3)
    BAT.sampler(td::_TestDensityStruct) = BAT.sampler(MvNormal(ones(3), PDMat(Matrix{Float64}(I,3,3))))

    @testset "rand" begin
        td = _TestDensityStruct()
        @test rand(MersenneTwister(7002), sampler(td)) ≈ [-2.415270938, 0.7070171342, 1.0224848653]
    end

    # ToDo: Add more tests
end


text:
module Generate

using LightGraphs
using MetaGraphs
using Random
import Base:occursin
using StatsBase

include("../utils/filters.jl")
include("../utils/maximalcliques.jl")
include("../utils/misc.jl")
include("../utils/chordalextension.jl")
include("../utils/constructgraph.jl")
include("../utils/graph.jl")
include("../utils/merge.jl")
include("filteredges.jl")
include("addedges.jl")
include("./Decompositions/Cholesky.jl")
include("./Decompositions/EG.jl")
include("./Decompositions/MD.jl")
include("./Decompositions/LEX_M.jl")
include("./Decompositions/MCS_M.jl")
include("./Decompositions/LC.jl")
include("./Decompositions/Trig.jl")
include("./decomposition.jl")
include("./Merge/merge.jl")
include("./Merge/molzahn.jl")
include("./Merge/sliwak.jl")
include("./merge.jl")
include("./combine.jl")

end # module

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
include("PsPot_HGH.jl")
include("ylmr2.jl")
include("../common/PWGrid_v02.jl")
include("../common/calc_strfact_v2.jl")

function test_main()

  #
  # Setup plane wave basis
  #
  Ns = [32, 32, 32]
  LatVecs = 18.0*diagm(ones(3))
  pw = PWGrid( Ns, LatVecs )
  Ω = pw.Ω

  # use G-vectors for wavefunction
  # CAUTION: this should be depend on k-points
  idx_gw2r = pw.gvecw.idx_gw2r
  G = pw.gvec.G[:,idx_gw2r]
  Gm = sqrt( pw.gvec.G2[idx_gw2r] )

  Ngwx = pw.gvecw.Ngwx
  Ng = pw.gvec.Ng

  #
  # Setup structure and pseudopotentials
  #
  const Nspecies = 2
  const Natoms   = 3
  atmsymb = ["Cu", "Mo"] # unique list of atomic symbols
  atm2species = [1, 2, 1]  # mapping from atom to species

  # not really used
  Zv = zeros(Float64,Nspecies)
  Zv = [11.0, 14.0]

  Xpos = reshape( [0.0, 0.0, 0.0,  # Cu
                   1.5, 0.0, 0.0,  # Mo
                   0.0, 1.5, 0.0], # Cu
                   (3,Natoms) )

  # Calculate structure factor
  Sf = calc_strfact( Xpos, Nspecies, atm2species, pw.gvec.G ) # Sf[Ng,Nspecies]

  psp = Array{PsPot_HGH,1}(Nspecies)

  psp[1] = PsPot_HGH(1, "Cu", "LDA_HGH/29cu.11.hgh")
  psp[2] = PsPot_HGH(2, "Mo", "LDA_HGH/42mo.14.hgh")

  # for multispecies system this should be max over all species
  lmaxkb = 0
  for isp = 1:Nspecies
    if lmaxkb < psp[isp].lmax
      lmaxkb = psp[isp].lmax
    end
  end
  @printf("lmaxkb = %d\n", lmaxkb)

  ylm = ylmr2( lmaxkb, G )

  # from init_us_1

  NH = zeros(Int,Nspecies)

  for isp = 1:Nspecies
    ih = 0
    println("snprj = ", psp[isp].snprj)
    for ib = 1:psp[isp].snprj
      l = psp[isp].lll[ib]
      println("l = ", l)
      for m = 1:(2*l+1)
        #nhtol[ih,isp] = l
        #nhtolm[ih,isp] = l*l + m
        #indv[ih,isp] = ib
        ih = ih + 1
      end
    end
    NH[isp] = ih
    @printf("isp = %d, ih = %d\n", isp, ih)
  end

  Nhmax = maximum(NH)

  println("Nhmax = ", Nhmax)

  Nkb = 0
  for ia = 1:Natoms
    isp = atm2species[ia]
    Nkb = Nkb + NH[isp]
  end
  println("Nkb   = ", Nkb)

  nhtol  = zeros(Int,Nhmax,Nspecies)
  nhtolm = zeros(Int,Nhmax,Nspecies)
  indv   = zeros(Int,Nhmax,Nspecies)

  for isp = 1:Nspecies
    ih = 1
    for ib = 1:psp[isp].snprj
      l = psp[isp].lll[ib]
      for m = 1:(2*l+1)
        nhtol[ih,isp] = l
        nhtolm[ih,isp] = l*l + m
        indv[ih,isp] = ib
        ih = ih + 1
      end
    end
  end

  vkb1 = zeros(Float64,Ngwx,Nhmax)
  Vkb = zeros(Complex128,Ngwx,Nkb)

  jkb = 0

  for isp = 1:Nspecies

    snprj = psp[isp].snprj
    @printf("\nSpecies = %d\n", isp)
    for ib = 1:snprj
      l = psp[isp].lll[ib]
      iproj = psp[isp].ipr[ib]
      # remember that Gm (the magnitudes) are used, not G
      Vprj = eval_HGH_proj_G( psp[isp], l, iproj, Gm, Ω )
      @printf("l = %d, iproj = %d done\n", l, iproj)

      for ih = 1:NH[isp]
        if ib == indv[ih,isp]
          lm = nhtolm[ih,isp]
          @printf("Same ib = %d, lm = %d\n", ib, lm)
          for ig=1:Ngwx
            vkb1[ig,ih] = ylm[ig,lm]*Vprj[ig]
          end
          # XXXXXX  need to zero out vkb1 ???
        end
      end
    end

    for ia = 1:Natoms
      if atm2species[ia] == isp
        # arg = dot( kpts[:], atpos[:,ia] ) )  # needed for kpts other than gamma (0,0,0)
        phase = 1.   # special case for gamma point sampling only
        for ih = 1:NH[isp]
          jkb = jkb + 1
          pre = (-1.0*im)^nhtol[ih,isp]*phase
          for ig = 1:Ngwx
            Vkb[ig,jkb] = vkb1[ig,ih]*Sf[idx_gw2r[ig],isp]*pre
          end
          # need to zero out multiple k-points
        end
      end
    end

  end # loop over species

end

test_main()


text:
# EXCLUDE FROM TESTING
# NOTE: This is temporarily not working until we merge the DMSTAG routines with the new Clang branch
#
# This implements src/snes/examples/tutorials/ex2.c from PETSc using the PETSc.jl package, using SNES
#
# This is the same as SNES_ex2b.j, except that we show how automatic differentiation can be used to
# compute the jacobian.
#
# Newton method to solve u'' + u^{2} = f, sequentially.

using PETSc, MPI, LinearAlgebra, SparseArrays, Plots, ForwardDiff

if ~MPI.Initialized()
    MPI.Init()
end

PETSc.initialize()

```
    Computes initial guess
```
function FormInitialGuess!(dm,x)

    x_Local = PETSc.DMCreateLocalVector(dm);
    x_array = PETSc.DMStagVecGetArray(dm,x_Local);
    x_array[1:2,1:2] .= 1;
    PETSc.DMLocalToGlobal(dm, x_Local, PETSc.INSERT_VALUES, x)
end

```
    Computes the residual f, given solution vector x
```
function FormResidual!(f,x)
    n       =   length(x);
    xp      =   LinRange(0.0,1.0, n);
    F       =   6.0.*xp .+ (xp .+1.e-12).^6.0;      # define source term function

    dx      =   1.0/(n-1.0);
    f[1]    =   x[1] - 0.0;
    for i=2:n-1
        f[i] = (x[i-1] - 2.0*x[i] + x[i+1])/dx^2 + x[i]*x[i] - F[i]
    end
    f[n]    =   x[n] - 1.0;

end

```
    Wrapper which makes it easier to compute the jacobian using automatic differntiation
```
function  ForwardDiff_routine(x)

    f   = zero(x)               # vector of zeros, of same type as e
    FormResidual!(f,x);

    return f;
end

```
    This copies a julia sparse matrix to PETSc MatSeqAIJ format
```
function Mat_JuliaToPETSc!(J::PETSc.MatSeqAIJ, J_julia::SparseMatrixCSC)

    for i = 1:size(J_julia,1)
        col = J_julia[i,:];
        row = ones(Int32,length(col.nzind))*i;
        for j=1:length(col.nzind)
            J[i, col.nzind[j]] = col.nzval[j];
        end
    end
    PETSc.assemble(J);  # finalize assembly

end

```
    Computes the jacobian, given solution vector x
```
function FormJacobian!(x, args...)

    J        =   args[1];        # preconditioner = args[2], in case we want it to be different from J

    # Use AD to compute jacobian; by transferring x into sparse, the output will be sparse
    J_julia  =   ForwardDiff.jacobian(ForwardDiff_routine,sparse(x));

    if typeof(J) <: PETSc.AbstractMat
        Mat_JuliaToPETSc!(J, J_julia);  # transfer julia sparse matrix 2 petsc
    else
        J .= J_julia;
    end
end

# ==========================================
# Main code

# Compute initial solution
nx   =   5;
nz   =   5;
x0   =   -5;
xend =    5;
z0   = - 10;
zend =    0;

dm  = PETSc.DMStagCreate2d(MPI.COMM_SELF,PETSc.DM_BOUNDARY_NONE,PETSc.DM_BOUNDARY_NONE,nx,nz,1,1,0,0,1,PETSc.DMSTAG_STENCIL_BOX,2);
PETSc.DMStagSetUniformCoordinatesExplicit(dm, x0, xend, z0, zend);
x   =   PETSc.DMCreateGlobalVector(dm);

#FormInitialGuess!(dm,x);

# Compute initial jacobian using a julia structure to obtain the nonzero structure
# Note that we can also obtain this structure in a different manner
#Jstruct  = zeros(n,n);
#FormJacobian!(x, Jstruct);                              # jacobian in julia form
#Jsp      =   sparse(Float64.(abs.(Jstruct) .> 0))       # sparse julia, with 1.0 in nonzero spots
#PJ       =   PETSc.MatSeqAIJ(Jsp);                      # transfer to PETSc (initialize matrix with correct nonzero pattern)

# Setup SNES
#x_s = PETSc.VecSeq(x);                  # solution vector
#res = PETSc.VecSeq(zeros(size(x)));     # residual vector

#S = PETSc.SNES{Float64}(MPI.COMM_SELF;
#        snes_rtol=1e-12,
#        snes_monitor=true,
#        snes_converged_reason=true);
#PETSc.setfunction!(S, FormResidual!, res)
#PETSc.setjacobian!(S, FormJacobian!, PJ, PJ)

# solve
#PETSc.solve!(x_s, S);

# Extract & plot solution
#x_sol = x_s.array;                  # convert solution to julia format
#FormResidual!(res.array,x_sol)      # just for checking, compute residual
#@show norm(res.array)

#plot(LinRange(0,1,n),x_sol,xlabel="width",ylabel="solution")

#PETSc.finalize()


text:
@test fire(dfsm, "can", "dance") == true
@test fire(dfsm, "can", "hop") == false
@test_throws ErrorException fire(dfsm, "skip")
@test_throws ErrorException fire(dfsm, "startup")

fire(dfsm, "dance")

@test dfsm.current == "first"
@test fire(dfsm, "can", "dance") == false
@test fire(dfsm, "can", "hop") == true
@test_throws ErrorException fire(dfsm, "skip")
@test_throws ErrorException fire(dfsm, "dance")

fire(dfsm, "hop")

@test dfsm.current == "second"

fire(dfsm, "skip")

@test dfsm.current == "third"

fire(dfsm, "sleep")

@test dfsm.current == "third"

@test fire(fsm, "is", "first") == true
@test fire(fsm, "can", "hop") == true
@test fire(fsm, "can", "jump") == false
@test fire(fsm, "finished") == false
@test_throws ErrorException fire(fsm, "jump")
@test_throws ErrorException fire(fsm, "skip")
@test_throws ErrorException fire(fsm, "startup")
@test_throws ErrorException fire(fsm, "humblebrag")

fire(fsm, "hop")

@test fsm.current == "second"
@test fire(fsm, "is", "first") == false
@test fire(fsm, "is", "second") == true
@test fire(fsm, "can", "hop") == false
@test fire(fsm, "can", "skip") == true
@test fire(fsm, "finished") == false
@test_throws ErrorException fire(fsm, "jump")
@test_throws ErrorException fire(fsm, "hop")
@test_throws ErrorException fire(fsm, "startup")
@test_throws ErrorException fire(fsm, "humblebrag")

fire(fsm, "skip")

@test fsm.current == "third"
@test fire(fsm, "is", "third") == true
@test fire(fsm, "can", "jump") == true
@test fire(fsm, "can", "skip") == false
@test fire(fsm, "finished") == false
@test_throws ErrorException fire(fsm, "hop")
@test_throws ErrorException fire(fsm, "skip")
@test_throws ErrorException fire(fsm, "startup")
@test_throws ErrorException fire(fsm, "humblebrag")

fire(fsm, "jump")

@test fsm.current == "fourth"
@test fire(fsm, "is", "fourth") == true
@test fire(fsm, "can", "hop") == false
@test fire(fsm, "can", "jump") == false
@test fire(fsm, "finished") == true
@test_throws ErrorException fire(fsm, "jump")
@test_throws ErrorException fire(fsm, "skip")
@test_throws ErrorException fire(fsm, "hop")
@test_throws ErrorException fire(fsm, "startup")
@test_throws ErrorException fire(fsm, "humblebrag")

fire(mfsm, "hop")

@test mfsm.current == "second"

fire(mfsm, "skip")

@test mfsm.current == "third"

fire(mfsm, "sleep")

@test mfsm.current == "third"

fire(mfsm, "hop")

@test mfsm.current == "second"


text:
module test_bukdu_routing

using Test

function penetrate_segments(segments)
    vals = Expr[]
    path_params = Expr[]
    for seg in segments
        if startswith(seg, ':')
            name = seg[2:end]
            mangled = Symbol(name, :_)
            pair = :(Pair($name, String(first(typeof($mangled).parameters))))
            push!(path_params, pair)
            expr = :($mangled::Any)
        else
            expr = :(::Val{Symbol($seg)})
        end
        push!(vals, expr)
    end
    return (vals, path_params)
end

struct Conn
end

abstract type ApplicationController end

struct WelcomeController <: ApplicationController
    conn::Conn
end

struct Route
    C::Type{<:ApplicationController}
    action
    path_params::Dict{String, String}
end

index(::WelcomeController) = "hello"

C = WelcomeController
action = index

url = "/a/:b/:c"
segments = split(url, '/')
(vals, path_params) = penetrate_segments(segments)
method = "GET"
@eval route(::Val{Symbol($method)}, $(vals...)) = Route($C, $action, Dict($(path_params...)))

target_path = "/a/25/36"
segments = split(target_path, '/')
vals = [Val(Symbol(seg)) for seg in segments]
r = route(Val(Symbol(method)), vals...)

@test r.C == WelcomeController
@test r.action == index
@test r.path_params == Dict("c"=>"36","b"=>"25")

end # module test_bukdu_routing

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module DocCheckTests

using Test

using Markdown
using Documenter.DocChecks: linkcheck
using Documenter.Documents

@testset "DocChecks" begin
    # The linkcheck tests are currently not reliable on CI, so they are disabled.
    @testset "linkcheck" begin
        src = md"""
            [HTTP (HTTP/1.1) success](http://www.google.com)
            [HTTPS (HTTP/2) success](https://www.google.com)
            [FTP success](ftp://ftp.iana.org/tz/data/etcetera)
            [FTP (no proto) success](ftp.iana.org/tz/data/etcetera)
            [Redirect success](google.com)
            [HEAD fail GET success](https://codecov.io/gh/invenia/LibPQ.jl)
            """

        Documents.walk(Dict{Symbol, Any}(), src) do block
            doc = Documents.Document(; linkcheck=true)
            result = linkcheck(block, doc)
            @test_skip doc.internal.errors == Set{Symbol}()
            result
        end

        src = Markdown.parse("[FILE failure](file://$(@__FILE__))")
        doc = Documents.Document(; linkcheck=true)
        Documents.walk(Dict{Symbol, Any}(), src) do block
            linkcheck(block, doc)
        end
        @test_skip doc.internal.errors == Set{Symbol}([:linkcheck])
    end
end

end


text:
using DataStructures
using Test
using ChemistryFeaturization
include("../../src/featurize.jl")
using .featurize: weave_featurize

mol = chem.MolFromSmiles("C1=CC=CC=C1")  # benzene
carbon = get(mol.GetAtoms(),0)
bond = get(mol.GetBonds(),0)

@testset "weave featurization" begin
    weave_mol = weave_featurize(["C1=CC=CC=C1"])[1]
    @test sum(weave_mol.atom_features[end,:]) == sum([sum(v) for v in values(smiles_atom_features(carbon))])
    correct_pair_feat_sum = sum([sum(v) for v in values(smiles_bond_features(bond))]) + 2
    count = 0
    for i in 1:size(weave_mol.pair_features,1)
        if sum(weave_mol.pair_features[i,:]) == correct_pair_feat_sum
            count += 1
        end
    end
    @test count == 12
    @test weave_mol.num_atoms == 6
    @test weave_mol.atom_feat_num == 75

    @test_logs (:warn,"No bonds found in C! Pair feature matrix filled with zeros!") weave_featurize(["C"])[1]

end


text:
"""
Dynamic path planning problem of hovercraft trying to minimize thrust usage
to hit all the waypoints at the alotted target times.
"""

using InfiniteOpt, JuMP, Ipopt, PyPlot

# Waypoints
xw = [1 4 6 1; 1 3 0 1] # positions
tw = [0, 25, 50, 60]    # times

# Initialize the model
m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

# Set the parameters and variables
@infinite_parameter(m, t in [0, 60], num_supports = 61)
@infinite_variable(m, x[1:2](t), start = 1) # position
@infinite_variable(m, v[1:2](t), start = 0) # velocity
@infinite_variable(m, u[1:2](t), start = 0) # thruster input

# Specify the objective
@objective(m, Min, integral(u[1]^2 + u[2]^2, t))

# Set the initial conditions
@BDconstraint(m, initial_velocity[i = 1:2](t == 0), v[i] == 0)

# Define the point physics
@constraint(m, [i = 1:2], deriv(x[i], t) == v[i])
@constraint(m, [i = 1:2], deriv(v[i], t) == u[i])

# Hit all the waypoints
@BDconstraint(m, [i = 1:2, j = eachindex(tw)](t == tw[j]), x[i] == xw[i, j])

# Optimize the model
optimize!(m)

# Get the results
if has_values(m)
    x_opt = value.(x)
    v_opt = value.(v)
    u_opt = value.(u)
    opt_supports = supports.(x)
    obj_opt = objective_value(m)
end

# Plot the results
figure()
scatter(xw[1,:], xw[2,:], label = "Waypoints")
plot(x_opt[1], x_opt[2], label = "Trajectory", color = "C1")
xlabel(L"$x_1$")
ylabel(L"$x_2$")
legend()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct SoftwareApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_apiV1FleetSoftwareGet(_api::SoftwareApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Nothing, "/api/v1/fleet/software", ["bearerAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""Software

Params:

Return: Nothing
"""
function apiV1FleetSoftwareGet(_api::SoftwareApi; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetSoftwareGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetSoftwareGet(_api::SoftwareApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetSoftwareGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export apiV1FleetSoftwareGet


text:
#clearconsole();
using LightGraphs
using GraphPlot
using MetaGraphs
using GraphRecipes#
using Colors
using ColorSchemes
using Plots
#gr()
#Plots.GRBackend()
pyplot()
#plotly()
# metaGrapher
function metaGrapher(mat)
    mg = MetaDiGraph(size(mat,1))
    set_props!.([mg],vertices(mg),[Dict(:marque => false, :parent=>-1, :distSource => Inf)])
    for i in 1:size(mat,1), j in 1:size(mat,2)
        if mat[i,j] != 0
                add_edge!(mg,i,j,Dict(:weight => mat[i,j],:flow => -1))
        end
    end
    return mg
end
# afficheGraph
function afficheGraph(Mat)
    g = metaGrapher(Mat)
    afficheGraph(g)
end

# returns a random graph adjacency matrix with n vertex
function randomGraphAdj(n)
    mat = rand(1:15*n,n,n)
    for i in 1:n
        mat[i,i] = 0
    end
    return mat
end

# pathfinder
function pathFinder(g,v)
   # imagine that 1 is the source and that we have the following parents
    N=size(g)[1]
    #path= Array{Int}(undef,N)
    path= Array{Int}(undef,N)
    global i=1;
    path[1]=v;

     while(path[i] != -1 )
        path[i+1]=get_prop(g,path[i],:parent)
        i+=1;
     end
     #path=[1,2,3,8];
     if(i<=2)
        pathRet=[];
     else
        pathRet=path[i-1:-1:1];
     end
     return pathRet;
end
##
#Test for PathFinder
g = smallgraph(:petersen)
g = metaGrapher(adjacency_matrix(g))
# imagine that 1 is the source and that we have the following parents
set_prop!(g,8,:parent,3)
set_prop!(g,3,:parent,2)
set_prop!(g,2,:parent,1)
# following the parents we have 1-> 2-> 3-> 8, so we should have [1,2,3,8]
@assert pathFinder(g,8) == [1,2,3,8]
# if the vertex is the source, its path should be empty
@assert pathFinder(g,1) == []
## dijkstra
function dijkstra(Mat,s::Int)
    g1=Mat;
    g = metaGrapher(Mat)
    N=size(g)[1]
    V=vertices(g)
    # We initialize the source with a distance of     0 to the source
    set_prop!(g,s,:distSource,0.0)
    #up=(zeros(1,10))
    for cout=1:N
            # Pick the minimum distance vertex from
            # the set of vertices not yet processed.
            # u is always equal to src in first iteration
            u = minDistance((get_prop.([g],1:N,:distSource)),  get_prop.([g],1:N,:marque))
            println("minDist-->",u)
            set_prop!(g,u,:marque,true)
            for v = 1:N
                if (g1[u,v] > 0.0 && get_prop(g,v,:marque)==false && get_prop(g,v,:distSource) > Float64(get_prop(g,u,:distSource) + g1[u,v]))
                    set_prop!(g,v,:distSource,Float64((get_prop(g,u,:distSource) + g1[u,v])))
                    println(Float64(get_prop(g,u,:distSource) + g1[u,v]))
                end
            end
      end
       #println(pathFinder(g,3))
       nei=[];
       resultH=Vector{Vector{Int64}}([])
       for ip=1:N
           zk=-1;
            ip1=-1;
           while zk!=s
              if(zk==-1)
               ip1=ip;
               set_prop!.([g],1:N,:parent,-1)
               nei=LightGraphs.outneighbors(g,ip1)
               println("neigbour1-->",nei)
              end
              val1=findmin(get_prop.([g],nei,:distSource))
              println(get_prop.([g],nei,:distSource))
              println(val1)
              set_prop!(g,ip1,:parent,nei[val1[2]])
              zk=nei[val1[2]]
              ip1=zk;
              nei=LightGraphs.outneighbors(g,zk)
              println("neigbour2-->",nei)
           end
            if(s!=ip)
                Result=pathFinder(g,ip)
                println(Result)
            else
                Result=[];
                println(Result)
            end

            push!(resultH,Result)
            resultH
       end
          print(resultH)
          return resultH
end

function minDistance(dist,sptSet)
    # Initialize min value
    N=size(dist,1)
    min = Inf;
    min_index= undef;
    for v=1:N
        if (sptSet[v] == false && dist[v] < min)
            min = dist[v];
            min_index = v;
        end
    end
    println(min_index)
    return min_index;
end

#Test Dijkistra (1)
# Petersen graph
g = [ 0  1  0  0  1  1  0  0  0  0
     1  0  1  0  0  0  1  0  0  0
     0  1  0  1  0  0  0  1  0  0
     0  0  1  0  1  0  0  0  1  0
     1  0  0  1  0  0  0  0  0  1
     1  0  0  0  0  0  0  1  1  0
     0  1  0  0  0  0  0  0  1  1
     0  0  1  0  0  1  0  0  0  1
     0  0  0  1  0  1  1  0  0  0
     0  0  0  0  1  0  1  1  0  0]
## to visualize this graph do displayGraph (g)
@assert dijkstra(g,3) == [[3, 2, 1],[3, 2], [], [3, 4],
 [3, 4, 5], [3, 8, 6], [3, 2, 7],
 [3, 8], [3, 4, 9], [3, 8, 10]]
 ## Visualize the GraphPlot
x=[2, 8, 6, 16, 10, 9, 4, 8, 15, 14];
y=[4, 6, 4, 8, 8, 6, 10, 6, 12, 7];
 nodelabel=1:10
 n =10
 edgelabel_dict = Dict()
 edgelabel_dict2 = Dict()
 cs1 =ColorScheme(range(colorant"red", colorant"green", length=100))
 edgelabel_mat = Array{String}(undef, n, n)
 for i in 1:n
     for j in 1:n
         edgelabel_mat[i, j] = edgelabel_dict[(i, j)] = string("edge ", i, " to ", j)
         edgelabel_dict2[(i, j)]=cs1[i];
     end
 end
 GraphRecipes.graphplot(g)
#GraphRecipes.graphplot(g,x=x,y=y,nodesize=3,edgecolor = edgelabel_dict2, names=nodelabel,curvature_scalar=0.1,arrow=arrow(:closed, :head, 1, 1), shorten=0.023)
# Assign x,y to Graph(V,E)
ac=1;
x1=x;
y1=y;
dest=1;
 anim= @animate for i1=1:60
     global x,y,dest
     if(mod(i1,10)==0)
         dest+=1
         ac=1
     end
     ss=dijkstra(g,dest);
     global ac
     #ac=i1;
     for i in 1:n
         for j2 in 1:n
             edgelabel_mat[i, j2] = edgelabel_dict[(i, j2)] = string("edge ", i, " to ", j2);
             edgelabel_dict2[(i, j2)]=nothing;
         end
     end
 if(!isempty(ss[ac]))
     #nodecolor[dest]=colorant"red"
     #nodecolor[ss[ac][end]]=colorant"green"
     spz=reverse(ss[ac]);
     global ab=1;
     while((ab+1)<=size(spz)[1])
         print((spz[ab],spz[ab+1]))
         sp1=[spz[ab], spz[ab+1]]
         #edgelabel_dict2[(spz[ab],spz[ab+1])]=:red;
         edgelabel_dict2[(sort(sp1)[1],sort(sp1)[2])]=:red;
         ab+=1;
     end

     membership = [1,1,1,1,1,1,1,1,1,1]
     nodecolor2 = range(colorant"lightyellow", stop=colorant"yellow", length=n)
     nodecolor2[dest]=colorant"red"
     nodecolor2[spz[1]]=colorant"green"

     nodeshape2=[:hexagon]
     nodefillc1 = nodecolor2#[membership]
     nodeshape1 = nodeshape2[membership]

     node_weights = 1:n
    # nodecolor = range(colorant"yellow", stop=colorant"red", length=n)

     x=x.+rand(10)
     y=y.+rand(10)
    # print(x)
     (GraphRecipes.graphplot(g,x=x,y=y,linewidth=1, self_edge_size=0.25,  nodestrokewidth=1.5,
              markerstrokestyle=:solid,
              markerstrokealpha=1,
              markerstrokecolor=:black,markersize=2,nodeshape=nodeshape1,fontsize=10,nodesize=1,edgecolor = edgelabel_dict2,nodecolor=nodefillc1, names=nodelabel,curvature_scalar=0.08,arrow=arrow(:closed, :head,1, 1), shorten=0.003));
              title!("Dijkstra Routing with Animation(JULIA)")
 end
     ac+=1;
end
gif(anim, "anim_fps15.gif", fps = 5)
mp4(anim, "movie.mp4", fps=4)

max_stars_repo_name: amburosesekar/DijkstraRoutingAnimnate

max_issues_repo_name: amburosesekar/DijkstraRoutingAnimnate

max_forks_repo_name: amburosesekar/DijkstraRoutingAnimnate


text:
using ESGIWeatherLogistics2021: ESGIWeatherLogistics2021, Data
using Plots
using Dates
using DataFrames

const obs = Data.get_obs()
const obsdaily = ESGIWeatherLogistics2021.aggregate_daily_obs(obs)
const prices = Data.get_prices()

function plot_year(obsdaily, year; area, crop)
    obsidx = (obsdaily.area .== area) .& (Dates.year.(obsdaily.date) .== year)
    pricesidx = (prices.crop .== crop) .& (Dates.year.(prices.date) .== year)
    plot(obsdaily.date[obsidx], obsdaily.temperature_mean[obsidx])
    plot!(obsdaily.date[obsidx], obsdaily.temperature_minimum[obsidx])
    plot!(obsdaily.date[obsidx], obsdaily.temperature_maximum[obsidx])
    plot!(prices.date[pricesidx], 10 .* prices.price[pricesidx])
end


text:
import Schedulers

let flags = map(lowercase ∘ strip, split(get(ENV, "SCHEDULERS_JL_TEST_ENABLE", ""), ","))
    enable_full_debugging = "full" in flags
    enable_recording = enable_full_debugging || "recording" in flags
    enable_logging = enable_full_debugging || "logging" in flags
    enable_debugging = enable_full_debugging || "debugging" in flags
    if enable_recording
        @info "Enable recording"
        Schedulers.Internal.Debug.enable()
    end
    if enable_logging
        @info "Enable logging"
        Schedulers.Internal.Debug.enable_logging()
    end
    if enable_debugging
        @info "Enable debugging"
        Schedulers.Internal.enable_debug()
    end
end

using TestFunctionRunner
TestFunctionRunner.@run(paths = ["../benchmark/SchedulersBenchmarks"])


text:
"""
    @forward MyType.prop method1, method2, ...

Extend the provided methods by forwarding the property `prop` of `MyType` instances.
This will give, for a given `method`:
```julia
method(x::MyType, args...; kwargs...) = method(x.prop, args...; kwargs...)
```

"""
macro forward(ex, fs)
    T, prop = @match ex begin
        :($T.$prop) => (T, prop)
        _ => error("Invalid expression $ex, expected <Type>.<prop>")
    end

    fs = @match fs begin
        :(($(fs...),)) => fs
        :($mod.$method) => [fs]
        ::Symbol => [fs]
        _ => error("Expected a method or a tuple of methods, got $fs")
    end

    defs = map(fs) do f
        esc(:($f(x::$T, args...; kwargs...) = $f(x.$prop, args...; kwargs...)))
    end

    Expr(:block, defs...)
end


text:
# This file includes code that was formerly a part of Julia.
# License is MIT: LICENSE.md

using ModuleInterfaceTools

using Random

@api test StrBase

# Should test GenericString also, once overthing else is working
const UnicodeStringTypes = (String, UTF8Str, )
    # (String, UTF16Str, UTF32Str, UniStr, UTF8Str)
const ASCIIStringTypes = (String, UTF8Str, ASCIIStr, LatinStr)
    #    (UnicodeStringTypes..., ASCIIStr, LatinStr, UCS2Str)

function cvtchar(T, ch)
    try
        T(ch)
    catch
        Text4Chr(ch)
    end
end

include("basic.jl")

@testset "Invalid sequences" begin include("invalid.jl")  end
@testset "Valid sequences"   begin include("valid.jl")    end
@testset "Bounds Errors"     begin include("bounds.jl")   end
@testset "UTF-16 tests"      begin include("utf16.jl")    end
@testset "UTF-32 tests"      begin include("utf32.jl")    end
@testset "Conversion errors" begin include("convert.jl")  end
@testset "Pointer functions" begin include("pointer.jl")  end
@testset "Search functions"  begin include("search.jl")   end
@testset "SubStrings"        begin include("substr.jl")   end
@testset "Utility functions" begin include("util.jl")     end
@testset "IO functions"      begin include("io.jl")       end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Blosc2"
version = v"2.0.4"

# Collection of sources required to build Blosc2
sources = [
    GitSource("https://github.com/Blosc/c-blosc2.git", "38da659cd5862dccbb0457266fe752cfd6541d4f"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/c-blosc2/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_TESTS=OFF \
    -DBUILD_BENCHMARKS=OFF \
    -DBUILD_EXAMPLES=OFF \
    -DBUILD_STATIC=OFF \
    -DCMAKE_SHARED_LIBRARY_LINK_C_FLAGS="" \
    -DPREFER_EXTERNAL_ZLIB=ON \
    -DPREFER_EXTERNAL_ZSTD=ON \
    -DPREFER_EXTERNAL_LZ4=ON \
    ..
make -j${nproc}
make install
install_license ../LICENSES/*.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libblosc2", :libblosc2),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("Zstd_jll"),
    Dependency("Lz4_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5.2")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#declaration/assignment, declaration is optional
x::Int32 = 1
#datatypes are inferred dynamically, but can also be set thru convert functions and datatype literals
x = 1 #x is inferred as Int, which is Int32 for 32-bit machines, Int64 for 64-bit machines
#variable reference
julia>x
1

x = int8(1) #x is of type Int8
x = 1.0 #x is Float64
x = y = 1 #assign both x and y to 1
global x = 1 #assigns 1 to global variable x (used inside scope blocks)
local x = 1 #assigns 1 to local variable x (used inside scope blocks)
x = 'a' #x is a 'Char' type, designated by single quotes
x = "a" #x is a 1-element string, designated by double quotes

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type HMCMetric end
struct DiagEuclideanMetric <: HMCMetric end
struct UnitEuclideanMetric <: HMCMetric end
struct DenseEuclideanMetric <: HMCMetric end

function AHMCMetric(metric::DiagEuclideanMetric, dim::Integer)
    return AdvancedHMC.DiagEuclideanMetric(dim)
end

function AHMCMetric(metric::UnitEuclideanMetric, dim::Integer)
    return AdvancedHMC.UnitEuclideanMetric(dim)
end

function AHMCMetric(metric::DenseEuclideanMetric, dim::Integer)
    return AdvancedHMC.DenseEuclideanMetric(dim)
end


text:
using .Curve: BN, FP, bn_read_bin!
using .Config: @EP, @EP2, @ID, SIGNATURE_SIZE, PUBLIC_KEY_SIZE, PRIVATE_KEY_SIZE, ORDER
using .Shamir: BarycentricWeightGenerator

abstract type AbstractPrivateKey end
struct PrivateKey <: AbstractPrivateKey
    sk::BN
    PrivateKey(sk) = PrivateKey(BN(sk))
    function PrivateKey(sk::BN)
        safesk = mod(sk, ORDER)
        safesk == sk || @warn "Unsafe private key detected"
        return new(safesk)
    end
    function PrivateKey(bytes::Vector{UInt8})
        length(bytes) >= PRIVATE_KEY_SIZE || @warn "Weak private key used"
        return PrivateKey(bn_read_bin!(BN(), bytes))
    end
end

abstract type AbstractPublicKey end
struct PublicKey <: AbstractPublicKey
    pk::@EP2
    PublicKey(pk::@EP2) = new(pk)
    PublicKey(sk::PrivateKey) = PublicKey(Util.genpk(sk.sk))
    function PublicKey(bytes::Vector{UInt8})
        length(bytes) == PUBLIC_KEY_SIZE || error("bytes array has wrong size: $(length(bytes))")
        return PublicKey(Util.decode(@EP2, bytes))
    end
end

abstract type AbstractIdentity end
struct Identity <: AbstractIdentity
    id::@ID
    Identity(id::@ID) = new((@ID)(id & typemax(@ID)))
    Identity(pk::PublicKey) = Identity((@ID)(pk.pk))
end

abstract type AbstractHash end
struct Hash <: AbstractHash
    hash::@EP
    Hash(hash::@EP) = new(hash)
    Hash(msg::Vector{UInt8}) = Hash(Util.hash2curve(msg))
    Hash(msg::String) = Hash(Vector{UInt8}(msg))
end

abstract type AbstractSignature end
struct Signature <: AbstractSignature
    sig::@EP
    Signature(sig::@EP) = new(sig)
    Signature(sk::PrivateKey, hash::Hash) = Signature(Util.sign(sk.sk, hash.hash))
    function Signature(bytes::Vector{UInt8})
        length(bytes) == SIGNATURE_SIZE || error("bytes array has wrong size: $(length(bytes))")
        return Signature(Util.decode(@EP, bytes))
    end
end

abstract type AbstractPrivateKeyPoly end
struct PrivateKeyPoly <: AbstractPrivateKeyPoly
    coeffs::Vector{BN}
    PrivateKeyPoly(degree::Int) = new([Util.gensk() for i in 0:degree])
end

abstract type AbstractPublicKeyPoly end
struct PublicKeyPoly <: AbstractPublicKeyPoly
    coeffs::Vector{@EP2}
    PublicKeyPoly(coeffs::Vector{@EP2}) = new(coeffs)
    PublicKeyPoly(skpoly::PrivateKeyPoly) = PublicKeyPoly([Util.genpk(sk) for sk in skpoly])
end

abstract type AbstractSignatureShares end
struct SignatureShares <: AbstractSignatureShares
    shares::Dict{@ID,@EP}
    weights::BarycentricWeightGenerator
    SignatureShares() = new(Dict{@ID,@EP}(), BarycentricWeightGenerator())
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# # Binary heap and Dijkstra's algorithm
#
# * 黒木玄
# * 2020-06-26～2020-06-28, 2020-08-23, 2021-10-14
#
# 2020-06-28: [CSV.jl v0.7](https://discourse.julialang.org/t/ann-csv-jl-0-7-release/42162)以上に対応した.
#
# 2020-08-23: Julia v1.6.0-DEVを使うようにし, @inbounds を付けた. ([gist](https://gist.github.com/genkuroki/4b8dd01ccc931abe16246cdd8a0d75af))
#
# 2021-10-14: Julia v1.7.0-rc1 (2021-09-12) を使うようにし, `isless`, `rowvals`, `nzrange` を使うようにした. 警告: `isless` を使うようにしたら, `Float64` 型成分の配列のヒープソートが遅くなってしまった.

# %% [markdown] toc=true
#
目次

#

1  MyHeap型の定義

2  Heap sort

3  Dijkstra's algorithm

3.1  疎行列による重み付きグラフの取り扱い

3.2  Dijkstra's algorithm の実装

4  Tokyo_Edgelist example

4.1  read CSV

4.2  Benchmark

4.3  おまけ: 距離行列の成分をAny型にするとどうなるか

# %% [markdown]
# ## MyHeap型の定義

# %%
"""
    MyHeap{T}(value::T)

The heap structure realized on an array-like object value.
"""
struct MyHeap{T}
    value::T
end

parent_index(k) = k ÷ 2
left_index(k)   = 2k
right_index(k)  = 2k + 1

function heapifyup!(v, k)
    @inbounds while !isone(k)
        p = parent_index(k)
        isless(v[k], v[p]) || break
        v[k], v[p] = v[p], v[k]
        k = p
    end
end

function heapifydown!(v, k)
    N = length(v)
    @inbounds while true
        l = left_index(k)
        l > N && break
        r = right_index(k)
        if r > N
            isless(v[l], v[k]) || break
            v[k], v[l] = v[l], v[k]
            break
        end
        m =  isless(v[l], v[r]) ? l : r
        isless(v[m], v[k]) || break
        v[k], v[m] = v[m], v[k]
        k = m
    end
end

# function MyHeap(v)
#     N = length(v)
#     for k in 1:N
#         heapifyup!(v, k)
#     end
#     MyHeap{typeof(v)}(v)
# end

function MyHeap(v)
    N = length(v)
    for k in N:-1:1
        heapifydown!(v, k)
    end
    MyHeap{typeof(v)}(v)
end

MyHeap{T}() where T = MyHeap(T())
MyHeap() where T = MyHeap{Array{Float64, 1}}()

using Printf

function print_heap(v)
    N = length(v)
    iszero(N) && return
    L = floor(Int, log(2, N))
    for l in 0:L
        for i in 2^l:min(2^(l+1)-1, N)
            print(" "^((6 ÷ 2)*(2^(L-l)-1)))
            @printf("%6s", v[i])
            print(" "^((6 ÷ 2)*(2^(L-l)-1)))
        end
        println()
    end
end

print_heap(H::MyHeap) = print_heap(H.value)

Base.isempty(H::MyHeap) = isempty(H.value)
Base.length(H::MyHeap) = length(H.value)
Base.getindex(H::MyHeap, k) = H.value[k]

function Base.setindex!(H::MyHeap, val, k)
    v = H.value
    @inbounds v[k] = val
    @inbounds if !isone(k) && isless(v[k], v[parent_index(k)])
        heapifyup!(v, k)
        return
    end
    heapifydown!(v, k)
    val
end

"""
    push!(H::MyHeap, val)

Add val to H.
"""
function Base.push!(H::MyHeap, val)
    v = H.value
    N = length(v)
    push!(v, val)
    heapifyup!(v, N+1)
    H
end

"""
    popindex!(H::MyHeap, k)

Get and remove the k-th value of H.
"""
function popindex!(H::MyHeap, k)
    v = H.value
    @inbounds v[k], v[end] = v[end], v[k]
    val = pop!(v)
    heapifydown!(v, k)
    val
end

"""
    popfirst!(H::MyHeap)

Get and remove the first value of H.
"""
Base.popfirst!(H::MyHeap) = popindex!(H, 1)

# %%
# test of MyHeap

@show H = MyHeap(rand(1:100, 27))
print_heap(H)

# %%
# test of setindex!

H = MyHeap(rand(0:99, 31))
print_heap(H)
println()
@show H[31]
println()

v = copy(H.value)
v[31] = -99
print_heap(v)
println()
@show v[31]
println()

H[31] = -99
print_heap(H)

println()
println()
println()
println()

H = MyHeap(rand(0:99, 31))
print_heap(H)
println()
@show H[3]
println()

v = copy(H.value)
v[3] = 999
print_heap(v)
println()
@show v[3]
println()

H[3] = 999
print_heap(H)

# %%
# test of push!

H = MyHeap(rand(0:99, 20))
print_heap(H)
println()

w = copy(H.value)
push!(w, 5)
print_heap(w)
println()

push!(H, 5)
print_heap(H)

# %%
# test of popfirst!

H = MyHeap(rand(0:99, 12))
print_heap(H)
println()

for _ in 1:length(H)
    @show popfirst!(H)
    print_heap(H)
    println()
end

# %% [markdown]
# ## Heap sort

# %% [markdown]
# 上で定義したヒープを用いたソートと `SortingAlgorithms.jl` のヒープ・ソートを比較してみる.

# %%
function my_heap_sort(v)
    H = MyHeap(v)
    w = similar(v)
    @inbounds for i in 1:length(v)
        w[i] = popfirst!(H)
    end
    w
end

using SortingAlgorithms

v = rand(1:100, 10)
@show v
@show sort(v; alg=HeapSort)
@show my_heap_sort(v)
println()

X = rand(1:10^7, 10^7)
@time Y = sort(X; alg=HeapSort)
v = copy(X)
@time w = my_heap_sort(v)
w == Y

# %% [markdown]
# ## Dijkstra's algorithm

# %% [markdown]
# ### 疎行列による重み付きグラフの取り扱い

# %%
using SparseArrays
using LinearAlgebra

"""
    MyWeightedGraph{Tv, Ti}(W::SparseMatrixCSC{S, T})

Weighted graph structure given by the weight matrix W.
Assume that `W` is symmetric and off-diagonal.
"""
struct MyWeightedGraph{Tv, Ti}
    W::SparseMatrixCSC{Tv, Ti}
end

function MyWeightedGraph(W::AbstractMatrix{Tv}) where Tv
    W = sparse(Symmetric(sparse(W)))
    MyWeightedGraph{Tv, Int}(W)
end

const MyWGraph = MyWeightedGraph

weight_matrix(g::MyWeightedGraph) = g.W
nvertices(g::MyWeightedGraph) = size(g.W, 1)
vertices(g::MyWeightedGraph) = 1:nvertices(g)

nzptr(W::SparseArrays.AbstractSparseMatrixCSC, k) = view(rowvals(W), nzrange(W, k))
neighborhood(g::MyWeightedGraph, k) = nzptr(g.W, k)

# %% [markdown]
# 以下のセルで `SimpleWeightedGraphs.jl`, `LightGraphs.jl`, `GraphPlot` を読み込むが, 重み付きグラフのプロットのみに使用され, 最短経路の計算には使用されない.

# %%
using SimpleWeightedGraphs: SimpleWeightedGraphs
using LightGraphs: LightGraphs
using GraphPlot

function weighted_graph_plot(g; kwargs...)
    simplified_g = LightGraphs.SimpleGraph(g)
    nodelabel = 1:LightGraphs.nv(simplified_g)
    edgelabel = [LightGraphs.weights(g)[e.src, e.dst] for e in LightGraphs.edges(simplified_g)]
    gplot(simplified_g; nodelabel=nodelabel, edgelabel=edgelabel, kwargs...)
end

# %%
i = [1,  1, 2, 3, 4, 5, 6, 6,  7, 8]
j = [2,  5, 3, 4, 5, 6, 7, 9,  8, 9]
d = [1, 10, 2, 3, 4, 5, 6, 20, 7, 8]
W = sparse(i, j, d, 10, 10)
W = Symmetric(W)
W = sparse(W)
g = MyWGraph(W)
swg = SimpleWeightedGraphs.WGraph(W)

SimpleWeightedGraphs.weights(swg) |> Matrix |> display
weighted_graph_plot(swg; layout=circular_layout)

# %%
W = weight_matrix(g)
for i in 1:nvertices(g)
    for j in neighborhood(g, i)
        @show i, j, W[i, j]
    end
end

# %% [markdown]
# ### Dijkstra's algorithm の実装

# %% [markdown]
# 位置を表す key と出発点からの距離を表す val の組のヒープを扱うことによって, ヒープ中のデータの中で出発点からの距離が最小の位置が自動的に得られる.

# %%
struct KeyVal{Tv, Ti}
    key::Ti
    val::Tv
end

Base.isless(x::KeyVal, y::KeyVal) = isless(x.val, y.val)

# %%
# test of MyHeap of KeyVal

v = [KeyVal(k, rand(0:0.1:1)) for k in 1:10]
vv = v .|> x -> (x.key, x.val)
@show vv

H = MyHeap(v)
v = H.value

keys = [v[k].key for k in 1:10]
print_heap(keys)
println()
vals = [v[k].val for k in 1:10]
print_heap(vals)

H

# %% [markdown]
# Dijkstra's algorithm の解説については以下を参照せよ.
#
# * https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
# * https://www.youtube.com/results?search_query=dijkstra%27s+algorithm

# %%
"""
    DijkstraResult{S,T}(parent::S, distance::T)

The result of Dijkstra's algorithm searching the shortest path.
"""
struct DijkstraResult{S,T}
    parent::S
    distance::T
end

"""
    dijkstra(g::MyWeightedGraph{Tv, Ti}, start)

returns the result of Dijkstra's algorithm searching the shortest path from `start`.
"""
function dijkstra(g::MyWeightedGraph{Tv, Ti}, start) where {Tv, Ti}
    W = weight_matrix(g)
    N = nvertices(g)

    distance = fill(typemax(Tv), N)
    parent = zeros(Ti, N)
    visited = zeros(Bool, N)

    distance[start] = zero(Tv)
    visited[start] = true
    H = MyHeap([KeyVal(start, distance[start])])

    @inbounds while !isempty(H)
        i = popfirst!(H).key
        d_parent = distance[i]
        for j in neighborhood(g, i)
            d = d_parent + W[i, j]
            d ≥ typemax(Tv) && continue
            if !visited[j]
                visited[j] = true
                distance[j] = d
                parent[j] = i
                push!(H, KeyVal(j, d))
            elseif d < distance[j]
                distance[j] = d
                parent[j] = i
                push!(H, KeyVal(j, d))
            end
        end
    end

    DijkstraResult(parent, distance)
end

"""
    shortest_path(r::DijkstraResult, k)

returns a shortest path from `start` to `k`.
"""
function shortest_path(r::DijkstraResult, k)
    parent = r.parent
    distance = r.distance
    p = eltype(parent)[]
    distance[k] ≥ typemax(eltype(distance)) && return p
    push!(p, k)
    k = parent[k]
    while !iszero(k)
        pushfirst!(p, k)
        k = parent[k]
    end
    p
end

"""
    shortest_distance(r::DijkstraResult, k)

returns the shortest distance from `start` to `k`.
"""
shortest_distance(r::DijkstraResult, k) = r.distance[k]

# %%
i = [1,  1, 2, 3, 4, 5, 6, 6,  7, 8]
j = [2,  5, 3, 4, 5, 6, 7, 9,  8, 9]
d = [1, 10, 2, 3, 4, 5, 6, 20, 7, 8]
W = sparse(i, j, d, 10, 10)
W = Symmetric(W)
W = sparse(W)
g = MyWGraph(W)
swg = SimpleWeightedGraphs.WGraph(W)

weighted_graph_plot(swg; layout=circular_layout)

# %%
r = dijkstra(g, 1)

# %%
[(from=1, to=k, distance=shortest_distance(r, k), path=shortest_path(r, k)) for k in 1:10]

# %% [markdown]
# ## Tokyo_Edgelist example

# %% [markdown]
# ### read CSV

# %%
@time using CSV

# %%
# download("https://ndownloader.figshare.com/files/3663336", "Tokyo.zip")
# run(`unzip -o Tokyo.zip Tokyo_Edgelist.csv`)

# %%
@time csv = CSV.File("Tokyo_Edgelist.csv") # warmup
@show csv.names
@show csv.types
csv.columns

# %% [markdown]
# ### Benchmark

# %%
@time begin
    @time csv = CSV.File("Tokyo_Edgelist.csv")
    @time start_node = csv.START_NODE
    @time end_node = csv.END_NODE
    @time edge_length = csv.LENGTH
    @time D = sparse(start_node, end_node, edge_length)
    @time weighted_graph = MyWGraph(D)
    println()
    print("Dijkstra's algorithm:")
    @time dijkstra_result = dijkstra(weighted_graph, 1)
    println()
    @time dijkstra_path = shortest_path(dijkstra_result, 2539)
    println()
    @show dijkstra_path
    println()
end # warmup

# %%
@time begin
    @time csv = CSV.File("Tokyo_Edgelist.csv")
    @time start_node = csv.START_NODE
    @time end_node = csv.END_NODE
    @time edge_length = csv.LENGTH
    @time D = sparse(start_node, end_node, edge_length)
    @time weighted_graph = MyWGraph(D)
    println()
    print("Dijkstra's algorithm:")
    @time dijkstra_result = dijkstra(weighted_graph, 1)
    println()
    @time dijkstra_path = shortest_path(dijkstra_result, 2539)
    println()
    @show dijkstra_path
    println()
end

# %%
# Comparison between the above and the reki2000_route result.

reki2000_path = vec([1 110273 110421 110498 110574 110658 110706 110746 110865 110998 111102 111186 111371 111680 111808 111820 112054 112134 112217 112452 112503 112530 112546 112774 113010 113100 113710 113830 113926 113999 114010 114092 114308 114458 114704 114804 114969 115015 115091 115303 115352 115645 115667 115733 115771 115859 116007 116227 116606 116605 116815 117103 117447 117610 117607 117860 118052 426 510 567 617 630 666 749 833 939 956 1029 1127 1214 1350 1483 1521 1678 1742 1752 1813 1960 2013 2119 2232 2351 2408 2434 2464 2469 2539])

dijkstra_path == reki2000_path

# %% [markdown]
# See https://qiita.com/reki2000/items/55ef54b96b26d80ad694 for the reki2000_route result.

# %%
[(from=1, to=k, distance=dijkstra_result.distance[k],
        path=shortest_path(dijkstra_result, k)) for k in 1:20]

# %%
D

# %% [markdown]
# ### おまけ: 距離行列の成分をAny型にするとどうなるか

# %%
using BenchmarkTools

Base.typemax(::Type{Any}) = typemax(Float64)
Base.zero(::Type{Any}) = zero(Float64)

D = sparse(start_node, end_node, edge_length)
D_Any = SparseMatrixCSC{Any, Int}(D)
weighted_graph = MyWGraph(D)
weighted_graph_Any = MyWGraph(D_Any)
dijkstra(weighted_graph_Any, 1)

@show typeof(weighted_graph)
@show typeof(weighted_graph_Any)
println()

print("Float64:")
@btime dijkstra(weighted_graph, 1)
print("Any:    ")
@btime dijkstra(weighted_graph_Any, 1)

# %%
D_Any

# %%

max_stars_repo_path: 0022/Binary heap and Dijkstra's algorithm.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public