text:
using Bukdu # ApplicationController Conn HTML Router CLI render routes get post
using Bukdu.HTML5.Form # change
                       # form_for label_for
                       # text_area text_input radio_button checkbox
                       # submit
import Documenter.Utilities.DOM: @tags

struct FormController <: ApplicationController
    conn::Conn
end

struct User
end

function layout(body)
    """

  <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />

Bukdu.HTML5.Form

$body

"""
end

global changeset = Changeset(User, (name="", famous=false, season="summer", intro=""))

function index(c::FormController)
    global changeset
    @tags div
    form1 = form_for(changeset, (FormController, post_result), method=post, multipart=true) do f
        div(
            div.([
                text_area(f, :name, placeholder="Name", rows="3", cols="50"),
                label_for(checkbox(f, :famous), "Famous"),
                submit("Submit"),
                " multipart/form-data",
            ])
        )
    end
    form2 = form_for(changeset, (FormController, post_result), method=post, multipart=false) do f
        div(
            div.([
                text_area(f, :intro, placeholder="Intro", rows="5", cols="50"),
                [
                    label_for(radio_button(f, :season, "winter"), "Winter"),
                    label_for(radio_button(f, :season, "spring"), "Spring"),
                    label_for(radio_button(f, :season, "summer"), "Summer"),
                    label_for(radio_button(f, :season, "autumn"), "Autumn"),
                ],
                submit("Submit"),
                " application/x-www-form-urlencoded",
            ])
        )
    end
    body = div(form1, form2)
    render(HTML, layout(body))
end

function post_result(c::FormController)
    global changeset
    @tags div a h3 li p
    result = change(changeset, c.params)
    if !isempty(result.changes)
        changeset.changes = merge(changeset.changes, result.changes)
    end
    body = div(
        h3( a[:href => "/"]("back") ),
        p(isempty(result.changes) ?
            "nothing's changed" :
            "✅ changed"),
        li( "result: ", string(result) ),
        li( "global changeset: ", string(changeset) ),
    )
    render(HTML, layout(body))
end

if PROGRAM_FILE == basename(@__FILE__)

routes() do
    get("/", FormController, index)
    post("/post_result", FormController, post_result)
end

Bukdu.start(8080)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
#!/usr/bin/env julia

println("\e[34mInitialising\e[0m")

# Function definitions
# -----------------------------------------------------------------------------

# Initialise output store
function init_out()
  return Dict{String,Union{Nothing,String,Array{String,1}}}(
    "id" => nothing,
    "text" => Array{String,1}()
  )
end

# Prepare the text for the XML output
function xml_text(content::Array{String,1})::String
  local xml_start::Array{String,1} = [
    """<?xml version="1.0" encoding="UTF-8"?>""",
    """<hmdb xmlns="http://www.hmdb.ca">"""
  ]
  local xml_end::String = "</hmdb>"
  return "$(join(xml_start, "\n"))\n$(join(content, "\n"))\n$xml_end\n"
end

# Write the data out
function write_xml(
    out::Dict{String,Union{Nothing,String,Array{String,1}}},
    dir_out::String
  )::Bool
  local id::String
  id = isnothing(out["id"]) ? "IDNA_$(lpad(string(idna), 7, "0"))" : out["id"]
  open("$dir_out/$(out["id"]).xml", "w") do file
    write(file, xml_text(out["text"]))
  end
  return true
end

# Operations
# -----------------------------------------------------------------------------

cd("../hmdb_data/")
if !isdir("./all_split")
  mkdir("./all_split")
end

# Define paths to use
path = Dict{Symbol,String}(
  :xml => "./hmdb_metabolites.xml",
  :dir_out => "./all_split"
)

# DEBUG
# global stop = 0

println("\e[34mRunning operations\e[0m")

global idna = 0
global out = init_out()

for line in eachline(path[:xml])
  # Ignore matching
  if occursin(r"^\<(\/)?(hmdb|\?xml)", line)
    continue
  end
  push!(out["text"], line)
  # While no ID, check if current row is ID
  if isnothing(out["id"])
    m = match(r"\<accession\>(?<id>HMDB\d+)\<\/accession\>", line)
    if !isnothing(m)
      out["id"] = string(m[:id])
    end
  end
  # If end of metabolite data, write out and reset
  if occursin(r"^\<\/metabolite\>", line)
    write_xml(out, path[:dir_out])
    if isnothing(out["id"])
      global idna += 1
    end
    global out = init_out()
    # DEBUG
    # global stop += 1
    # if stop > 9
    #   break
    # end
  end
end

println("\e[32mOperations complete\e[0m")

ext: jl


text:
const ST = Union{Int32, Int64, Float32, Float64}

for (A,F) in ((:besselj, :arb_hypgeom_bessel_j), (:bessely, :arb_hypgeom_bessel_y),
              (:besseli, :arb_hypgeom_bessel_i), (:besselk, :arb_hypgeom_bessel_k))
    @eval begin
        function ($A)(nu::ArbReal{P}, x::ArbReal{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong),
                                      z, nu, x, prec)
        return z
        end
    end
end

for (A,F) in ((:besselj, :acb_hypgeom_bessel_j), (:bessely, :acb_hypgeom_bessel_y),
              (:besseli, :acb_hypgeom_bessel_i), (:besselk, :acb_hypgeom_bessel_k))
    @eval begin
        function ($A)(nu::ArbComplex{P}, x::ArbComplex{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong),
                                      z, nu, x, prec)
            return z
        end
    end
end

besselj(nu::ArbFloat{P}, x::ArbFloat{P}, prec::Int=P) where {P} = ArbFloat{P}(besselj(ArbReal{P}(nu), ArbReal{P}(x)), prec)
bessely(nu::ArbFloat{P}, x::ArbFloat{P}, prec::Int=P) where {P} = ArbFloat{P}(bessely(ArbReal{P}(nu), ArbReal{P}(x)), prec)
besseli(nu::ArbFloat{P}, x::ArbFloat{P}, prec::Int=P) where {P} = ArbFloat{P}(besseli(ArbReal{P}(nu), ArbReal{P}(x)), prec)
besselk(nu::ArbFloat{P}, x::ArbFloat{P}, prec::Int=P) where {P} = ArbFloat{P}(besselk(ArbReal{P}(nu), ArbReal{P}(x)), prec)

besselj(nu::S, x::ArbFloat{P}, prec::Int=P) where {P, S<:ST} = ArbFloat{P}(besselj(ArbReal{P}(nu), ArbReal{P}(x)), prec)
bessely(nu::S, x::ArbFloat{P}, prec::Int=P) where {P, S<:ST} = ArbFloat{P}(bessely(ArbReal{P}(nu), ArbReal{P}(x)), prec)
besseli(nu::S, x::ArbFloat{P}, prec::Int=P) where {P, S<:ST} = ArbFloat{P}(besseli(ArbReal{P}(nu), ArbReal{P}(x)), prec)
besselk(nu::S, x::ArbFloat{P}, prec::Int=P) where {P, S<:ST} = ArbFloat{P}(besselk(ArbReal{P}(nu), ArbReal{P}(x)), prec)

besselj(nu::S, x::ArbReal{P}, prec::Int=P) where {P, S<:ST} = besselj(ArbReal{P}(nu), x, prec)
bessely(nu::S, x::ArbReal{P}, prec::Int=P) where {P, S<:ST} = bessely(ArbReal{P}(nu), x, prec)
besseli(nu::S, x::ArbReal{P}, prec::Int=P) where {P, S<:ST} = besseli(ArbReal{P}(nu), x, prec)
besselk(nu::S, x::ArbReal{P}, prec::Int=P) where {P, S<:ST} = besselk(ArbReal{P}(nu), x, prec)

besselj0(x::ArbFloat{P}, prec::Int=P) where {P} = besselj(zero(ArbReal{P}), ArbReal{P}(x), prec)
besselj1(x::ArbFloat{P}, prec::Int=P) where {P} = besselj(one(ArbReal{P}), ArbReal{P}(x), prec)
bessely0(x::ArbFloat{P}, prec::Int=P) where {P} = bessely(zero(ArbReal{P}), ArbReal{P}(x), prec)
bessely1(x::ArbFloat{P}, prec::Int=P) where {P} = bessely(one(ArbReal{P}), ArbReal{P}(x), prec)

besselj0(x::ArbReal{P}, prec::Int=P) where {P} = besselj(zero(ArbReal{P}), x, prec)
besselj1(x::ArbReal{P}, prec::Int=P) where {P} = besselj(one(ArbReal{P}), x, prec)
bessely0(x::ArbReal{P}, prec::Int=P) where {P} = bessely(zero(ArbReal{P}), x, prec)
bessely1(x::ArbReal{P}, prec::Int=P) where {P} = bessely(one(ArbReal{P}), x, prec)

besselj0(x::ArbComplex{P}, prec::Int=P) where {P} = besselj(zero(ArbComplex{P}), x, prec)
besselj1(x::ArbComplex{P}, prec::Int=P) where {P} = besselj(one(ArbComplex{P}), x, prec)
bessely0(x::ArbComplex{P}, prec::Int=P) where {P} = bessely(zero(ArbComplex{P}), x, prec)
bessely1(x::ArbComplex{P}, prec::Int=P) where {P} = bessely(one(ArbComplex{P}), x, prec)

besselj(nu::S, x::ArbComplex{P}, prec::Int=P) where {P, S<:ST} = besselj(ArbComplex{P}(nu), x, prec)
bessely(nu::S, x::ArbComplex{P}, prec::Int=P) where {P, S<:ST} = bessely(ArbComplex{P}(nu), x, prec)
besseli(nu::S, x::ArbComplex{P}, prec::Int=P) where {P, S<:ST} = besseli(ArbComplex{P}(nu), x, prec)
besselk(nu::S, x::ArbComplex{P}, prec::Int=P) where {P, S<:ST} = besselk(ArbComplex{P}(nu), x, prec)

function airyai(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{ArbReal}, Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Clong),
              z, C_NULL, C_NULL, C_NULL, x, prec)
    return z
end

function airyaiprime(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{ArbReal}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Clong),
              C_NULL, z, C_NULL, C_NULL, x, prec)
    return z
end

function airybi(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Ref{Cvoid}, Ref{ArbReal}, Clong),
              C_NULL, C_NULL, z, C_NULL, x, prec)
    return z
end

function airybiprime(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Ref{ArbReal}, Clong),
              C_NULL, C_NULL, C_NULL, z, x, prec)
    return z
end

for F in (:airyai, :airyaiprime, :airybi, :airybiprime)
    @eval begin
        function $F(x::ArbFloat{P}, prec::Int=P) where {P}
            r = ArbReal{P}(x)
            r = airy(r, prec)
            z = ArbFloat{P}(r)
            return r
        end
    end
end

function airyai(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{ArbComplex}, Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Clong),
              z, C_NULL, C_NULL, C_NULL, x, prec)
    return z
end

function airyaiprime(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{ArbComplex}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Clong),
              C_NULL, z, C_NULL, C_NULL, x, prec)
    return z
end

function airybi(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Ref{Cvoid}, Ref{ArbComplex}, Clong),
              C_NULL, C_NULL, z, C_NULL, x, prec)
    return z
end

function airybiprime(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Ref{ArbComplex}, Clong),
              C_NULL, C_NULL, C_NULL, z, x, prec)
    return z
end


text:
# ellipse to corpus callosum like shape
# initial and final landmark positions observed

using Revise
using BridgeLandmarks
const BL = BridgeLandmarks
using Distributions
using LinearAlgebra
using Random
using FileIO

# using RCall
# using Random
# using DataFrames
# using DelimitedFiles
# using CSV
# using StaticArrays
#using JLD2

Random.seed!(9)

workdir = @__DIR__
cd(workdir)
outdir = workdir

dat = load("data_exp-dgs.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]

p_ms = Pars_ms(δmom=0.01/n,cinit = 0.2, γinit = 2.0/√n, ρinit = 0.95)
p_ahs = Pars_ahs(δmom=0.1,cinit = 0.02, γinit = 0.2, ρinit = 0.95, stdev=0.75, db=[2.5, 2.5])

ainit = mean(norm.(diff(xobs0)))
ups = [:innov, :mala_mom]
ups = [:innov, :rmrw_mom, :parameter]

@time landmarksmatching(xobs0,xobsT; ITER=100,pars=p_ms, outdir=outdir, ainit=ainit, updatescheme=ups)
#@time landmarksmatching(xobs0,xobsT; ITER=500,pars=p_ahs, outdir=outdir, ainit=ainit, updatescheme=ups)

plotlandmarksmatching(outdir)


text:
# VCF Record
# ==========

type VCFRecord
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # indexes
    chrom::UnitRange{Int}
    pos::UnitRange{Int}
    id::Vector{UnitRange{Int}}
    ref::UnitRange{Int}
    alt::Vector{UnitRange{Int}}
    qual::UnitRange{Int}
    filter::Vector{UnitRange{Int}}
    infokey::Vector{UnitRange{Int}}
    format::Vector{UnitRange{Int}}
    genotype::Vector{Vector{UnitRange{Int}}}
end

"""
    VCFRecord()

Create an unfilled `VCFRecord` object.
"""
function VCFRecord()
    return VCFRecord(
        # data and filled
        UInt8[], 1:0,
        # chrom-alt
        1:0, 1:0, UnitRange{Int}[], 1:0, UnitRange{Int}[],
        # qual-genotype
        1:0, UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[])
end

"""
    VCFRecord(data::Vector{UInt8})

Create a `VCFRecord` object from `data` containing a VCF record.
This function verifies the format and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new `VCFRecord` object.
"""
function VCFRecord(data::Vector{UInt8})
    return convert(VCFRecord, data)
end

function Base.convert(::Type{VCFRecord}, data::Vector{UInt8})
    record = VCFRecord(
        # data and filled
        data, 1:0,
        # chrom-alt
        1:0, 1:0, UnitRange{Int}[], 1:0, UnitRange{Int}[],
        # qual-genotype
        1:0, UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[])
    index!(record)
    return record
end

"""
    VCFRecord(str::AbstractString)

Create a `VCFRecord` object from `str` containing a VCF record.
This function verifies the format and indexes fields for accessors.
"""
function VCFRecord(str::AbstractString)
    return convert(VCFRecord, str)
end

function Base.convert(::Type{VCFRecord}, str::AbstractString)
    return VCFRecord(convert(Vector{UInt8}, str))
end

function initialize!(record::VCFRecord)
    record.filled = 1:0
    record.chrom = 1:0
    record.pos = 1:0
    empty!(record.id)
    record.ref = 1:0
    empty!(record.alt)
    record.qual = 1:0
    empty!(record.filter)
    empty!(record.infokey)
    empty!(record.format)
    empty!(record.genotype)
    return record
end

function isfilled(record::VCFRecord)
    return !isempty(record.filled)
end

function datarange(record::VCFRecord)
    return record.filled
end

function Base.:(==)(record1::VCFRecord, record2::VCFRecord)
    if isfilled(record1) == isfilled(record2) == true
        r1 = datarange(record1)
        r2 = datarange(record2)
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function checkfilled(record::VCFRecord)
    if !isfilled(record)
        throw(ArgumentError("unfilled VCF record"))
    end
end

function VCFRecord(base::VCFRecord;
                   chromosome=nothing, position=nothing, identifier=nothing,
                   reference=nothing, alternate=nothing, quality=nothing,
                   filter=nothing, information=nothing, genotype=nothing)
    checkfilled(base)
    buf = IOBuffer()

    if chromosome == nothing
        write(buf, base.data[base.chrom])
    else
        print(buf, string(chromosome))
    end

    print(buf, '\t')
    if position == nothing
        write(buf, base.data[base.pos])
    else
        print(buf, convert(Int, position))
    end

    print(buf, '\t')
    if identifier == nothing
        if isempty(base.id)
            print(buf, '.')
        else
            for (i, r) in enumerate(base.id)
                if i != 1
                    print(buf, ';')
                end
                write(buf, base.data[r])
            end
        end
    else
        if !isa(identifier, Vector)
            identifier = [identifier]
        end
        if isempty(identifier)
            print(buf, '.')
        else
            for (i, x) in enumerate(identifier)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(x))
            end
        end
    end

    print(buf, '\t')
    if reference == nothing
        write(buf, base.data[base.ref])
    else
        print(buf, string(reference))
    end

    print(buf, '\t')
    if alternate == nothing
        if isempty(base.alt)
            print(buf, '.')
        else
            for (i, r) in enumerate(base.alt)
                if i != 1
                    print(buf, ';')
                end
                write(buf, base.data[r])
            end
        end
    else
        if !isa(alternate, Vector)
            alternate = [alternate]
        end
        if isempty(alternate)
            print(buf, '.')
        else
            for (i, x) in enumerate(alternate)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(x))
            end
        end
    end

    print(buf, '\t')
    if quality == nothing
        write(buf, base.data[base.qual])
    else
        print(buf, convert(Float64, quality))
    end

    print(buf, '\t')
    if filter == nothing
        if isempty(base.filter)
            print(buf, '.')
        else
            for (i, r) in enumerate(base.filter)
                if i != 1
                    print(buf, ';')
                end
                write(buf, base.data[r])
            end
        end
    else
        if !isa(filter, Vector)
            filter = [filter]
        end
        if isempty(filter)
            print(buf, '.')
        else
            for (i, x) in enumerate(filter)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(x))
            end
        end
    end

    print(buf, '\t')
    if information == nothing
        if isempty(base.infokey)
            print(buf, '.')
        else
            write(buf, base.data[first(base.infokey[1]):last(infovalrange(base, endof(base.infokey)))])
        end
    else
        if !isa(information, Associative)
            throw(ArgumentError("information must be an associative object"))
        elseif isempty(information)
            print(buf, '.')
        else
            for (i, (key, val)) in enumerate(information)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(key))
                if val != nothing
                    print(buf, '=', vcfformat(val))
                end
            end
        end
    end

    print(buf, '\t')
    if genotype == nothing
        if isempty(base.format)
            print(buf, '.')
        else
            write(buf, base.data[first(base.format[1]):last(base.format[end])])
        end
        if !isempty(base.genotype)
            for indiv in base.genotype
                print(buf, '\t')
                for (i, r) in enumerate(indiv)
                    if i != 1
                        print(buf, ':')
                    end
                    write(buf, base.data[r])
                end
            end
        end
    else
        if !isa(genotype, Vector)
            genotype = [genotype]
        end
        if isempty(genotype)
            print(buf, '.')
        else
            allkeys = String[]
            for indiv in genotype
                if !isa(indiv, Associative)
                    throw(ArgumentError("individual must be an associative"))
                end
                append!(allkeys, keys(indiv))
            end
            allkeys = sort!(unique(allkeys))
            if !isempty(allkeys)
                join(buf, allkeys, ':')
                for indiv in genotype
                    print(buf, '\t')
                    for (i, key) in enumerate(allkeys)
                        if i != 1
                            print(buf, ':')
                        end
                        print(buf, vcfformat(get(indiv, key, '.')))
                    end
                end
            end
        end
    end

    return VCFRecord(takebuf_array(buf))
end

function vcfformat(val)
    return string(val)
end

function vcfformat(val::Vector)
    return join(map(vcfformat, val), ',')
end

function Base.copy(record::VCFRecord)
    return VCFRecord(
        copy(record.data),
        record.filled,
        record.chrom,
        record.pos,
        copy(record.id),
        record.ref,
        copy(record.alt),
        record.qual,
        copy(record.filter),
        copy(record.infokey),
        copy(record.format),
        deepcopy(record.genotype))
end

function Base.write(io::IO, record::VCFRecord)
    checkfilled(record)
    return write(io, record.data)
end

# Accessor functions
# ------------------

"""
    chromosome(record::VCFRecord)::String

Get the chromosome name of `record`.
"""
function chromosome(record::VCFRecord)::String
    checkfilled(record)
    if ismissing(record, record.chrom)
        missingerror(:chromosome)
    end
    return String(record.data[record.chrom])
end

function haschromosome(record::VCFRecord)
    return isfilled(record) && !ismissing(record, record.chrom)
end

"""
    leftposition(record::VCFRecord)::Int

Get the reference position of `record`.
"""
function Bio.leftposition(record::VCFRecord)::Int
    checkfilled(record)
    if ismissing(record, record.pos)
        missingerror(:leftposition)
    end
    # TODO: no-copy accessor
    return parse(Int, String(record.data[record.pos]))
end

function Bio.hasleftposition(record::VCFRecord)
    return isfilled(record) && !ismissing(record, record.pos)
end

"""
    identifier(record::VCFRecord)::Vector{String}

Get the identifiers of `record`.
"""
function identifier(record::VCFRecord)::Vector{String}
    checkfilled(record)
    if isempty(record.id)
        missingerror(:identifier)
    end
    return [String(record.data[r]) for r in record.id]
end

function hasidentifier(record::VCFRecord)
    return isfilled(record) && !isempty(record.id)
end

"""
    reference(record::VCFRecord)::String

Get the reference bases of `record`.
"""
function reference(record::VCFRecord)::String
    checkfilled(record)
    if ismissing(record, record.ref)
        missingerror(:reference)
    end
    return String(record.data[record.ref])
end

function hasreference(record::VCFRecord)
    return isfilled(record) && !ismissing(record, record.ref)
end

"""
    alternate(record::VCFRecord)::Vector{String}

Get the alternate bases of `record`.
"""
function alternate(record::VCFRecord)::Vector{String}
    checkfilled(record)
    if isempty(record.alt)
        missingerror(:alternate)
    end
    return [String(record.data[r]) for r in record.alt]
end

function hasalternate(record::VCFRecord)
    return isfilled(record) && !isempty(record.alt)
end

"""
    quality(record::VCFRecord)::Float64

Get the quality score of `record`.
"""
function quality(record::VCFRecord)::Float64
    checkfilled(record)
    if ismissing(record, record.qual)
        missingerror(:quality)
    end
    # TODO: no-copy parse
    return parse(Float64, String(record.data[record.qual]))
end

function hasquality(record::VCFRecord)
    return isfilled(record) && !ismissing(record, record.qual)
end

"""
    filter_(record::VCFRecord)::Vector{String}

Get the filter status of `record`.
"""
function filter_(record::VCFRecord)::Vector{String}
    checkfilled(record)
    if isempty(record.filter)
        missingerror(:filter_)
    end
    return [String(record.data[r]) for r in record.filter]
end

function hasfilter_(record::VCFRecord)
    return isfilled(record) && !isempty(record.filter)
end

"""
    information(record::VCFRecord)::Vector{Pair{String,String}}

Get the additional information of `record`.
"""
function information(record::VCFRecord)::Vector{Pair{String,String}}
    checkfilled(record)
    if isempty(record.infokey)
        missingerror(:information)
    end
    ret = Pair{String,String}[]
    for (i, key) in enumerate(record.infokey)
        val = infovalrange(record, i)
        push!(ret, String(record.data[key]) => String(record.data[val]))
    end
    return ret
end

function hasinformation(record::VCFRecord)
    return isfilled(record) && !isempty(record.infokey)
end

"""
    information(record::VCFRecord, key::String)::String

Get the additional information of `record` with `key`.
Keys without corresponding values return an empty string.
"""
function information(record::VCFRecord, key::String)::String
    checkfilled(record)
    i = findinfokey(record, key)
    if i == 0
        throw(KeyError(key))
    end
    val = infovalrange(record, i)
    if isempty(val)
        return ""
    else
        return String(record.data[val])
    end
end

function hasinformation(record::VCFRecord, key::String)
    return isfilled(record) && findinfokey(key) > 0
end

function findinfokey(record::VCFRecord, key::String)
    for i in 1:endof(record.infokey)
        if isequaldata(key, record.data, record.infokey[i])
            return i
        end
    end
    return 0
end

"""
    infokeys(record::VCFRecord)::Vector{String}

Get the keys of the additional information of `record`.
This function returns an empty vector when the INFO field is missing.
"""
function infokeys(record::VCFRecord)::Vector{String}
    checkfilled(record)
    return [String(record.data[key]) for key in record.infokey]
end

# Returns the data range of the `i`-th value.
function infovalrange(record::VCFRecord, i::Int)
    checkfilled(record)
    data = record.data
    key = record.infokey[i]
    if last(key) + 1 ≤ endof(data) && data[last(key)+1] == UInt8('=')
        endpos = search(data, ';', last(key)+1)
        if endpos == 0
            endpos = search(data, '\t', last(key)+1)
            @assert endpos != 0
        end
        return last(key)+2:endpos-1
    else
        return last(key)+1:last(key)
    end
end

"""
    format(record::VCFRecord)::Vector{String}

Get the genotype format of `record`.
"""
function format(record::VCFRecord)::Vector{String}
    checkfilled(record)
    if isempty(record.format)
        missingerror(:format)
    end
    return [String(record.data[r]) for r in record.format]
end

function hasformat(record::VCFRecord)
    return isfilled(record) && !isempty(record.format)
end

"""
    genotype(record::VCFRecord)::Vector{Vector{String}}

Get the genotypes of `record`.
"""
function genotype(record::VCFRecord)
    checkfilled(record)
    ret = Vector{String}[]
    for i in 1:endof(record.genotype)
        push!(ret, genotype_impl(record, i, 1:endof(record.format)))
    end
    return ret
end

"""
    genotype(record::VCFRecord, index::Integer)::Vector{String}

Get the genotypes of the `index`-th individual in `record`.
This is effectively equivalent to `genotype(record)[index]` but more efficient.
"""
function genotype(record::VCFRecord, index::Integer)
    checkfilled(record)
    return genotype_impl(record, index, 1:endof(record.format))
end

"""
    genotype(record::VCFRecord, indexes, keys)

Get the genotypes in `record` that match `indexes` and `keys`.
`indexes` and `keys` can be either a scalar or a vector value.
Trailing fields that are dropped are filled with `"."`.
"""
function genotype(record::VCFRecord, index::Integer, key::String)::String
    checkfilled(record)
    k = findgenokey(record, key)
    if k == 0
        throw(KeyError(key))
    end
    return genotype_impl(record, index, k)
end

function genotype(record::VCFRecord, index::Integer, keys::AbstractVector{String})::Vector{String}
    checkfilled(record)
    return [genotype(record, index, key) for key in keys]
end

function genotype{T<:Integer}(record::VCFRecord, indexes::AbstractVector{T}, key::String)::Vector{String}
    checkfilled(record)
    k = findgenokey(record, key)
    if k == 0
        throw(KeyError(key))
    end
    return [genotype_impl(record, i, k) for i in indexes]
end

function genotype{T<:Integer}(record::VCFRecord, indexes::AbstractVector{T}, keys::AbstractVector{String})::Vector{Vector{String}}
    checkfilled(record)
    ks = Vector{Int}(length(keys))
    for i in 1:endof(keys)
        key = keys[i]
        k = findgenokey(record, key)
        if k == 0
            throw(KeyError(key))
        end
        ks[i] = k
    end
    return [genotype_impl(record, i, ks) for i in indexes]
end

function genotype(record::VCFRecord, ::Colon, key::String)::Vector{String}
    return genotype(record, 1:endof(record.genotype), key)
end

function findgenokey(record::VCFRecord, key::String)
    return findfirst(r -> isequaldata(key, record.data, r), record.format)
end

function genotype_impl(record::VCFRecord, index::Int, keys::AbstractVector{Int})
    return [genotype_impl(record, index, k) for k in keys]
end

function genotype_impl(record::VCFRecord, index::Int, key::Int)
    geno = record.genotype[index]
    if key > endof(geno)  # dropped field
        return "."
    else
        return String(record.data[geno[key]])
    end
end

function Base.show(io::IO, record::VCFRecord)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   chromosome: ", haschromosome(record) ? chromosome(record) : "<missing>")
        println(io, "     position: ", hasleftposition(record) ? leftposition(record) : "<missing>")
        println(io, "   identifier: ", hasidentifier(record) ? join(identifier(record), " ") : "<missing>")
        println(io, "    reference: ", hasreference(record) ? reference(record) : "<missing>")
        println(io, "    alternate: ", hasalternate(record) ? join(alternate(record), " ") : "<missing>")
        println(io, "      quality: ", hasquality(record) ? quality(record) : "<missing>")
        println(io, "       filter: ", hasfilter_(record) ? join(filter_(record), " ") : "<missing>")
          print(io, "  information: ")
        if hasinformation(record)
            for (key, val) in information(record)
                print(io, key)
                if !isempty(val)
                    print(io, '=', val)
                end
                print(io, ' ')
            end
        else
            print(io, "<missing>")
        end
        println(io)
          print(io, "       format: ", hasformat(record) ? join(format(record), " ") : "<missing>")
        if hasformat(record)
            println(io)
            print(io, "     genotype:")
            for i in 1:endof(record.genotype)
                print(io, " [$(i)] ", let x = genotype(record, i); isempty(x) ? "." : join(x, " "); end)
            end
        end
    else
        print(io, " <not filled>")
    end
end

function ismissing(record::VCFRecord, range::UnitRange{Int})
    return length(range) == 1 && record.data[first(range)] == UInt8('.')
end

# Check if `str == data[range]`
function isequaldata(str::String, data::Vector{UInt8}, range::UnitRange{Int})
    rlen = length(range)
    return rlen == sizeof(str) && memcmp(pointer(data, first(range)), pointer(str), rlen) == 0
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Void}, Ptr{Void}, Csize_t), p1, p2, n)
end


text:
"""
The main `Module` containing the necessary types of functions for integration
with a Kroki service.

Defines `Base.show` and corresponding `Base.showable` methods for different
output formats and [`Diagram`](@ref Kroki.Diagram) types, so they render in
their most optimal form in different environments (e.g. the documentation
system, Documenter output, Jupyter, etc.).
"""
module Kroki

using Base64: base64encode
using CodecZlib: ZlibCompressor, transcode
using HTTP: request
using HTTP.ExceptionRequest: StatusError

include("./kroki/documentation.jl")
using .Documentation
@setupDocstringMarkup()

include("./kroki/service.jl")
using .Service: ENDPOINT

"""
A representation of a diagram that can be rendered by a Kroki service.

# Examples

```
julia> Kroki.Diagram(:PlantUML, "Kroki -> Julia: Hello Julia!")
     ┌─────┐          ┌─────┐
     │Kroki│          │Julia│
     └──┬──┘          └──┬──┘
        │ Hello Julia!   │
        │───────────────>│
     ┌──┴──┐          ┌──┴──┐
     │Kroki│          │Julia│
     └─────┘          └─────┘
```
"""
struct Diagram
  "The textual specification of the diagram"
  specification::AbstractString

  """
  The type of diagram specification (e.g. ditaa, Mermaid, PlantUML, etc.). This
  value is case-insensitive.
  """
  type::Symbol

  Diagram(type::Symbol, specification::AbstractString) = new(specification, type)
end

"""
An `Exception` to be thrown when a [`Diagram`](@ref) representing an invalid
specification is passed to [`render`](@ref).
"""
struct InvalidDiagramSpecificationError <: Exception
  error::String
  cause::Diagram
end

Base.showerror(io::IO, error::InvalidDiagramSpecificationError) = print(
  io,
  """
  $(RenderErrorHeader(error))

  This is (likely) caused by an invalid diagram specification.
  """,
)

"""
An `Exception` to be thrown when a [`Diagram`](@ref) is [`render`](@ref)ed to
an unsupported or invalid output format.
"""
struct InvalidOutputFormatError <: Exception
  error::String
  cause::Diagram
end

Base.showerror(io::IO, error::InvalidOutputFormatError) = print(
  io,
  """
  $(RenderErrorHeader(error))

  This is (likely) caused by an invalid or unknown output format.
  """,
)

# Helper function to render common headers when showing render errors
function RenderErrorHeader(
  error::Union{InvalidDiagramSpecificationError, InvalidOutputFormatError},
)
  """
  The Kroki service responded with:
  $(error.error)

  In response to a '$(error.cause.type)' diagram with the specification:
  $(error.cause.specification)
  """
end

"""
Compresses a [`Diagram`](@ref)'s `specification` using
[zlib](https://zlib.net), turning the resulting bytes into a URL-safe Base64
encoded payload (i.e. replacing `+` by `-` and `/` by `_`) to be used in
communication with a Kroki service.

See the [Kroki documentation](https://docs.kroki.io/kroki/setup/encode-diagram)
for more information.
"""
UriSafeBase64Payload(diagram::Diagram) = foldl(
  replace,
  ['+' => '-', '/' => '_'];
  init = base64encode(transcode(ZlibCompressor, diagram.specification)),
)

# Rewrites generic `HTTP.ExceptionRequest.StatusError`s into more specific
# errors based on Kroki's response if possible
function RenderError(diagram::Diagram, exception::StatusError)
  # Both errors related to invalid diagram specifications and invalid or
  # unsupported output formats are denoted by 400 responses, so further
  # processing of the response is necessary
  service_response = String(exception.response.body)

  if occursin("Unsupported output format", service_response)
    InvalidOutputFormatError(service_response, diagram)
  elseif occursin("Syntax Error", service_response)
    InvalidDiagramSpecificationError(service_response, diagram)
  else
    exception
  end
end
RenderError(::Diagram, exception::Exception) = exception

"""
Renders a [`Diagram`](@ref) through a Kroki service to the specified output
format.

If the Kroki service responds with an error throws an
[`InvalidDiagramSpecificationError`](@ref) or
[`InvalidOutputFormatError`](@ref) if a know type of error occurs. Other errors
(e.g. `HTTP.ExceptionRequest.StatusError` for connection errors) are propagated
if they occur.
"""
render(diagram::Diagram, output_format::AbstractString) =
  try
    getfield(
      request(
        "GET",
        join(
          [
            ENDPOINT[],
            lowercase("$(diagram.type)"),
            output_format,
            UriSafeBase64Payload(diagram),
          ],
          '/',
        ),
      ),
      :body,
    )
  catch exception
    throw(RenderError(diagram, exception))
  end

"""
Some MIME types are not supported by all diagram types, this constant contains
all these limitations. The union of all values corresponds to all supported
[`Diagram`](@ref) `type`s.
"""
const LIMITED_DIAGRAM_SUPPORT = Dict{AbstractString, Tuple{Symbol, Vararg{Symbol}}}(
  "application/pdf" => (
    :blockdiag,
    :seqdiag,
    :actdiag,
    :nwdiag,
    :packetdiag,
    :rackdiag,
    :erd,
    :graphviz,
    :vega,
    :vegalite,
  ),
  "image/jpeg" => (:c4plantuml, :erd, :graphviz, :plantuml, :umlet),
  "image/png" => (
    :blockdiag,
    :seqdiag,
    :actdiag,
    :nwdiag,
    :packetdiag,
    :rackdiag,
    :c4plantuml,
    :ditaa,
    :erd,
    :graphviz,
    :plantuml,
    :umlet,
    :vega,
    :vegalite,
  ),
  # Although all diagram types support SVG, these _only_ support SVG so are
  # included separately
  "image/svg+xml" => (:mermaid, :nomnoml, :svgbob, :wavedrom),
  "text/plain" => (:c4plantuml, :plantuml),
)

# `Base.show` methods should only be defined for diagram types that actually
# support the desired output format. This would make sure incompatible formats
# are not accidentally rendered on compatible `AbstractDisplay`s causing
# [`InvalidOutputFormatError`](@ref)s. As the diagram type information is only
# available within [`Diagram`](@ref) instances, the `show` method is defined
# generically, but then restricted using `Base.showable` to only those types
# that actually support the format
Base.show(io::IO, ::MIME{Symbol("image/png")}, diagram::Diagram) =
  write(io, render(diagram, "png"))
Base.showable(::MIME{Symbol("image/png")}, diagram::Diagram) =
  diagram.type ∈ LIMITED_DIAGRAM_SUPPORT["image/png"]

# SVG output is supported by _all_ diagram types, so there's no additional
# checking for support. This makes sure SVG output also works for new diagram
# types if they get added to Kroki, but not yet to this package
Base.show(io::IO, ::MIME"image/svg+xml", diagram::Diagram) =
  write(io, render(diagram, "svg"))

# PlantUML is capable of rendering textual representations, all other diagram
# types are not
Base.show(io::IO, diagram::Diagram) =
  if endswith(lowercase("$(diagram.type)"), "plantuml")
    write(io, render(diagram, "utxt"))
  else
    write(io, diagram.specification)
  end

for diagram_type in map(
  # The union of the values of `LIMITED_DIAGRAM_SUPPORT` correspond to all
  # supported `Diagram` types. Converting the `Symbol`s to `String`s improves
  # readability of the `macro` bodies
  String,
  collect(Set(Iterators.flatten(values(LIMITED_DIAGRAM_SUPPORT)))),
)
  macro_name = Symbol("$(diagram_type)_str")
  macro_signature = Symbol("@$macro_name")

  docstring = "Shorthand for instantiating $diagram_type [`Diagram`](@ref)s."

  @eval begin
    export $macro_signature

    @doc $docstring macro $macro_name(specification::AbstractString)
      Diagram(Symbol($diagram_type), specification)
    end
  end
end

end


text:
let cmd = `$(Base.julia_cmd()) --depwarn=error --startup-file=no runtests_exec.jl`
    for test_nthreads in (1, 2, 4) # run once to try single-threaded mode, then try a couple times to trigger bad races
        new_env = copy(ENV)
        new_env["JULIA_NUM_THREADS"] = string(test_nthreads)
        println("\n# Threads = $test_nthreads")
        run(pipeline(setenv(cmd, new_env), stdout = stdout, stderr = stderr))
    end
end


text:
"""
    interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::Vector{Float64}, Zq::Vector{Float64})

Interpolates the values given in `v` at the points `r, z` and interpolates
them at the query points 'Rq, Zq' and returns `Vq`.

`v`, `r`, and `z` are vectors of size N×1.
`Vq`, `Rq`, and `Rz` are vectors of size M×1.

"""
function interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::Vector{Float64}, Zq::Vector{Float64})
    Vq = zeros(Float64, length(Rq)) # init result vector

    # loop over every query point
    for i in 1:length(Rq)
        rq = Rq[i]; zq = Zq[i]; # extract current point
        corners = _findcorners(r, z, rq, zq); # get corners of interpolation quad

        # check if all corners were found
        if sum(corners .== 0) == 0
            X = _X(r, z, corners); # get matrix representation of corners
            V = [v[corners[1]]; v[corners[2]]; v[corners[3]]; v[corners[4]]]; # get vector of corner values
            Vq[i] = _interppoint(X, V, rq, zq); # add interpolated value to vector
        else
            # if no interpolation quad could be found
            Vq[i] = NaN;
        end
    end

    return Vq
end
interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::LinRange{Float64}, Zq::LinRange{Float64}) = interpolation2Dsym(r,   z,   v, collect(Rq), collect(Zq))
interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::LinRange{Float64}, Zq::Vector{Float64})   = interpolation2Dsym(r,   z,   v, collect(Rq), Zq)
interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::Vector{Float64},   Zq::LinRange{Float64}) = interpolation2Dsym(r,   z,   v, Rq,          collect(Zq))
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::Vector{Float64},   Zq::Vector{Float64})   = interpolation2Dsym(p.r, p.z, v, Rq,          Zq)
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::LinRange{Float64}, Zq::LinRange{Float64}) = interpolation2Dsym(p.r, p.z, v, collect(Rq), collect(Zq))
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::LinRange{Float64}, Zq::Vector{Float64})   = interpolation2Dsym(p.r, p.z, v, collect(Rq), Zq)
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::Vector{Float64},   Zq::LinRange{Float64}) = interpolation2Dsym(p.r, p.z, v, Rq,          collect(Zq))

"""
Finds the closest corner points of the enclosing quad element.
"""
function _findcorners(r::Vector{Float64}, z::Vector{Float64}, rq::Float64, zq::Float64)
    d = sqrt.( (rq .- r).^2 + (zq .- z).^2 );
    corners = zeros(Int64, 4);
    for i in 1:length(d)
        r[i]<rq && z[i]<zq ? (corners[1]==0 || d[i]<d[corners[1]] ? corners[1] = i : continue) :
        r[i]>rq && z[i]<zq ? (corners[2]==0 || d[i]<d[corners[2]] ? corners[2] = i : continue) :
        r[i]>rq && z[i]>zq ? (corners[3]==0 || d[i]<d[corners[3]] ? corners[3] = i : continue) :
        r[i]<rq && z[i]>zq ? (corners[4]==0 || d[i]<d[corners[4]] ? corners[4] = i : continue) : continue
    end
    return corners
end

"""
Interpolates at the query point `rq, zq` given the matrix representation of the
corner coordinates `X` (16×2) and the vector representation of the corner values
`V` (4×1).
"""
function _interppoint(X::Matrix{Float64}, V::Vector{Float64}, rq::Float64, zq::Float64; λ=0.1, rtol=1e-4, N=250)
    ξi = [0.0; 0.0] # initial guess in natural coordinates

    # loop to find natural coordinates corresponding to [rq; zq]
    counter = 1
    while true
        counter += 1; if counter > N; break; end

        Ni = _N(ξi) # shape functions
        xi = Ni * X[1:2:end, 1] # coordinates of guess ξi
        dNdξi = _dNdξ(ξi) # derivative of shape functions wrt. natural coordinates
        dxdξi = dNdξi * X # derivative of coordinates wrt. natural coordinates
        dξdxi = inv(dxdξi) # derivative of natural coordinates wrt. coordinates

        # update natural coordinate guess
        ξi += dξdxi * ([rq; zq] - xi) * λ

        # break condition
        sqrt((rq-xi[1])^2 + (zq-xi[2])^2) < rtol ? break : continue
    end

    # calculate interpolated value with shape functions
    N = transpose(_N(ξi)[1, 1:2:end])
    return N * V
end

"""
Generates the matrix representation of the corner coordinates `X` (16×2).
"""
function _X(r::Vector{Float64}, z::Vector{Float64}, corners::Vector{Int64})
    return [
        r[corners[1]] 0;
        0 r[corners[1]];
        z[corners[1]] 0;
        0 z[corners[1]];

        r[corners[2]] 0;
        0 r[corners[2]];
        z[corners[2]] 0;
        0 z[corners[2]];

        r[corners[3]] 0;
        0 r[corners[3]];
        z[corners[3]] 0;
        0 z[corners[3]];

        r[corners[4]] 0;
        0 r[corners[4]];
        z[corners[4]] 0;
        0 z[corners[4]];
    ]
end

"""
Generates the shape functions for the mapping of local coordinates (r, z) to
natural coordinates (ξ₁, ξ₂) as a matrix `N` (2×8).
"""
function _N(ξ::Vector{Float64})
    N1 = 1/4 * (1 - ξ[1]) * (1 - ξ[2]);
    N2 = 1/4 * (1 + ξ[1]) * (1 - ξ[2]);
    N3 = 1/4 * (1 + ξ[1]) * (1 + ξ[2]);
    N4 = 1/4 * (1 - ξ[1]) * (1 + ξ[2]);
    return [N1 0 N2 0 N3 0 N4 0; 0 N1 0 N2 0 N3 0 N4]
end

"""
Generates the matrix representation of the gradients of the shape functions wrt.
the natural coordinates (ξ₁, ξ₂) as `dNdξ` (2×16).
"""
function _dNdξ(ξ::Vector{Float64})
    dN1dξ1 = - 1/4 * (1 - ξ[2]); dN1dξ2 = - 1/4 * (1 - ξ[1]);
    dN2dξ1 = + 1/4 * (1 - ξ[2]); dN2dξ2 = - 1/4 * (1 - ξ[1]);
    dN3dξ1 = + 1/4 * (1 - ξ[2]); dN3dξ2 = + 1/4 * (1 - ξ[1]);
    dN4dξ1 = - 1/4 * (1 - ξ[2]); dN4dξ2 = + 1/4 * (1 - ξ[1]);

    return [
        dN1dξ1 dN1dξ2 0 0 dN2dξ1 dN2dξ2 0 0 dN3dξ1 dN3dξ2 0 0 dN4dξ1 dN4dξ2 0 0;
        0 0 dN1dξ1 dN1dξ2 0 0 dN2dξ1 dN2dξ2 0 0 dN3dξ1 dN3dξ2 0 0 dN4dξ1 dN4dξ2
    ]
end


text:
module math

import builtin

self.epsilon = 0.000001
self.pi = 3.14159265359
self.e = 2.71828182845

def min(args) {
	m = args[0]
  for i in args {
    if (i < m) {
      m = i
    }
  }
  m
}

def max(args) {
	m = args[0]
  for i in args {
    if (i > m) {
      m = i
    }
  }
  m
}

def pow(num, power) {
  builtin.pow__(num, power)
}

; if 1 arg is specified, defaults to natural log.
def log(args) {
  builtin.log__(args)
}

def ln(x) {
  builtin.log__(x)
}

def log10(x) {
  log(10, x)
}

def sqrt(x) pow(x, 0.5)


text:
"""
	mask_elements(img, threshold, connectivity)
Remove every element that is below the `threshold` and
adjacent to more than the `connectivity` amount of other
elements similarly below the threshold. Returns a mask where
all elements are 1 for elements that should be kept based on
`threshold` and `connectivity` criteria.

# Arguments
- img: array which will be used to determine mask
- threshold: cutoff value for elements that should be
    removed
- connectivity: minimum number of adjacent elements that are
    required to remove an element
"""
function mask_elements(img, threshold, connectivity)
    mask = img .< threshold
    label = Images.label_components(mask)
    n = length(unique(label))
    indices = []
    for lbl in 0:(n - 1)
        lbl_mask = map(x -> x == lbl, label)
        num = count(lbl_mask)
        if num < connectivity
            rmv_indices = findall(x -> x == lbl, label)
            push!(indices, rmv_indices)
        end
    end
    for i in 1:length(indices)
        idx = indices[i]
        mask[idx] .= 0
    end
    return iszero.(mask)
end


text:
"""
    linearfilter(N::BinaryNetwork, α::Vector{Float64}=[0.25, 0.25, 0.25, 0.25])

Given a network `N` compute the linear filter scores according to Stock et al.
(2017). High scores for negative interactions indicate potential false negative
or missing interactions. Though this it returned as a probabilistic network,
score do not necessary convey a probabilistic interpretation.
"""
function linearfilter(N::T; α::Vector{Float64}=fill(0.25, 4)) where {T <: BinaryNetwork}
  @assert length(α) == 4
  @assert all(α .≥ 0.0)
  α = α./sum(α) # This ensures that α sums to 1.0

  # Get the size of the network
  n = richness(N; dims=1)
  m = richness(N; dims=2)

  # Get the in and out degree, as well as the total number of interactions
  k_out = degree(N; dims=1)
  k_in = degree(N; dims=2)
  A = links(N)

  # Prepare a return object
  return_type = T <: AbstractBipartiteNetwork ? BipartiteProbabilisticNetwork : UnipartiteProbabilisticNetwork
  F = return_type(
    zeros(Float64, size(N)),
    EcologicalNetworks.species_objects(N)...
  )

  # Get probabilities
  for i in 1:richness(N; dims=1)
    s1 = species(N; dims=1)[i]
    for j in 1:richness(N; dims=2)
      s2 = species(N; dims=2)[j]
      t_val = α[1]*N[i,j] + (α[2]/m)*k_out[s1] + (α[3]/n)*k_in[s2] + α[4]/(n*m)*A
      F[s1,s2] = max(min(t_val, 1.0), 0.0)
    end
  end

  return F
end

"""
    linearfilterzoo(N::BinaryNetwork, α::Vector{Float64}=[0.25, 0.25, 0.25, 0.25])

Compute the zero-one-out version of the linear filter (`linearfilter`), i.e.
the score for each interaction if that interaction would not occur in the network.
For example, if N[4, 6] = 1 (interaction between species 4 and 6), the result at
postion (4, 6) is the score of the filter using a network in which that interaction
did not occur. This function is useful for validating the filter whether it can
detect false negative (missing) interactions.
"""
function linearfilterzoo(N::T; α::Vector{Float64}=fill(0.25, 4)) where {T <: BinaryNetwork}
  F = linearfilter(N, α=α)  # depart from scores of the filter
  Fzoo = F  # update the values to zero-one-out scores, reference is to keep
            # the code readable

  # Get the size of the network
  n = richness(N; dims=1)
  m = richness(N; dims=2)
  # Get probabilities
  for s1 in species(N; dims=1)
    for s2 in species(N; dims=2)
      t_val = F[s1,s2] - (α[1] + (α[2]/m) + (α[3]/n) + α[4]/(n * m)) * N[s1,s2]
      Fzoo[s1,s2] = max(min(t_val, one(t_val)), zero(t_val))
    end
  end

  return Fzoo
end


text:
using SymbolicUtils: FnType

const IndexMap = Dict{Char,Char}(
            '0' => '₀',
            '1' => '₁',
            '2' => '₂',
            '3' => '₃',
            '4' => '₄',
            '5' => '₅',
            '6' => '₆',
            '7' => '₇',
            '8' => '₈',
            '9' => '₉')
function map_subscripts(indices)
    str = string(indices)
    join(IndexMap[c] for c in str)
end

"""
$(TYPEDEF)

A named variable which represents a numerical value. The variable is uniquely
identified by its `name`, and all variables with the same `name` are treated
as equal.

# Fields
$(FIELDS)

For example, the following code defines an independent variable `t`, a parameter
`α`, a function parameter `σ`, a variable `x`, which depends on `t`, a variable
`y` with no dependents, a variable `z`, which depends on `t`, `α`, and `x(t)`
and parameters `β₁` and `β₂`.

```julia
t = Num(Variable{ModelingToolkit.Parameter{Real}}(:t))  # independent variables are treated as known
α = Num(Variable{ModelingToolkit.Parameter{Real}}(:α))  # parameters are known
σ = Num(Variable{ModelingToolkit.FnType{Tuple{Any},Real}}(:σ)) # left uncalled, since it is used as a function
w = Num(Variable{ModelingToolkit.FnType{Tuple{Any},Real}}(:w)) # unknown, left uncalled
x = Num(Variable{ModelingToolkit.FnType{Tuple{Any},Real}}(:x))(t)  # unknown, depends on `t`
y = Num(Variable(:y))   # unknown, no dependents
z = Num(Variable{ModelingToolkit.FnType{NTuple{3,Any},Real}}(:z))(t, α, x)  # unknown, multiple arguments
β₁ = Num(Variable(:β, 1)) # with index 1
β₂ = Num(Variable(:β, 2)) # with index 2

expr = β₁ * x + y^α + σ(3) * (z - t) - β₂ * w(t - 1)
```
"""
struct Variable{T} <: Function
    """The variable's unique name."""
    name::Symbol
    Variable(name) = Sym{Real}(name)
    Variable{T}(name) where T = Sym{T}(name)
    function Variable{T}(name, indices...) where T
        var_name = Symbol("$(name)$(join(map_subscripts.(indices), "ˏ"))")
        Sym{T}(var_name)
    end
end

function Variable(name, indices...)
    var_name = Symbol("$(name)$(join(map_subscripts.(indices), "ˏ"))")
    Variable(var_name)
end

rename(x::Sym{T},name) where T = Sym{T}(name)
rename(x::Term, name) where T = x.op isa Sym ? rename(x.op, name)(x.args...) : error("can't rename $x to $name")

# Build variables more easily
function _parse_vars(macroname, type, x)
    ex = Expr(:block)
    var_names = Symbol[]
    # if parsing things in the form of
    # begin
    #     x
    #     y
    #     z
    # end
    x = x isa Tuple && first(x) isa Expr && first(x).head == :tuple ? first(x).args : x # tuple handling
    x = flatten_expr!(x)
    for _var in x
        iscall = isa(_var, Expr) && _var.head == :call
        isarray = isa(_var, Expr) && _var.head == :ref
        issym  = _var isa Symbol
        @assert iscall || isarray || issym "@$macroname expects a tuple of expressions or an expression of a tuple (`@$macroname x y z(t) v[1:3] w[1:2,1:4]` or `@$macroname x, y, z(t) v[1:3] w[1:2,1:4]`)"

        if iscall
            var_name, expr = _construct_vars(_var.args[1], type, _var.args[2:end])
        else
            var_name, expr = _construct_vars(_var, type, nothing)
        end
        push!(var_names, var_name)
        push!(ex.args, expr)
    end
    rhs = build_expr(:tuple, var_names)
    push!(ex.args, :(($(var_names...),) = $rhs))
    return ex
end

function _construct_vars(_var, type, call_args)
    issym  = _var isa Symbol
    isarray = isa(_var, Expr) && _var.head == :ref
    if isarray
        var_name = _var.args[1]
        indices = _var.args[2:end]
        expr = _construct_array_vars(var_name, type, call_args, indices...)
    else
        # Implicit 0-args call
        var_name = _var
        expr = _construct_var(var_name, type, call_args)
    end
    var_name, :($var_name = $expr)
end

function _construct_var(var_name, type, call_args)
    expr = if call_args === nothing
        :($Num($Sym{$type}($(Meta.quot(var_name)))))
    elseif !isempty(call_args) && call_args[end] == :..
        :($Num($Sym{$FnType{Tuple, $type}}($(Meta.quot(var_name))))) # XXX: using Num as output
    else
        :($Num($Sym{$FnType{NTuple{$(length(call_args)), Any}, $type}}($(Meta.quot(var_name)))($(map(x->:($value($x)), call_args)...))))
    end
end

function _construct_var(var_name, type, call_args, ind)
    # TODO: just use Sym here
    if call_args === nothing
        :($Num($Variable{$type}($(Meta.quot(var_name)), $ind...)))
    elseif !isempty(call_args) && call_args[end] == :..
        :($Num($Variable{$FnType{Tuple{Any}, $type}}($(Meta.quot(var_name)), $ind...))) # XXX: using Num as output
    else
        :($Num($Variable{$FnType{NTuple{$(length(call_args)), Any}, $type}}($(Meta.quot(var_name)), $ind...)($(map(x->:($value($x)), call_args)...))))
    end
end

function _construct_array_vars(var_name, type, call_args, indices...)
    :(map(Iterators.product($(indices...))) do ind
        $(_construct_var(var_name, type, call_args, :ind))
    end)
end

"""
$(SIGNATURES)

Define one or more unknown variables.

```julia
@parameters t α σ(..) β[1:2]
@variables w(..) x(t) y() z(t, α, x)

expr = β₁* x + y^α + σ(3) * (z - t) - β₂ * w(t - 1)
```

Note that `@parameters` and `@variables` implicitly add `()` to values that
are not given a call. The former specifies the values as known, while the
latter specifies it as unknown. `(..)` signifies that the value should be
left uncalled.

Sometimes it is convenient to define arrays of variables to model things like `x₁,…,x₃`.
The `@variables` and `@parameters` macros support this with the following syntax:

```julia
@variables x[1:3];
x

3-element Array{Operation,1}:
 x₁()
 x₂()
 x₃()

# support for arbitrary ranges and tensors
@variables y[2:3,1:5:6];
y

2×2 Array{Operation,2}:
    y₂̒₁() y₂̒₆()
    y₃̒₁() y₃̒₆()

# also works for dependent variables
@parameters t; @variables z[1:3](t);
z

3-element Array{Operation,1}:
 z₁(t())
 z₂(t())
 z₃(t())
```
"""
macro variables(xs...)
    esc(_parse_vars(:variables, Real, xs))
end

"""
$(TYPEDSIGNATURES)

Renames the variable `x` to have `name`.
"""
function rename(x::Variable,name::Symbol)
    Variable{symtype(x)}(name)
end

TreeViews.hastreeview(x::Variable) = true
function TreeViews.treelabel(io::IO,x::Variable,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Text(x.name))
end

"""
varmap_to_vars(varmap,varlist)

Takes a list of pairs of variables=>values and an ordered list of variables and
creates the array of values in the correct order
"""
function varmap_to_vars(varmap::AbstractArray{<:Pair},varlist)
    out = similar(varmap,typeof(last(first(varmap))))
    for (ivar, ival) in varmap
        j = findfirst(isequal(ivar),varlist)
        if isnothing(j)
            throw(ArgumentError("Value $(ivar) provided in map not found in $(varlist)"))
        elseif j > length(varmap)
            throw(ArgumentError("Missing value in $(varmap), need $(varlist)"))
        end
        out[j] = ival
    end

    # Make output match varmap in type and shape
    # Does things like MArray->SArray
    ArrayInterface.restructure(varmap,out)
end
varmap_to_vars(varmap::AbstractArray,varlist) = varmap
varmap_to_vars(varmap::DiffEqBase.NullParameters,varlist) = varmap
varmap_to_vars(varmap::Nothing,varlist) = varmap


text:
module Impute

using BSON
using CSV
using DataDeps
using Distances
using IterTools
using Missings
using NamedDims
using NearestNeighbors
using Random
using Statistics
using StatsBase
using TableOperations
using Tables: Tables, materializer, istable

using Base.Iterators
using LinearAlgebra
using LinearAlgebra: Diagonal

include("utils.jl")
include("declaremissings.jl")
include("imputors.jl")
include("filter.jl")
include("validators.jl")
include("chain.jl")
include("deprecated.jl")

include("functional.jl")
include("data.jl")

__init__() = register_datadep()

end  # module


text:
using Documenter, SixDOF

makedocs(;
    modules=[SixDOF],
    format=Documenter.HTML(),
    pages=[
        "Guide" => "guide.md",
        "Theory" => "theory.md",
    ],
    repo="https://github.com/byuflowlab/SixDOF.jl/blob/{commit}{path}#L{line}",
    sitename="SixDOF.jl",
    authors="Andrew Ning <aning@byu.edu>",
    # assets=String[],
)

deploydocs(
    repo = "github.com/byuflowlab/SixDOF.jl.git",
)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

"""
    get_beam_forces_vector_2d(X1, X2, qt, qn, f)

Function integrates forces vector for 6 DOF Euler-Bernoulli beam element in 2D.

X1 = beams left node coordinates
X2 = beams right node coordinates
qt = Tangential uniformly distributed load
qn = Normal uniformly distributed load
f = Point forces vector in global coordinates
"""
function get_beam_forces_vector_2d(X1, X2, qt, qn, f)
    le=norm(X2-X1)                      # Lenght of element
    a=atan((X2[2]-X1[2])/(X2[1]-X1[1])) # Rotation angle of the element
    nn=2                            # Number of nodes
    nd=3*1+3                           # Number of DOFs
    Gp=[-1/3*sqrt(5+2*sqrt(10/7)) -1/3*sqrt(5-2*sqrt(10/7)) 0 1/3*sqrt(5-2*sqrt(10/7)) 1/3*sqrt(5+2*sqrt(10/7))] # Five Gauss integration points
    w=[(322-13*sqrt(70))/900 (322+13*sqrt(70))/900 128/225 (322+13*sqrt(70))/900 (322-13*sqrt(70))/900] # Five Gauss integration weights
    # Rotation matrix
    B=[cos(a) sin(a) 0  0     0       0;
      -sin(a) cos(a) 0  0     0       0;
       0      0      1  0     0       0;
       0      0      0  cos(a) sin(a) 0;
       0      0      0 -sin(a) cos(a) 0;
       0      0      0  0      0      1]
    # Integration of the equivalent forces vector
    # For truss element tfq
    detJ=le/2
    function tfqint(wi::Float64, xi::Float64)
        N1s=1+(-1/2)*(1+xi)
        N2s=(1/2)*(1 + xi)
        Ns=[N1s N2s]
        return wi*Ns'
    end
    tfq=zeros(2,1)
    for i = 1:size(Gp,2)
        tfq +=tfqint(w[i],Gp[i])
    end
    tfq=qn*tfq*detJ
    # For 4 DOF beam element bfqe
    detJ=le/2
    function bfqint(wi::Float64, xi::Float64)
        N1=1/4*(1-xi)^2*(2+xi)
        N2=le/8*(1-xi)^2*(xi+1)
        N3=1/4*(1+xi)^2*(2-xi)
        N4=le/8*(1+xi)^2*(xi-1)
        N=[N1 N2 N3 N4]
        return wi*N'
    end
    bfq=zeros(4,1)
    for i = 1:size(Gp,2)
        bfq +=bfqint(w[i],Gp[i])
    end
    bfq=qt*bfq*detJ
    # Assembly of the 6 DOF beam element equivalent forces vector fqe
    fqe=zeros(6,1)
    fqe[1,1],fqe[4,1]=tfq[1,1],tfq[2,1]
    fqe[2,1],fqe[3,1],fqe[5,1],fqe[6,1]=bfq[1,1],bfq[2,1],bfq[3,1],bfq[4,1]
    # Rotation
    fqe=B'*fqe
    # Adding equivalent forces vector to point forces vector
    f +=fqe
    return f
end


text:
"""
    ceil_val(x)

Finds the ceiling of x as an Integer

# Example
```julia
ceil_val(1.3)   # 2.0
ceil_val(2.0)   # returns 2.0
ceil_val(-1.5)  #returns -1.0
````

# Reference
- https://en.wikipedia.org/wiki/Floor_and_ceiling_functions

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function ceil_val(x)
    return x -trunc(x) <= 0 ? trunc(x) : trunc(x)+1
end

"""
    floor_val(x)

Finds the floor of x as an Integer

# Example
```julia
floor_val(1.3)  # 1
floor_val(2.0)  # returns 2.0
floor_val(-1.7) # returns -2.0
````

# Reference
- https://en.wikipedia.org/wiki/Floor_and_ceiling_functions

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function floor_val(x)
    return x -trunc(x) >= 0 ? trunc(x) : trunc(x)-1
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia