text:
using DataStructures, AutoHashEquals, Base.Dates, Unitful, UnitfulAngles

import Base: push!, empty!, delete!, isempty, ==, in, show

export  VideoFile, Point, POI, Run, Repetition, Association,
        replace!, findVideoFiles, getVideoFiles, save,
        loadLogVideoFiles, loadPOIs, loadRuns, loadAssociation

const exiftool_base = joinpath(Pkg.dir("BeetleWay"), "deps", "src", "exiftool", "exiftool")
const exiftool = exiftool_base*(is_windows() ? ".exe" : "")
const exts = [".webm", ".mkv", ".flv", ".flv", ".vob", ".ogv", ".ogg", ".drc", ".mng", ".avi", ".mov", ".qt", ".wmv", ".yuv", ".rm", ".rmvb", ".asf", ".amv", ".mp4", ".m4p", ".m4v", ".mpg", ".mp2", ".mpeg", ".mpe", ".mpv", ".mpg", ".mpeg", ".m2v", ".m4v", ".svi", ".3gp", ".3g2", ".mxf", ".roq", ".nsv", ".flv", ".f4v", ".f4p", ".f4a", ".f4b", ".MTS", ".DS_Store"]

@auto_hash_equals immutable VideoFile
    file::String
    datetime::DateTime
end

function getDateTime(folder::String, file::String)::DateTime
    fullfile = joinpath(folder, file)
    dateTimeOriginal, createDate, modifyDate = strip.(split(readstring(`$exiftool -T -AllDates -n $fullfile`), '\t'))
    datetime = DateTime(now())
    for i in [dateTimeOriginal, createDate, modifyDate]
        m = match(r"^(\d\d\d\d):(\d\d):(\d\d) (\d\d):(\d\d):(\d\d).?(\d?\d?\d?)", i)
        m == nothing && continue
        ts = zeros(Int, 7)
        for (j, t) in enumerate(m.captures)
            if !isempty(t)
                ts[j] = parse(Int, t)
            else
                break
            end
        end
        isnull(Dates.validargs(DateTime, ts...)) || continue
        datetime = min(datetime, DateTime(ts...))
    end
    return datetime
end

@auto_hash_equals immutable Point
    file::String
    time::Second
end

Point(;file = "", time = Second(0)) = Point(file, time)
Point(f::String, t::Time) = Point(f, Second(t.instant))
Point(f::String, h::Int, m::Int, s::Int) = Point(f, Time(h, m, s))

@auto_hash_equals immutable POI
    name::String
    start::Point
    stop::Point
    label::String
    comment::String

    function POI(name, start, stop, label, comment)
        # @assert start.file != stop.file || start.time <= stop.time
        new(name, start, stop, label, comment)
    end
end

POI(;name = "", start = Point(), stop = Point(), label = "", comment = "") = POI(name, start, stop, label, comment)

# time differences
function duration(poi::POI, folder::String)::Int
    Δ = poi.stop.time.value - poi.start.time.value
    poi.start.file == poi.stop.file && return Δ
    fullfile = joinpath(folder, poi.start.file)
    t = round(Int, parse(readstring(`$exiftool -T -Duration -n $fullfile`)))
    return t + Δ
end

# time of day
function timeofday(poi::POI, folder::String)
    filescsv = joinpath(folder, "log", "files.csv")
    a, _ = readcsv(filescsv, String, header = true, quotes = true, comments = false)
    for i = 1:size(a,1)
        strip(a[i, 1]) == poi.start.file && return Time(DateTime(strip(a[i, 2])) + poi.start.time)
    end
end

@auto_hash_equals immutable Run
    metadata::Dict{Symbol, String}
    comment::String
end

Run(;metadata = Dict{Symbol, String}(), comment = "") = Run(metadata, comment)

@auto_hash_equals immutable Repetition
    run::Run
    repetition::Int
end

@auto_hash_equals immutable Association
    pois::OrderedSet{POI}
    runs::OrderedSet{Repetition}
    associations::Set{Tuple{POI, Repetition}}
end

Association() = Association(OrderedSet{POI}(), OrderedSet{Repetition}(), Set{Tuple{POI, Repetition}}())

# in

in(x::POI, a::Association) = x in a.pois
in(x::Repetition, a::Association) = x in a.runs
in(x::Tuple{POI, Repetition}, a::Association) = x in a.associations

# equal keys
==(a::Base.KeyIterator, b::Base.KeyIterator) = length(a)==length(b) && all(k->in(k,b), a)

# pushes

function run2repetition(xs::OrderedSet{Repetition}, r::Run)
    isempty(xs) || @assert keys(last(xs).run.metadata) == keys(r.metadata)
    Repetition(r, reduce((x, y) -> max(x, y.run.metadata == r.metadata ? y.repetition : 0), 0, xs) + 1)
end

push!(xs::OrderedSet{Repetition}, r::Run) = push!(xs, run2repetition(xs, r))

function push!(a::Association, t::POI)
    push!(a.pois, t)
    return a
end

function push!(a::Association, r::Run)
    push!(a.runs, r)
    return a
end

function push!(a::Association, x::Tuple{POI, Repetition})
    @assert first(x) in a.pois
    @assert last(x) in a.runs
    push!(a.associations, x)
    return a
end

# replace

## runs

function replace!(a::Association, o::Run, n::Run)
    # @assert o in a
    @assert keys(last(a.runs).run.metadata) == keys(o.metadata) == keys(n.metadata)

    runs = OrderedSet{Repetition}()
    associations = Set{Tuple{POI, Repetition}}()
    for r1 in a.runs
        if r1.run == o
            push!(runs, n)
        else
            push!(runs, r1.run)
        end
        r2 = last(runs)
        for (p, r) in a.associations
            if r1 == r
                push!(associations, (p, r2))
            end
        end
    end

    empty!(a.runs)
    push!(a.runs, runs...)
    isempty(associations) && return a
    empty!(a.associations)
    push!(a.associations, associations...)
    return a
end

function replace!(a::Association, o::Repetition, n::Run)
    @assert o in a
    @assert keys(last(a.runs).run.metadata) == keys(o.run.metadata) == keys(n.metadata)

    runs = OrderedSet{Repetition}()
    associations = Set{Tuple{POI, Repetition}}()
    for r1 in a.runs
        if r1 == o
            push!(runs, n)
        else
            push!(runs, r1.run)
        end
        r2 = last(runs)
        for (p, r) in a.associations
            if r1 == r
                push!(associations, (p, r2))
            end
        end
    end

    empty!(a.runs)
    push!(a.runs, runs...)
    isempty(associations) && return a
    empty!(a.associations)
    push!(a.associations, associations...)
    return a
end

## pois

replace!(xs::OrderedSet{POI}, o::POI, n::POI) = OrderedSet{POI}(x == o ? n : x for x in xs)
replace!(xs::Set{Tuple{POI, Repetition}}, o::POI, n::POI) = Set{Tuple{POI, Repetition}}(first(x) == o ? (n, last(x)) : x for x in xs)
function replace!(a::Association, o::POI, n::POI)
    o == n && return a
    @assert o in a
    pois = replace!(a.pois, o, n)
    empty!(a.pois)
    push!(a.pois, pois...)
    isempty(a.associations) && return a
    associations = replace!(a.associations, o, n)
    empty!(a.associations)
    push!(a.associations, associations...)
    return a
end

# deletes

function delete!(a::Association, r::Repetition)
    r in a || return a
    delete!(a.runs, r)
    filter!(x -> last(x) != r, a.associations)
    for x in a.runs
        if x.run.metadata == r.run.metadata && x.repetition > r.repetition
            replace!(a, x, x.run)
        end
    end
    return a
end

function delete!(a::Association, p::POI)
    p in a || return a
    delete!(a.pois, p)
    filter!(x -> first(x) != p, a.associations)
    return a
end

function delete!(a::Association, x::Tuple{POI, Repetition})
    x in a || return a
    @assert first(x) in a.pois
    @assert last(x) in a.runs
    delete!(a.associations, x)
    return a
end

# saves

function prep_file(folder::String, what::String)::String
    folder = joinpath(folder, "log")
    isdir(folder) || mkdir(folder)
    return joinpath(folder, "$what.csv")
end

function save(folder::String, x::OrderedSet{VideoFile})
    file = prep_file(folder, "files")
    #isempty(x) && rm(file, force=true)
    n = length(x)
    a = Matrix{String}(n + 1,2)
    a[1,:] .= ["file", "date and time"]
    for (i, v) in enumerate(x)
        a[i + 1, :] .= [v.file, string(v.datetime)]
    end
    writecsv(file, a)
end

function save(folder::String, x::OrderedSet{POI})
    file = prep_file(folder, "pois")
    n = length(x)
    a = Matrix{String}(n + 1,7)
    a[1,:] .= ["name", "start file", "start time (seconds)", "stop file", "stop time (seconds)", "label", "comments"]
    for (i, t) in enumerate(x)
        a[i + 1, :] .= [t.name, t.start.file, string(t.start.time.value), t.stop.file, string(t.stop.time.value), t.label, t.comment]
    end
    a .= strip.(a)
    writecsv(file, a)
end

function save(folder::String, x::OrderedSet{Repetition})
    ks = keys(first(x).run.metadata)
    @assert reduce((x, y) -> x && keys(y.run.metadata) == ks, true, x)
    file = prep_file(folder, "runs")
    ks = sort(collect(ks))
    header = string.(ks)
    push!(header, "Comment", "Repetition")
    n = length(x)
    a = Matrix{String}(n + 1, length(header))
    a[1,:] .= header
    for (i, r) in enumerate(x)
        for (j, k) in enumerate(ks)
            a[i + 1, j] = r.run.metadata[k]
        end
        a[i + 1, end - 1] = r.run.comment
        a[i + 1, end] = string(r.repetition)
    end
    a .= strip.(a)
    writecsv(file, a)
end

function save(folder::String, a::Association)
    save(folder, a.pois)
    save(folder, a.runs)
    file = prep_file(folder, "associations")
    open(file, "w") do o
        println(o, "POI number, run number")
        for (t, r) in a.associations
            ti = findfirst(a.pois, t)
            ri = findfirst(a.runs, r)
            println(o, ti, ",", ri)
        end
    end
end

# loads

function loadLogVideoFiles(folder::String)::Dict{String, DateTime}
    filescsv = joinpath(folder, "log", "files.csv")
    vfs = Dict{String, DateTime}()
    if isfile(filescsv)
        a, _ = readcsv(filescsv, String, header = true, quotes = true, comments = false)
        a .= strip.(a)
        @assert allunique(a[:,1])
        nrow, ncol = size(a)
        @assert ncol == 2
        for i = 1:nrow
            vfs[a[i, 1]] = DateTime(a[i, 2])
        end
    end
    return vfs
end

function findVideoFiles!(log::Dict{String, DateTime}, folder::String)
    for (root, dir, files) in walkdir(folder)
        for file in files
            file[1] == '.' && continue
            last(splitext(file)) in exts || continue
            fname = relpath(joinpath(root, file), folder)
            if !haskey(log, fname)
                log[fname] = getDateTime(folder, fname)
            end
        end
    end
end

function getVideoFiles(folder::String)::OrderedSet{VideoFile}
    log = loadLogVideoFiles(folder)
    findVideoFiles!(log, folder)
    o = sort([(k, v) for (k, v) in log], by = last)
    return OrderedSet{VideoFile}(VideoFile(k, v) for (k,v) in o)
end

function loadPOIs(folder::String)::OrderedSet{POI}
    filescsv = joinpath(folder, "log", "pois.csv")
    tgs = OrderedSet{POI}()
    if isfile(filescsv)
        a, _ = readcsv(filescsv, String, header = true, quotes = true, comments = false)
        a .= strip.(a)
        nrow, ncol = size(a)
        @assert ncol == 7
        for i = 1:nrow
            tg = POI(a[i, 1], Point(a[i, 2], Second(parse(Int, a[i, 3]))), Point(a[i, 4], Second(parse(Int, a[i, 5]))), a[i, 6], a[i, 7])
            @assert !(tg in tgs)
            push!(tgs, tg)
        end
    end
    return tgs
end

function getmetadata(folder)
    metadata = Dict{Symbol, String}()
    file = joinpath(folder, "metadata", "run.csv")
    if isfile(file)
        b = readcsv(file, comments = false)
        for i = 1:size(b,1)
            metadata[Symbol(strip(b[i, 1]))] = length(b[i,:]) < 2 ? "" : strip(b[i, 2])
        end
    end
    return metadata
end

function loadRuns(folder::String)::OrderedSet{Repetition}
    filescsv = joinpath(folder, "log", "runs.csv")
    rs = OrderedSet{Repetition}()
    if isfile(filescsv)
        a, ks = readcsv(filescsv, String, header = true, quotes = true, comments = false)
        ks = Symbol.(strip.(vec(ks)))
        a .= strip.(a)
        metadata = getmetadata(folder)
        for (k, v) in metadata
            if !(k in ks)
                unshift!(ks, k)
                a = [repmat([v], size(a, 1)) a]
            end
        end
        nks = length(ks)
        nrow, ncol = size(a)
        @assert nks == ncol > 2
        for i = 1:nrow
            metadata = Dict{Symbol, String}()
            for j = 1:ncol - 2
                metadata[ks[j]] = a[i, j]
            end
            comment = a[i, ncol - 1]
            repetition = parse(Int, a[i, ncol])
            r = Repetition(Run(metadata, comment), repetition)
            @assert !(r in rs)
            push!(rs, r)
        end
    end
    return rs
end

function loadAssociation(folder::String)::Association
    ts = loadPOIs(folder)
    rs = loadRuns(folder)
    filescsv = joinpath(folder, "log", "associations.csv")
    as = Set{Tuple{POI, Repetition}}()
    if isfile(filescsv)
        a, ks = readcsv(filescsv, Int, header = true, comments = false)
        nrow, ncol = size(a)
        @assert ncol == 2
        for i = 1:nrow
            push!(as, (ts[a[i,1]], rs[a[i, 2]]))
        end
    end
    return Association(ts, rs, as)
end

# empty

function empty!(a::Association)
    empty!(a.pois)
    empty!(a.runs)
    empty!(a.associations)
    return a
end

isempty(a::Association) = isempty(a.pois) && isempty(a.runs) && isempty(a.associations)

function delete_empty_metadata!(a::Association)
    for k in keys(a.runs[1].run.metadata)
        if all(isempty(r.run.metadata[k]) for r in a.runs)
            for r in a.runs
                delete!(r.run.metadata, k)
            end
        end
    end
end


text:
using Documenter, Chess, Chess.Book, Chess.PGN, Chess.UCI

makedocs(;
    modules = [Chess, Chess.Book, Chess.PGN, Chess.UCI],
    format = Documenter.HTML(assets = String[]),
    pages = [
        "Home" => "index.md",
        "Tutorial" => "manual.md",
        "API Reference" => "api.md",
        "Index" => "api-index.md",
    ],
    repo = "https://github.com/romstad/Chess.jl/blob/{commit}{path}#L{line}",
    sitename = "Chess.jl",
    authors = "Tord Romstad",
)

deploydocs(repo = "github.com/romstad/Chess.jl.git")

ext: jl


text:
"""
$(SIGNATURES)

Build `EnsembleProblem` object for different open-system models. For `:lindblad` and `:ame`, it build the ensemble for the quantum trajectories method. For `:stochastic`, it builds the ensemble for the stochastic Schrödinger equation. For `:redfield`, it builds the ensemble to infuse stochastic nosie into the Redfield equation.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total annealing time.
- `type::Symbol`: type of the ensemble to build. Available options are `:lindblad`, `:ame`, `:stochastic` and `:redfield`.
- `tspan = (0, tf)`: time interval to solve the dynamics.
- `initializer = DEFAULT_INITIALIZER`: initializer for the ensemble problem. Currently it is only supported by the `:stochastic` ensemble.
- `output_func`: The function determines what is saved from the solution to the output array. Defaults to saving the solution itself. The output is (out,rerun) where out is the output and rerun is a boolean which designates whether to rerun. It is part of the `DifferentialEquations`'s parallel interface.
- `reduction`: This function determines how to reduce the data in each batch. Defaults to appending the data from the batches. The second part of the output determines whether the simulation has converged. If true, the simulation will exit early. By default, this is always false. It is part of the `DifferentialEquations`'s parallel interface.
- `kwargs`: other keyword arguments supported by the specific solver or by `DifferentialEquations`.
...
"""
function build_ensembles(
    A::Annealing,
    tf::Real,
    type::Symbol;
    tspan = (0.0, tf),
    output_func = (sol, i) -> (sol, false),
    reduction = (u, data, I) -> (append!(u, data), false),
    initializer = DEFAULT_INITIALIZER,
    kwargs...,
)
    if type == :stochastic
        res = build_ensemble_stochastic(
            A,
            tf,
            output_func,
            reduction;
            tspan = tspan,
            initializer = initializer,
            kwargs...,
        )
    elseif type == :ame
        res = build_ensemble_ame(
            A,
            tf,
            output_func,
            reduction;
            tspan = tspan,
            initializer = initializer,
            kwargs...,
        )
    elseif type == :lindblad
        res = build_ensemble_lindblad(
            A,
            tf,
            output_func,
            reduction;
            tspan = tspan,
            kwargs...,
        )
    else
        error("Ensemble problem of type $type is not supported.")
    end
    res
end

function build_ensembles(
    A::Annealing,
    tf::Real,
    unitary,
    type::Symbol;
    vectorize::Bool = false,
    output_func = (sol, i) -> (sol, false),
    reduction = (u, data, I) -> (append!(u, data), false),
    int_atol = 1e-8,
    int_rtol = 1e-6,
    Ta = tf,
    initializer = DEFAULT_INITIALIZER,
    kwargs...,
)
    if type == :redfield
        build_ensemble_redfield(
            A,
            tf,
            unitary,
            output_func,
            reduction;
            vectorize = vectorize,
            int_atol = int_atol,
            int_rtol = int_rtol,
            Ta = Ta,
            initializer = DEFAULT_INITIALIZER,
            kwargs...,
        )
    else
        error("Ensemble problem of type $type is not implemented.")
    end
end

DEFAULT_PROB_FUNC(prob, i, repeat) = (prob)

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
#=
s2Cmds.jl
    Written by Matthew White
    11/16/2019
=#
using DelimitedFiles

function saveMat(name,mat)
    # Save the output matrix
    exts = ["dat","txt"] # valid file extensions
    while true
        if name == ""
            name = ask("File Name: ")
        end
        if lowercase(name) == "exit" || lowercase(name) == ""
            break
        end
        parts = split(name,".")
        if occursin("/",name) # can change in the future to look for available folders
            errorHandler(name,"$name must not be inside a folder.")

        elseif length(parts) > 2 # extra periods
            errorHandler(name,"$name is an invalid file name.")

        elseif sum(parts[end] .== exts) == 0 # not a valid extention
            errorHandler(name,"$(parts[end]) is an invalid extension.")

        else
            writedlm(name,mat)
            break
        end
        name = ""
    end
end

function runSim(commandDict,surf,curfield,nsteps,delvort,sim)
    # Run the specificed simulation, returns mat
    sim = lowercase(sim)
    simList = ["lve","ldvm"]
    while true
        if sim == "lve"
            println()
            println(" Running LVE Simulation...")
            writeInterval = floor(nsteps/commandDict.s2["wcount"])
            frames, IFR_field, mat = LVE(surf,curfield,nsteps,commandDict.s1["dt"],commandDict.s2["wcount"],"UI Window",1000,commandDict.s2["anim"],30,commandDict.s2["wflag"],writeInterval)
            println("\t Done!")
            println()
            return mat
        elseif sim == "ldvm"
            println()
            println(" Running LDVM Simulation...")
            writeInterval = commandDict.s1["runtime"]/commandDict.s2["wcount"]
            commandDict.s2["wflag"] ? wflag = 1 : wflag = 0
            commandDict.s2["sflag"] ? sflag = 1 : sflag = 0
            mat, surf, curfield = ldvm(surf, curfield, nsteps, commandDict.s1["dt"],sflag, wflag, writeInterval, delvort)
            println("\t Done!")
            println()
            return mat
        end

        sim = ask("Simulation Type: ")
        if sim == ""
            break
        elseif sum(simList .== sim) == 0 # if not on simulation list
            errorHandler(simList,"lsit")
        end
    end
end

function plotResults(mat,subCmd)
    # Plot anything from sim results
    ploty = []
    plotCnt = []
    plotx = []
    matCol = ["t","alpha","h","u","lesp","cl","cd","cm"]
    matIdx = DataStructures.SortedDict{String,Any}(matCol[i] => i for i = 1:length(matCol))
    if haskey(matIdx,lowercase(subCmd)) # subcmd is plot type (infers only one plot to be calculated)
        plotCnt = 1
        ploty = [subCmd]
    else
        plotCnt = tryparse(Int64,subCmd)
        if plotCnt == nothing # subcmd is number of plots
            plotCnt = []
        end
    end
    # Gather plot info
    while true
        if plotCnt == [] # number of plots
            plotCnt = svPrompt("i","Number of plots: ")
            if plotCnt == nothing
                return nothing
            end
        end
        if plotCnt > 7 # too many plots
            plotCnt = []
            errorHandler(plotCnt,"$plotCnt exceeds the maximum number of plots.")
            continue
        else
            break
        end
    end
    if ploty == [] # determine y axis(es)
        for i = 1:plotCnt
            while true
                ploty = [ploty;svPrompt("s","Plot $i y-axis: ")]
                if ploty[i] == nothing
                    return nothing
                end
                if !haskey(matIdx,lowercase(ploty[i])) # not a vaild y axis
                    errorHandler(matCol,"list")
                    ploty = ploty[1:end-1] # reset ploty
                    continue
                else
                    break
                end
            end
        end
    end
    while true # determine x axis
        plotx = svPrompt("s","X axis: ")
        if plotx == nothing
            return nothing
        end
        if !haskey(matIdx,lowercase(plotx)) # not a vaild x axis
            errorHandler(matCol,"list")
            continue
        else
            break
        end
    end
    plotMat(mat,plotCnt,plotx,ploty,matIdx)

    if ynAsk("Save figure?")
        fn = svPrompt("s","File name: ")
        if occursin(".",fn)
            temp = split(fn,".") # pyplot only supports .png
            fn = temp[1]
        end
        fn = "$fn.png"
        savefig(fn)
    end
end

function plotTime(tt  = "")
    # Asks for a specific time to plot step file results
    ttparse = tryparse(Float64, tt)
    if tt == "all"
        return tt
    elseif tt == "" || ttparse == nothing
        if ynAsk("Plot all step files?")
            return "all"
        else
            return svPrompt("f","Time to plot: ")
        end
    else
        return ttparse
    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

table_header = [html_thead(html_tr([html_th("First Name"), html_th("Last Name")]))];

row1 = html_tr([html_td("Arthur"), html_td("Dent")]);
row2 = html_tr([html_td("Ford"), html_td("Prefect")]);
row3 = html_tr([html_td("Zaphod"), html_td("Beeblebrox")]);
row4 = html_tr([html_td("Trillian"), html_td("Astra")]);

table_body = [html_tbody([row1, row2, row3, row4])];

table = dbc_table([table_header; table_body], bordered = true);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
@testset "Utils" begin
    @testset "Caching" begin
        m = Model()
        @variable(m, x)
        c = @constraint(m, x <= 8.)
        @test Lagrangian.sense(c) == :le
        @test Lagrangian.getslack(c) == AffExpr(x - 8)
        m.linconstr[c.idx].lb = 8.
        @test Lagrangian.sense(c) == :eq
        @test Lagrangian.getslack(c) == AffExpr(x - 8)
        m.linconstr[c.idx].ub = Inf
        @test Lagrangian.sense(c) == :ge
        @test Lagrangian.getslack(c) == AffExpr(x - 8)
        m.linconstr[c.idx].ub = 10.
        @test_throws Exception Lagrangian.sense(c)
        @test Lagrangian.getslack(c) == AffExpr(x - 10)
    end
    @testset "Satisfied" begin
        m = Model()
        x = @variable(m, [i=1:2])
        m.colVal[x[1].col] = 1.
        m.colVal[x[2].col] = 2.
        @test Lagrangian.issatisfied( 1e-9, :le, 1e-8) == true
        @test Lagrangian.issatisfied(-1e-9, :le, 1e-8) == true
        @test Lagrangian.issatisfied(  -1., :le, 1e-8) == true
        @test Lagrangian.issatisfied(   1., :le, 1e-8) == false
        @test Lagrangian.issatisfied( 1e-9, :eq, 1e-8) == true
        @test Lagrangian.issatisfied(-1e-9, :eq, 1e-8) == true
        @test Lagrangian.issatisfied(  -1., :eq, 1e-8) == false
        @test Lagrangian.issatisfied(   1., :eq, 1e-8) == false
        @test Lagrangian.issatisfied( 1e-9, :ge, 1e-8) == true
        @test Lagrangian.issatisfied(-1e-9, :ge, 1e-8) == true
        @test Lagrangian.issatisfied(  -1., :ge, 1e-8) == false
        @test Lagrangian.issatisfied(   1., :ge, 1e-8) == true
        @test_throws Exception Lagrangian.issatisfied(AffExpr(), :x)
    end
    @testset "LevelMethod" begin
        @test_throws Exception levelmethod = LevelMethod(initialbound=0.0, level=2, quadsolver=GLPKSolverMIP())
        # Need to specify quadsolver
        @test_throws Exception levelmethod = LevelMethod(initialbound=0.0)
    end
    @testset "SubgradientDescent" begin
        @testset "CorrectSign" begin
            pi = [1.; -1.; 1.; -1.; 0.]
            Lagrangian.correctsign!(pi, [:le; :ge; :ge; :le; :anything])
            @test isapprox(norm(pi - [1.; -1.; -1.; 1.; 0.]), 0., atol=1e-9)
        end
    end
    @testset "Tolerance" begin
        tol = Absolute(0.1)
        @test isclose(100.0, 101.0, tol) == false
        @test isclose(0.01, 0.012, tol)  == true
        tol = Relative(0.1)
        @test isclose(100.0, 101.0, tol) == true
        @test isclose(-100.0, -109.9, tol) == true
        @test isclose(0.01, 0.012, tol)  == false
        @test isclose(1e-9, 1.1e-9, tol) == true
        tol = Unit(0.1)
        @test isclose(100.0, 101.0, tol) == true
        @test isclose(0.01, 0.012, tol)  == true
        @test isclose(0.01, 1.012, tol)  == false
    end
end

@testset "Lagrangian" begin
    @testset "Relax and Recover" begin
        m=Model(solver=GLPKSolverMIP())
        @variable(m, -1 <= x[1:2] <=2)
        @constraints(m, begin
            c1, x[1] <= 0.
            c2, x[1] >= 0.
            c3, x[2] == 0.
        end)
        @objective(m, :Max, sum(x))
        solve(m)

        subgradient = LevelMethod(initialbound=-5.0, quadsolver=IpoptSolver(print_level=0))
        relaxedbounds = [10.; -10.; 10.]
        l  = LinearProgramData(m.obj,
            [c1; c2; c3],
            relaxedbounds,
            method=subgradient)
        @testset "Relax" begin
            Lagrangian.relaxandcache!(l, m)
            @test l.old_bound ≈ zeros(3)
            @test m.linconstr[c1.idx].lb == -Inf
            @test m.linconstr[c1.idx].ub ≈ 10.
            @test m.linconstr[c2.idx].lb ≈ -10.
            @test m.linconstr[c2.idx].ub == Inf
            @test m.linconstr[c3.idx].lb == -Inf
            @test m.linconstr[c3.idx].ub ≈ 10.
        end
        @testset "Recover" begin
            Lagrangian.recover!(l, m, -ones(3))
            @test m.linconstr[c1.idx].lb == -Inf
            @test m.linconstr[c1.idx].ub ≈ 0.
            @test m.linconstr[c2.idx].lb ≈ 0.
            @test m.linconstr[c2.idx].ub == Inf
            @test m.linconstr[c3.idx].lb == 0.
            @test m.linconstr[c3.idx].ub ≈ 0.
            @test getdual.([c1; c2; c3]) ≈ ones(3)
        end
    end
end

const examples_dir = joinpath(dirname(dirname(@__FILE__)), "examples")

@testset "Lagrangian Examples" begin
    for example in [
            "small_MIP.jl",
            "rcsp.jl"
        ]
        @testset "$example" begin
            include(joinpath(examples_dir, "Lagrangian", example))
        end
    end
end


text:
using Test
using MARY_fRG.Drawers, MARY_fRG.Triangulated, MARY_fRG.Basics

@testset "六边型的切割" begin
    hexa1 = EqHexagon(Point2D(0., 0.), 2)
    ltris, ladjs = split_hexagon(hexa1, 10)
    #plt = draw_figure(
    #    nothing, hexa1.edges, nothing
    #)
    #@savepng plt "sh1"
    #
    tri2idx = Dict(((tri, idx) for (idx, tri) in enumerate(ltris)))
    #pts = vcat([ltris[1].center; ladjs[1][2].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[1]], tri2idx[ladjs[1][2]]])
    #draw_lines!(plt, [ltris[1].edges[2]])
    #
    #pts = vcat([ltris[22].center; ladjs[22][2].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[22]], tri2idx[ladjs[22][2]]])
    #draw_lines!(plt, [ltris[22].edges[2]])
    #
    #pts = vcat([ltris[21].center; ladjs[21][1].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[21]], tri2idx[ladjs[21][1]]])
    #draw_lines!(plt, [ltris[21].edges[1]])
    #
    #pts = vcat([ltris[44].center; ladjs[44][1].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[44]], tri2idx[ladjs[44][1]]])
    #draw_lines!(plt, [ltris[44].edges[1]])
    #
    #
    #pts = vcat([ltris[56].center; ladjs[56][1].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[56]], tri2idx[ladjs[56][1]]])
    #draw_lines!(plt, [ltris[56].edges[1]])
    @test tri2idx[ladjs[56][1]] == 57
    #
    #
    #pts = vcat([ltris[556].center; ladjs[556][1].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[556]], tri2idx[ladjs[556][1]]])
    #draw_lines!(plt, [ltris[556].edges[1]])
    @test ismissing(ladjs[556][2])
    #
    #pts = vcat([ltris[600].center; ladjs[600][1].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[600]], tri2idx[ladjs[600][1]]])
    #draw_lines!(plt, [ltris[600].edges[1]])
    @test ismissing(ladjs[600][2])
    #
    #pts = vcat([ltris[580].center; ladjs[580][2].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[580]], tri2idx[ladjs[580][2]]])
    #draw_lines!(plt, [ltris[580].edges[2]])
    @test ismissing(ladjs[580][1])
    #
    #
    #pts = vcat([ltris[532].center; ladjs[532][2].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[532]], tri2idx[ladjs[532][2]]])
    #draw_lines!(plt, [ltris[532].edges[2]])
    @test ismissing(ladjs[532][1])
    #
    #pts = vcat([ltris[521].center; ladjs[521][2].center])
    #draw_points!(plt, pts;
    #text=[tri2idx[ltris[521]], tri2idx[ladjs[521][2]]])
    #draw_lines!(plt, [ltris[521].edges[2]])
    @test tri2idx[ladjs[521][2]] == 522
    #
    #@savepng plt "sh2"
    #
end


text:
module JTetra

importall jInv.Mesh
using jInv.Utils

export TetraMesh

type TetraMesh <: AbstractMesh

	Points::Array{Float64}      # Points of tesselation
  Tetras::Array{Int64}        # Connectivity map

	nc::Int  # number of cells
	nf::Int  # number of faces
	ne::Int  # number of edges
	nn::Int  # number of nodes

	Div::SparseMatrixCSC
	Curl::SparseMatrixCSC
	Grad::SparseMatrixCSC

	V::Array{Float64,1} # cell volume
	F::Array{Float64,1} # face area
	L::Array{Float64,1} # edge lengths

	faces::Array # FaceNumbering
	edges::Array # EdgeNumbering

end # type TetraMesh

include("getTetraMesh.jl")
include("getTetraMeshNumbering.jl")
include("findRow.jl")
include("getNodalMassMatrix.jl")
include("getEdgeMassMatrix.jl")
include("getFaceMassMatrix.jl")
include("getNodalGradientMatrix.jl")
include("getCurlMatrix.jl")
include("getDivergenceMatrix.jl")
include("getSizes.jl")
include("getEdgeIntegralOfPolygonalChain.jl")

import Base.clear!
function clear!(M::TetraMesh)
	M.Div  = clear!(M.Div)
	M.Curl = clear!(M.Curl)
	M.Grad = clear!(M.Grad)
end # function clear

export clear!

end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
using DiffEqNoiseProcess, SciMLBase
W = WienerProcess(0.0,1.0,1.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
sol = solve(prob; dt=0.1)
sol = solve(ensemble_prob, EnsembleSerial(); trajectories=100, dt = 0.1)


text:
using XMLRPC
using Base.Test

# Use the Odoo demo instance for testing

v = XMLRPCProxy("http://demo.odoo.com/start")
v.url == "http://demo.odoo.com/start"

res = v["start"]()

@test typeof(res) == Dict{Any, Any}

url = res["host"]
pw = res["password"]
db = res["database"]
un = res["user"]

common = XMLRPCProxy(url*"/xmlrpc/2/common")

uid = common["authenticate"](db, un, pw, [])

models = XMLRPCProxy(url*"/xmlrpc/2/object")

b = models["execute_kw"](db, uid, pw,
    "res.partner", "check_access_rights",
    ["read"], Dict("raise_exception"=> false))

@test b == true

c = models["execute_kw"](db, uid, pw,
    "res.partner", "search",
    Any[Any[Any["is_company", "=", true], Any["customer", "=", true]]])

@test typeof(c) == Vector{Any}
@test typeof(c[1]) == Int32


text:
# A* shortest-path search with two generalizations over a_star_spath:
# 1. The graph is implicit; edges are not explicitly stored but instead generated
# during vertex expansion by a user-defined visitor function
# 2. The include_vertex! visitor function can also ADD new vertices, i.e. the size
# of the graph can increase during the search
#
# We do not bother with deleting vertices as it is typically not worth the effort
# for implicit search

"""
Data structure for heap entry associated with each expanded vertex.

Attributes:
    - `v_idx::Int64` The integer index of the vertex being expanded
    - `gvalue::D` The cost-to-come of the vertex from the source
    - `fvalue::D` The cost-to-come + heuristic cost-to-go to the goal
"""
struct AStarHEntry{D <: Number}
    v_idx::Int64
    gvalue::D
    fvalue::D
end

Base.isless(e1::AStarHEntry, e2::AStarHEntry) = (e1.fvalue, e1.gvalue) < (e2.fvalue, e2.gvalue) # Lexicographic ordering

"""
Compared with a_star_spath, the attributes of AStarStates are all Dicts instead of
fixed-length Vectors, as arbitrarily many vertices can be added by the visitors.
"""
@with_kw mutable struct AStarStates{D<:Number}
    parent_indices::Dict{Int64,Int64} = Dict{Int64,Int64}()
    dists::Dict{Int64,D} = Dict{Int64,D}()
    colormap::Dict{Int64,Int64} = Dict{Int64,Int64}()
    heap::MutableBinaryMinHeap{AStarHEntry{D}} = MutableBinaryMinHeap{AStarHEntry{D}}()
    hmap::Dict{Int64,Int64} = Dict{Int64,Int64}()
end

function set_source!(state::AStarStates{D}, s::Int64) where {D <: Number, V}
    state.parent_indices[s] = s
    state.dists[s] = zero(D)
    state.colormap[s] = 2
end

"""
Execute expand operation on the chosen node, while implicitly generating its neighbors based on a
visitor method, and populating `neighbors` with the neighboring vertices.
"""
function process_neighbors_implicit!(
    state::AStarStates{D},
    graph::AbstractGraph{V},
    edge_wt_fn::Function,
    neighbors::Vector{Int64},
    u::Int64, du::D, visitor::AbstractDijkstraVisitor,
    heuristic::Function) where {V, D <: Number}

    dv = zero(D)

    for iv in neighbors

        # Default color 0
        v_color = get(state.colormap, iv, 0)

        if v_color == 0

            # Inserting for the first time
            state.dists[iv] = dv = du + edge_wt_fn(graph.vertices[u], graph.vertices[iv])
            state.parent_indices[iv] = u
            state.colormap[iv] = 1
            Graphs.discover_vertex!(visitor, graph.vertices[u], graph.vertices[iv], dv)
            state.hmap[iv] = push!(state.heap, AStarHEntry(iv, dv, dv + heuristic(graph.vertices[iv])))

        elseif v_color == 1

            dv = du + edge_wt_fn(graph.vertices[u], graph.vertices[iv])

            # Update cost-to-come if cheaper from current parent
            if dv < state.dists[iv]

                state.dists[iv] = dv
                state.parent_indices[iv] = u

                Graphs.update_vertex!(visitor, graph.vertices[u], graph.vertices[iv], dv)
                update!(state.heap, state.hmap[iv], AStarHEntry(iv, dv, dv + heuristic(graph.vertices[iv])))
            end
        end
    end
end

"""
Run A* search on the given list-of-indices graph. The graph visitor must generate the neighbors of the expanded
vertex just-in-time, based on the user-defined external logic.

Attributes:
    - `graph::AbstractGraph{V} where {V}`
    - `edge_wt_fn::F1 where {F1 <: Function}` Maps (u,v) to the edge weight of the u -> v edge
    - `source::Int64`
    - `visitor::AbstractDijkstraVisitor` The graph visitor that implements include_vertex!, which is called when
                                         a vertex is expanded
    - `heuristic::F2 where {F2 <: Function}` Maps u to h(u); an admissible heuristic for the cost-to-go
    - `::Type{D} = Float64` The type of the edge weight value

Returns:
    - `::AStarStates` The result of the A* search
"""
function a_star_implicit_shortest_path!(
    graph::AbstractGraph{V},                # the graph
    edge_wt_fn::F1, # distances associated with edges
    source::Int64,             # the source index
    visitor::AbstractDijkstraVisitor,# visitor object
    heuristic::F2,      # Heuristic function for vertices
    ::Type{D} = Float64) where {V, D <: Number, F1 <: Function, F2 <: Function}

    state = AStarStates{D}()
    set_source!(state, source)
    source_entry = AStarHEntry(source, zero(D), heuristic(graph.vertices[source]))
    state.hmap[source] = push!(state.heap, source_entry)

    while ~(isempty(state.heap))

        # pick next vertex to include
        entry = pop!(state.heap)
        ui = entry.v_idx
        du = entry.gvalue

        state.colormap[ui] = 2

        # Will be populated by include_vertex!
        nbrs = Vector{Int64}(undef, 0)

        if Graphs.include_vertex!(visitor, graph.vertices[state.parent_indices[ui]], graph.vertices[ui], du, nbrs) == false
            return state
        end

        # process u's neighbors
        process_neighbors_implicit!(state, graph, edge_wt_fn, nbrs, ui, du, visitor, heuristic)
        Graphs.close_vertex!(visitor, graph.vertices[ui])
    end

    return state
end

"""
Extract the sequence of vertex indices on the shortest path using the parent pointer dictionary.

Arguments:
    - `parent_indices::Dict{Int64, Int64}` The dictionary mapping vertex indices to their parents.
    - `graph::AbstractGraph{V}` The graph (as a vector of indices)
    - `source_idx::Int64`
    - `target_idx::Int64`

Returns:
    - `::Vector{Int64}` The list of indices on the shortest path in order.
"""
function shortest_path_indices(parent_indices::Dict{Int64, Int64}, graph::AbstractGraph{V},
                               source_idx::Int64, target_idx::Int64) where {V}

    @assert haskey(parent_indices, target_idx) == true "Target has no parent!"

    sp = [target_idx]

    # Walk back from target to source using parent indices
    curr_idx = target_idx
    while curr_idx != source_idx
        curr_idx = parent_indices[curr_idx]
        pushfirst!(sp, curr_idx)
    end

    return sp
end


text:
function initial_data_inner_loop!(u0, f0, nodes, method, mesh, i)
    if method == :average
        u0[:,i] .= num_integrate(f0,nodes[i], nodes[i+1])/cell_volume(mesh,i)
    elseif method == :eval_in_centers
        u0[:,i] .= f0((nodes[i]+nodes[i+1])/2)
    else
        error("invalid initial data processing method: ", method)
    end
    nothing
end

function getInitialState(mesh, f0; method =:average, use_threads = false, MType = Float64)
    N = getncells(mesh)
    numvars = size(f0(getnodecoords(mesh, 1)[1]),1)
    u0 = MMatrix{numvars,N,MType}(undef)
    nodes = get_nodes_matrix(mesh)
    if use_threads
        Threads.@threads for i in 1:getncells(mesh)
            initial_data_inner_loop!(u0, f0, nodes, method, mesh, i)
        end
    else
        for i in 1:getncells(mesh)
            initial_data_inner_loop!(u0, f0, nodes, method, mesh, i)
        end
    end
    return numvars > 1 ? u0 : u0[1,:]
end

function getSemiDiscretization(f,alg::AbstractFVAlgorithm,
    mesh,dbcs; Df = nothing, numvars = 1,
    Type = Float64, use_threads = false, dt = 0.0)
    probtype = numvars > 1 ? GeneralProblem() : ScalarProblem()
    internal_mesh = mesh_setup(mesh, dbcs,probtype)
    fluxes = numvars > 1 ? MMatrix{numvars,getnnodes(mesh),Type}(undef) : MVector{getnnodes(mesh),Type}(undef)
    Flux = numvars > 1 ? flux_function(f, Df) : scalar_flux_function(f, Df)
    FVIntegrator(alg,internal_mesh,Flux,numvars, fluxes, dt, use_threads)
end

function fv_solution(sol::AbstractODESolution{T,N}, mesh; vars::Int = 1) where {T,N}
    FVSolution{T,N,typeof(mesh), vars}(sol, mesh)
  end


text:
using FileIO, CSV, ImageAxes, Compose

"""
    readSEManticsImage(fn::AbstractString)

Read a SEMantics PNG image and the create an `ImageAxes` (`AxisArray`) which associates
`:x`, `:y` coordinates with the image axes according to the X, Y and FOV in the
images.txt file in the same directory.
"""
function readSEManticsImage(fn::AbstractString)
    img, sp = load(fn), splitpath(fn)
    img_txt = joinpath(sp[1:end-1]..., "images.txt")
    if isfile(img_txt)
        imgs = CSV.File(img_txt, delim = "\t") |> DataFrame
        (nm, _) = splitext(sp[end])
        rg = Regex(replace(nm[1:end-3], "_" => "[.\\\\:]", "["=>"\\[", "]"=>"\\]"))
        mimgs = filter(r -> !isnothing(match(rg, r[:Name])), imgs)
        if nrow(mimgs) > 0
            # Since an image can be written more than once, take the last...
            x, y, fov = mimgs[end, :X], mimgs[end, :Y], mimgs[end, :FOV]
            pix = fov / size(img, 2) # square pixels determined by horizontal FOV
            img = AxisArray(
                img, #
                Axis{:y}(StepRangeLen(y, pix, size(img,1), size(img,1)÷2)),
                Axis{:x}(StepRangeLen(x, -pix, size(img,2), size(img,2)÷2))
            )

        else
            @info "There is no meta-data in image.txt for images matching $rg."
        end
    else
        @info "There is no meta-data file images.txt for images matching $rg."
    end
    return img
end

function _ns_closest(srl::AbstractRange, val)
    idx = searchsortedfirst(srl, val, rev = step(srl) < 0.0)
    ip = idx + (step(srl) < 0.0 ? -1 : 1)
    return (val < min(srl[1], srl[end])) || (val > max(srl[1], srl[end])) ? nothing : #
           (abs(srl[idx] - val) < abs(srl[ip] - val) ? idx : ip)
end

function _ns_closest(img::AxisArray, coord)
    return (
        _ns_closest(ImageAxes.axisvalues(img)[2], coord[1]),
        _ns_closest(ImageAxes.axisvalues(img)[1], coord[2]),
    )
end

function build_scale_bar(img::AxisArray; marker = HSVA(60, 1, 1, 0.3))
    ax = AxisArrays.axisvalues(img)[1]
    fov = abs(ax[end] - ax[1]) / 2.0
    sc = 10.0^floor(log10(fov))
    l = (fov / sc >= 5.0 ? 5.0 : (fov / sc >= 2.5 ? 2.5 : 1.0)) * sc
    dp = l / abs(ax[end] - ax[1])
    xo, yo = 0.02, 0.02
    return Tuple[
        (context(), Compose.line([(1.0 - 2.0 * xo - dp, yo + 0.07), (1.0 - 2.0 * xo, yo + 0.07)]), stroke(marker), linewidth(1.0)),
        (context(), Compose.text(1.0 - 2.0 * xo - dp / 2.0, yo + 0.045, "$(round(Int,1000.0*l)) μm", hcenter, vbottom), fill(marker), stroke("transparent"), fontsize(8.0)),
        (context(), Compose.rectangle(1.0 - 3.0 * xo - dp, yo, dp + 2.0 * xo, yo + 0.07), stroke("transparent"), fill(RGBA(0.1, 0.1, 0.1, 0.3)))
    ]
end

function build_acquisition_points(
    img::AxisArray,
    coords::AbstractArray{<:AbstractDict{Symbol, <:AbstractFloat}};
    label = true,
    marker = HSVA(60, 1, 1, 0.3)
)
    x, y = Float64[], Float64[]
    for coord in coords
        pt = _ns_closest(img, (coord[:X], coord[:Y]))
        if !any(isnothing.(pt))
            push!(x, pt[1] / size(img, 1))
            push!(y, pt[2] / size(img, 2))
        end
    end
    #  Place circles at acquisition point
    items = Tuple[(context(), Compose.circle(x, y, [0.005]), fill(marker), stroke("transparent"))]
    if label
        fx(x) = x < 0.9 ? x + 0.01 : x - 0.01
        fh(x) = x < 0.9 ? hleft : hright
        push!(items, (context(), Compose.text(fx.(x), y, repr.(eachindex(coords)), fh.(x), [ vcenter]), fill(marker), stroke("transparent")))
    end
    return items
end

function build_thumbnails(
    img::AxisArray,
    thumbnails::Union{Nothing, AbstractArray{<:AxisArray}}=nothing;
    label = true,
    marker = HSVA(60, 1, 1, 0.3)
)
    ya, xa = axisvalues(img)[1:2]
    rx, ry, rw, rh = Float64[], Float64[], Float64[], Float64[]
    for thumb in thumbnails
        yt, xt = axisvalues(thumb)[1:2]
        y1o = (yt[1]-ya[1])/(ya[end]-ya[1])
        yeo = (yt[end]-ya[1])/(ya[end]-ya[1])
        x1o = (xt[1]-xa[1])/(xa[end]-xa[1])
        xeo = (xt[end]-xa[1])/(xa[end]-xa[1])
        if (max(x1o,xeo)>0.0) && (min(x1o,xeo)<1.0) &&
            (max(y1o,yeo)>0.0) && (min(y1o,yeo)<1.0)
            push!(rx, x1o)
            push!(ry, y1o)
            push!(rw, xeo-x1o)
            push!(rh, yeo-y1o)
        end
    end
    #  Place rectangles at image points
    items = Tuple[ (context(), Compose.rectangle(rx, ry, rw, rh), fill(RGBA(0.1, 0.1, 0.1, 0.3)), stroke(marker)) ]
    if label
        push!(items, (context(), Compose.text(rx .+ 0.001, ry .+ 0.005, repr.(eachindex(thumbnails)), [ hleft ], [ vtop ]), fill(marker), stroke("transparent")))
    end
    return items
end

"""
    annotate(
        img::AxisArray;
        scale::Bool = true,
        coords::Union{Nothing, AbstractArray{<:AbstractDict{Symbol, <:AbstractFloat}}} = nothing,
        spectra::Union{Nothing, <:AbstractArray{<:Spectrum}} = nothing,
        thumbnails::Union{Nothing, AbstractArray{<:AxisArray}}=nothing,
        labelcoords::Bool = true,
        labelthumbnails::Bool = true,
        marker::Colorant = HSVA(60, 1, 1, 0.3)
    )

Add annotations like scale-bars, acquisition points, image areas to an image
and display the result.

  * `coords` and `spectra` display as circles (with/without numeric labels).
  * `thumbnails` display as rectangles overlaying the image
  * `scale` displays as a scale-bar in the upper-right corner
  * `marker` is the color of the annotations which is by-default slightly transparent
"""
function annotate(
    img::AxisArray;
    scale::Bool = true,
    coords::Union{Nothing, AbstractArray{<:AbstractDict{Symbol, <:AbstractFloat}}} = nothing,
    spectra::Union{Nothing, <:AbstractArray{<:Spectrum}} = nothing,
    thumbnails::Union{Nothing, AbstractArray{<:AxisArray}}=nothing,
    labelcoords::Bool = true,
    labelthumbnails::Bool = true,
    marker::Colorant = HSVA(60, 1, 1, 0.3)
)
    items = Tuple[]
    if scale
        append!(items, build_scale_bar(img, marker=marker))
    end
    allcoords = AbstractDict{Symbol, <:AbstractFloat}[]
    if !isnothing(coords)
        append!(allcoords, coords)
    end
    if !isnothing(spectra)
        foreach(sp->push!(allcoords, sp[:StagePosition]), spectra)
    end
    if length(allcoords)>0
        append!(items, build_acquisition_points(img, allcoords; label=labelcoords, marker=marker))
    end
    if !isnothing(thumbnails)
        append!(items, build_thumbnails(img, thumbnails; label=labelthumbnails, marker=marker))
    end
    io = IOBuffer()
    ImageIO.save(Stream{format"PNG"}(io), img)
    return compose(
        context(0.0, 0.0, 1.0, 1.0),
        items...,
        ( context(), Compose.bitmap("image/png", take!(io), 0.0, 0.0, 1.0, 1.0) ),
    )
end

"""
    shannon_entropy(img::AbstractArray{Gray{N0f8}})

Computes the log2-entropy of the data in the image.
The entropy(...) in Images.jl 24.1 is buggy and is removed in 25.0
"""
function shannon_entropy(img::AbstractArray{Gray{N0f8}})
    # c = @MVector zeros(Int(typemax(UInt8))+1) # Reduces allocations to zero..
    c = zeros(Int(typemax(UInt8))+1)
    for v in img
        @inbounds c[reinterpret(UInt8, v.val)+1] += 1.0/length(img)
    end
    l2(v) = v==0 ? 0.0 : log2(v)
    return -sum(v->v*l2(v), c)
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
#! /usr/bin/env julia

# Copyright 2017 John Hanley.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# The software is provided "AS IS", without warranty of any kind, express or
# implied, including but not limited to the warranties of merchantability,
# fitness for a particular purpose and noninfringement. In no event shall
# the authors or copyright holders be liable for any claim, damages or
# other liability, whether in an action of contract, tort or otherwise,
# arising from, out of or in connection with the software or the use or
# other dealings in the software.

# usage:
#     src/install_modules.jl

# Pkg.update()
Pkg.add("SQLite")
run(`sudo apt install -y zlib1g-dev`)


text:
include("init.jl")

@testset "Broker" begin
    broker = createBroker(1)
    initstring = "--broker="
    identifier = h.helicsBrokerGetIdentifier(broker)
    initstring = initstring * identifier
    initstring = initstring * " --broker_address "
    address = h.helicsBrokerGetAddress(broker)
    initstring = initstring * address
    @test initstring == "--broker=mainbroker --broker_address tcp://127.0.0.1:23404"
    destroyBroker(broker)

end

@testset "MessageFilter Registration" begin

    broker = createBroker(2)

    fFed, ffedinfo = createMessageFederate(1, "filter")
    mFed, mfedinfo = createMessageFederate(1, "message")

    h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter1")
    f2 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter2")
    h.helicsFederateRegisterEndpoint(fFed, "fout", "")
    h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM,  "filter0/fout")
    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    filter_name = h.helicsFilterGetName(f1)
    @test filter_name == "Testfilter/filter1"

    filter_name = h.helicsFilterGetName(f2)
    @test filter_name == "Testfilter/filter2"

    # filter_target = h.helicsFilterGetTarget(f2)
    # @test filter_target == "port2"

    h.helicsFederateFinalize(mFed)
    h.helicsFederateFinalize(fFed)

    destroyFederate(fFed, ffedinfo)
    destroyFederate(mFed, mfedinfo)
    sleep(1.0)

    destroyBroker(broker)
end

@testset "MessageFilter Info" begin

    broker = createBroker(2)

    fFed, ffedinfo = createMessageFederate(1, "filter")
    mFed, mfedinfo = createMessageFederate(1, "message")

    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    h.helicsEndpointSetInfo(p1, "p1_test")
    h.helicsEndpointSetInfo(p2, "p2_test")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    h.helicsFilterSetInfo(f1, "f1_test")

    f2 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter2")
    h.helicsFilterAddSourceTarget(f2, "port2")
    h.helicsFilterSetInfo(f2, "f2_test")

    ep1 = h.helicsFederateRegisterEndpoint(fFed, "fout", "")
    h.helicsEndpointSetInfo(ep1, "ep1_test")

    f3 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM,  "c4")
    h.helicsFilterAddSourceTarget(f3, "filter0/fout")
    h.helicsFilterSetInfo(f3, "f3_test")

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    filter_name = h.helicsFilterGetName(f1)
    @test filter_name == "Testfilter/filter1"

    filter_name = h.helicsFilterGetName(f2)
    @test filter_name == "Testfilter/filter2"

    @test h.helicsEndpointGetInfo(p1) == "p1_test"
    @test h.helicsEndpointGetInfo(p2) == "p2_test"
    @test h.helicsEndpointGetInfo(ep1) == "ep1_test"

    @test h.helicsFilterGetInfo(f1) == "f1_test"
    @test h.helicsFilterGetInfo(f2) == "f2_test"
    @test h.helicsFilterGetInfo(f3) == "f3_test"

    h.helicsFederateFinalize(mFed)
    h.helicsFederateFinalize(fFed)

    destroyFederate(fFed, ffedinfo)
    destroyFederate(mFed, mfedinfo)
    sleep(1.0)

    destroyBroker(broker)
end

@testset "MessageFilter Function" begin
    broker = createBroker(2)

    fFed, ffedinfo = createMessageFederate(1, "filter")
    mFed, mfedinfo = createMessageFederate(1, "message")

    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "random")

    f1 = h.helicsFederateRegisterGlobalFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    f2 = h.helicsFederateRegisterGlobalFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter2")
    h.helicsFilterAddSourceTarget(f2, "port1")
    h.helicsFederateRegisterEndpoint(fFed,"fout","")
    f3 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_RANDOM_DELAY, "filter3")
    h.helicsFilterAddSourceTarget(f3,"filter/fout")

    h.helicsFilterSet(f2, "delay", 2.5)
    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)
    state = h.helicsFederateGetState(fFed)
    @test state == 2
    data = "hello world"

    filt_key = h.helicsFilterGetName(f1)
    @test filt_key == "filter1"

    filt_key = h.helicsFilterGetName(f2)
    @test filt_key == "filter2"

    h.helicsEndpointSendMessageRaw(p1, "port2", data)
    h.helicsFederateRequestTimeAsync(mFed, 1.0)
    grantedtime = h.helicsFederateRequestTime(fFed, 1.0)
    @test grantedtime == 1.0
    grantedtime = h.helicsFederateRequestTimeComplete(mFed)
    @test grantedtime == 1.0
    res = h.helicsFederateHasMessage(mFed)
    @test res == 0
    res = h.helicsEndpointHasMessage(p2)
    @test res == 0
    #grantedtime = h.helicsFederateRequestTime(fFed, 3.0)
    #@test res==h.helics_true

    h.helicsFederateFinalize(mFed)
    h.helicsFederateFinalize(fFed)
    #f2 = h.helicsFederateRegisterDestinationFilter(fFed, h.helics_custom_filter, "filter2", "port2")
    #ep1 = h.helicsFederateRegisterEndpoint(fFed, "fout", "")
    #f3 = h.helicsFederateRegisterSourceFilter(fFed, h.helics_custom_filter, "", "filter0/fout")

    destroyFederate(fFed, ffedinfo)
    destroyFederate(mFed, mfedinfo)
    sleep(1.0)
    destroyBroker(broker)

end


text:
# NOTE: Add staggered grids
struct Grid
    xmax::Float64
    ymax::Float64

    nx::Int64
    ny::Int64

    Δx::Float64
    Δy::Float64

    xRange::Array{Float64,1}
    yRange::Array{Float64,1}

    function Grid(xmax, ymax, nx, ny)
        Δx = 2*xmax/(nx-1)
        Δy = 2*ymax/(ny-1)

        xRange = (-xmax:Δx:xmax)
        yRange = (-ymax:Δy:ymax)

        return new(xmax, ymax, nx, ny, Δx, Δy, xRange, yRange)
    end
end


text:
module Configuration

Settings = Dict(
    :api_access_control_allow_origin => "*",
    :dev => false,
    :html_base_filename => "html/base",
    :html_favicon_name => "html/favicon",
    :log_filename => "log/dancejl",
    :routes_filename => "routes",
    :server_host => "",
    :server_port => ""
)

"""
    is_prod(file_path)

Check if is Production environment
"""
function is_prod() :: Bool
    global Settings
    return !parse(Bool, Settings[:dev])
end

"""
    populate(file_path)

- Read settings/Global.jl from project
- Populate/update Settings dict accordingly
"""
function populate(file_path) :: Bool
    is_success:Bool = false

    ### Parse project settings ###
    function parse() :: Nothing
        function _parse_file(path) :: Nothing
            open(path) do file
                for line in eachline(file)
                    if length(line) > 0
                        # Trim left space
                        line = lstrip(line)

                        if startswith(line, "include(")
                            _parse_file(joinpath(
                                file_path, "settings",
                                replace(replace(line, "include(\"" => ""), "\")" => "")
                            ))
                        elseif line[1]==':'
                            line_split::Array{String} = split(line, "=")
                            key::String = lstrip(rstrip(line_split[1]), ':')
                            value::Union{Float64, Int64, String, SubString{String}} = lstrip(line_split[2])

                            # Remove trailing comment
                            if occursin(" #", value)
                                value = rstrip(split(value, " #")[1])
                            end

                            # Remove enclosing quotation marks if is String value
                            if value[1]=='"'
                                value = lstrip(rstrip(value, '"'), '"')
                            end

                            # Convert SubString{String} from to String
                            value = string(value)

                            # Convert to Float/Int if possible
                            if tryparse(Float64, value) isa Number
                                if tryparse(Int64, value) !== nothing
                                    value = tryparse(Int64, value)
                                else
                                    value = parse(Float64, value)
                                end
                            end

                            # Populate/update Settings dict
                            if haskey(Settings, Symbol(key))
                                global Settings
                                Settings[Symbol(key)] = value
                            else
                                global Settings
                                push!(Settings, Symbol(key) => value)
                            end
                        end
                    end
                end
            end
        end

        return _parse_file(joinpath(file_path, "settings/Global.jl"))
    end

    parse()
    is_success = true

    ### Set default server values for Dev environment ###
    if !is_prod()
        if Settings[:server_host]==""
            global Settings
            Settings[:server_host] = "127.0.0.1"
        end
        if Settings[:server_port] ==""
            global Settings
            Settings[:server_port] = 8000
        end
    end

    ### Server values cannot be blank for Prod environment, as no default ###
    for key in [:server_host, :server_port]
        if Settings[key] == ""
            @error "Please set valid $key value, cannot be blank"
            is_success = false
        end
    end

    return is_success
end

end


text:
using CUDA, Distributed
ngpu=length(devices())
addprocs(ngpu-1; exeflags="--project")

@everywhere using CUDA

# assign devices
asyncmap((zip(workers(), collect(devices())[2:ngpu]))) do (p, d)
    remotecall_wait(p) do
        @info "Worker $p uses $d"
        device!(d)
    end
end


text:
# reversing

# the kernel works by treating the array as 1d. after reversing by dimension x an element at
# pos [i1, i2, i3, ... , i{x},            ..., i{n}] will be at
# pos [i1, i2, i3, ... , d{x} - i{x} + 1, ..., i{n}] where d{x} is the size of dimension x

# out-of-place version, copying a single value per thread from input to output
function _reverse(input::AnyCuArray{T, N}, output::AnyCuArray{T, N};
                  dims::Integer=1) where {T, N}
    @assert size(input) == size(output)
    shape = [size(input)...]
    numelemsinprevdims = prod(shape[1:dims-1])
    numelemsincurrdim = shape[dims]

    function kernel(input::AbstractArray{T, N}, output::AbstractArray{T, N}) where {T, N}
        offset_in = blockDim().x * (blockIdx().x - 1)

        index_in = offset_in + threadIdx().x

        if index_in <= length(input)
            element = @inbounds input[index_in]

            # the index of an element in the original array along dimension that we will flip
            #assume(numelemsinprevdims > 0)
            #assume(numelemsincurrdim > 0)
            ik = ((cld(index_in, numelemsinprevdims) - 1) % numelemsincurrdim) + 1

            index_out = index_in + (numelemsincurrdim - 2ik + 1) * numelemsinprevdims

            @inbounds output[index_out] = element
        end

        return
    end

    nthreads = 256
    nblocks = cld(prod(shape), nthreads)
    shmem = nthreads * sizeof(T)

    @cuda threads=nthreads blocks=nblocks kernel(input, output)
end

# in-place version, swapping two elements on half the number of threads
function _reverse(data::AnyCuArray{T, N}; dims::Integer=1) where {T, N}
    shape = [size(data)...]
    numelemsinprevdims = prod(shape[1:dims-1])
    numelemsincurrdim = shape[dims]

    function kernel(data::AbstractArray{T, N}) where {T, N}
        offset_in = blockDim().x * (blockIdx().x - 1)

        index_in = offset_in + threadIdx().x

        # the index of an element in the original array along dimension that we will flip
        #assume(numelemsinprevdims > 0)
        #assume(numelemsincurrdim > 0)
        ik = ((cld(index_in, numelemsinprevdims) - 1) % numelemsincurrdim) + 1

        index_out = index_in + (numelemsincurrdim - 2ik + 1) * numelemsinprevdims

        if index_in <= length(data) && index_in < index_out
            @inbounds begin
                temp = data[index_out]
                data[index_out] = data[index_in]
                data[index_in] = temp
            end
        end

        return
    end

    # NOTE: we launch twice the number of threads, which is wasteful, but the ND index
    #       calculations don't allow using only the first half of the threads
    #       (e.g. [1 2 3; 4 5 6] where threads 1 and 2 swap respectively (1,2) and (2,1)).

    nthreads = 256
    nblocks = cld(prod(shape), nthreads)
    shmem = nthreads * sizeof(T)

    @cuda threads=nthreads blocks=nblocks kernel(data)
end

# n-dimensional API

# in-place
function Base.reverse!(data::AnyCuArray{T, N}; dims::Integer) where {T, N}
    if !(1 ≤ dims ≤ length(size(data)))
      ArgumentError("dimension $dims is not 1 ≤ $dims ≤ $length(size(input))")
    end

    _reverse(data; dims=dims)

    return data
end

# out-of-place
function Base.reverse(input::AnyCuArray{T, N}; dims::Integer) where {T, N}
    if !(1 ≤ dims ≤ length(size(input)))
      ArgumentError("dimension $dims is not 1 ≤ $dims ≤ $length(size(input))")
    end

    output = similar(input)
    _reverse(input, output; dims=dims)

    return output
end

# 1-dimensional API

# in-place
Base.@propagate_inbounds function Base.reverse!(data::AnyCuVector{T}, start, stop=length(data)) where {T}
    _reverse(view(data, start:stop))
    return data
end

Base.reverse(data::AnyCuVector{T}) where {T} = @inbounds reverse(data, 1, length(data))

# out-of-place
Base.@propagate_inbounds function Base.reverse(input::AnyCuVector{T}, start, stop=length(input)) where {T}
    output = similar(input)

    start > 1 && copyto!(output, 1, input, 1, start-1)
    _reverse(view(input, start:stop), view(output, start:stop))
    stop < length(input) && copyto!(output, stop+1, input, stop+1)

    return output
end

Base.reverse!(data::AnyCuVector{T}) where {T} = @inbounds reverse!(data, 1, length(data))


text:
# SimpleMeshFunction2.jl
# Created for the new simpleMesh type

@doc """
### createMesh

It creates a simple mesh with triangle elements in a rectangular domain in 2D. In This code The horizontal axis is referred to as X-axis and the vertical axis is referred to as the Y-axis.

**Inputs**

*  `lengthx` : Length on the X-axis
*  `lengthy` : Length on the Y-axis
*  `nedx`    : Number of element edges along X-axis
*  `nedy`    : Number of element edges along Y-axis
*  `nnpe`    : Number of nodes per element edge (same for all edges)
*  `numDofPerNode` : Number of degrees of freedom per node

**Outputs**

*  `xnodes` : Number of nodes along X-axis in one row
*  `ynodes` : Number of nodes along Y-axis in one column
*  `nnp`    : Number of nodal points
*  `nel`    : Number of elements mesh
*  `dofs`   : Element information matrix

We assume a counter-clockwise numbering scheme for the the elements and edges with the numbering starting from the base of the geometry.
"""->

function createMesh(lengthx,lengthy,nedx,nedy,nnpe,numDofPerNode)
  xnodes = (nnpe-1)*nedx + 1 # Total number of nodes on all adjacent element edges along X-axis
  ynodes = (nnpe-1)*nedy + 1 # Total number of nodes on all adjacent element edges along Y-axis
  nel = 2*nedx*nedy # number of elements
  nsd = 2  # number of spatial dimensions

  if nnpe == 2
    nen = 3 # Number of element nodes
    dofs = Array{Int}(numDofPerNode, nen, nel)
    nnp = xnodes*ynodes # number of nodal points
    m = 0
    for i = 1:2:nel
      for j = 1:numDofPerNode
        dofs[j,1,i] = m + j
      end
      #
      for j = 1:numDofPerNode
        dofs[j,2,i] = m + numDofPerNode + j
      end
      #m += numDofPerNode*xnodes
      for j = 1:numDofPerNode
        dofs[j,3,i] = m + numDofPerNode*xnodes + j
      end
      m += numDofPerNode
      for k = 1:nedy
        if i == 2*nedx*k - 1
          m +=numDofPerNode
        end
      end
    end
    m = numDofPerNode
    for i = 2:2:nel
      for j = 1:numDofPerNode
        dofs[j,1,i] = m + numDofPerNode*xnodes + j
      end
      for j = 1:numDofPerNode
        dofs[j,2,i] = m + numDofPerNode*(xnodes-1) + j
      end
      for j = 1:numDofPerNode
        dofs[j,3,i] = m + j
      end
      m += numDofPerNode
      for k = 1:nedy
        if i == 2*nedx*k
          m +=numDofPerNode
        end
      end
    end

    return dofs, nnp, nel, nen
  elseif nnpe == 3
    nen = 7 # Number of element nodes
    nnp = xnodes*ynodes + nel;
    dofs = Array{Int}(numDofPerNode, nen, nel)
    m = 0
    for j = 1:nedy
      for i = 1:2:2*nedx
        index = i+2*(j-1)*nedx
        for k = 1:numDofPerNode
          dofs[k,:,index] = [m+k;
                             m+2*numDofPerNode + k;
                             m+2*xnodes*numDofPerNode + k;
                             m+numDofPerNode + k;
                             m+(xnodes+1)*numDofPerNode + k;
                             m+xnodes*numDofPerNode + k;0];
        end
        m = m+2*numDofPerNode;
        if index == 2*nedx*j - 1
          m = m + (xnodes+1)*numDofPerNode;
        end
      end
    end
    m = 2*numDofPerNode;
    for j = 1:nedy
      for i = 2:2:2*nedx
        index = i + 2*(j-1)*nedx;
        for k = 1:numDofPerNode
          dofs[k,:,index] = [m+2*xnodes*numDofPerNode+k;
                             m+(2*xnodes-2)*numDofPerNode+k;
                             m+k;
                             m+((2*xnodes)-1)*numDofPerNode+k;
                             m+(xnodes-1)*numDofPerNode+k;
                             m+xnodes*numDofPerNode+k;0];
        end
        m = m+2*numDofPerNode;
        if index == 2*nedx*j
          m = m + (xnodes+1)*numDofPerNode;
        end
      end
    end
    m = xnodes*ynodes*numDofPerNode;
    for j = 1:nedy
      for i = 1:2*nedx
        index = i + 2*(j-1)*nedx;
        for k = 1:numDofPerNode
          dofs[k,7,index] = m + k;
        end
        m = m+numDofPerNode;
      end
    end

    return dofs, nnp, nel, nen
  elseif nnpe == 4
    nen = 12 # Number of element nodes
    nnp = xnodes*(nedy+1) + (nnpe-2)*nedy*(2*nedx+1) + 3*nel;
    dofs = Array{Int}(numDofPerNode, nen, nel)
    m = 0
    nmen = 2*nedx +1; # Number of midpoint element nodes in a row
    for j = 1:nedy
      k = 2;
      for i = 1:2:2*nedx
        index = i+2*(j-1)*nedx;
        intv = (j-1)*(xnodes+2*nmen);     # Intermediate variable
        for l = 1:numDofPerNode
          dofs[l,:,index] = [m+l;
                             m+(nnpe-1)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen)*numDofPerNode+l;
                             m+numDofPerNode+l;
                             m+2*numDofPerNode+l;
                             (xnodes+k+intv-1)*numDofPerNode+l;
                             (xnodes+nmen+k+intv-1)*numDofPerNode+l;
                             (xnodes+(k-1)+intv-1)*numDofPerNode+l;
                             (xnodes+nmen+(k-1)+intv-1)*numDofPerNode+l;0;0;0];
        end
        m = m+(nnpe-1)*numDofPerNode;
        k = k+2;
        if index == 2*nedx*j - 1
          m = m+((nnpe-2)*(2*nedx+1)+1)*numDofPerNode;
        end
      end
    end
    m = 3*numDofPerNode;
    for j = 1:nedy
      k = 2;
      for i = 2:2:2*nedx
        index = i + 2*(j-1)*nedx;
        intv = (j-1)*(xnodes+2*nmen);     # Intermediate variable
        for l = 1:numDofPerNode
          dofs[l,:,index] = [m+(xnodes+(nnpe-2)*nmen)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen-3)*numDofPerNode+l;
                             m+l;
                             m+(xnodes+(nnpe-2)*nmen-1)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen-2)*numDofPerNode+l;
                             (xnodes+nmen+k+intv-1)*numDofPerNode+l;
                             (xnodes+k+intv-1)*numDofPerNode+l;
                             (xnodes+nmen+k+intv)*numDofPerNode+l;
                             (xnodes+k+intv)*numDofPerNode+l;0;0;0];
        end
        m = m+(nnpe-1)*numDofPerNode;
        k = k+2;
        if index == 2*nedx*j
          m = m+((nnpe-2)*(2*nedx+1)+1)*numDofPerNode;
        end
      end
    end

    m = (xnodes*(nedy+1) + (nnpe-2)*nedy*(2*nedx+1))*numDofPerNode;
    for j = 1:nedy
      for i = 1:2*nedx
        index = i + 2*(j-1)*nedx;
        for l = 1:numDofPerNode
          dofs[l,10,index] = m+l;
          dofs[l,11,index] = m+numDofPerNode+l;
          dofs[l,12,index] = m+2*numDofPerNode+l;
        end
        m =  m+3*numDofPerNode;
      end
    end

  return dofs, nnp, nel, nen
 elseif nnpe == 5
    nen = 18 # Number of element nodes
    nnp = xnodes*(nedy+1) + (nnpe-2)*nedy*(2*nedx+1) + 6*nel;
    dofs = Array{Int}(numDofPerNode, nen, nel)
    m = 0
    nmen = 2*nedx +1; # Number of midpoint element nodes in a row
    for j = 1:nedy
      k = 2;
      for i = 1:2:2*nedx
        index = i + 2*(j-1)*nedx;
        intv = (j-1)*(xnodes+(nnpe-2)*nmen);     # Intermediate variable
        for l = 1:numDofPerNode
          dofs[l,:,index] = [m+l;
                             m+(nnpe-1)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen)*numDofPerNode+l;
                             m+numDofPerNode+l;
                             m+2*numDofPerNode+l;
                             m+3*numDofPerNode+l;
                             (xnodes+k+intv-1)*numDofPerNode+l;
                             (xnodes+nmen+k+intv-1)*numDofPerNode+l;
                             (xnodes+(2*nmen)+k+intv-1)*numDofPerNode+l;
                             (xnodes+k-2+intv)*numDofPerNode+l;
                             (xnodes+nmen+k-2+intv)*numDofPerNode+l;
                             (xnodes+(2*nmen)+k-2+intv)*numDofPerNode+l;
                             0;0;0;0;0;0];
        end
        m = m+(nnpe-1)*numDofPerNode;
        k = k+2;
        if index == 2*nedx*j - 1
          m = m+((nnpe-2)*(2*nedx+1)+1)*numDofPerNode;
        end
      end
    end
    m = 4*numDofPerNode;
    for j = 1:nedy
      k = 2;
      for i = 2:2:2*nedx
        index = i + 2*(j-1)*nedx;
        intv = (j-1)*(xnodes+(nnpe-2)*nmen);     # Intermediate variable
        for l = 1:numDofPerNode
          dofs[l,:,index] = [m+(xnodes+(nnpe-2)*nmen)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen-4)*numDofPerNode+l;
                             m+l;
                             m+(xnodes+(nnpe-2)*nmen-1)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen-2)*numDofPerNode+l;
                             m+(xnodes+(nnpe-2)*nmen-3)*numDofPerNode+l;
                             (xnodes+(2*nmen)+k+intv-1)*numDofPerNode+l;
                             (xnodes+nmen+k+intv-1)*numDofPerNode+l;
                             (xnodes+k+intv-1)*numDofPerNode+l;
                             (xnodes+(2*nmen)+k+intv)*numDofPerNode+l;
                             (xnodes+nmen+k+intv)*numDofPerNode+l;
                             (xnodes+k+intv)*numDofPerNode+l;
                             0;0;0;0;0;0];
        end
        m = m+(nnpe-1)*numDofPerNode;
        k = k+2;
        if index == 2*nedx*j
          m = m+((nnpe-2)*(2*nedx+1)+1)*numDofPerNode;
        end
      end
    end

    m = (xnodes*(nedy+1) + (nnpe-2)*nedy*(2*nedx+1))*numDofPerNode;
    for j = 1:nedy
      for i = 1:2*nedx
        index = i + 2*(j-1)*nedx;
        for l = 1:numDofPerNode
          dofs[l,13,index] = m+l;
          dofs[l,14,index] = m+2*numDofPerNode+l;
          dofs[l,15,index] = m+4*numDofPerNode+l;
          dofs[l,16,index] = m+3*numDofPerNode+l;
          dofs[l,17,index] = m+5*numDofPerNode+l;
          dofs[l,18,index] = m+numDofPerNode+l;
        end
        m = m+6*numDofPerNode;
      end
    end
    return dofs, nnp, nel, nen
  end # ends the if-else statement
end # Ends the function

@doc """
### vtxLocation

It creates 3D array storing the X,Y & Z coordinates of the vertices of every element in addition to calculating the element edge length along the X and Y direction

**Inputs**

*  `lengthx` : Length on the X-axis
*  `lengthy` : Length on the Y-axis
*  `nedx`    : Number of element edges along X-axis
*  `nedy`    : Number of element edges along Y-axis
*  `nel`     : Number of elements in the mesh

**Outputs**

*  `vtxCoord` : 3D array storing the X,Y & Z coordinates of the vertices. First index stores the coordinates, 2nd index has the local vertex numbering, 3rd index is the element number
*  `elem_edge_lengthx : Length of element edge along the X-axis
*  `elem_edge_lengthy : Length of element edge along the Y-axis

"""->

function vtxLocation(lengthx, lengthy, nedx, nedy, nel)
  elemEdgeLengthx = lengthx/nedx;
  elemEdgeLengthy = lengthy/nedy;
  vtxCoord = zeros(Float64, 3, 3, nel) # Dim1: xyz coord, Dim2: Local vtx numbering, Dim3 Element no
  xvtx = nedx + 1 # Total number of nodes on all adjacent element edges along X-axis
  yvtx = nedy + 1 # Total number of nodes on all adjacent element edges along Y-axis
  m = 0
  n = 0
  for j = 1:nedy
    for i = 1:2:2*nedx
      index = i + 2*(j-1)*nedx;
      vtxCoord[1,1,index] = m*elemEdgeLengthx
      vtxCoord[2,1,index] = n*elemEdgeLengthy
      vtxCoord[1,2,index] = (m+1)*elemEdgeLengthx
      vtxCoord[2,2,index] = n*elemEdgeLengthy
      vtxCoord[1,3,index] = m*elemEdgeLengthx
      vtxCoord[2,3,index] = (n+1)*elemEdgeLengthy
      m += 1
      if index == 2*nedx*j - 1
        m = 0
        n += 1
      end
    end
  end
  m = 1
  n = 1
  for j = 1:nedy
    for i = 2:2:2*nedx
      index = i + 2*(j-1)*nedx;
      vtxCoord[1,1,index] = m*elemEdgeLengthx
      vtxCoord[2,1,index] = n*elemEdgeLengthy
      vtxCoord[1,2,index] = (m-1)*elemEdgeLengthx
      vtxCoord[2,2,index] = n*elemEdgeLengthy
      vtxCoord[1,3,index] = m*elemEdgeLengthx
      vtxCoord[2,3,index] = (n-1)*elemEdgeLengthy
      m += 1
      if index == 2*nedx*j
        m = 1
        n += 1
      end
    end
  end

  return vtxCoord, elemEdgeLengthx, elemEdgeLengthy
end # Ends the function

@doc """
### boundaryNodeInfo

It generates information about nodes on the boundary for implementing boundary conditions.

**Inputs**

*  `nedx`   : Number of element edges along X-axis
*  `nedy`   : Number of element edges along Y-axis
*  `nnpe`   : Number of nodes per element edge (same for all edges)

**Outputs**

*  `NodeEdgex` : A 2*xnodes matrix that stores nodes on the boundaries along the X-axis. First row stores the nodes on edge 1 and the scond row stores nodes on edge 3 of the entire geometry. The nodes are stored from left to right.
*  `NodeEdgey` : A 2*ynodes matrix that stores nodes on the boundaries along the Y-axis.First row stores nodes on edge 2 and second row stores nodes on edge 4 of the entire geometry. The nodes are stored from bottom to top.

"""->

function boundaryNodeInfo(nedx,nedy,nnpe)
  if nnpe == 2
    xnodes::Int = nedx + 1 # Total number of nodes on all adjacent element edges along X-axis
    ynodes::Int = nedy + 1 # Total number of nodes on all adjacent element edges along Y-axis
    NodeEdgex = zeros(2,xnodes)
    NodeEdgey = zeros(2,ynodes)
    for i = 1:xnodes
      NodeEdgex[1,i] = i
      NodeEdgex[2,i] = xnodes*nedy + i
    end
    for j = 1:ynodes
      NodeEdgey[1,j] = j*xnodes
      NodeEdgey[2,j] = 1+(j-1)*xnodes
    end
    return NodeEdgex, NodeEdgey

  elseif nnpe == 3
    xnodes = (nedx*2)+1; # Total number of nodes on all adjacent element edges along X-axis
    ynodes = (nedy*2)+1; # Total number of nodes on all adjacent element edges along Y-axis
    NodeEdgex = zeros(2,xnodes)
    NodeEdgey = zeros(2,ynodes)
    for i = 1:xnodes
      NodeEdgex[1,i] = i;
      NodeEdgex[2,i] = xnodes*nedy + (nnpe-2)*(2*nedx+1)*nedy + i;
    end
    for j = 1:ynodes
      NodeEdgey[1,j] = j*xnodes
      NodeEdgey[2,j] = 1+(j-1)*xnodes
    end
    return NodeEdgex, NodeEdgey

  elseif nnpe == 4
    xnodes = (nnpe-1)*nedx + 1; # Total number of nodes on all adjacent element edges along X-axis
    ynodes = (nnpe-1)*nedy + 1; # Total number of nodes on all adjacent element edges along Y-axis
    NodeEdgex = zeros(2,xnodes)
    NodeEdgey = zeros(2,ynodes)
    intv = 2*nedx + 1 # Intermediate variable which stores the value of number of midpoint nodes in a row
    for i = 1:xnodes
      NodeEdgex[1,i] = i;
      NodeEdgex[2,i] = xnodes*nedy + (nnpe-2)*(intv)*nedy + i;
    end
    m = 1;
    for j = 1:nedy
      NodeEdgey[1,m] = xnodes + (j-1)*(xnodes+(nnpe-2)*(intv));
      NodeEdgey[1,m+1] = (xnodes + intv) + (j-1)*(xnodes+(nnpe-2)*(intv));
      NodeEdgey[1,m+2] = (xnodes + 2*intv) + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m] = 1 + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m+1] = 1 + xnodes + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m+2] = 1 + xnodes + intv  + (j-1)*(xnodes+(nnpe-2)*intv);
      if j == nedy
        NodeEdgey[1,m+3] = 2*(xnodes+intv) + (j-1)*(xnodes+(nnpe-2)*intv);
        NodeEdgey[2,m+3] = 1 + xnodes + 2*intv + (j-1)*(xnodes+(nnpe-2)*intv);
      end
      m += 3;
    end
    return NodeEdgex, NodeEdgey

  elseif nnpe == 5
    xnodes = (nnpe-1)*nedx + 1; # Total number of nodes on all adjacent element edges along X-axis
    ynodes = (nnpe-1)*nedy + 1; # Total number of nodes on all adjacent element edges along Y-axis
    NodeEdgex = zeros(2,xnodes)
    NodeEdgey = zeros(2,ynodes)
    intv = 2*nedx + 1
    for i = 1:xnodes
      NodeEdgex[1,i] = i;
      NodeEdgex[2,i] = xnodes*nedy + (nnpe-2)*(2*nedx+1)*nedy + i;
    end
    m = 1;
    for j = 1:nedy
      NodeEdgey[1,m] = xnodes + (j-1)*(xnodes+(nnpe-2)*(intv));
      NodeEdgey[1,m+1] = (xnodes + intv) + (j-1)*(xnodes+(nnpe-2)*(intv));
      NodeEdgey[1,m+2] = (xnodes + 2*intv) + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[1,m+3] = (xnodes + 3*intv) + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m] = 1 + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m+1] = 1 + xnodes + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m+2] = 1 + xnodes + intv  + (j-1)*(xnodes+(nnpe-2)*intv);
      NodeEdgey[2,m+3] = 1 + xnodes + 2*intv + (j-1)*(xnodes+(nnpe-2)*intv);
      if j == nedy
        NodeEdgey[1,m+4] = 2*xnodes + 3*intv + (j-1)*(xnodes+(nnpe-2)*intv);
        NodeEdgey[2,m+4] = 1 + xnodes + 3*intv + (j-1)*(xnodes+(nnpe-2)*intv);
      end
      m += 4;
    end
    return NodeEdgex, NodeEdgey
  end # Ends the if statement
end # Ends function boundaryNodeInfo

@doc """
### boundaryEdges

This function identifies the edges and elements on the oundary. It is of the type Boundary

**Inputs**

*  `nedx`   : Number of element edges along X-axis
*  `nedy`   : Number of element edges along Y-axis
*  `nel`    : Number of elements in the mesh

**Outputs**

*  `boundaries` : Array containing the edges and elements on the boundary
*  `nedges`     : Total number of edges in the mesh

"""->

function boundaryEdges(nedx,nedy,nel)
  nedges = nedx*(nedy+1) + nedx*nedy + (nedx+1)*nedy;
  boundaries = Array{Boundary}(2*(nedx+nedy))
  m = 1
  for i = 1:nedx
    boundaries[i] = Boundary(m,1)
    m += 2
  end
  m = nel - 2*nedx+2
  for i = 1:nedx
    boundaries[nedx+i] = Boundary(m,1)
    m += 2
  end
  m = 1
  for i = 1:nedy
    boundaries[2*nedx+i] = Boundary(m,3)
    m += 2*nedx
  end
  m = 2*nedx
  for i = 1:nedy
    boundaries[2*nedx+nedy+i] = Boundary(m,3)
    m += 2*nedx
  end

  # println(boundaries)
  return boundaries, nedges
end

@doc """
### interiorEdgeinfo

This function identifies all the element interfaces (edges in 2D and faces in 3D) in the grid.
It is of the type Interface.

**Inputs**

*  `nedx`   : Number of element edges along X-axis
*  `nedy`   : Number of element edges along Y-axis
*  `nel`    : Number of elements in the mesh

**Outputs**

* `interfaces` : Array of type Interface containing the interior edge info

"""->

function interiorEdgeinfo(nedx,nedy,nel)
  nedges = nedx*(nedy+1) + nedx*nedy + (nedx+1)*nedy;
  nIntFace = 0
  interfaces = Array{Interface}(nIntFace)
  for j = 1:nedy
    for i = 1:2:2*nedx
      index = i + 2*(j-1)*nedx;
      if index == 1 + 2*(j-1)*nedx
        if index == 1
          nIntFace += 1
          interfaces = resize!(interfaces,nIntFace)
          interfaces[nIntFace] = Interface(index,index+1,2,2)
        else
          nIntFace += 2
          interfaces = resize!(interfaces,nIntFace)
          interfaces[nIntFace-1] = Interface(index,index-2*nedx+1,1,1)
          interfaces[nIntFace] = Interface(index,index+1,2,2)
        end
      elseif index < 2*nedx
        nIntFace += 2
        interfaces = resize!(interfaces,nIntFace)
        interfaces[nIntFace-1] = Interface(index,index-1,3,3)
        interfaces[nIntFace] = Interface(index,index+1,2,2)
      else
        nIntFace += 3
        interfaces = resize!(interfaces,nIntFace)
        interfaces[nIntFace-2] = Interface(index,index-2*nedx+1,1,1)
        interfaces[nIntFace-1] = Interface(index,index-1,3,3)
        interfaces[nIntFace] = Interface(index,index+1,2,2)
      end
    end
  end

  return interfaces
end

@doc """
### boundaryEdgeInfo

This function identifies the mesh edges on the boundary of the geometry.

**Inputs**

*  `nedx`   : Number of element edges along X-axis
*  `nedy`   : Number of element edges along Y-axis

**Outputs**

*  `HBedges` : A 2*nedx array which stores the edge numbers on the boundary along the X-axis. The first row stores the edges along the base of the rectangle. The second row stores edges on the top of the rectangle.
*  `VBedges` : A 2*nedy array which stores the edge numbers on the boundary along the Y-axis. The first row stores the vertical edges along the right side of the rectangle. The second row stores the edges on the left side of the rectangle.
*  `nedges`  : Total number of mesh edges in the geometry.

"""->

function boundaryEdgeInfo(nedx,nedy)
  # Boundary Edges
  nedges = nedx*(nedy+1) + nedx*nedy + (nedx+1)*nedy;
  HBedges = zeros(2,nedx);  # Horizontal Boundary edges
  VBedges = zeros(2,nedy);  # Vertical Boundary edges
  for i = 1:nedx
    HBedges[1,i] = i;
    HBedges[2,i] = nedx*nedy + i;
  end
  for i = 1:nedy
    VBedges[1,i] = nedx*(nedy+1) + nedx*nedy + 1 + (i-1)*(nedx+1);
    VBedges[2,i] = nedx*(nedy+1) + nedx*nedy + i*(nedx+1);
  end
  return HBedges, VBedges, nedges
end   # Ends function boundaryEdgeInfo


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 = "VOLK"
version = v"2.5.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/gnuradio/volk/releases/download/v$(version)/volk-$(version).tar.gz", "d9183b9f86a32cdbb8698cbbeb15de574962c05200ccf445c1058629073521f8"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/volk-*

if [[ ${target} == *-freebsd* ]]; then
    #this is not in 0.6.0 release, but has been added on master, can probably remove after next release
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/add-freebsd-macros.patch"
elif [[ ${target} == x86_64-w64-mingw* ]]; then
    #disable avx512 on x86_64-mingw to avoid "Error: invalid register for .seh_savexmm". This is not needed on i686-mingw for some reason?
    #copied from https://github.com/xianyi/OpenBLAS/issues/1801
    export CFLAGS="${CFLAGS} -fno-asynchronous-unwind-tables"
fi

mkdir build
cd build/

pip install mako

cmake .. \
-DCMAKE_INSTALL_PREFIX=${prefix} \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DCROSSCOMPILE_MULTILIB=true \
-DENABLE_TESTING=OFF \
-DENABLE_MODTOOL=OFF

make -j${nproc}
make install

"""

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

#volk by default builds with -D_GLIBCXX_USE_CXX11_ABI=1 see: https://github.com/gnuradio/volk/blob/2ff7b768a4b2174379a1f2e9214051677082f14b/CMakeLists.txt#L53
filter!(x -> cxxstring_abi(x) != "cxx03", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libvolk", :libvolk),
    ExecutableProduct("volk-config-info", :volk_config_info),
    ExecutableProduct("volk_profile", :volk_profile),
    ExecutableProduct("list_cpu_features", :list_cpu_features)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("boost_jll"; compat="=1.76.0")
    HostBuildDependency(PackageSpec(name="ORC_jll", uuid="fb41591b-4dee-5dae-bf56-d83afd04fbc0"))
]

# Build the tarballs, and possibly a `build.jl` as well.
#need at least gcc5 for Wincompatible-pointer-type flag for cpu_features, needs at least gcc6 to avoid std::inf errors on musl see: https://github.com/gnuradio/volk/issues/375
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using BenchmarkTools
using MarchingCubes
using PlyIO
using Test

@testset "march" begin
    mc = MarchingCubes.scenario()

    march(mc)
    @test length(mc.vertices) == 11_333
    @test length(mc.triangles) == 22_620
    MarchingCubes.output(PlyIO, mc, tempname())

    march_legacy(mc)
    @test length(mc.vertices) == 11_333
    @test length(mc.triangles) == 22_732
    MarchingCubes.output(PlyIO, mc, tempname())
end

@testset "normalize" begin
    nx, ny, nz = 10, 20, 30
    start_x, stop_x = -0.1, 0.1
    start_y, stop_y = -10, 10
    start_z, stop_z = 100, 200
    mc = MC(
        rand(nx, ny, nz) .- 0.5,
        Int;
        x = collect(Float64, range(start_x, stop_x, length = nx)),
        y = collect(Float64, range(start_y, stop_y, length = ny)),
        z = collect(Float64, range(start_z, stop_z, length = nz)),
    )
    for callable ∈ (march, march_legacy)
        callable(mc)
        @test all(map(v -> start_x ≤ v[1] ≤ stop_x, mc.vertices))
        @test all(map(v -> start_y ≤ v[2] ≤ stop_y, mc.vertices))
        @test all(map(v -> start_z ≤ v[3] ≤ stop_z, mc.vertices))
    end
end

@testset "plane" begin
    mc = MarchingCubes.scenario(case = :plane)
    bytes = @allocated march(mc)
    @test bytes == 0
    @test length(mc.vertices) == 7_038
    @test length(mc.triangles) == 13_720
end

@testset "sphere" begin
    mc = MarchingCubes.scenario(case = :sphere)
    bytes = @allocated march(mc)
    @test bytes == 0
    @test length(mc.vertices) == 792
    @test length(mc.triangles) == 1572
end

@testset "hyperboloid" begin
    mc = MarchingCubes.scenario(case = :hyperboloid)
    bytes = @allocated march(mc)
    @test bytes == 0
    @test length(mc.vertices) == 12_297
    @test length(mc.triangles) == 24_118
end

@testset "isovalue" begin
    dat = Float32[(x - 3)^2 + (y - 3)^2 + (z - 3)^2 for x ∈ 1:5, y ∈ 1:5, z ∈ 1:5]

    m1 = MC(dat)
    march(m1, 5.0)

    m2 = MC(dat .- 5.0)
    march(m2)

    @test all(m1.triangles .≈ m2.triangles)
    @test all(m1.vertices .≈ m2.vertices)
    @test all(m1.normals .≈ m2.normals)
end