text:
using DataFrames
using CSV

include("coordinates.jl")

# All relevant filenames
const LESS_FILE = joinpath("data","MELBOURNE_HOUSE_PRICES_LESS.csv")
const LESS_FILE_PROCESSED = joinpath("data", "MELBOURNE_HOUSE_PRICES_LESS_PROCESSED.csv")
const FULL_FILE = joinpath("data", "Melbourne_housing_FULL_processed.csv")
const FULL_FILE_PROCESSED = joinpath("data","Melbourne_housing_FULL_processed_2.csv")

# Set to true if current processed files need to be overwritten
const OVERWRITE_PROCESSED = true

# Set to true if we must overwrite the coordinates
const OVERWRITE_COORDINATES = false

# Gonna put this in a .env later - if this was a real repo, I'd do it earlier
const API_KEY = "AIzaSyAnxxxFGwpJKRaGGAv0De9hNL0OqlHT6B0"

"""
    fix_coordinate_typos(df::AbstractDataFrame)::DataFrame

Returns a copy of a `DataFrame` with fixed typos in Latitude and Longitude columns
"""
fix_coordinate_typos(df::AbstractDataFrame)::DataFrame = if ("Lattitude" in names(df) && "Longtitude" in names(df)) rename(df, Dict(:Lattitude => "Latitude", :Longtitude => "Longitude")) else df end

"""
    uppercasewordfirst(name::AbstractString)::String

Returns copy of a `String` where only the first letter of each word is capitalised.

# Examples
```jldoctest
julia>uppercasewordfirst("fOO BAr")
Foo Bar
```
"""
uppercasewordfirst(name::AbstractString)::String = join(uppercasefirst.(lowercase.(split(name, " "))), " ")

"""
    consolidate_string_case(df::AbstractDataFrame, cols=[:Suburb, :Address, :Regionname, :CouncilArea])::DataFrame

Returns a copy of the given `DataFrame`, where the (`String``) data in each of the given columns is changed as per the `uppercasewordfirst` function.
"""
consolidate_string_case(df::AbstractDataFrame, cols=[:Suburb, :Address, :Regionname, :CouncilArea])::DataFrame = transform(df, [col => ByRow(x -> uppercasewordfirst(x)) => col for col in cols]...)

replace_na_w_missing(df::AbstractDataFrame, cols=[:Distance, :Postcode, :CouncilArea, :Regionname]) = transform(df, [col => ByRow(x -> if (x |> uppercase == "#N/A") missing else x end) => col for col in cols]...)

df_less = DataFrame(CSV.File(LESS_FILE)) |> consolidate_string_case

df_full = DataFrame(CSV.File(FULL_FILE)) |> consolidate_string_case |> fix_coordinate_typos |> replace_na_w_missing

x = combine(groupby(df_less, [:Suburb]), nrow => :count)

y = filter(row -> row.count == 1, x)

@debug y

x_1 = combine(groupby(df_less, [:Regionname]), nrow => :count)
y_1 = filter(row -> row.count == 1, x_1)

@debug y_1

"""
    find_missing_subset(dataframe::AbstractDataFrame, columns=:)::DataFrame

Returns a copy of the given `DataFrame` that has `missing` records in any of the passed columns (all columns by default)
"""
find_missing_subset(dataframe::AbstractDataFrame, columns=:)::DataFrame = dataframe[.!completecases(dataframe, columns), :]

"""
    list_missing_data(dataframe::AbstractDataFrame, columns=:)

If `JULIA_DEBUG`, prints some information about missing data in the given columns. Otherwise, it does nothing.
"""
function list_missing_data(dataframe::AbstractDataFrame, columns=:)
    info = []
    for col in propertynames(dataframe[!, columns])
        n_rows = nrow(find_missing_subset(dataframe, col))
        if n_rows > 0
            push!(info, "$(col) has $(n_rows) missing rows.")
        end
    end
    return info
end

@info "Missing LESS data: \n\t$(join(list_missing_data(df_less), "\n\t"))"

@info "Missing FULL data: \n\t$(join(list_missing_data(df_full), "\n\t"))"

@debug "Missing latitudes/longitudes:" find_missing_subset(df_full, [:Latitude, :Longitude])

"""
    get_block_of_missing_coords(df::AbstractDataFrame, to_edit::Bool, block_size=100)::DataFrame

Specialised function. Returns a `DataFrame` that contains the last `block_size` rows of `df` that are missing their `Latitude` and `Longitude` if `to_edit` is `true`; otherwise, returns a `DataFrame` that contains all records with missing coordinates save for the last `block_size` rows.
"""
function get_block_of_missing_coords(df::AbstractDataFrame, to_edit::Bool, block_size=100)::DataFrame
    missing_df = df[.!completecases(df, [:Latitude, :Longitude]), :]
    return missing_df[vcat(max(nrow(missing_df) - block_size, 0) |> trues, min(block_size, nrow(missing_df)) |> falses) .⊻ to_edit, :]
end

# Handle file renaming, creating

OVERWRITE_PROCESSED && @warn "Overwriting existing files"

no_file_message(type::String) = "Could not find $(type) file. Creating new."

!isfile(FULL_FILE_PROCESSED) && @warn no_file_message("FULL")
(OVERWRITE_PROCESSED || !isfile(FULL_FILE_PROCESSED)) && CSV.write(FULL_FILE_PROCESSED, df_full)

!isfile(LESS_FILE_PROCESSED) && @warn no_file_message("LESS")
(OVERWRITE_PROCESSED || !isfile(LESS_FILE_PROCESSED)) && CSV.write(LESS_FILE_PROCESSED, df_less)

# Loop through our missing coordinates and fill them in chunk-by-chunk. Writes to the file and reads from it once each loop.
OVERWRITE_COORDINATES && @warn "About to make a whole lotta API calls, watch out"
while OVERWRITE_COORDINATES
    # Read a fresh DataFrame from the file (slow, but saves on API calls in case of error)
    df_filled = DataFrame(CSV.File(FULL_FILE_PROCESSED))

    # Populate the missing set we won't touch, and the one we will
    missing_coords_yesedit = get_block_of_missing_coords(df_filled, true)
    missing_coords_noedit = get_block_of_missing_coords(df_filled, false)

    if nrow(missing_coords_yesedit) == 0
        @info "No more coordinates to fill."
        break
    end

    @info "Editing $(nrow(missing_coords_yesedit)) rows."
    @info "Leaving behind $(nrow(missing_coords_noedit)) rows."

    # Populate the missing coordinates
    transform!(missing_coords_yesedit, [:Suburb, :Address] => ByRow((subr, addr) -> get_coordinates(addr, subr, API_KEY)) => [:Latitude, :Longitude])

    len_filled = nrow(df_filled)

    # Re-arrange and fill our file-read DataFrame with our new data
    filled_df = append!(append!(dropmissing(df_filled, [:Latitude, :Longitude]) |> allowmissing!,  missing_coords_yesedit), missing_coords_noedit)

    row_diff = len_filled - nrow(filled_df)

    @info "Difference in rows: $(row_diff)"
    if row_diff != 0
        # Likely that something's gone wrong, will need a human eye
        @error "Output `DataFrame` differs from input `DataFrame`. Exiting."
        break
    end

    # Overwrite the file
    CSV.write(FULL_FILE_PROCESSED, filled_df)
end


text:
module MortIntl

using Colors, CSV, DataFrames, JSON, PGFPlotsX
srcpath = @__DIR__
mainpath = normpath(@__DIR__, "..")
datapath = normpath(mainpath, "data")
conf = JSON.parsefile(normpath(datapath, "mortintl.json"))
sexes = ["m", "f"]
agest = [0; 0:4; 5:5:95]
ageend = ["ω"; 0:4; 9:5:94; "ω"]

function hmddf_ctrysex(ctry, sex)
        iso = conf["countries"]["$ctry"]["iso3166"]
        hmddf(normpath(datapath, "$(iso)_$(sexes[sex])ltper_1x1.txt"))
end

function hmddf(path)
        of = tempname()
        run(pipeline(`sed 's/110+/110/' $path`, stdout = of))
        CSV.File(of, delim=' ', header=3, ignorerepeated=true) |> DataFrame
end

save_ctry_caprop(cc, path) = cc[:propframe] |>
CSV.write(normpath(path, fname(cc[:ctry], cc[:ca1], cc[:ca2])))

fname(ctry, ca1, ca2) = "$(ca1)$(ca2)$(ctry).csv"

function ctry_caprop(ctry, ca1, ca2, loadpath, savepath)
        ca1d = conf["causes"][ca1]
        ca2d = conf["causes"][ca2]
        loadfullpath = normpath(loadpath, fname(ctry, ca1, ca2))
        if (loadpath != "" && ispath(loadfullpath))
                df = CSV.File(loadfullpath, types = Dict(:li => String)) |> DataFrame
        else
                df = caprop(ctry, ca1d["causeexpr"], ca2d["causeexpr"])
        end
        cc = Dict(:ca1 => ca1, :ca2 => ca2, :ctry => ctry, :propframe => df)
        if savepath != ""
                save_ctry_caprop(cc, savepath)
        end
        cc
end

function caprop(ctry, ca1es, ca2es)
        df = DataFrame()
        for (li, ca1e) in ca1es
                of = tempname()
                le = conf["listexpr"][li]
                ca2e = ca2es[li]
                run(pipeline(Cmd(`$srcpath/propyrs_ctry.sh $ctry "$ca1e" "$ca2e" "$le"`, dir = srcpath),
                             stdout = of))
                dfli = CSV.File(of) |> DataFrame
                df = vcat(df, dfli)
        end
        sort(df, [:yr, :age])
end

function caprop_eplot(ctries, sex, ca1, ca2, caage, eage, lang, loadpath, savepath)
        ca1lab = conf["causes"][ca1]["alias"][lang]
        ca2lab = conf["causes"][ca2]["alias"][lang]
        sexlab = conf["sexes"]["$sex"]["alias"][lang]
        dlab = conf["deaths"]["alias"][lang]
        elab = uppercasefirst(conf["expectancy"]["alias"][lang])
        if lang == "sv"
                pgfnf = @pgf {"/pgf/number format/use comma"}
        else
                pgfnf = @pgf {}
        end
        p = @pgf Axis({xlabel = "e($(eage))",
                      ylabel = "$dlab $ca1lab/$ca2lab $(agest[caage])–$(ageend[caage])",
                      "yticklabel style"=pgfnf,
                      "legend style"={"font=\\tiny"},
                      legend_pos="outer north east",
                      title = "$elab vs $dlab $sexlab", xmajorgrids, ymajorgrids})
        plotcolors = distinguishable_colors(length(ctries)+1, [RGB(1,1,1)])[2:end]
        for (i, ctry) in enumerate(ctries)
                caprop_eplot_ctry(p, plotcolors[i], ctry, sex, ca1, ca2,
                                  caage, eage, loadpath, savepath)
        end
        p
end

function licmp(li1, li2)
        if li1 == li2
                return "{}"
        else
                return li1
        end
end

function caprop_eplot_ctry(p, pcol, ctry, sex, ca1, ca2, caage, eage, loadpath, savepath)
        propf = ctry_caprop(ctry, ca1, ca2, loadpath, savepath)[:propframe]
        propf_sex_caage = propf[(propf[!, :sex].==sex) .& (propf[!, :age].==caage), :]
        ef = hmddf_ctrysex(ctry, sex)
        ef_age = ef[ef[!, :Age].==eage, :]
        pef = innerjoin(ef_age, propf_sex_caage, on = [:Year=>:yr])
        liprev = vcat("", pef[!, :li][1:length(pef[!, :li])-1])
        listart = map(licmp, pef[!, :li], liprev)
        iso = conf["countries"]["$(ctry)"]["iso3166"]
        syr = pef[!, :Year][1]
        eyr = pef[!, :Year][end]
        @pgf push!(p, PlotInc({"mark=+", "nodes near coords",
                               "point meta=explicit symbolic", color = pcol},
                              Table({meta = "meta"}, ["x" => pef[!, :ex], "y" => pef[!, :rat],
                                     "meta" => listart])),
                   LegendEntry("$(iso) $(syr)–$(eyr)"))
end

end

ext: jl


text:
using Images,FileIO,JSON

function Load(imgPath)
    img = load(imgPath)
    data = Dict()
    for p in img
        c = convert(RGB,p)
        data[c] = c
    end
    return collect(values(data))
end

function SaveJson(fileName,dat)
    f = open(fileName,"w")
    json = JSON.json(dat)
    write(f,json)
    close(f)
end

function CvtColor(img,palette)
    img_c = img
    for index = 1:length(img)
        c = img_c[index]
        src = convert(RGB,c)
        minDst = 256*3
        minClr = RGB()
        for p in palette
            dis = abs(src-p)
            if dis<minDst
                minDst = dis
                minClr = p
            end
        end
        img_c[index] = RGBA(minClr,alpha(c))
    end
    return img_c
end

println("start")

dat = Load("colored.png")
#SaveJson("t.json",dat)

img = load("Sample.png")
img = CvtColor(img,dat)
save(File(format"PNG","t.png"),img)

println("end")

max_stars_repo_name: fenixnix/GoAscii

max_issues_repo_name: fenixnix/GoAscii

max_forks_repo_name: fenixnix/GoAscii


text:
module UtilityModels
    using Parameters, Distributions
    import Distributions: mean, var, std
    export UtilityModel, TAX, ProspectTheory, ExpectedUtility, Gamble
    export mean, var , std
    abstract type UtilityModel end

    include("Gamble.jl")
    include("UtilityModel.jl")
    include("ProspectTheory.jl")
    include("ExpectedUtility.jl")
    include("TAX.jl")
end


text:
import Base: getindex, setindex!

function dynamic_eq(def)
  name, val = def.args
  if isexpr(name, :(::))
    name, T = name.args
    :(const $(esc(name)) = Ref{$(esc(T))}($(esc(val))))
  else
    :(const $(esc(name)) = Ref($(esc(val))))
  end
end

function dynamic_let(ex)
  # bindings = [:(bind($(esc(b.args[1])), $(esc(b.args[2])), t)) for b in block(ex.args[1]).args]
  bs = [(esc(b.args[1]), esc(b.args[2])) for b in block(ex.args[1]).args]
  xs = [gensym() for n = 1:length(bs)]
  save = [:($(esc(x)) = $b[]) for (x, (b, v)) in zip(xs, bs)]
  set = [:($b[] = $v) for (x, (b, v)) in zip(xs, bs)]
  unset = [:($b[] = $(esc(x))) for (x, (b, v)) in zip(xs, bs)]
  :(let $(xs...)
    try
      $(save...)
      $(set...)
      $(esc(ex.args[2]))
    finally
      $(unset...)
    end
  end)
end

macro dynamic(def)
  isexpr(def, :(=)) ? dynamic_eq(def) :
  isexpr(def, :let) ? dynamic_let(def) :
  error("Unsupported @dynamic expression")
end


text:
"""
    The number of edges in each "gadget" per original edge
"""
const GADGET_EDGE_MULTIPLIER = 5
const GADGET_VTX_MULTIPLIER = 2

"""
    It is necessary to convert the env graph to an undirected graph because the
    gadget is based on undirected edges. Self-edges also need to be removed, as
    these are already accounted for in the gadget.
"""
function convert_env_graph_to_undirected(G)
    env_graph = Graph(G)
    # for v in vertices(env_graph)
    #     rem_edge!(env_graph,v,v)
    # end
    # return env_graph
end

"""
    GadgetGraph

Represents a time-extended graph useful for MILP formulations

"""
@with_kw_noshow struct GadgetGraph <: AbstractGraph{Int}
    G::DiGraph          = DiGraph()
    vtxs::Vector{Int}   = Vector{Int}()
    tvec::Vector{Int}   = Vector{Int}()
end
Base.zero(graph::GadgetGraph) = GadgetGraph()
LightGraphs.edges(graph::GadgetGraph) = edges(graph.G)
LightGraphs.is_directed(graph::GadgetGraph) = true
function LightGraphs.add_vertex!(graph::GadgetGraph,t=-1,v=-1)
    if add_vertex!(graph.G)
        push!(graph.vtxs,v)
        push!(graph.tvec,t)
        return true
    end
    return false
end
for op in [
    :edgetype,:has_edge,:has_vertex,:inneighbors,:ne,:nv,:outneighbors,
    :vertices,:add_edge!
    ]
    @eval LightGraphs.$op(g::GadgetGraph,args...) = $op(g.G,args...)
end
get_vtx_from_var(G::GadgetGraph,v) = get(G.vtxs,v,-1)
get_t_from_var(G::GadgetGraph,v) = get(G.t_vec,v,-1)

"""
    add_movement_vtx!

Flags a vertex of the gadget graph as corresponding to a non-"wait" edge
"""
add_movement_vtx!(graph::GadgetGraph,t) = add_vertex!(graph,t,0)
is_movement_vtx(graph::GadgetGraph,v) = get_vtx_from_var(graph,v) == 0

"""
    Return a source map such that source_map[v][t] points to the corresponding
    vertex in the gadget graph.
"""
function get_source_map(G::GadgetGraph,env_graph,T_MAX)
    # source_map = zeros(Int,nv(env_graph),T_MAX)
    source_map = [Dict{Int,Int}() for v in 1:nv(env_graph)]
    for (idx,v,t) in zip(vertices(G),G.vtxs,G.tvec)
        if has_vertex(env_graph,v)
            source_map[v][t] = idx
        end
    end
    return source_map
end

function add_vertex_layer!(G::GadgetGraph,incoming::Dict{Int,Vector{Int}},t)
    outgoing = Dict{Int,Int}()
    for v in sort(collect(keys(incoming)))
        add_vertex!(G,t,v)
        outgoing[v] = nv(G)
        for v2 in incoming[v]
            add_edge!(G,v2,nv(G))
        end
    end
    return outgoing
end
function add_gadget_layer!(G::GadgetGraph,edge_list,t,vtx_map=Dict{Int,Int}())
    outgoing = Dict{Int,Vector{Int}}()
    for e in edge_list
        if !haskey(vtx_map,e.src) && !haskey(vtx_map,e.dst)
            continue
        end
        if e.src == e.dst
            # shortcut
            add_vertex!(G,t)
            add_edge!(G,vtx_map[e.src],nv(G))
            push!(get!(outgoing,e.dst,Int[]),nv(G))
        else
            if !haskey(vtx_map,e.src)
                e = reverse(e)
            end
            # Gadget
            add_movement_vtx!(G,t)
            v1 = nv(G)
            add_edge!(G,vtx_map[e.src],v1)
            if haskey(vtx_map,e.dst)
                add_movement_vtx!(G,t)
                v2 = nv(G)
                add_edge!(G,vtx_map[e.dst],v2)
            end
            add_movement_vtx!(G,t)
            v3 = nv(G)
            add_edge!(G,v1,v3)
            if haskey(vtx_map,e.dst)
                add_edge!(G,v2,v3)
                add_movement_vtx!(G,t)
                v4 = nv(G)
                add_edge!(G,v3,v4)
                push!(get!(outgoing,e.src,Int[]),v4)
            end
            add_movement_vtx!(G,t)
            v5 = nv(G)
            add_edge!(G,v3,v5)
            push!(get!(outgoing,e.dst,Int[]),v5)
        end
    end
    return outgoing
end

function construct_gadget_graph(env_graph,T_MAX,t0=0,cap=true)
    G = GadgetGraph()
    incoming = Dict(v=>Int[] for v in vertices(env_graph))
    edge_list = collect(edges(env_graph))
    for t in t0:T_MAX-1
        outgoing = add_vertex_layer!(G,incoming,t)
        incoming = add_gadget_layer!(G,edge_list,t,outgoing)
    end
    if cap
        add_vertex_layer!(G,incoming,T_MAX)
    end
    return G
end

function add_point_constraint!(model,source_map,flow,v,t=0)
    idx = source_map[v][t]
    @constraint(model, flow[idx] == 1)
    return model
end
function add_continuity_constraints!(model,G,flow)
    for v in vertices(G)
        @constraint(model,1 >= sum(map(vp->flow[vp],outneighbors(G,v))))
        @constraint(model,1 >= sum(map(vp->flow[vp],inneighbors(G,v))))
        if outdegree(G,v) > 0
            @constraint(model,
                flow[v] <= sum(map(vp->flow[vp],outneighbors(G,v)))
                )
        end
        if indegree(G,v) > 0
            @constraint(model,
                flow[v] <= sum(map(vp->flow[vp],inneighbors(G,v)))
                )
        end
    end
    return model
end
function add_carrying_constraints!(model,G,robot_flow,object_flow)
    for v in vertices(G)
        if is_movement_vtx(G,v)
            # object may not traverse an edge unless it is being carried
            @constraint(model, object_flow[v] - robot_flow[v] <= 0)
        end
    end
    return model
end
function add_precedence_constraints!(model,source_map,inflow,outflow,goal,start,Δt=0)
    for t in 0:size(source_map,2)-Δt
        vtx1 = source_map[goal][t]
        vtx2 = source_map[start][t+1+Δt]
        # if flow1 is not at vtx1, flow2 must still be at vtx2
        @constraint(model,outflow[vtx2] + inflow[vtx1] >= 1)
    end
end
function add_makespan_constraint!(model,T,flow,source_map,goal)
    for (t,v) in source_map[goal]
        @constraint(model,T >= t*(1-flow[v]))
    end
    return model
end
function add_total_flow_constraint!(model,source_map,flow,n,t=0)
    @constraint(model,sum(map(d->flow[d[t]],source_map)) == n)
end

struct PCTAPF_MILP
    model::JuMP.Model
    G::GadgetGraph
    source_map::Vector{Dict{Int,Int}}
    robot_flow::Vector{VariableRef}
    object_flows::Dict{ObjectID,Vector{VariableRef}}
end
JuMP.optimize!(m::PCTAPF_MILP) = JuMP.optimize!(m.model)
extract_robot_flow(m::PCTAPF_MILP) = Int.(round.(value.(m.robot_flow)))
extract_object_flows(m::PCTAPF_MILP) = Dict(id=>Int.(round.(value.(var))) for (id,var) in m.object_flows)
function extract_flow_path(m::PCTAPF_MILP,flow,v0,t0)
    vtx = m.source_map[v0][0]
    path = Int[vtx]
    while outdegree(m.G,vtx) > 0
        for vp in outneighbors(m.G,vtx)
            if flow[vp] == 1
                vtx = vp
                break
            end
        end
        push!(path,vtx)
    end
    return path
end
function extract_true_path(m::PCTAPF_MILP,flow_path)
    path = Int[]
    for vtx in flow_path
        if get_vtx_from_var(m.G,vtx) > 0
            push!(path,get_vtx_from_var(m.G,vtx))
        end
    end
    return path
end
function extract_robot_paths(prob::PC_TAPF,m::PCTAPF_MILP)
    sched = get_schedule(get_env(prob))
    robot_ICs = get_robot_ICs(sched)
    flow = Int.(round.(value.(m.robot_flow)))
    paths = Dict{RobotID,Vector{Int}}()
    for (id,pred) in robot_ICs
        v0 = get_id(get_initial_location_id(pred))
        t0 = get_t0(get_env(prob),id)
        flow_path = extract_flow_path(m,flow,v0,t0)
        paths[id] = extract_true_path(m,flow_path)
    end
    return paths
end
function extract_object_paths(prob::PC_TAPF,m::PCTAPF_MILP)
    sched = get_schedule(get_env(prob))
    object_ICs = get_object_ICs(sched)
    paths = Dict{ObjectID,Vector{Int}}()
    for (id,pred) in object_ICs
        flow = Int.(round.(value.(m.object_flows[id])))
        v0 = get_id(get_initial_location_id(pred))
        t0 = get_t0(get_env(prob),id)
        flow_path = extract_flow_path(m,flow,v0,t0)
        paths[id] = extract_true_path(m,flow_path)
    end
    return paths
end

"""
    formulate_big_milp

Formulate a PCTAPF problem as a giant network flow MILP.
"""
function formulate_big_milp(prob::PC_TAPF,T_MAX,model = JuMP.Model())
    sched = get_schedule(get_env(prob))
    robot_ICs = get_robot_ICs(sched)
    object_ICs = get_object_ICs(sched)
    env_graph = convert_env_graph_to_undirected(get_graph(get_env(prob)).graph)
    G = construct_gadget_graph(env_graph,T_MAX) # for robot flow
    source_map = get_source_map(G,env_graph,T_MAX)
    # Robot flows
    @variable(model,robot_flow[1:nv(G)], binary=true)
    # Object flows
    object_flows = Dict(id=>@variable(model,[1:nv(G)], binary=true) for (id,pred) in object_ICs)
    # Flow constraints
    add_continuity_constraints!(model,G,robot_flow)
    # initial constraints
    for (id,pred) in robot_ICs
        v0 = get_id(get_initial_location_id(pred))
        t0 = get_t0(get_env(prob),id)
        add_point_constraint!(model,source_map,robot_flow,v0,t0)
    end
    add_total_flow_constraint!(model,source_map,robot_flow,length(robot_ICs))
    # makespan
    @variable(model,T >= 0)
    for (id,pred) in object_ICs
        object_flow = object_flows[id]
        add_total_flow_constraint!(model,source_map,object_flow,1)
        # initial location
        v0 = get_id(get_initial_location_id(pred))
        t0 = get_t0(get_env(prob),id)
        add_point_constraint!(model,source_map,object_flow,v0,t0)
        add_continuity_constraints!(model,G,object_flow)
        add_carrying_constraints!(model,G,robot_flow,object_flow)
    end
    # precedence constraints
    for op in values(get_operations(sched))
        Δt = duration(op)
        for (id1,o2) in postconditions(op)
            outflow = object_flows[id2]
            start = get_id(get_initial_location_id(o2))
            @assert get_id(get_destination_location_id(object_ICs[id2])) == start
            for (id1,o1) in preconditions(op)
                # makespans
                # id1 = get_object_id(o1)
                inflow = object_flows[id1]
                goal = get_id(get_destination_location_id(o1))
                # @assert get_id(get_destination_location_id(get_object_FCs(sched)[id1])) == goal
                add_precedence_constraints!(model,source_map,inflow,outflow,goal,start,Δt)
            end
        end
    end
    # Objective
    for op in values(get_operations(sched))
        for (id,o) in preconditions(op)
            # makespans
            # id = get_object_id(o)
            object_flow = object_flows[id]
            goal = get_id(get_destination_location_id(o))
            @show o, goal
            add_point_constraint!(model,source_map,object_flow,goal,T_MAX)
            add_makespan_constraint!(model,T,object_flow,source_map,goal)
        end
    end
    @objective(model,Min,T)
    return PCTAPF_MILP(model, G, source_map, robot_flow, object_flows)
end


text:
#__precompile__(false)

module Fastnet

using StaticArrays,Random,DataFrames,Formatting,UnicodePlots

export FastNet,FastSim,LinkType,
    adjacent,
    node_f, node,           # Node functions
    nodestate_f, nodestate,
    nodestate_f!, nodestate!,
    makenode_f!, makenode!,
    makenodes_f!, makenodes!,
    destroynode_f!, destroynode!,
    countnodes_f, countnodes,
    nodecounts_f, nodecounts,
    firstlinkin_f, firstlinkin,
    firstlinkout_f, firstlinkout,
    degree_f, degree,
    indegree_f, indegree,
    outdegree_f, outdegree,
    nodeexists_f,nodeexists,
    randomnode_f,randomnode,

    link_f, link,                   # link functions
    linkstate_f, linkstate,
    makelink_f!, makelink!,
    destroylink_f!, destroylink!,
    linksrc_f, linksrc,
    linkdst_f, linkdst,
    nextlinkout_f, nextlinkout,
    nextlinkin_f, nextlinkin,
    countlinks_f, countlinks,
    linkcounts_f, linkcounts,
    countlinksinstate_f, countlinksinstate,
    linkexists_f,linkexists,
    randomlink_f,randomlink,

    nullgraph!,             # networktools
    randomgraph!,
    emptygraph!,
    configmodel!,
    regulargraph!,
    rectlattice!,
    adjacency!,

    showlinks,shownodes,     # Debug Functions
    healthcheck,

    listnodes,               # Analaysis functions
    listlinks,
    listnodesinstate,
    listneighbors,
    listnodestates,
    results,
    savelinklist,
    savenodeinfo,
    degreedist,
    quickplotnodes,
    quickplotlinks,

    runsim!,                  # Simulation
    simstep!

include("helpers.jl")
include("structure.jl")
include("internal.jl")
include("access.jl")
include("node.jl")
include("link.jl")
include("safenode.jl")
include("safelink.jl")
include("simulation.jl")
include("nettools.jl")
include("debug.jl")
include("prettyprint.jl")

end # module


text:
module DirectTrajectoryOptimization

using LinearAlgebra
using SparseArrays
using Symbolics, IfElse
using Parameters
using Ipopt
using MathOptInterface
const MOI = MathOptInterface

include("costs.jl")
include("constraints.jl")
include("bounds.jl")
include("general_constraint.jl")
include("dynamics.jl")
include("options.jl")
include("data.jl")
include("solver.jl")
include("moi.jl")
include("utils.jl")

# objective
export Cost

# constraints
export Bound, Bounds, Constraint, Constraints, GeneralConstraint

# dynamics
export Dynamics

# solver
export Solver, Options, initialize_states!, initialize_controls!, solve!, get_trajectory

# utils
export linear_interpolation

end # module


text:
function calc_Vandermonde_dot_dot_dot_dot(x,ni,nj)

    V=zeros(BigFloat, ni, nj)
    for i=1:ni
        for j=5:nj
            V[i,j] = (BigInt(j)-BigInt(1))*(BigInt(j)-BigInt(2)*(BigInt(j)-BigInt(3))*(BigInt(j)-BigInt(4)))*x[i]^(BigInt(j)-BigInt(5))
        end
    end

    return V

end


text:
# Fractional Laplacians - fourier
# Author: Ziyi Yin, ziyi.yin@gatech.edu
# Date: July 2020

export laplacian_mask, apply_mask_fourier, laplacian_operator

###############################################################################
# Fractional laplacian implemented by 2D Discrete Fourier Transform

function laplacian_mask(n::Tuple, d::Tuple, order::Number)

	dk = 1/d[1]/n[1]
	mask = ones(Float32,n)
	center_x = (n[1]+1)/2
	center_z = (n[2]+1)/2

	x = reshape(collect(1:n[1]), :, 1)
	z = reshape(collect(1:n[2]), 1, :)

	f = dk*sqrt.((x.-center_x).^2f0.+(z.-center_z).^2f0)

	if order >= 0	# differentiation
		mask = (2*pi*f).^order
	else # integration - warning: null space of low frequency
		#delta = dk/10f0
		#mask = (2*pi*(f.+delta)).^order	# damped by delta
		center_x = Int(round(center_x))
		center_z = Int(round(center_z))
		mask[[1:center_x-1;center_x+1:end],:] = (2*pi*f[[1:center_x-1;center_x+1:end],:]).^order
		mask[center_x,[1:center_z-1;center_z+1:end]] = (2*pi*f[center_x,[1:center_z-1;center_z+1:end]]).^order
		mask[center_x,center_z] = 0.25*(mask[center_x+1,center_z]+mask[center_x-1,center_z]+mask[center_x,center_z+1]+mask[center_x,center_z-1])
	end

	return convert(Array{Float32,2},mask)
end

function apply_mask_fourier(q::Array,mask::Array,nsrc::Int)
	n = size(mask)
	q_t = reshape(q,n[1],n[2],nsrc)
	q_ans = deepcopy(q_t)
	for i = 1:nsrc
		q_f = fftshift(fft(q_t[:,:,i]))
		q_new = mask.*q_f
		q_ans[:,:,i] = real.(ifft(ifftshift(q_new)))
	end
	return reshape(q_ans,size(q))
end

function apply_mask_fourier(q::PhysicalParameter,mask::Array,nsrc::Int)
	return apply_mask_fourier(vec(q.data),mask,nsrc)
end

function apply_mask_fourier(q::judiWeights,mask::Array,nsrc::Int)
	n = size(mask)
	q_ans = deepcopy(q)
	for i = 1:nsrc
		q_t = q.weights[i]
		q_f = fftshift(fft(q_t))
		q_new = mask.*q_f
		q_ans.weights[i] = real.(ifft(ifftshift(q_new)))
	end
	return q_ans
end

function laplacian_operator(model::Model,order::Number,info::Info)
	n = model.n
	#if order >= 0
	if true
		P = joLinearFunctionFwd_T(prod(n), prod(n),
	                             v -> apply_mask_fourier(v,laplacian_mask(model,order),model,info),
	                             w -> apply_mask_fourier(w,laplacian_mask(model,order),model,info),
								 Float32,Float32,name="Fractional laplacian operator")
	else
		P1 = joLinearFunctionFwd_T(prod(n), prod(n),
								 v -> apply_mask_fourier(v,laplacian_mask(model,-order),model,info),
								 w -> apply_mask_fourier(w,laplacian_mask(model,-order),model,info),
								 Float32,Float32,name="Fractional laplacian operator")

		P = joLinearFunctionFwd_T(prod(n), prod(n),
								 v -> lsqr(P1,v,damp=1,verbose=true),
								 w -> lsqr(P1,w,damp=1,verbose=true),
								 Float32,Float32,name="Fractional laplacian operator")
	end

	return P
end


text:
# using Test
# using BioXP
using JSON

@testset "read master from file" begin
    # reactions = readmaster("../../ecg/test/userdata/kegg/master.json")
    reactions = readmaster("data/master-subset-3rxns.json")
    expected_reactions = ["R10514", "R09735", "R08827"]
    @test typeof(reactions) == Reactions
    @test Set(keys(reactions)) == Set(expected_reactions)

    @testset "correct reaction data" begin
        expected_left = [
            ["C00024","C20638"],
            ["C19848","C06232"],
            ["C00355","C08538"]]
        expected_right = [
            ["C00010","C20672"],
            ["C18237","C00020","C00001"],
            ["C08543","C00001"]]
        for (i,r) in enumerate(expected_reactions)
            @test reactions[r].left == expected_left[i]
            @test reactions[r].right == expected_right[i]
        end
    end
end

@testset "read master from open json" begin
    master = JSON.parsefile("data/master-subset-3rxns.json")
    @test typeof(master) <: Dict

    reactions = readmaster(master)
    expected_reactions = ["R10514", "R09735", "R08827"]
    @test typeof(reactions) == Reactions
    @test Set(keys(reactions)) == Set(expected_reactions)

    @testset "correct reaction data" begin
        expected_left = [
            ["C00024","C20638"],
            ["C19848","C06232"],
            ["C00355","C08538"]]
        expected_right = [
            ["C00010","C20672"],
            ["C18237","C00020","C00001"],
            ["C08543","C00001"]]
        for (i,r) in enumerate(expected_reactions)
            @test reactions[r].left == expected_left[i]
            @test reactions[r].right == expected_right[i]
        end
    end
end

@testset "read compounds from dir" begin
    compounds = readcompounds("data/compound")

    @testset "spot check compound properties C00379" begin
        c = compounds["C00379"]
        @test c.entry_id == "C00379"
        @test c.name == "Xylitol"
        @test c.formula == "C5H12O5"
        @test c.exact_mass == 152.0685
        @test c.mol_weight == 152.1458
        @test Set(c.reactions) == Set([
            "R01431",
            "R01896",
            "R01904",
            "R02136",
            "R05831",
            "R07152",
            "R09477"])

        expected_pathways = Dict("map00040"=> "Pentose and glucuronate interconversions",
        "map01100"=> "Metabolic pathways",
        "map02010"=> "ABC transporters")
        @test Set(keys(expected_pathways)) == Set(keys(c.pathways))
        for k in keys(expected_pathways)
            @test c.pathways[k] == expected_pathways[k]
        end

        expected_enzymes = [
            "1.1.1.9",
            "1.1.1.10",
            "1.1.1.14",
            "1.1.1.15",
            "1.1.1.21",
            "1.1.1.307",
            "1.1.3.41",
            "2.7.1.122"]
        @test Set(keys(expected_enzymes)) == Set(keys(c.enzymes))
    end

    @testset "spot check compound properties C00380" begin
        c = compounds["C00380"]
        @test c.entry_id == "C00380"
        @test c.name == "Cytosine"
        @test c.formula == "C4H5N3O"
        @test c.exact_mass == 111.0433
        @test c.mol_weight == 111.102
        @test Set(c.reactions) == Set([
            "R00510",
            "R00974",
            "R02137",
            "R02296",
            "R10837"])

        expected_pathways = Dict("map00240"=>  "Pyrimidine metabolism",
        "map01100"=> "Metabolic pathways")
        @test Set(keys(expected_pathways)) == Set(keys(c.pathways))
        for k in keys(expected_pathways)
            @test c.pathways[k] == expected_pathways[k]
        end

        expected_enzymes = [
            "2.4.2.2",
            "2.4.2.57",
            "3.2.2.8",
            "3.2.2.10",
            "3.5.4.1"]
        @test Set(keys(expected_enzymes)) == Set(keys(c.enzymes))
    end

end
@testset "read reaction IDs" begin
    @testset "from dir" begin
        reactions = readids("data/taxon_reactions/1234567892.json")
        @test Set(reactions) == Set([
            "R01773",
            "R01775",
            "R02003",
            "R05556"])

    end

    @testset "from vector" begin
        reactions = readids([
            "R01773",
            "R01775",
            "R02003",
            "R05556"])

        @test Set(reactions) == Set([
            "R01773",
            "R01775",
            "R02003",
            "R05556"])

    end
end

@testset "read seed and target IDs" begin
    @testset "from dir" begin
        seeds = readkeyedids("data/seeds/seeds.json")

        @test Set(seeds["Goldford-SAFR-032"]) == Set(["C00001","C00011","C00288","C00283","C00014","C00697","C00058","C00033"])
        @test Set(seeds["Enceladus_20-SAFR-032"]) == Set(["C00001","C00011","C00237","C00282","C00067","C00132","C06548","C00469","C00283","C00014","C00697","C01326","C01438","C01548","C06547","C11505","C20783","C01407"])
        @test Set(seeds["Enceladus_20-SAFR-032_P"]) == Set(["C00001","C00011","C00237","C00282","C00067","C00132","C06548","C00469","C00283","C00014","C00697","C01326","C01438","C01548","C06547","C11505","C20783","C01407","C00009"])
        @test Set(seeds["targets_Freilich09"]) == Set(["C00001","C00002","C00003","C00004","C00005","C00006","C00008","C00015","C00016","C00020","C00024","C00025","C00035","C00037","C00041","C00043","C00044","C00047","C00049","C00054","C00055","C00062","C00063","C00064","C00065","C00073","C00075","C00078","C00079","C00082","C00097","C00105","C00112","C00116","C00123","C00131","C00135","C00144","C00148","C00152","C00183","C00188","C00234","C00239","C00249","C00255","C00286","C00350","C00360","C00362","C00364","C00399","C00407","C00458","C00459","C00641","C00748","C01050","C05764","C05890","C05894","C05899","C05980","C06040","C15672","C16221"])
    end

    @testset "from dict" begin
        seeds = readkeyedids(Dict("Goldford-SAFR-032"=>["C00001","C00011","C00288","C00283","C00014","C00697","C00058","C00033"],
                 "Enceladus_20-SAFR-032"=>["C00001","C00011","C00237","C00282","C00067","C00132","C06548","C00469","C00283","C00014","C00697","C01326","C01438","C01548","C06547","C11505","C20783","C01407"]))

        @test Set(seeds["Goldford-SAFR-032"]) == Set(["C00001","C00011","C00288","C00283","C00014","C00697","C00058","C00033"])
        @test Set(seeds["Enceladus_20-SAFR-032"]) == Set(["C00001","C00011","C00237","C00282","C00067","C00132","C06548","C00469","C00283","C00014","C00697","C01326","C01438","C01548","C06547","C11505","C20783","C01407"])

    end

    @testset "from vector" begin
        seeds = readids(["C00001","C00011","C00288","C00283","C00014","C00697","C00058","C00033"])

        @test Set(seeds) == Set(["C00001","C00011","C00288","C00283","C00014","C00697","C00058","C00033"])
    end
end

@testset "symbolize keys" begin
    testdicts = [Dict("key1"=>1,"key2"=>2), Dict("key1"=>"1","key2"=>"2")]

    for d in testdicts
        s =  BioXP.symbolizekeys(d)
        for k in keys(s)
            @test typeof(k) == Symbol
        end
    end

end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
### DIC ###

# Batch DIC
function deviance(param::T, loglike::Function) where T
  return -2 * loglike(param)
end

function DIC(params::Vector{T}, loglike::Function, paramMean::Function) where T
  D(param) = deviance(param, loglike)

  Dbar = mean(D.(params))
  paramBar = paramMean(params)
  pD = Dbar - D(paramBar)

  return Dbar + pD
end

# Streaming DIC
"""
Methods to implement:
- updateParams(d::DICstream{T}, param::P) where {T, P}, to add param to d.paramSum
- paramMeanCompute(d::DICstream{T}) where T, to compute the mean of the parameters
- loglike(param::P) computes the loglikelihood based on the parameters (param)
- convert(state::State) convets the current state to param::Param
"""
mutable struct DICstream{T}
  paramSum::T
  Dsum::Float64
  counter::Int

  DICstream{T}(tmp) where T = new(tmp, 0, 0)
end

"""
convert: converts state to param
"""
function updateDIC(x::DICstream{T}, state::S, updateParams::Function,
                   loglike::Function, convert::Function) where {T, S}
  @assert x.counter >= 0

  param = convert(state)

  # update Dsum
  x.Dsum += deviance(param, loglike)

  # update paramSum
  if x.counter == 0
    x.paramSum = deepcopy(param)
  else
    updateParams(x, param)
  end

  # update counter
  x.counter += 1
end

"""
DIC is a measure of fit. Higher is better.

DIC = Dmean + pD, where
- Dmean is a measure of fit. Higher is better.
- pD is a measure of (non) complexity. Higher is better.
"""
function computeDIC(x::DICstream{T}, loglike::Function,
                    paramMeanCompute::Function;
                    return_Dmean_pD::Bool=false) where {T}
  Dmean = x.Dsum / x.counter
  paramMean = paramMeanCompute(x)
  pD = Dmean - deviance(paramMean, loglike)

  if return_Dmean_pD
    return Dmean, pD
  else
    return Dmean + pD
  end
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
######################################################################
# trajectory_manager.jl
#
# handles the trajectory manager
######################################################################

using StatsBase: Weights, sample

abstract type Dynamics end
abstract type Initializer end
abstract type Descender end
abstract type IntegrationScheme end

type TrajectoryManager

	# needed for all trajectories
	N::Int
	h::Float64
	x0::Vector{Float64}

	# Cost functions
	q::Float64
	R::Matrix{Float64}
	Qn::Matrix{Float64}
	Rn::Matrix{Float64}
	barrier_cost::Float64

	initializer::Initializer
	descender::Descender
	dynamics::Dynamics
	int_scheme::IntegrationScheme

	function TrajectoryManager(x0::Vector{Float64}, h::Real, N::Int, i::Initializer=RandomInitializer())
		tm = new()

		# needed for all trajectories
		tm.N = N
		tm.h = h
		tm.x0 = deepcopy(x0)

		# needed for ergodic trajectories
		tm.Qn = eye(2)
		tm.q = 1.0
		tm.R = 0.01 * eye(2)
		tm.Rn = eye(2)
		tm.barrier_cost = 0.
		tm.initializer = i
		tm.descender = ArmijoLineSearch()

		# dynamics stuff
		tm.dynamics = LinearDynamics(eye(2), tm.h*eye(2))
		tm.int_scheme = ForwardEuler()

		return tm
	end
end

#typealias VTM Vector{TrajectoryManager}
const VTM = Vector{TrajectoryManager}

# computes controls from a trajectory
# TODO: really, this is a general tool useful for other code
#  it should go somewhere else
export compute_controls
function compute_controls(xd::VVF, h::Float64)
	N = length(xd) - 1
	ud = Array(Vector{Float64}, N)
	for n = 1:N
		ud[n] = (xd[n+1] - xd[n]) / h
	end

	return ud
end

# creates a sample_trajectory
function sample_trajectory(em::ErgodicManager, tm::TrajectoryManager)
	return initialize(SampleInitializer(), em, tm)
end


text:
using PtFEM
using Test

# write your own tests here

examples = [
  "PtFEM/4 Static Equilibrium/Ex41.1.jl",
  "PtFEM/4 Static Equilibrium/Ex41.2.jl",
  "PtFEM/4 Static Equilibrium/Ex42.2.jl",
  "PtFEM/4 Static Equilibrium/Ex43.1b.jl",
  "PtFEM/4 Static Equilibrium/Ex44.2.jl",
  "PtFEM/4 Static Equilibrium/Ex45.2.jl",
  "PtFEM/4 Static Equilibrium/Ex46.2.jl",
  "PtFEM/4 Static Equilibrium/Ex47.1.jl",
  "PtFEM/5 Elastic Solids/Ex51.1.jl",
  "PtFEM/5 Elastic Solids/Ex51.2.jl",
  "PtFEM/5 Elastic Solids/Ex51.3.jl",
  "PtFEM/5 Elastic Solids/Ex51.4.jl",
  "PtFEM/5 Elastic Solids/Ex51.5.jl",
  "PtFEM/5 Elastic Solids/Ex52.1.jl",
  "PtFEM/5 Elastic Solids/Ex53.1.jl",
  "PtFEM/5 Elastic Solids/Ex54.1.jl",
  "PtFEM/5 Elastic Solids/Ex54.2.jl",
  "PtFEM/5 Elastic Solids/Ex55.1.jl",
  "PtFEM/5 Elastic Solids/Ex56.1.jl",
  #"PtFEM/5 Elastic Solids/Ex57.1.jl",        # Takes a while ...
  #"PtFEM/5 Elastic Solids/Ex57.1b.jl",       # Requires too much memory
  "PtFEM/6 Material Nonlinearity/Ex61.1.jl",
  "PtFEM/6 Material Nonlinearity/Ex62.1.jl",
  "Nodal Formulas/figure-23.jl",
  "Nodal Formulas/figure-23b.jl",
  "Nodal Formulas/figure-24.jl",
]

println("\nRunning examples:\n")

if VERSION.minor < 7
  println("\nRunning Julia 0.5 examples\n")
  for ex in examples
    println("\n\nRunning $(ex).\n\n")
    include(joinpath(dirname(@__FILE__), "..", "examples", ex))
  end
end


text:
import Base: sum, prod, maximum, minimum, sum!, prod!, maximum!, minimum!
for (fname, _fname, fname!, mpiop) in [(:sum, :_sum, :sum!, :(MPI.SUM)), (:prod, :_prod, :prod!, :(MPI.PROD)), (:maximum,:_maximum, :maximum!, :(MPI.MAX)), (:minimum, :_minimum, :minimum!, :(MPI.MIN))]

    @eval begin
        function ($_fname)(f::Function, a::MPIArray{T,N,A}, ::Colon) where {T,N,A}
            partial = ($fname)(f, a.localarray)
            sync()
            MPI.Allreduce(partial, $mpiop, MPI.COMM_WORLD)
        end

        function ($fname!)(r::MPIArray{T,N,A}, arr::MPIArray{T,N,A}) where {T,N,A}
            for (dr, darr) in zip(size(r), size(arr))
                @assert dr == darr || dr == 1
            end
            @assert size(r)[end] == size(arr)[end]
            ($fname!)(r.localarray, arr.localarray)
            r
        end

        function ($fname!)(r::AbstractArray{T,N}, arr::MPIArray{T,N,A}) where {T,N,A}
            @assert size(r)[end] == 1
            ($fname!)(r, arr.localarray)
            sync()
            MPI.Allreduce!(r, $mpiop, MPI.COMM_WORLD)
            r
        end

        function ($_fname)(a::MPIArray{T,N,A}, dims::NTuple{N2, Integer}) where {T,N,A,N2}
            outsize = map(x -> x in dims ? 1 : size(a, x), 1:N)
            if N in dims
                out = A{T}(undef, outsize...)
            else
                out = similar(a, outsize...)
            end
            ($fname!)(out, a)
            out
        end

        ($_fname)(f::Function, a::MPIArray{T,N,A}, dims::Integer) where {T,N,A} = ($_fname)(f, a, (dims,))
        ($fname)(f::Function, a::MPIArray{T,N,A}; dims=:) where {T,N,A} = ($_fname)(f, a, dims)
        ($_fname)(a::MPIArray{T,N,A}, ::Colon) where {T,N,A} = ($_fname)(identity, a, :)
        ($_fname)(a::MPIArray{T,N,A}, dims::Integer) where {T,N,A} = ($_fname)(a, (dims,))
        ($fname)(a::MPIArray{T,N,A}; dims=:) where {T,N,A} = ($_fname)(a, dims)
    end
end


text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "make full adjacency matrix" begin
    try
        i2f = Dict([(1, "#ab"), (2, "abc"), (3, "bc#"), (4, "#bc"), (5, "ab#")])
        JudiLing.make_full_adjacency_matrix(i2f)
        @test true
    catch e
        @test e == false
    end
end

@testset "make combined adjacency matrix" begin
    try
        latin_full = DataFrame(CSV.File(joinpath(
            @__DIR__,
            "data",
            "latin_mini.csv",
        )))

        latin_train = latin_full[1:3, :]
        latin_val = latin_full[10:15, :]

        A = JudiLing.make_combined_adjacency_matrix(
            latin_train,
            latin_val,
            grams = 3,
            target_col = :Word,
            tokenized = false,
            keep_sep = false,
        )

        @test A[1, 2] == 1
        @test A[2, 3] == 1
        @test A[3, 4] == 1
        @test A[4, 5] == 1
        @test A[2, 6] == 1
        @test A[6, 7] == 1
    catch e
        @test false
    end
end


text:
function bbsession(bbscripts::Vector{T};
    livepull::Bool = error("livepull must be true or false")) where T<:Val
  S = startbbsession()
  try
    for v ∈ bbscripts
      bbscript(S,v, livepull=livepull)
    end
  finally
    endbbsession(S) #always close the session when done
  end
end

#helper methods to make the above easier to call
bbsession(;kws...) = bbsession(Vector{Val}(); kws...)
bbsession(v::Val;kws...) = bbsession([v]; kws...)
bbsession(bbscripts::Vector{Symbol};kws...) = bbsession([Val{s}() for s ∈ bbscripts]; kws...)
bbsession(s::Symbol;kws...) = bbsession([s]; kws...)

startbbsession() = BLPData.Session(client_mode=BLPData.BLPAPI_CLIENTMODE_SAPI)
endbbsession(S) = BLPData.stop(S)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#EasyPlot core type, constant & function definitions
#-------------------------------------------------------------------------------

#==Types & constants
===============================================================================#
#Dispatchable symbol (minimize namespace pollution):
struct DS{T}; end
DS(v::Symbol) = DS{v}()

"`Optional{T} = Union{Nothing, T}`"
const Optional{T} = Union{Nothing, T}

struct Default; end
const default = Default()

struct NoOverwrite; end #Alt?: NoChange, KeepValue
const nooverwrite = NoOverwrite()

#==Accessors
===============================================================================#
"""`NoOverwrite(v)`

Test whether a value is of type "NoOverwrite"."""
NoOverwrite(v) = false
NoOverwrite(v::NoOverwrite) = true

#==Helper functions
===============================================================================#
SI(v; ndigits=3) = NumericIO.formatted(v, :SI, ndigits=ndigits)

#==Documentation
===============================================================================#
@doc """`Default`

Empty type whose singleton instance `default` represents a default value.
""" Default

@doc """`default`

Singleton instance of type `Default` representing a default value.
""" default

@doc """`NoOverwrite`

Empty type whose singleton instance `nooverwrite` indicates an attribute should
not be overwritten.
""" NoOverwrite

@doc """`nooverwrite`

Singleton instance of type `NoOverwrite`.
Identifies an attribute that should not be overwritten.
""" nooverwrite

#Last line


text:
using Pkg
Pkg.add("Latexify")
using Latexify

Pkg.add("Plots")
using Plots
Pkg.add("PyPlot")
pyplot()

# Julia code for linear Stommel problem solution in "F. X. GIRALDO AND M. RESTELLI"
# High-order semi-implicit time-integrators for a triangular discontinuous Galerkin
# oceanic shallow water model

L=1e6 # Basin size
τ₀=0.2
β=1e-11
f₀=1e-4
γ=1e-6
g=10
H=1000
ρ=1000

# Derived terms
y_m=L/2
ϕ_B=g*H

ETMP01=:( ((β/γ)^2+4*(π/L)^2)^0.5 )
λ₁_def=:( λ₁ = (-β/γ+$ETMP01)/2 )
λ₂_def=:( λ₂ = (-β/γ-$ETMP01)/2 )
display( latexify(λ₁_def) )
display( latexify(λ₂_def) )

C₀_def=:( C₀ = τ₀/(γ*ρ*H)*(L/π) )
C₁_def=:( C₁ =  C₀*( (1-exp(λ₂*L))  / (exp(λ₂*L) - exp(λ₁*L)) ) )
C₂_def=:( C₂ = -C₀*( (1-exp(λ₁*L))  / (exp(λ₂*L) - exp(λ₁*L)) ) )
display( latexify(C₀_def) )
display( latexify(C₁_def) )
display( latexify(C₂_def) )

ψ_def=:( ψ = (C₀ + C₁*exp(λ₁*x) + C₂*exp(λ₂*x))*sin((π*y)/L) )
display( latexify(ψ_def) )

display( λ₁_def )
display( λ₂_def )
display( C₀_def )
display( C₁_def )
display( C₂_def )
display( ψ_def  )

λ₁ = (-β / γ + ((β / γ) ^ 2 + 4 * (π / L) ^ 2) ^ 0.5) / 2
λ₂ = (-β / γ - ((β / γ) ^ 2 + 4 * (π / L) ^ 2) ^ 0.5) / 2
C₀ = (τ₀ / (γ * ρ * H)) * (L / π)
C₁ = C₀ * ((1 - exp(λ₂ * L)) / (exp(λ₂ * L) - exp(λ₁ * L)))
C₂ = -C₀ * ((1 - exp(λ₁ * L)) / (exp(λ₂ * L) - exp(λ₁ * L)))
ψ(x,y) = (C₀ + C₁ * exp(λ₁ * x) + C₂ * exp(λ₂ * x)) * sin((π * y) / L)

nx=ny=50;
x=collect(range(0.,L,step=L/nx))
y=collect(range(0.,L,step=L/ny))
ψval=zeros(nx+1,ny+1)
for j=1:ny
 for i=1:nx
  ψval[i,j]=ψ(x[i],y[j])
 end
end
ψval

p1 = contour(x, y, ψval')
plot(p1)

ETMP01=:( C₀*β*L/π*cos(π*y/L) +fofy*ψ )
ETMP02=:( C₁/λ₁*exp(λ₁*x)+C₂/λ₂*exp(λ₂*x) )
display( latexify(ETMP01) )
display( latexify(ETMP02) )
ϕ_s_def=:( ϕ_s = $ETMP01 + γ*π/L*cos(π*y/L)*$ETMP02 )
display( latexify(ϕ_s_def) )

ϕ_s_func(x,y,fofy)=(((C₀ * β * L) / π) * cos((π * y) / L) + fofy * ψval[x,y]) + ((γ * π) / L) * cos((π * y) / L) * ((C₁ / λ₁) * exp(λ₁ * x) + (C₂ / λ₂) * exp(λ₂ * x))

ϕ_s_val=zeros(nx+1,ny+1)
for j=1:ny+1
 for i=1:nx+1
  ϕ_s_val[i,j]=ϕ_s_func( i,j,f₀+β*y[j] )
 end
end

p1 = contour(x, y, ϕ_s_val')
plot(p1)

max_stars_repo_name: christophernhill/climahacks

max_issues_repo_name: christophernhill/climahacks

max_forks_repo_name: christophernhill/climahacks


text:
module Points

export Point

using ..Elements
import ..Motions: induce_velocity, mutually_induce_velocity!, self_induce_velocity!, advect

#== Type definition ==#

"""
    Point{T} <: Elements.Element

An immutable structure representing a point source/vortex

## Fields

- `z::ComplexF64`: position
- `S::T`: strength/circulation
"""
struct Point{T <: Number} <: Element
    z::ComplexF64
    S::T
    Point{T}(z, s::Real) where T <: Complex = new(z, im*s)
    Point{T}(z, s)       where T <: Complex = new(z, s)
    Point{T}(z, s::Real) where T <: Real = new(z, s)
end
Elements.kind(::Point) = Singleton
Elements.kind(::Type{Point{T}}) where T = Singleton

#== Methods to be extended ==#

Elements.position(p::Point) = p.z

Elements.streamfunction(z::ComplexF64, p::Point) = real(-0.5p.S*log(z - p.z)/π)

Elements.complexpotential(z::ComplexF64, p::Point) = -0.5im*p.S*log(z - p.z)/π

cauchy_kernel(z) = z != zero(z) ? 0.5im/(π*conj(z)) : zero(z)

function induce_velocity(z::ComplexF64, p::Point, t)
    p.S'*cauchy_kernel(z - p.z)
end

function mutually_induce_velocity!(ws₁, ws₂,
                                   points₁::Vector{Point{T₁}},
                                   points₂::Vector{Point{T₂}}, t) where {T₁, T₂}
    for (s, source) in enumerate(points₁)
        for (t, target) in enumerate(points₂)
            K = cauchy_kernel(target.z - source.z)
            ws₂[t] += source.S'*K
            ws₁[s] -= target.S'*K
        end
    end
    nothing
end

function self_induce_velocity!(ws, points::Vector{Point{T}}, t) where T
    N = length(points)

    for s in 1:N, t in s+1:N
        K = cauchy_kernel(points[t].z - points[s].z)
        ws[t] += points[s].S'*K
        ws[s] -= points[t].S'*K
    end
    ws
end

function advect(p::Point{T}, w::ComplexF64, Δt::Float64) where T
    Point{T}(p.z + w*Δt, p.S)
end

#function Base.show(io::IO, p::Point)
#    print(io, "Vortex.Point(z = $(round(p.z, digits=3)), Γ = $(round(p.Γ, digits=3)))")
#end

end


text:
#!/usr/bin/env julia

include("PassiveCable.jl")

using JSON
using SIUnits.ShortUnits
using PassiveCable

# This should run the same effective model
# as rallpack1, but with differing
# electrical parameters (see below).

function run_cable(x_prop, ts)
    # Physical properties:

    # f is a fudge factor. rM needs to be the same
    # the same as in arbor, where we cannot yet set
    # the membrane conductance parameter. Scaling
    # other parameters proportionally, however,
    # gives the same dynamics.

    f = 0.1/4

    diam =   1.0µm        # cable diameter
    L    =   1.0mm        # cable length
    I    =   0.1nA     /f # current injection
    rL   =   1.0Ω*m    *f # bulk resistivity
    erev = -65.0mV        # (passive) reversal potential
    rM   =   4Ω*m^2    *f # membrane resistivity
    cM   =   0.01F/m^2 /f # membrane specific capacitance

    # convert to linear resistivity, length and time constants
    area = pi*diam^2/4
    r    = rL/area

    lambda = sqrt(diam/4 * rM/rL)
    tau = cM*rM

    # compute solutions
    tol = 1e-8mV
    return [cable(L*x_prop, t, L, lambda, tau, r, erev, -I, tol=tol) for t in ts]
end

function run_rallpack1(x_prop, ts)
    return [rallpack1(0.001*x_prop, t/s)*V for t in ts]
end

# Generate traces at x=0, x=0.3L, x=L

ts = collect(0s: 0.025ms: 250ms)
trace = Dict(
    :name => "membrane voltage",
    :sim => "numeric",
    :model => "rallpack1",
    :units => "mV",
    :data => Dict(
        :time => map(t->t/ms, ts),
        Symbol("cable.x0.0") => map(v->v/mV, run_cable(0, ts)),
        Symbol("cable.x0.3") => map(v->v/mV, run_cable(0.3, ts)),
        Symbol("cable.x1.0") => map(v->v/mV, run_cable(1.0, ts))
    )
)

println(JSON.json([trace]))

max_stars_repo_name: kabicm/arbor

max_issues_repo_name: kabicm/arbor

max_forks_repo_name: kabicm/arbor


text:
using CategoricalArrays, DataFrames
include("experiment_model.jl")

function summarize_data(datapath)
    data = DataFrame!(CSV.File(datapath))
    recode_raw_survey_data!(data)
    println(describe(data, :mean, :std; cols=["Q1", "Q2", "Q3", "Q4"]))
    println(length(data.Q1))
end

function recode_raw_survey_data!(df)
    recode!(df.Q1, 1=>5.0, 2=> 15.0, 3=> 25.0, 4=> 35.0, 5=> 45.0, 6=>55.0,
                    7=> 65.0, 8=>75.0, 9=>85.0, 10=> 95.0, 11=> 125.0, 12=> 150.0)
    recode!(df.Q2, 1=> 21.0, 2=> 29.0, 3=> 39.0, 4=> 49.0, 5=> 59.0, 6=> 69.0, 7=> 79.0, 8=> 85.0)
    recode!(df.Q3, 1=> 10.0, 2=> 12.0, 3=> 13.0, 4=> 14.0, 5=> 16.0, 6=> 18.0, 7=> 20.0)
    recode!(df.Q4, 1=> 1.0, 2=>0.0)
end

function recode_perturb!(df, perturb)
    recode!(df.perturb1, 2=> 1.0, 1=> -1.0)
    df.perturb1 = df.perturb1.*perturb[1]
    recode!(df.perturb2, 2=> 1.0, 1=> -1.0)
    df.perturb2 = df.perturb2.*perturb[2]
    recode!(df.perturb3, 2=> 1.0, 1=> -1.0)
    df.perturb3 = df.perturb3.*perturb[3]
end

function check_bonus_payments(df, beta)
    dscale=1
    perturb = [0, beta[2:4]...]./3
    xmeans = [1, mean(df.Q2), mean(df.Q3), mean(df.Q4)]
    println("perturbs: ", perturb)
    println("mean prediction: ", sum(beta.*xmeans)*dscale)
    println("mean prediction at max perturb: ", sum((beta.+perturb).*xmeans)*dscale)
    println("mean prediction at min perturb: ", sum((beta.-perturb).*xmeans)*dscale)
    println("maximum possible prediction: ", sum((beta.+ perturb).*[1, 85, 20, 1])*dscale)
    println("minimum possible prediction: ", sum((beta.-perturb).*[1, 21, 10, 0])*dscale)
    println("minimum no perturb: ", sum((beta).*[1, 21, 10, 0])*dscale)
    println("maximum no perturb: ", sum(beta.*[1, 85, 20, 1])*dscale)
end

function run_ols(datapath)
    data = DataFrame!(CSV.File(datapath))
    recode_raw_survey_data!(data)
    ols = lm(@formula(Q1 ~ Q2 +Q3 +Q4), data)
    println(ols)
    return coef(ols), data
end

function run_gradient(datapath, β, α = [0.5, 1/80000.0, 1/150.0, 1/3.0])
    perturb = β[2:length(β)]./3
    data = DataFrame!(CSV.File(datapath))
    recode_raw_survey_data!(data)
    recode_perturb!(data, perturb)
    data_in = ExperimentData(data.Q1, Matrix(data[!, ["Q2", "Q3", "Q4"]]), Matrix(data[!, ["perturb1", "perturb2", "perturb3"]]))
    β_new = update_β_robust(data_in, β, α)
    return β_new
end

function o_mse_calc(files, beta)
    mses = zeros(length(files))
    n  = zeros(length(files))
    for (i, f) in zip(1:length(files), files)
        data = DataFrame!(CSV.File(f))
        recode_raw_survey_data!(data)
        x = Matrix([ones(size(data)[1]) data[!, ["Q2", "Q3", "Q4"]]])
        y = data.Q1
        yhat = x*beta
        mses[i] =  1 - mean((yhat .- y).^2)/var(y)
        n[i] = length(yhat)
    end
    println(mses)
    return sum(mses.*n)/sum(n)
end

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
function readarray(fileName)
    f=open(fileName);
    arr=Float64[];
    a = readline(f);    # Get the first line
    while a != ""
        push!(arr,parsefloat(Float64,a));
        a=readline(f);  # Get subsequent lines
    end
    return arr;
end


text: include("./notebook.jl")


text:
####################################
# Torneo binario
###################################

"""
    TorneoBinario(generación, tamaño)
Devuelve `tamaño` índices  ganadores de cada enfrentamiento.
"""
function TorneoBinario(evaluación_generación, tamaño)
    índices_ganadores = zeros(Int, tamaño)
    for n in 1:tamaño
        # Índices a enfrentar
        i = rand(1:tamaño)
        j = rand(1:tamaño)
        while j==i
            j = rand(1:tamaño)
        end
        # guardamos el índice del mejor
        índices_ganadores[n] = evaluación_generación[i] < evaluación_generación[j] ? j : i
    end

    return índices_ganadores
end

"""
    TorneoBinario(generación, tamaño)
Toma dos elementos aleatorios y los enfrenta,
devuelve el mejor candidato
"""
function TorneoBinarioEntreDos(evaluación_generación, tamaño)
    # Índices a enfrentar
    i = rand(1:tamaño)
    j = rand(1:tamaño)
    while j==i
        j = rand(1:tamaño)
    end
    # guardamos el índice del mejor
    índice_ganador = evaluación_generación[i] < evaluación_generación[j] ? j : i

    return índice_ganador
end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ac65c4de-58f3-11eb-00c2-f9f83a4c01a2
using PlutoUI

# ╔═╡ 5580aee0-6165-11eb-2967-ff278167c8d5
md"""
#### Random hiring
Captain Bayes: **To become a member of my crew please show how well you can produce random numbers. Proof your randomness!**

👉 Enter **200 digits** ´0-9´ in the textbox below!

$$G(x,y)=\frac{1}{2\pi \sigma^2}e^{\frac{-(x^2+y^2)}{2\sigma^2}}$$
"""

# ╔═╡ be49b1d0-58f3-11eb-1770-23439de5d5b9
@bind random_string TextField((40, 5); default="")

# ╔═╡ 57cfdfc0-6149-11eb-19bd-5f941bb1d8a5
@bind test Slider(0:10)

# ╔═╡ 2e10cb10-61a6-11eb-2f51-6b30066a00e6
HTML("

")

# ╔═╡ 4a826c10-6163-11eb-0077-fb171c3a7379
begin
	function arra(inp)
	try
			if length(inp) == 0
				[]
			else
				parse.(Int,split(inp, ""))
			end
	catch
		0
	end
	end
end

# ╔═╡ 6650f3d2-6163-11eb-27e8-d7db2c942e54
arr = arra(random_string)

# ╔═╡ e76f7270-58f3-11eb-37b2-fd410fcf6082

# ╔═╡ 697e0370-58f5-11eb-3da4-510b6c3d0118
length(arr) + test

# ╔═╡ dba94ef0-58f5-11eb-1857-cf90de0c4e14
prob = rand()

# ╔═╡ d37c5710-6163-11eb-36c6-dde1fbcc2099

if arr == 0

	md"""
!!! danger "Wrong characters!"
    **Please enter digits**: 0, 1, 2, 3, 4, 5, 6, 7 ,8 ,9
	"""

	else
		if isempty(arr)
	Markdown.MD(Markdown.Admonition("warning", "Hiring task!", [md"Please enter 200 digits below in the textbox."]))
elseif length(arr) == 200
	#!!! correct "Completed!"

	#Captain Bayes asserts your random series stems from a good flat random generator with a probability of **$(prob)** %
	Markdown.MD(Markdown.Admonition("correct", "Completed!", [md"Captain Bayes asserts your random series stems from a good flat random generator with a probability of **$(prob)** %"]))
else
		Markdown.MD(Markdown.Admonition("warning", "Hiring task!", [md"Please enter $(200-length(arr)) digits below in the textbox."]))
end

end

# ╔═╡ Cell order:
# ╟─ac65c4de-58f3-11eb-00c2-f9f83a4c01a2
# ╟─5580aee0-6165-11eb-2967-ff278167c8d5
# ╟─d37c5710-6163-11eb-36c6-dde1fbcc2099
# ╟─be49b1d0-58f3-11eb-1770-23439de5d5b9
# ╟─57cfdfc0-6149-11eb-19bd-5f941bb1d8a5
# ╟─2e10cb10-61a6-11eb-2f51-6b30066a00e6
# ╟─4a826c10-6163-11eb-0077-fb171c3a7379
# ╟─6650f3d2-6163-11eb-27e8-d7db2c942e54
# ╟─e76f7270-58f3-11eb-37b2-fd410fcf6082
# ╟─697e0370-58f5-11eb-3da4-510b6c3d0118
# ╟─dba94ef0-58f5-11eb-1857-cf90de0c4e14

max_stars_repo_name: fonsp/bayes-imoox

max_issues_repo_name: fonsp/bayes-imoox

max_forks_repo_name: fonsp/bayes-imoox


text:
using MPI
using ClimateMachine
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods.FVReconstructions: FVConstant, FVLinear
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using LinearAlgebra
using Printf
using Test
import ClimateMachine.VTK: writevtk, writepvtu
import ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using Dates

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output = parse(Bool, lowercase(get(ENV, "JULIA_CLIMA_OUTPUT", "false")))

include("advection_diffusion_model.jl")

struct Pseudo1D{ns, α, β, μ, δ} <: AdvectionDiffusionProblem end

function init_velocity_diffusion!(
    ::Pseudo1D{ns, α, β},
    aux::Vars,
    geom::LocalGeometry,
) where {ns, α, β}
    # Direction of flow is ns[i] with magnitude α
    aux.advection.u = hcat(ntuple(i -> α * ns[i], Val(length(ns)))...)

    # Diffusion of strength β in the ns[i] direction
    aux.diffusion.D = hcat(ntuple(i -> β * ns[i] * ns[i]', Val(length(ns)))...)
end

function gaussian(x, t, α, β, μ, δ)
    exp(-(x - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
end
function ∇gaussian(n, x, t, α, β, μ, δ)
    -2n * (x - μ - α * t) / (4 * β * (δ + t)) *
    exp(-(x - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
end
function initial_condition!(
    ::Pseudo1D{ns, α, β, μ, δ},
    state,
    aux,
    localgeo,
    t,
) where {ns, α, β, μ, δ}
    ρ = ntuple(Val(length(ns))) do i
        ξn = dot(ns[i], localgeo.coord)
        gaussian(ξn, t, α, β, μ, δ)
    end
    state.ρ = length(ns) == 1 ? ρ[1] : ρ
end

inhomogeneous_data!(::Val{0}, P::Pseudo1D, x...) = initial_condition!(P, x...)

function inhomogeneous_data!(
    ::Val{1},
    ::Pseudo1D{ns, α, β, μ, δ},
    ∇state,
    aux,
    x,
    t,
) where {ns, α, β, μ, δ}
    ∇state.ρ = hcat(ntuple(Val(length(ns))) do i
        ξn = dot(ns[i], x)
        ∇gaussian(ns[i], ξn, t, α, β, μ, δ)
    end...)
end

function do_output(mpicomm, vtkdir, vtkstep, dgfvm, Q, Qe, model, testname)
    ## Name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
    exactnames = statenames .* "_exact"

    writevtk(filename, Q, dgfvm, statenames, Qe, exactnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## Name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## Name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(
            pvtuprefix,
            prefixes,
            (statenames..., exactnames...),
            eltype(Q),
        )

        @info "Done writing VTK: $pvtuprefix"
    end
end

function test_run(
    mpicomm,
    dim,
    fvmethod,
    polynomialorders,
    level,
    ArrayType,
    FT,
    vtkdir,
    direction,
)

    n_hd =
        dim == 2 ? SVector{3, FT}(1, 0, 0) :
        SVector{3, FT}(1 / sqrt(2), 1 / sqrt(2), 0)

    n_vd = dim == 2 ? SVector{3, FT}(0, 1, 0) : SVector{3, FT}(0, 0, 1)

    n_dg =
        dim == 2 ? SVector{3, FT}(1 / sqrt(2), 1 / sqrt(2), 0) :
        SVector{3, FT}(1 / sqrt(3), 1 / sqrt(3), 1 / sqrt(3))
    connectivity = dim == 2 ? :face : :full
    if direction isa EveryDirection
        ns = (n_hd, n_vd, n_dg)
    elseif direction isa HorizontalDirection
        ns = (n_hd,)
    elseif direction isa VerticalDirection
        ns = (n_vd,)
    end

    α = FT(1)
    β = FT(1 // 100)
    μ = FT(-1 // 2)
    δ = FT(1 // 10)

    # Grid/topology information
    base_num_elem = 4
    Ne = 2^(level - 1) * base_num_elem
    N = polynomialorders
    L = ntuple(j -> FT(j == dim ? 1 : N[1]) / 4, dim)
    brickrange = ntuple(j -> range(-L[j]; length = Ne + 1, stop = L[j]), dim)
    periodicity = ntuple(j -> false, dim)
    bc = ntuple(j -> (1, 2), dim)

    topl = StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = periodicity,
        boundary = bc,
        connectivity = connectivity,
    )

    dt = (α / 4) * L[1] / (Ne * polynomialorders[1]^2)
    timeend = 1
    outputtime = timeend / 10
    @info "time step" dt

    @info @sprintf """Test parameters:
    FVM Reconstruction          = %s
    ArrayType                   = %s
    FloatType                   = %s
    Dimension                   = %s
    Direction                   = %s
    Horizontal polynomial order = %s
    Vertical polynomial order   = %s
      """ fvmethod ArrayType FT dim direction polynomialorders[1] polynomialorders[end]

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorders,
    )

    bcs = (InhomogeneousBC{0}(), InhomogeneousBC{1}())
    # Model being tested
    model = AdvectionDiffusion{dim}(
        Pseudo1D{ns, α, β, μ, δ}(),
        bcs,
        num_equations = length(ns),
    )

    # Main DG discretization
    dgfvm = DGFVModel(
        model,
        grid,
        fvmethod,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        direction = direction,
    )

    # Initialize all relevant state arrays and create solvers
    Q = init_ode_state(dgfvm, FT(0))

    eng0 = norm(Q, dims = (1, 3))
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0[1]

    solver = LSRK54CarpenterKennedy(dgfvm, Q; dt = dt, t0 = 0)

    # Set up the information callback
    starttime = Ref(Dates.now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = Dates.now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %.16e
                runtime = %s
                norm(Q) = %.16e""",
                gettime(solver),
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end
    callbacks = (cbinfo,)
    if ~isnothing(vtkdir)
        # Create vtk dir
        mkpath(vtkdir)

        vtkstep = 0
        # Output initial step
        do_output(
            mpicomm,
            vtkdir,
            vtkstep,
            dgfvm,
            Q,
            Q,
            model,
            "advection_diffusion",
        )

        # Setup the output callback
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            Qe = init_ode_state(dgfvm, gettime(solver))
            do_output(
                mpicomm,
                vtkdir,
                vtkstep,
                dgfvm,
                Q,
                Qe,
                model,
                "advection_diffusion",
            )
        end
        callbacks = (callbacks..., cbvtk)
    end

    solve!(Q, solver; timeend = timeend, callbacks = callbacks)

    # Reference solution
    engf = norm(Q, dims = (1, 3))
    Q_ref = init_ode_state(dgfvm, FT(timeend))

    engfe = norm(Q_ref, dims = (1, 3))
    errf = norm(Q_ref .- Q, dims = (1, 3))

    metrics = @. (engf, engf / eng0, engf - eng0, errf, errf / engfe)

    @info begin
        j = 1
        msg = "Finished\n"
        if direction isa HorizontalDirection || direction isa EveryDirection
            msg *= @sprintf """
            Horizontal field:
              norm(Q)                 = %.16e
              norm(Q) / norm(Q₀)      = %.16e
              norm(Q) - norm(Q₀)      = %.16e
              norm(Q - Qe)            = %.16e
              norm(Q - Qe) / norm(Qe) = %.16e
            """ ntuple(f -> metrics[f][j], 5)...
            j += 1
        end

        if direction isa VerticalDirection || direction isa EveryDirection
            msg *= @sprintf """
            Vertical field:
              norm(Q)                 = %.16e
              norm(Q) / norm(Q₀)      = %.16e
              norm(Q) - norm(Q₀)      = %.16e
              norm(Q - Qe)            = %.16e
              norm(Q - Qe) / norm(Qe) = %.16e
            """ ntuple(f -> metrics[f][j], 5)...
            j += 1
        end

        if direction isa EveryDirection
            msg *= @sprintf """
            Diagonal field:
              norm(Q)                 = %.16e
              norm(Q) / norm(Q₀)      = %.16e
              norm(Q) - norm(Q₀)      = %.16e
              norm(Q - Qe)            = %.16e
              norm(Q - Qe) / norm(Qe) = %.16e
            """ ntuple(f -> metrics[f][j], 5)...
        end
        msg
    end

    return Tuple(errf)
end

"""
    main()

Run this test problem
"""
function main()

    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD

    # Dictionary keys: dim, level, and FT
    expected_result = Dict()
    expected_result[2, 1, Float32, FVConstant()] =
        (2.3391991853713989e-02, 5.0738707184791565e-02, 4.1857458651065826e-02)
    expected_result[2, 2, Float32, FVConstant()] =
        (2.0331495907157660e-03, 3.3669386059045792e-02, 2.3904174566268921e-02)
    expected_result[2, 3, Float32, FVConstant()] =
        (2.6557327146292664e-05, 2.0002065226435661e-02, 1.2790882028639317e-02)

    expected_result[2, 1, Float32, FVLinear()] =
        (2.3391995579004288e-02, 4.9536284059286118e-02, 3.5911198705434799e-02)
    expected_result[2, 2, Float32, FVLinear()] =
        (2.0331430714577436e-03, 2.2621221840381622e-02, 1.5531544573605061e-02)
    expected_result[2, 3, Float32, FVLinear()] =
        (2.6568108296487480e-05, 7.5304862111806870e-03, 6.4526572823524475e-03)

    expected_result[3, 1, Float32, FVConstant()] =
        (8.7377587333321571e-03, 7.1755334734916687e-02, 4.3733172118663788e-02)
    expected_result[3, 2, Float32, FVConstant()] =
        (6.2517996411770582e-04, 4.7615684568881989e-02, 2.3986011743545532e-02)
    expected_result[3, 3, Float32, FVConstant()] =
        (3.5004395613214001e-05, 2.8287241235375404e-02, 1.2639496475458145e-02)

    expected_result[3, 1, Float32, FVLinear()] =
        (8.7377615272998810e-03, 7.0054858922958374e-02, 3.7571556866168976e-02)
    expected_result[3, 2, Float32, FVLinear()] =
        (6.2518107006326318e-04, 3.1991228461265564e-02, 1.6405479982495308e-02)
    expected_result[3, 3, Float32, FVLinear()] =
        (3.5004966775886714e-05, 1.0649790056049824e-02, 7.1499347686767578e-03)

    expected_result[2, 1, Float64, FVConstant()] =
        (2.3391871809628567e-02, 5.0738761087541523e-02, 4.1857480220018319e-02)
    expected_result[2, 2, Float64, FVConstant()] =
        (2.0332783913617892e-03, 3.3669399006499491e-02, 2.3904204301839819e-02)
    expected_result[2, 3, Float64, FVConstant()] =
        (2.6572168347086839e-05, 2.0002110124502395e-02, 1.2790993871268752e-02)
    expected_result[2, 4, Float64, FVConstant()] =
        (1.9890000550154039e-07, 1.0929144069594809e-02, 6.5110938897763263e-03)

    expected_result[2, 1, Float64, FVLinear()] =
        (2.3391871809628550e-02, 4.9536356061135857e-02, 3.5911213637569099e-02)
    expected_result[2, 2, Float64, FVLinear()] =
        (2.0332783913618278e-03, 2.2621252826152356e-02, 1.5531565558700888e-02)
    expected_result[2, 3, Float64, FVLinear()] =
        (2.6572168347111800e-05, 7.5305291439555161e-03, 6.4526740563561544e-03)
    expected_result[2, 4, Float64, FVLinear()] =
        (1.9890000549995218e-07, 2.5302226025389427e-03, 2.7792905025059711e-03)

    expected_result[3, 1, Float64, FVConstant()] =
        (8.7378337431297422e-03, 7.1755444068009461e-02, 4.3733196512722658e-02)
    expected_result[3, 2, Float64, FVConstant()] =
        (6.2510740807095622e-04, 4.7615720711942776e-02, 2.3986017606198881e-02)
    expected_result[3, 3, Float64, FVConstant()] =
        (3.4995405318038341e-05, 2.8287255414151377e-02, 1.2639742577376042e-02)
    expected_result[3, 4, Float64, FVConstant()] =
        (1.4362091045094841e-06, 1.5456143768350493e-02, 6.3677406803847331e-03)

    expected_result[3, 1, Float64, FVLinear()] =
        (8.7378337431297439e-03, 7.0054986572200995e-02, 3.7571599264936015e-02)
    expected_result[3, 2, Float64, FVLinear()] =
        (6.2510740807095286e-04, 3.1991282544615307e-02, 1.6405489038457288e-02)
    expected_result[3, 3, Float64, FVLinear()] =
        (3.4995405318035868e-05, 1.0649776447227678e-02, 7.1502921502446283e-03)
    expected_result[3, 4, Float64, FVLinear()] =
        (1.4362091045110612e-06, 3.5782751203335653e-03, 3.1186151510318319e-03)

    @testset "Variable degree DG: advection diffusion model" begin
        for FT in (Float32, Float64)
            for direction in
                (EveryDirection(), HorizontalDirection(), VerticalDirection())
                numlevels =
                    integration_testing ||
                    ClimateMachine.Settings.integration_testing ?
                    (FT == Float64 ? 4 : 3) : 1
                for dim in 2:3
                    for fvmethod in (FVConstant(), FVLinear(), FVLinear{3}())
                        polynomialorders = (4, 0)
                        result = Dict()
                        for level in 1:numlevels
                            vtkdir =
                                output ?
                                "vtk_advection" *
                                "_poly$(polynomialorders)" *
                                "_dim$(dim)_$(ArrayType)_$(FT)" *
                                "_fvmethod$(fvmethod)" *
                                "_level$(level)" :
                                nothing
                            result[level] = test_run(
                                mpicomm,
                                dim,
                                fvmethod,
                                polynomialorders,
                                level,
                                ArrayType,
                                FT,
                                vtkdir,
                                direction,
                            )
                            fv_key =
                                fvmethod isa FVLinear ? FVLinear() : fvmethod
                            if direction isa EveryDirection
                                @test all(
                                    result[level] .≈
                                    FT.(expected_result[
                                        dim,
                                        level,
                                        FT,
                                        fv_key,
                                    ]),
                                )
                            elseif direction isa HorizontalDirection
                                @test result[level][1] .≈ FT.(expected_result[
                                    dim,
                                    level,
                                    FT,
                                    fv_key,
                                ])[1]
                            elseif direction isa VerticalDirection
                                @test result[level][1] .≈ FT.(expected_result[
                                    dim,
                                    level,
                                    FT,
                                    fv_key,
                                ])[2]
                            end
                        end
                        @info begin
                            msg = ""
                            for l in 1:(numlevels - 1)
                                rate = @. log2(result[l]) - log2(result[l + 1])
                                msg *= @sprintf("\n  rates for level %d", l)
                                j = 1
                                if direction isa HorizontalDirection ||
                                   direction isa EveryDirection
                                    msg *=
                                        @sprintf(", Horizontal = %e", rate[j])
                                    j += 1
                                end
                                if direction isa VerticalDirection ||
                                   direction isa EveryDirection
                                    msg *= @sprintf(", Vertical = %e", rate[j])
                                    j += 1
                                end
                                if direction isa EveryDirection
                                    msg *= @sprintf(", Diagonal = %e", rate[j])
                                end
                                msg *= "\n"
                            end
                            msg
                        end
                    end
                end
            end
        end
    end
end

main()


text:
# WaveletLayerState
immutable WaveletLayerState{B<:ScatteredBlob} <: AbstractScatteredLayerState
    blobs::Vector{B}
    layer::WaveletLayer
end

function Mocha.setup(
        backend::Mocha.Backend,
        layer::WaveletLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    pathkey = layer.bank.behavior.pathkey
    chromakey = prepend(:χ, layer.bank.behavior.pathkey)
    chromarange = chromakey => 0:1:(layer.bank.spec.n_filters_per_octave-1)
    octavekey = prepend(:j, layer.bank.behavior.pathkey)
    innodes = inputs[1].nodes
    outnodes = DataStructures.SortedDict{Path,AbstractNode,
        Base.Order.ForwardOrdering}()
    for inpath in keys(innodes)
        innode = innodes[inpath]
        inranges = innode.ranges
        inpathkeys = [ pair.first for pair in collect(inranges) ]
        subscripts = find(pathkey .== inpathkeys)
        insizes = collect(size(innode.data))
        if isa(innode, RealFourierNode)
            insizes[subscripts] = 2 * (insizes[subscripts] - 1)
        end
        outranges = [inranges..., chromarange]
        outsizes = [insizes ; layer.bank.spec.n_filters_per_octave]
        for j in layer.bank.behavior.j_range
            ψ_log2_sampling = layer.bank.behavior.ψ_log2_samplings[1+j]
            for subscript in subscripts
                outsizes[subscript] =
                    insizes[subscript] >> (-ψ_log2_sampling)
                inrange = inranges[subscript].second
                outranges[subscript] = outranges[subscript].first =>
                    (inrange.start):(1<<-ψ_log2_sampling):(inrange.stop)
            end
            outdata = zeros(eltype(innode.data), tuple(outsizes...))
            inds = [fill(Colon(), ndims(innode.data)) ; 0]
            for χ in 0:(layer.bank.spec.n_filters_per_octave-1)
                ψ = layer.bank.ψs[1, 1+χ, 1+j]
                inds[end] = 1 + χ
                transform!(view(outdata, inds...), ψ, innode, subscripts...)
            end
            outpath = Path(octavekey => j)
            outnodes[outpath] = Node(outdata, tuple(outranges...))
        end
    end
    WaveletLayerState([ScatteredBlob(outnodes)], layer)
end


text:
using Graphs
using SparseArrays

function dijkstra_digraph(mat, nrows, ncols, startindex, endindex)
    dist = spzeros(Int, length(mat), length(mat))
    g = SimpleDiGraph(length(mat))
    for i in 1:length(mat)
        if i % nrows != 0
            add_edge!(g, i, i + 1)
            dist[i, i + 1] = mat[i + 1]
            add_edge!(g, i + 1, i)
            dist[i + 1, i] = mat[i]
        end
        if i <= length(mat) - nrows
            add_edge!(g, i, i + nrows)
            dist[i, i + nrows] = mat[i + nrows]
            add_edge!(g, i + nrows, i)
            dist[i + nrows, i] = mat[i]
        end
    end
    return dijkstra_shortest_paths(g, startindex, dist).dists[endindex]
end

function day15()
    part = [0, 0]
    lines = filter(!isempty, strip.(readlines("AoCdata/AoC_2021_day15.txt")))
    nrows, ncols = length(lines), length(first(lines))
    mat = zeros(Int, nrows, ncols)
    for i in 1:nrows
        mat[i, :] = parse.(Int, collect(lines[i]))
    end
    part[1] = dijkstra_digraph(mat, nrows, ncols, 1, length(mat))

    mat2 = zeros(Int, ncols * 5, nrows * 5)
    for c in CartesianIndices(mat2)
        mat2[c] = mod1(
            mat[mod1(c[1], nrows), mod1(c[2], ncols)] +
            (c[1] - 1) ÷ nrows +
            (c[2] - 1) ÷ ncols,
            9,
        )
    end
    part[2] = dijkstra_digraph(mat2, nrows * 5, ncols * 5, 1, length(mat2))

    return part
end

part = day15()
println("Part 1: ", part[1])
println("Part 2: ", part[2])


text:
using DelaunayTriangulations

@testset "Geometry Predicates" begin
    @testset "Orient" begin
        @testset "Orient Fast" begin
            @test orient_fast((0.1, 0.1), (0.2, 0.2), (0.3, 0.3)) == 0.0
            @test orient_fast((0.0, 3.14159e23), (0.0, 42.4242e11), (0.0, 1.234567e2)) == 0.0
            @test orient_fast((-1.0e-7, -1.0e-7), (2.0, 2.0), (3.0e7, 3.0e7)) == 0.0
            @test orient_fast((Float32(0.1), Float32(0.1)), (Float32(0.2), Float32(0.2)), (Float32(0.3), Float32(0.3))) == Float32(0.0)
            @test orient_fast((Float16(0.1), Float16(0.1)), (Float16(0.2), Float16(0.2)), (Float16(0.3), Float16(0.3))) == Float16(0.0)
            @test orient_fast((0.0, 0.0), (0.1, 0.1), (0.3, 0.2)) < 0.0
            @test orient_fast((0.0, 0.0), (0.1, 0.1), (0.1, 0.2)) > 0.0
            # Case where orient fast fails
            @test orient_fast((0.0, 0.0), (0.1, 0.1), (nextfloat(0.4), 0.4)) > 0.0 # Should be < 0.0
        end
        @testset "Orient Exact" begin
            @test orient_exact((0.1, 0.1), (0.2, 0.2), (0.3, 0.3)) == 0.0
            @test orient_exact((0.0, 3.14159e23), (0.0, 42.4242e11), (0.0, 1.234567e2)) == 0.0
            @test orient_exact((-1.0e-7, -1.0e-7), (2.0, 2.0), (3.0e7, 3.0e7)) == 0.0
            @test orient_exact((Float32(0.1), Float32(0.1)), (Float32(0.2), Float32(0.2)), (Float32(0.3), Float32(0.3))) == Float32(0.0)
            @test orient_exact((Float16(0.1), Float16(0.1)), (Float16(0.2), Float16(0.2)), (Float16(0.3), Float16(0.3))) == Float16(0.0)
            @test orient_exact((0.0, 0.0), (0.1, 0.1), (0.3, 0.2)) < 0.0
            @test orient_exact((0.0, 0.0), (0.1, 0.1), (0.1, 0.2)) > 0.0
            # Case where orient fast fails
            @test orient_exact((0.0, 0.0), (0.1, 0.1), (nextfloat(0.4), 0.4)) < 0.0 # Should be < 0.0
        end
        @testset "Orient" begin
            @test orient((0.1, 0.1), (0.2, 0.2), (0.3, 0.3)) == 0.0
            @test orient((0.0, 3.14159e23), (0.0, 42.4242e11), (0.0, 1.234567e2)) == 0.0
            @test orient((-1.0e-7, -1.0e-7), (2.0, 2.0), (3.0e7, 3.0e7)) == 0.0
            @test orient((Float32(0.1), Float32(0.1)), (Float32(0.2), Float32(0.2)), (Float32(0.3), Float32(0.3))) == Float32(0.0)
            @test orient((Float16(0.1), Float16(0.1)), (Float16(0.2), Float16(0.2)), (Float16(0.3), Float16(0.3))) == Float16(0.0)
            @test orient((0.0, 0.0), (0.1, 0.1), (0.3, 0.2)) < 0.0
            @test orient((0.0, 0.0), (0.1, 0.1), (0.1, 0.2)) > 0.0
            # Case where orient fast fails
            @test orient((0.0, 0.0), (0.1, 0.1), (nextfloat(0.4), 0.4)) < 0.0 # Should be < 0.0
        end
    end

    @testset "Incircle" begin
        @testset "Incircle Fast" begin
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, 0.0)) > 0.0
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (1.0, 0.0)) < 0.0
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, -0.5)) == 0.0
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.500000001, 0.0)) < 0.0
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.499999999, 0.0)) > 0.0
            @test incircle_fast((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.0))) > Float32(0.0)
            @test incircle_fast((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(1.0), Float32(0.0))) < Float32(0.0)
            @test incircle_fast((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(0.0), -Float32(0.5))) == Float32(0.0)
            @test incircle_fast((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.0))) > Float16(0.0)
            @test incircle_fast((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(1.0), Float16(0.0))) < Float16(0.0)
            @test incircle_fast((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(0.0), -Float16(0.5))) == Float16(0.0)
            # Degenerate circle
            incircle_fast((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, 0.0)) == 0.0
            incircle_fast((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, 0.1)) < 0.0
            incircle_fast((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, -0.1)) > 0.0
            # Case where incircle fast fails
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, prevfloat(-0.5))) == 0.0 # should be < 0.0
            @test incircle_fast((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, nextfloat(-0.5))) == 0.0 # should be > 0.0
        end
        @testset "Incircle Exact" begin
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, 0.0)) > 0.0
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (1.0, 0.0)) < 0.0
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, -0.5)) == 0.0
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.500000001, 0.0)) < 0.0
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.499999999, 0.0)) > 0.0
            @test incircle_exact((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.0))) > Float32(0.0)
            @test incircle_exact((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(1.0), Float32(0.0))) < Float32(0.0)
            @test incircle_exact((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(0.0), -Float32(0.5))) == Float32(0.0)
            @test incircle_exact((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.0))) > Float16(0.0)
            @test incircle_exact((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(1.0), Float16(0.0))) < Float16(0.0)
            @test incircle_exact((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(0.0), -Float16(0.5))) == Float16(0.0)
            # Degenerate circle
            incircle_exact((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, 0.0)) == 0.0
            incircle_exact((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, 0.1)) < 0.0
            incircle_exact((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, -0.1)) > 0.0
            # Case where incircle fast fails
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, prevfloat(-0.5))) < 0.0 # should be < 0.0
            @test incircle_exact((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, nextfloat(-0.5))) > 0.0 # should be > 0.0
        end
        @testset "Incircle" begin
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, 0.0)) > 0.0
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (1.0, 0.0)) < 0.0
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, -0.5)) == 0.0
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.500000001, 0.0)) < 0.0
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.499999999, 0.0)) > 0.0
            @test incircle((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.0))) > Float32(0.0)
            @test incircle((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(1.0), Float32(0.0))) < Float32(0.0)
            @test incircle((Float32(0.5), Float32(0.0)), (Float32(0.0), Float32(0.5)), (-Float32(0.5), Float32(0.0)), (Float32(0.0), -Float32(0.5))) == Float32(0.0)
            @test incircle((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.0))) > Float16(0.0)
            @test incircle((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(1.0), Float16(0.0))) < Float16(0.0)
            @test incircle((Float16(0.5), Float16(0.0)), (Float16(0.0), Float16(0.5)), (-Float16(0.5), Float16(0.0)), (Float16(0.0), -Float16(0.5))) == Float16(0.0)
            # Degenerate circle
            incircle((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, 0.0)) == 0.0
            incircle((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, 0.1)) < 0.0
            incircle((0.0, 0.0), (0.5, 0.0), (0.7, 0.0), (0.8, -0.1)) > 0.0
            # Case where incircle fast fails
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, prevfloat(-0.5))) < 0.0
            @test incircle((0.5, 0.0), (0.0, 0.5), (-0.5, 0.0), (0.0, nextfloat(-0.5))) > 0.0
        end
        @testset "Ghost Incircle" begin
            @testset "Upper Half Plane" begin
                @test incircle((0.42, floatmin()), (0.0, 0.0), (0.5, 0.0)) > 0.0
                @test incircle((0.42, -floatmin()), (0.0, 0.0), (0.5, 0.0)) < 0.0
                @test incircle((floatmin(), 0.42), (0.0, 0.0), (0.0, 0.5)) < 0.0
                @test incircle((-floatmin(), 0.42), (0.0, 0.0), (0.0, 0.5)) > 0.0
            end
            @testset "In Segment" begin
                @test incircle((-floatmin(), 0.0), (0.0, 0.0), (0.5, 0.0)) == 0.0
                @test incircle((0.0, 0.0), (0.0, 0.0), (0.5, 0.0)) == 0.0
                @test incircle((floatmin(), 0.0), (0.0, 0.0), (0.5, 0.0)) == 1.0
                @test incircle((prevfloat(0.5), 0.0), (0.0, 0.0), (0.5, 0.0)) == 1.0
                @test incircle((0.5, 0.0), (0.0, 0.0), (0.5, 0.0)) == 0.0
                @test incircle((nextfloat(0.5), 0.0), (0.0, 0.0), (0.5, 0.0)) == 0.0

                @test incircle((0.0, -floatmin()), (0.0, 0.0), (0.0, 0.5)) == 0.0
                @test incircle((0.0, 0.0), (0.0, 0.0), (0.0, 0.5)) == 0.0
                @test incircle((0.0, floatmin()), (0.0, 0.0), (0.0, 0.5)) == 1.0
                @test incircle((0.0, prevfloat(0.5)), (0.0, 0.0), (0.0, 0.5)) == 1.0
                @test incircle((0.0, 0.5), (0.0, 0.0), (0.0, 0.5)) == 0.0
                @test incircle((0.0, nextfloat(0.5)), (0.0, 0.0), (0.0, 0.5)) == 0.0
            end
        end
    end
end

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
using .TimeZones: ZonedDateTime

# Promote Rule for ZonedDateTime and InfExtendedTime
Base.promote_rule(::Type{InfExtendedTime{T}}, ::Type{S}) where {T<:TimeType, S<:ZonedDateTime} = InfExtendedTime(promote_type(T, S))